Python scipy.sparse.linalg.eigsh() Examples

def solve_modal(model,k:int):
    """
    Solve eigen mode of the MDOF system
    
    params:
        model: FEModel.
        k: number of modes to extract.
    """
    K_,M_=model.K_,model.M_
    if k>model.DOF:
        logger.info('Warning: the modal number to extract is larger than the system DOFs, only %d modes are available'%model.DOF)
        k=model.DOF
    omega2s,modes = sl.eigsh(K_,k,M_,sigma=0,which='LM')
    delta = modes/np.sum(modes,axis=0)
    model.is_solved=True
    model.mode_=delta
    model.omega_=np.sqrt(omega2s).reshape((k,1)) 

def test_eigs_truss():
    """Eigenvalues of a bar"""
    nnodes = 513
    
    x = np.linspace(0, np.pi, nnodes)
    nodes = np.zeros((nnodes, 3))
    nodes[:, 0] = range(nnodes)
    nodes[:, 1] = x
    cons = np.zeros((nnodes, 2))
    cons[:, 1] = -1
    cons[0, 0] = -1 
    cons[-1, 0] = -1
    mats = np.array([[1.0, 1.0, 1.0]])
    elements = np.zeros((nnodes - 1, 5 ), dtype=int)
    elements[:, 0] = range(nnodes - 1)
    elements[:, 1] = 6
    elements[:, 3] = range(nnodes - 1)
    elements[:, 4] = range(1, nnodes)
    
    assem_op, bc_array, neq = ass.DME(cons, elements)
    stiff, mass = ass.assembler(elements, mats, nodes, neq, assem_op)
    
    vals, _ = eigsh(stiff, M=mass, which="SM")
    assert np.allclose(vals, np.linspace(1, 6, 6)**2, rtol=1e-2) 

def manifold_harmonics(verts, tris, K, scaled=True, return_D=False, return_eigenvalues=False):
    Q, vertex_area = compute_mesh_laplacian(
        verts, tris, 'cotangent', 
        return_vertex_area=True, area_type='lumped_mass'
    )
    if scaled:
        D = sparse.spdiags(vertex_area, 0, len(verts), len(verts))
    else:
        D = sparse.spdiags(np.ones_like(vertex_area), 0, len(verts), len(verts))

    try:
        lambda_dense, Phi_dense = eigsh(-Q, M=D, k=K, sigma=0)
    except RuntimeError, e:
        if e.message == 'Factor is exactly singular':
            logging.warn("factor is singular, trying some regularization and cholmod")
            chol_solve = factorized(-Q + sparse.eye(Q.shape[0]) * 1.e-9)
            OPinv = sparse.linalg.LinearOperator(Q.shape, matvec=chol_solve)
            lambda_dense, Phi_dense = eigsh(-Q, M=D, k=K, sigma=0, OPinv=OPinv)
        else:
            raise e 

def _method_1(data, num_pcs=None):
    """Compute OPCA when num_observations > num_dimensions."""
    data = np.nan_to_num(data - nanmean(data, axis=0))
    T = data.shape[0]
    corr_offset = np.dot(data[1:].T, data[:-1])
    corr_offset += corr_offset.T
    if num_pcs is None:
        eivals, eivects = eigh(corr_offset)
    else:
        eivals, eivects = eigsh(corr_offset, num_pcs, which='LA')
    eivals = np.real(eivals)
    eivects = np.real(eivects)
    idx = np.argsort(-eivals)  # sort the eigenvectors and eigenvalues
    eivals = old_div(eivals[idx], (2. * (T - 1)))
    eivects = eivects[:, idx]
    return eivals, eivects, np.dot(data, eivects) 

def _cov_eigen(self, X):
        """
        Perform direct computation of covariance matrix eigen{values,vectors},
        given a scipy.sparse matrix.

        Parameters
        ----------
        X : WRITEME

        Returns
        -------
        WRITEME
        """

        v, W = eigen_symmetric(X.T.dot(X) / X.shape[0], k=self.num_components)

        # The resulting components are in *ascending* order of eigenvalue, and
        # W contains eigenvectors in its *columns*, so we simply reverse both.
        return v[::-1], W[:, ::-1] 

