code: def eigendecomp(G, d):
      N = G.shape[0]
      lamda, u = np.linalg.eig(G)
      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)