Python scipy.sparse.linalg.eigsh() Examples

The following are 30 code examples of scipy.sparse.linalg.eigsh(). You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may also want to check out all available functions/classes of the module scipy.sparse.linalg , or try the search function .
Example #1
Source File: dynamic.py    From StructEngPy with MIT License 6 votes vote down vote up
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)) 
Example #2
Source File: test_integration.py    From SolidsPy with MIT License 6 votes vote down vote up
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) 
Example #3
Source File: cmm.py    From cmm with GNU General Public License v2.0 6 votes vote down vote up
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 
Example #4
Source File: oPCA.py    From sima with GNU General Public License v2.0 6 votes vote down vote up
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) 
Example #5
Source File: pca.py    From TextDetector with GNU General Public License v3.0 6 votes vote down vote up
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] 
Example #6
Source File: ldos.py    From quantum-honeycomp with GNU General Public License v3.0 6 votes vote down vote up
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 
Example #7
Source File: ldos.py    From quantum-honeycomp with GNU General Public License v3.0 6 votes vote down vote up
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 
Example #8
Source File: utils_physics_test.py    From mpnum with BSD 3-Clause "New" or "Revised" License 6 votes vote down vote up
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 
Example #9
Source File: sf.py    From karateclub with GNU General Public License v3.0 6 votes vote down vote up
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 
Example #10
Source File: topology.py    From quantum-honeycomp with GNU General Public License v3.0 6 votes vote down vote up
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) 
Example #11
Source File: arpack.py    From twitter-stock-recommendation with MIT License 5 votes vote down vote up
def eigsh(A, *args, **kwargs):
    return _eigsh(A, *args, **kwargs) 
Example #12
Source File: diffussion.py    From manifold-diffusion with MIT License 5 votes vote down vote up
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 
Example #13
Source File: bicluster.py    From twitter-stock-recommendation with MIT License 5 votes vote down vote up
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 
Example #14
Source File: test_integration.py    From SolidsPy with MIT License 5 votes vote down vote up
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) 
Example #15
Source File: algebraicconnectivity.py    From qgisSpaceSyntaxToolkit with GNU General Public License v3.0 5 votes vote down vote up
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 
Example #16
Source File: lap.py    From BioNEV with MIT License 5 votes vote down vote up
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 
Example #17
Source File: normcut.py    From sima with GNU General Public License v2.0 5 votes vote down vote up
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 
Example #18
Source File: bandstructure.py    From quantum-honeycomp with GNU General Public License v3.0 5 votes vote down vote up
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 
Example #19
Source File: algebraicconnectivity.py    From aws-kube-codesuite with Apache License 2.0 5 votes vote down vote up
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 
Example #20
Source File: algebraicconnectivity.py    From Carnets with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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 
Example #21
Source File: spectrum.py    From quantum-honeycomp with GNU General Public License v3.0 5 votes vote down vote up
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 
Example #22
Source File: spectrum.py    From quantum-honeycomp with GNU General Public License v3.0 5 votes vote down vote up
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 
Example #23
Source File: density.py    From quantum-honeycomp with GNU General Public License v3.0 5 votes vote down vote up
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 
Example #24
Source File: gap.py    From quantum-honeycomp with GNU General Public License v3.0 5 votes vote down vote up
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.]) 
Example #25
Source File: algebra.py    From quantum-honeycomp with GNU General Public License v3.0 5 votes vote down vote up
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 
Example #26
Source File: util.py    From cesi with Apache License 2.0 5 votes vote down vote up
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 
Example #27
Source File: util.py    From Graph-WaveNet with MIT License 5 votes vote down vote up
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() 
Example #28
Source File: velocity_pseudotime.py    From scvelo with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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] 
Example #29
Source File: linalg_test.py    From mpnum with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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 
Example #30
Source File: linalg_test.py    From mpnum with BSD 3-Clause "New" or "Revised" License 5 votes vote down vote up
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)