def ldos0d_wf(h,e=0.0,delta=0.01,num_wf = 10,robust=False,tol=0):
  """Calculates the local density of states of a hamiltonian and
     writes it in file, using arpack"""
  if h.dimensionality==0:  # only for 0d
    intra = csc_matrix(h.intra) # matrix
  else: raise # not implemented...
  if robust: # go to the imaginary axis for stability
    eig,eigvec = slg.eigs(intra,k=int(num_wf),which="LM",
                        sigma=e+1j*delta,tol=tol) 
    eig = eig.real # real part only
  else: # Hermitic Hamiltonian
    eig,eigvec = slg.eigsh(intra,k=int(num_wf),which="LM",sigma=e,tol=tol) 
  d = np.array([0.0 for i in range(intra.shape[0])]) # initialize
  for (v,ie) in zip(eigvec.transpose(),eig): # loop over wavefunctions
    v2 = (np.conjugate(v)*v).real # square of wavefunction
    fac = delta/((e-ie)**2 + delta**2) # factor to create a delta
    d += fac*v2 # add contribution
#  d /= num_wf # normalize
  d /= np.pi # normalize
  d = spatial_dos(h,d) # resum if necessary
  g = h.geometry  # store geometry
  write_ldos(g.x,g.y,d,z=g.z) # write in file 

def ldos_arpack(intra,num_wf=10,robust=False,tol=0,e=0.0,delta=0.01):
  """Use arpack to calculate hte local density of states at a certain energy"""
  if robust: # go to the imaginary axis for stability
    eig,eigvec = slg.eigs(intra,k=int(num_wf),which="LM",
                        sigma=e+1j*delta,tol=tol) 
    eig = eig.real # real part only
  else: # Hermitic Hamiltonian
    eig,eigvec = slg.eigsh(intra,k=int(num_wf),which="LM",sigma=e,tol=tol) 
  d = np.array([0.0 for i in range(intra.shape[0])]) # initialize
  for (v,ie) in zip(eigvec.transpose(),eig): # loop over wavefunctions
    v2 = (np.conjugate(v)*v).real # square of wavefunction
    fac = delta/((e-ie)**2 + delta**2) # factor to create a delta
    d += fac*v2 # add contribution
#  d /= num_wf # normalize
  d /= np.pi # normalize
  return d 

def test_cXY_E0(nr_sites, gamma, rgen, ldim=2):
    if sys.version_info[:2] == (3, 3) and gamma == -0.5:
        # Skip this test on Python 3.3 because it fails on Travis (but
        # only for Python 3.3). eigsh() fails with:
        # scipy.sparse.linalg.eigen.arpack.arpack.ArpackNoConvergence:
        # ARPACK error -1: No convergence (xxx iterations, 0/1
        # eigenvectors converged) [ARPACK error -14: ZNAUPD did not
        # find any eigenvalues to sufficient accuracy.]
        pt.skip("Test fails on Travis for unknown reason")
        return

    # Verify that the analytical solution of the ground state energy
    # matches the numerical value from eigsh()
    E0 = physics.cXY_E0(nr_sites, gamma)
    H = physics.sparse_cH(physics.cXY_local_terms(nr_sites, gamma))
    # Fix start vector for eigsh()
    v0 = rgen.randn(ldim**nr_sites) + 1j * rgen.randn(ldim**nr_sites)
    ev = eigsh(H, k=1, which='SR', v0=v0, return_eigenvectors=False).min()
    assert abs(E0 - ev) <= 1e-13 

def _calculate_sf(self, graph):
        """
        Calculating the features of a graph.

        Arg types:
            * **graph** *(NetworkX graph)* - A graph to be embedded.

        Return types:
            * **embedding** *(Numpy array)* - The embedding of a single graph.
        """
        number_of_nodes = graph.number_of_nodes()
        L_tilde = nx.normalized_laplacian_matrix(graph, nodelist=range(number_of_nodes))
        if number_of_nodes <= self.dimensions:
            embedding = eigsh(L_tilde, k=number_of_nodes-1, which='LM',
                              ncv=10*self.dimensions, return_eigenvectors=False)

            shape_diff = self.dimensions - embedding.shape[0] - 1
            embedding = np.pad(embedding, (1, shape_diff), 'constant', constant_values=0)
        else:
            embedding = eigsh(L_tilde, k=self.dimensions, which='LM',
                              ncv=10*self.dimensions, return_eigenvectors=False)
        return embedding 

