code: def compute_covariance(L_aug, Y, k, p):
      
      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
docs: 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