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)