def occupied_states(hkgen,k,window=None,max_waves=None):
  """ Returns the WF of the occupied states in a 2d hamiltonian"""
  hk = hkgen(k) # get hamiltonian
  if max_waves is None: es,wfs = algebra.eigh(hk) # diagonalize all waves
  else:  es,wfs = slg.eigsh(csc_matrix(hk),k=max_waves,which="SA",
                      sigma=0.0,tol=arpack_tol,maxiter=arpack_maxiter)
  wfs = np.conjugate(wfs.transpose()) # wavefunctions
  occwf = []
  for (ie,iw) in zip(es,wfs):  # loop over states
    if window is None: # no energy window
      if ie < 0:  # if below fermi
        occwf.append(iw)  # add to the list
    else: # energy window provided
      if -abs(window)< ie < 0:  # between energy window and fermi
        occwf.append(iw)  # add to the list
  return np.array(occwf) 

def eigsh(A, *args, **kwargs):
    return _eigsh(A, *args, **kwargs) 

def fsr_rankR(qsims, Wn, alpha = 0.99, R = 2000):
    vals, vecs = s_linalg.eigsh(Wn, k = R)
    p2 = diags((1.0 - alpha) / (1.0 - alpha*vals))
    vc = csr_matrix(vecs)
    p3 =  vc.dot(p2)
    vc_norm =  (vc.multiply(vc)).sum(axis = 0)
    out_sims = []
    for i in range(qsims.shape[0]):
        qsims_sparse = csr_matrix(qsims[i:i+1,:])
        p1 =(vc.T).dot(qsims_sparse.T)
        diff_sim = csr_matrix(p3)*csr_matrix(p1)
        out_sims.append(diff_sim.todense().reshape(-1,1))
    out_sims = np.concatenate(out_sims, axis = 1)
    ranks = np.argsort(-out_sims, axis = 0)
    return ranks 

def _svd(self, array, n_components, n_discard):
        """Returns first `n_components` left and right singular
        vectors u and v, discarding the first `n_discard`.

        """
        if self.svd_method == 'randomized':
            kwargs = {}
            if self.n_svd_vecs is not None:
                kwargs['n_oversamples'] = self.n_svd_vecs
            u, _, vt = randomized_svd(array, n_components,
                                      random_state=self.random_state,
                                      **kwargs)

        elif self.svd_method == 'arpack':
            u, _, vt = svds(array, k=n_components, ncv=self.n_svd_vecs)
            if np.any(np.isnan(vt)):
                # some eigenvalues of A * A.T are negative, causing
                # sqrt() to be np.nan. This causes some vectors in vt
                # to be np.nan.
                A = safe_sparse_dot(array.T, array)
                random_state = check_random_state(self.random_state)
                # initialize with [-1,1] as in ARPACK
                v0 = random_state.uniform(-1, 1, A.shape[0])
                _, v = eigsh(A, ncv=self.n_svd_vecs, v0=v0)
                vt = v.T
            if np.any(np.isnan(u)):
                A = safe_sparse_dot(array, array.T)
                random_state = check_random_state(self.random_state)
                # initialize with [-1,1] as in ARPACK
                v0 = random_state.uniform(-1, 1, A.shape[0])
                _, u = eigsh(A, ncv=self.n_svd_vecs, v0=v0)

        assert_all_finite(u)
        assert_all_finite(vt)
        u = u[:, n_discard:]
        vt = vt[n_discard:]
        return u, vt.T 

