def compute_covariance(L_aug, Y, k, p):
    """Given label matrix L_aug and labels Y, compute the covariance.

    Args:
        L: (np.array {0,1}) [n, d] The augmented (indicator) label matrix
        Y: (np.array int) [n] The true labels in {1,...,k}
        k: (int) Cardinality
        p: (np.array float) [k] The class balance
    """
    n, d = L_aug.shape
    assert Y.shape[0] == n
    mu = compute_mu(L_aug, Y, k, p)
    return (L_aug.T @ L_aug) / n - mu @ np.diag(p) @ mu.T
Given label matrix L_aug and labels Y, compute the covariance.

    Args:
        L: (np.array {0,1}) [n, d] The augmented (indicator) label matrix
        Y: (np.array int) [n] The true labels in {1,...,k}
        k: (int) Cardinality
        p: (np.array float) [k] The class balance