docs: def eigendecomp(G, d):
      """
      Computes sorted eigendecomposition of G.
  
      :param G:  Matrix 
      :param d:  rank
  
      :return factor: vector of square root d of eigenvalues (biggest to smallest). 
      :return u: matrix with colums equal to the normalized eigenvectors corresponding to sorted eigenvalues.
      """
      N = G.shape[0]
      lamda, u = np.linalg.eig(G)
      # test decomposition of G.
      #G_hat = np.dot(np.dot(u, np.diag(lamda)), u.T)
      #assert np.linalg.norm(G_hat - G) < 1e-10, 'decomposition not exact: err {}'.format(np.linalg.norm(G_hat - G))
  
      # sort the eigenvalues in decreasing order
      lamda = np.real(lamda)
      indices = np.argsort(lamda)[::-1]
      lamda_sorted = lamda[indices]
      assert (lamda_sorted[
              :d] > -1e-10).all(), "{} not all positive!".format(lamda_sorted[:d])
  
      u = u[:, indices]
      factor = np.empty((N,), dtype=lamda.dtype)
      np.sqrt(lamda_sorted[:d], out=factor[0:d])
      factor[d:] = 0.0
      return factor, np.real(u)