def test_eigs_beam():
    """Eigenvalues of a cantilever beam"""
    
    nnodes = 10

    x = np.linspace(0, np.pi, nnodes)
    nodes = np.zeros((nnodes, 3))
    nodes[:, 0] = range(nnodes)
    nodes[:, 1] = x
    cons = np.zeros((nnodes, 3))
    cons[0, :] = -1
    cons[:, 0] = -1
    mats = np.array([[1.0, 1.0, 1.0, 1.0]])
    elements = np.zeros((nnodes - 1, 5 ), dtype=int)
    elements[:, 0] = range(nnodes - 1)
    elements[:, 1] = 7
    elements[:, 3] = range(nnodes - 1)
    elements[:, 4] = range(1, nnodes)
    
    assem_op, bc_array, neq = ass.DME(cons, elements, ndof_node=3)
    stiff, mass = ass.assembler(elements, mats, nodes, neq, assem_op)
    
    vals, _ = eigsh(stiff, M=mass, which="SM")
    vals_analytic = np.array([0.596864162694467, 1.49417561427335,
                              2.50024694616670,  3.49998931984744,
                              4.50000046151508, 5.49999998005609])
    assert np.allclose(vals**0.25, vals_analytic, rtol=1e-2) 

def _get_fiedler_func(method):
    """Return a function that solves the Fiedler eigenvalue problem.
    """
    match = _tracemin_method.match(method)
    if match:
        method = match.group(1)
        def find_fiedler(L, x, normalized, tol):
            q = 2 if method == 'pcg' else min(4, L.shape[0] - 1)
            X = asmatrix(normal(size=(q, L.shape[0]))).T
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
            return sigma[0], X[:, 0]
    elif method == 'lanczos' or method == 'lobpcg':
        def find_fiedler(L, x, normalized, tol):
            L = csc_matrix(L, dtype=float)
            n = L.shape[0]
            if normalized:
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
                L = D * L * D
            if method == 'lanczos' or n < 10:
                # Avoid LOBPCG when n < 10 due to
                # https://github.com/scipy/scipy/issues/3592
                # https://github.com/scipy/scipy/pull/3594
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
                                 return_eigenvectors=True)
                return sigma[1], X[:, 1]
            else:
                X = asarray(asmatrix(x).T)
                M = spdiags(1. / L.diagonal(), [0], n, n)
                Y = ones(n)
                if normalized:
                    Y /= D.diagonal()
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
                                  maxiter=n, largest=False)
                return sigma[0], X[:, 0]
    else:
        raise nx.NetworkXError("unknown method '%s'." % method)

    return find_fiedler 

def get_train(self):
        lap_mat = self.getLap()
        print('finish getLap...')
        w, vec = eigsh(lap_mat, k=self.rep_size)
        print('finish eigh(lap_mat)...')
        # start = 0
        # for i in range(self.node_size):
        #     if w[i] > 1e-10:
        #         start = i
        #         break
        # vec = vec[:, start:start+self.rep_size]

        return vec 

def normcut_vectors(affinity_matrix, k):
    """Return the normalized cut vectors.

    These vectors satisfy are the largest eigenvalues of $D^{-1/2} W D^{-1/2}$,
    or equivalently the smallest of $D^{-1/2} (D - W) D^{-1/2}$. The first
    eigenvalue should be constant in all entries, but the second eigenvalue
    can be used to determine the normalized cut. See Shi & Malik, 2000.

    Parameters
    ----------
    A : array
        The affinity matrix containing the weight between graph nodes.
    k : int
        The number of cut eigenvectors to return.

    Returns
    -------
    array
        The normcut vectors.  Shape (num_nodes, k).

    """

    node_degrees = np.array(affinity_matrix.sum(axis=0)).flatten()
    transformation_matrix = diags(np.sqrt(old_div(1., node_degrees)), 0)
    normalized_affinity_matrix = transformation_matrix * affinity_matrix * \
        transformation_matrix
    _, vects = eigsh(normalized_affinity_matrix, k + 1, sigma=1.001,
                     which='LM')  # Get the largest eigenvalues.
    cuts = transformation_matrix * vects
    return cuts 

def smalleig(m,numw=10,evecs=False):
  """
  Return the smallest eigenvalues using arpack
  """
  tol = arpack_tol
  eig,eigvec = slg.eigsh(m,k=numw,which="LM",sigma=0.0,
                                  tol=tol,maxiter=arpack_maxiter)
  if evecs:  return eig,eigvec.transpose()  # return eigenvectors
  else:  return eig  # return eigenvalues 

def _get_fiedler_func(method):
    """Return a function that solves the Fiedler eigenvalue problem.
    """
    match = _tracemin_method.match(method)
    if match:
        method = match.group(1)

        def find_fiedler(L, x, normalized, tol):
            q = 2 if method == 'pcg' else min(4, L.shape[0] - 1)
            X = asmatrix(normal(size=(q, L.shape[0]))).T
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
            return sigma[0], X[:, 0]
    elif method == 'lanczos' or method == 'lobpcg':
        def find_fiedler(L, x, normalized, tol):
            L = csc_matrix(L, dtype=float)
            n = L.shape[0]
            if normalized:
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
                L = D * L * D
            if method == 'lanczos' or n < 10:
                # Avoid LOBPCG when n < 10 due to
                # https://github.com/scipy/scipy/issues/3592
                # https://github.com/scipy/scipy/pull/3594
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
                                 return_eigenvectors=True)
                return sigma[1], X[:, 1]
            else:
                X = asarray(asmatrix(x).T)
                M = spdiags(1. / L.diagonal(), [0], n, n)
                Y = ones(n)
                if normalized:
                    Y /= D.diagonal()
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
                                  maxiter=n, largest=False)
                return sigma[0], X[:, 0]
    else:
        raise nx.NetworkXError("unknown method '%s'." % method)

    return find_fiedler 

def _get_fiedler_func(method):
    """Returns a function that solves the Fiedler eigenvalue problem.
    """
    if method == "tracemin":  # old style keyword <v2.1
        method = "tracemin_pcg"
    if method in ("tracemin_pcg", "tracemin_chol", "tracemin_lu"):
        def find_fiedler(L, x, normalized, tol, seed):
            q = 1 if method == 'tracemin_pcg' else min(4, L.shape[0] - 1)
            X = asmatrix(seed.normal(size=(q, L.shape[0]))).T
            sigma, X = _tracemin_fiedler(L, X, normalized, tol, method)
            return sigma[0], X[:, 0]
    elif method == 'lanczos' or method == 'lobpcg':
        def find_fiedler(L, x, normalized, tol, seed):
            L = csc_matrix(L, dtype=float)
            n = L.shape[0]
            if normalized:
                D = spdiags(1. / sqrt(L.diagonal()), [0], n, n, format='csc')
                L = D * L * D
            if method == 'lanczos' or n < 10:
                # Avoid LOBPCG when n < 10 due to
                # https://github.com/scipy/scipy/issues/3592
                # https://github.com/scipy/scipy/pull/3594
                sigma, X = eigsh(L, 2, which='SM', tol=tol,
                                 return_eigenvectors=True)
                return sigma[1], X[:, 1]
            else:
                X = asarray(asmatrix(x).T)
                M = spdiags(1. / L.diagonal(), [0], n, n)
                Y = ones(n)
                if normalized:
                    Y /= D.diagonal()
                sigma, X = lobpcg(L, X, M=M, Y=asmatrix(Y).T, tol=tol,
                                  maxiter=n, largest=False)
                return sigma[0], X[:, 0]
    else:
        raise nx.NetworkXError("unknown method '%s'." % method)

    return find_fiedler 

def ev2d(h,nk=50,nsuper=1,reciprocal=False,
               operator=None,k0=[0.,0.],kreverse=False):
  """ Calculate the expectation value of a certain operator"""
  if h.dimensionality!=2: raise  # continue if two dimensional
  hk_gen = h.get_hk_gen() # gets the function to generate h(k)
  kxs = np.linspace(-nsuper,nsuper,nk,endpoint=True)+k0[0]  # generate kx
  kys = np.linspace(-nsuper,nsuper,nk,endpoint=True)+k0[1]  # generate ky
  if kreverse: kxs,kys = -kxs,-kys
  kdos = [] # empty list
  kxout = []
  kyout = []
  if reciprocal: R = h.geometry.get_k2K() # get matrix
  else:  R = np.matrix(np.identity(3)) # get identity
  # setup the operator
  operator = operator2list(operator) # convert into a list
  fo = open("EV2D.OUT","w") # open file
  for x in kxs:
    for y in kxs:
      print("Doing",x,y)
      r = np.matrix([x,y,0.]).T # real space vectors
      k = np.array((R*r).T)[0] # change of basis
      hk = hk_gen(k) # get hamiltonian
      if not h.is_sparse: evals,waves = lg.eigh(hk) # eigenvalues
      else: evals,waves = slg.eigsh(hk,k=max(nindex)*2,sigma=0.0,
             tol=arpack_tol,which="LM") # eigenvalues
      waves = waves.transpose() # transpose
      eneg,wfneg = [],[] # negative
      for (e,w) in zip(evals,waves): # loop
        if e<0: # negative
          eneg.append(e)
          wfneg.append(w)
      fo.write(str(x)+"     "+str(y)+"   ") # write k-point
      for op in operator: # loop over operators
          c = sum([braket_wAw(w,op) for w in wfneg]).real # expectation value
          fo.write(str(c)+"  ") # write in file
      fo.write("\n") # write in file
  fo.close() # close file 

def total_energy(h,nk=10,nbands=None,use_kpm=False,random=False,
        kp=None,mode="mesh",tol=1e-1):
  """Return the total energy"""
  h.turn_dense()
  if h.is_sparse and not use_kpm: 
    print("Sparse Hamiltonian but no bands given, taking 20")
    nbands=20
  f = h.get_hk_gen() # get generator
  etot = 0.0 # initialize
  iv = 0
  def enek(k):
    """Compute energy in this kpoint"""
    hk = f(k)  # kdependent hamiltonian
    if use_kpm: # Kernel polynomial method
      return kpm.total_energy(hk,scale=10.,ntries=20,npol=100) # using KPM
    else: # conventional diagonalization
      if nbands is None: vv = lg.eigvalsh(hk) # diagonalize k hamiltonian
      else: 
          vv,aa = slg.eigsh(hk,k=4*nbands,which="LM",sigma=0.0) 
          vv = -np.sort(-(vv[vv<0.0])) # negative eigenvalues
          vv = vv[0:nbands] # get the negative eigenvlaues closest to EF
      return np.sum(vv[vv<0.0]) # sum energies below fermi energy
  # compute energy using different modes
  if mode=="mesh":
    from .klist import kmesh
    kp = kmesh(h.dimensionality,nk=nk)
    etot = np.mean(parallel.pcall(enek,kp)) # compute total energy
  elif mode=="random":
    kp = [np.random.random(3) for i in range(nk)] # random points
    etot = np.mean(parallel.pcall(enek,kp)) # compute total eenrgy
  elif mode=="integrate":
    from scipy import integrate
    if h.dimensionality==1: # one dimensional
        etot = integrate.quad(enek,-1.,1.,epsabs=tol,epsrel=tol)[0]
    elif h.dimensionality==2: # two dimensional
        etot = integrate.dblquad(lambda x,y: enek([x,y]),-1.,1.,-1.,1.,
                epsabs=tol,epsrel=tol)[0]
    else: raise
  else: raise
  return etot 

def diagonalize(intra,n=20,e=0.0,mode="arpack"):
  if mode=="arpack":
    eig,eigvec = slg.eigsh(csc_matrix(intra),k=n,which="LM",sigma=e,
                                  tol=arpack_tol,maxiter=arpack_maxiter)
  else:
    eig,eigvec = lg.eigh(csc_matrix(intra).todense())
  return eig,np.transpose(eigvec) # return 

def raw_gap(h,kpgen,sparse=True,nk=100):
  hk_gen = h.get_hk_gen() # get hamiltonian generator
  ks = np.linspace(0.,1.,nk)
  etot = [] # total eigenvalues
  for k in ks:
    kp = kpgen(k)
    hk = hk_gen(kp) # generate hamiltonian
    if sparse: 
      es,ew = lgs.eigsh(csc_matrix(hk),k=4,which="LM",sigma=0.0)
    else:
      es = lg.eigvalsh(hk) # get eigenvalues
    etot.append(es)
  etot = np.array(etot)
  return min(etot[etot>0.]) 

def smalleig(m,numw=10,evecs=False,tol=1e-7):
  """
  Return the smallest eigenvalues using arpack
  """
  m = csc_matrix(m) # sparse matrix
  eig,eigvec = slg.eigsh(m,k=numw,which="LM",sigma=0.0,
                                  tol=tol)
  if evecs:  return eig,eigvec.transpose()  # return eigenvectors
  else:  return eig  # return eigenvalues 

def init_nvecs(xs, ys, sz, rank, with_T=False):
	from scipy.sparse.linalg import eigsh

	T = to_tensor(xs, ys, sz)
	T = [Tk.tocsr() for Tk in T]
	S = sum([T[k] + T[k].T for k in range(len(T))])
	_, E = eigsh(sp.csr_matrix(S), rank)
	if not with_T:
		return E
	else:
		return E, T 

def calculate_scaled_laplacian(adj_mx, lambda_max=2, undirected=True):
    if undirected:
        adj_mx = np.maximum.reduce([adj_mx, adj_mx.T])
    L = calculate_normalized_laplacian(adj_mx)
    if lambda_max is None:
        lambda_max, _ = linalg.eigsh(L, 1, which='LM')
        lambda_max = lambda_max[0]
    L = sp.csr_matrix(L)
    M, _ = L.shape
    I = sp.identity(M, format='csr', dtype=L.dtype)
    L = (2 / lambda_max * L) - I
    return L.astype(np.float32).todense() 

def compute_eigen(self, n_comps=10, sym=None, sort="decrease"):
        if self._transitions_sym is None:
            raise ValueError("Run `.compute_transitions` first.")
        n_comps = min(self._transitions_sym.shape[0] - 1, n_comps)
        evals, evecs = linalg.eigsh(self._transitions_sym, k=n_comps, which="LM")
        self._eigen_values = evals[::-1]
        self._eigen_basis = evecs[:, ::-1] 

def test_eig_cXY_groundstate(nr_sites, gamma, rank, tol, rgen, local_dim=2):
    # Verify that linalg.eig() finds the correct ground state energy
    # of the cyclic XY model
    E0 = physics.cXY_E0(nr_sites, gamma)
    mpo = physics.mpo_cH(physics.cXY_local_terms(nr_sites, gamma))
    eigs = ft.partial(eigsh, k=1, which='SA', tol=1e-6)
    E0_mp, mineig_eigvec_mp = mpnum.linalg.eig(
        mpo, startvec_rank=rank, randstate=rgen, var_sites=2, num_sweeps=3,
        eigs=eigs)
    print(abs(E0_mp - E0))
    assert abs(E0_mp - E0) <= tol 

def test_eig_sum(nr_sites, local_dim, rank, rgen):
    # Need at least three sites for var_sites = 2
    if nr_sites < 3:
        pt.skip("Nothing to test")
        return
    rank = max(1, rank // 2)
    mpo = factory.random_mpo(nr_sites, local_dim, rank, randstate=rgen,
                             hermitian=True, normalized=True)
    mpo.canonicalize()
    mps = factory.random_mpa(nr_sites, local_dim, rank, randstate=rgen,
                             dtype=np.complex_, normalized=True)
    mpas = [mpo, mps]

    vec = mps.to_array().ravel()
    op = mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
    op += np.outer(vec, vec.conj())
    eigvals, eigvec = np.linalg.eigh(op)

    # Eigenvals should be real for a hermitian matrix
    assert (np.abs(eigvals.imag) < 1e-10).all(), str(eigvals.imag)
    mineig_pos = eigvals.argmin()
    mineig, mineig_eigvec = eigvals[mineig_pos], eigvec[:, mineig_pos]
    mineig_mp, mineig_eigvec_mp = mp.eig_sum(
        mpas, num_sweeps=5, startvec_rank=5 * rank, randstate=rgen,
        eigs=ft.partial(eigsh, k=1, which='SA', tol=1e-6),
        var_sites=2)
    mineig_eigvec_mp = mineig_eigvec_mp.to_array().flatten()

    overlap = np.inner(mineig_eigvec.conj(), mineig_eigvec_mp)
    assert_almost_equal(mineig_mp, mineig)
    assert_almost_equal(abs(overlap), 1)