Python scipy.sparse.linalg 模块，eigsh() 实例源码

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 eigen_mode(model:Model.fem_model,n):
    """
    Solve the eigen mode of the MDOF system\n
    n: number of modes to extract.
    Optional:\n
    sparse: True if system use a sparse K\n
    return: w, freq, t, mode
    """
    w2,mode=sl.eigsh(A=model.K,M=model.M,k=n,which='LM')
    freq=[]
    w=[]
    T=[]

    for i in range(n):
        w.append(np.sqrt(w2[i]))
        T.append(2*np.pi/w[-1])
        freq.append(1/T[-1])
    return w,freq,T,mode

def nvecs(X, n, rank, do_flipsign=True, dtype=np.float):
    """
    Eigendecomposition of mode-n unfolding of a tensor
    """
    Xn = X.unfold(n)
    if issparse_mat(Xn):
        Xn = csr_matrix(Xn, dtype=dtype)
        Y = Xn.dot(Xn.T)
        _, U = eigsh(Y, rank, which='LM')
    else:
        Y = Xn.dot(Xn.T)
        N = Y.shape[0]
        _, U = eigh(Y, eigvals=(N - rank, N - 1))
        #_, U = eigsh(Y, rank, which='LM')
    # reverse order of eigenvectors such that eigenvalues are decreasing
    U = array(U[:, ::-1])
    # flip sign
    if do_flipsign:
        U = flipsign(U)
    return U

def test_eig(nr_sites, local_dim, rank, which, var_sites, rgen, request):
    if nr_sites <= var_sites:
        pt.skip("Nothing to test")
        return  # No local optimization can be defined
    if not (_pytest_want_long(request) or
            (nr_sites, local_dim, rank, var_sites, which) in {
                (3, 2, 4, 1, 'SA'), (4, 3, 5, 1, 'LM'), (5, 2, 1, 2, 'LA'),
                (6, 2, 4, 2, 'SA'),
            }):
        pt.skip("Should only be run in long tests")
    # With startvec_rank = 2 * rank and this seed, eig() gets
    # stuck in a local minimum. With startvec_rank = 3 * rank,
    # it does not.
    mpo = factory.random_mpo(nr_sites, local_dim, rank, randstate=rgen,
                             hermitian=True, normalized=True)
    mpo.canonicalize()
    op = mpo.to_array_global().reshape((local_dim**nr_sites,) * 2)
    v0 = factory._zrandn([local_dim**nr_sites], rgen)
    eigval, eigvec = eigsh(op, k=1, which=which, v0=v0)
    eigval, eigvec = eigval[0], eigvec[:, 0]

    eig_rank = (4 - var_sites) * rank
    eigval_mp, eigvec_mp = mp.eig(
        mpo, num_sweeps=5, var_sites=1, startvec_rank=eig_rank, randstate=rgen,
        eigs=ft.partial(eigsh, k=1, which=which, tol=1e-6, maxiter=250))
    eigvec_mp = eigvec_mp.to_array().flatten()

    overlap = np.inner(eigvec.conj(), eigvec_mp)
    assert_almost_equal(eigval, eigval_mp, decimal=14)
    assert_almost_equal(1, abs(overlap), decimal=14)

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)

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 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 estimation(self, y, ds):
        """ Estimate KCCA.

        Parameters
        ----------
        y : (number of samples, dimension)-ndarray
             One row of y corresponds to one sample.
        ds : int vector
             Dimensions of the individual subspaces in y; ds[i] = i^th
             subspace dimension.

        Returns
        -------
        i : float
            Estimated value of KCCA.

        References
        ----------
        Francis Bach, Michael I. Jordan. Learning graphical models with
        Mercer kernels. International Conference on Neural Information
        Processing Systems (NIPS), pages 1033-1040, 2002.

        Examples
        --------
        i = co.estimation(y,ds)

        """

        # verification:
        self.verification_compatible_subspace_dimensions(y, ds)

        num_of_samples = y.shape[0]
        tol = num_of_samples * self.eta

        r = compute_matrix_r_kcca_kgv(y, ds, self.kernel, tol, self.kappa)
        eig_min = eigsh(r, k=1, which='SM')[0][0]
        i = -log(eig_min) / 2

        return i

def LoadKM():
    """ Loads m and k matrices from a full file """
    import numpy as np

    # Create file reader object
    fobj = pyansys.FullReader(fullfile)
    dofref, k, m = fobj.LoadKM(utri=False)

    # print results
    ndim = k.shape[0]
    print('Loaded {:d} x {:d} mass and stiffness matrices'.format(ndim, ndim))
    print('\t k has {:d} entries'.format(k.indices.size))
    print('\t m has {:d} entries'.format(m.indices.size))

    # compute natural frequencies if installed
    try:
        from scipy.sparse import linalg
    except BaseException:
        return

    # removing constrained nodes (this is not efficient)
    free = np.logical_not(fobj.const).nonzero()[0]
    k = k[free][:, free]
    m = m[free][:, free]

    import numpy as np
    # Solve
    w, v = linalg.eigsh(k, k=20, M=m, sigma=10000)

    # System natural frequencies
    f = np.real(w)**0.5 / (2 * np.pi)

    print('First four natural frequencies:')
    for i in range(4):
        print('{:.3f} Hz'.format(f[i]))

def eigsh(A,max_try=6,return_eigenvectors=True,**karg):
    '''
    Find the eigenvalues and eigenvectors of the real symmetric square matrix or complex hermitian matrix A.
    This is a wrapper for scipy.sparse.linalg.eigsh to handle the exceptions it raises.

    Parameters
    ----------
    A : An NxN matrix, array, sparse matrix, or LinearOperator
        The matrix whose eigenvalues and eigenvectors is to be computed.
    max_try : integer, optional
        The maximum number of tries to do the computation when the computed eigenvalues/eigenvectors do not converge.
    return_eigenvectors : logical, optional
        True for returning the eigenvectors and False for not.
    karg : dict
        Please refer to https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.eigsh.html for details.
    '''
    if A.shape==(1,1):
        assert 'M' not in karg and 'sigma' not in karg and karg.get('k',1)==1
        if return_eigenvectors:
            return A.dot(np.ones(1)).reshape(-1),np.ones((1,1),dtype=A.dtype)
        else:
            return A.dot(np.ones(1)).reshape(-1)
    else:
        ntry=1
        while True:
            try:
                result=pl.eigsh(A,return_eigenvectors=return_eigenvectors,**karg)
                break
            except pl.ArpackNoConvergence as err:
                if ntry<max_try:
                    ntry+=1
                else:
                    raise err
        return result

def test_lanczos():
    print 'test_lanczos'
    N,Nv,Ne,Niter=4000,15,1,750

    np.random.seed(1)
    m=np.random.random((N,N))+1j*np.random.random((N,N))
    m=m+m.T.conjugate()
    v0=np.random.random((Nv,N))
    exacteigs=sl.eigh(m,eigvals_only=True)[:Ne] if N<1000 else eigsh(m,k=Ne,return_eigenvectors=False,which='SA')[::-1]
    print 'Exact eigs:',exacteigs
    print

    a=Lanczos(m,deepcopy(v0),maxiter=Niter,keepstate=True)
    for i in xrange(Niter): a.iter()
    print 'diff from Hermitian: %s.'%sl.norm(a.T-a.T.T.conjugate())
    h=np.zeros((Niter,Niter),dtype=m.dtype)
    for i in xrange(Niter):
        hi=a.matrix.dot(a.vectors[i]).conjugate()
        for j in xrange(Niter):
            h[i,j]=hi.dot(a.vectors[j])
    print 'diff from h: %s.'%sl.norm(a.T-h)
    V=np.asarray(a.vectors)
    print 'vecotrs input diff: %s.'%sl.norm(V[0:Nv,:].T.dot(a.P)-v0.T)
    print 'vectors orthonormal diff: %s.'%sl.norm(np.identity(Niter)-V.dot(V.T.conjugate()))
    Leigs=a.eigs()[:Ne]
    print 'Lanczos eigs:',Leigs
    print 'eigenvalues diff: %s.'%sl.norm(exacteigs-Leigs)
    print

def solve(_m2, _m3, _n,  _k):

    vals, vecs = la.eigsh(_m2, _k)

    vals = vals + 0j
    W = np.array(vecs, dtype=np.complex128)
    for i in range(_k):
        for j in range(_n):
            W[j, i] /= np.sqrt(vals[i])

    T = np.dot(W.T, _m3.dot(W))

    vals, vecs = np.linalg.eig(T)

    v = np.dot(W, np.linalg.solve(np.dot(W.T, W), vecs))

    s = []
    for i in range(_k):
        s.append(sum(v[:, i]))

    for i in range(_k):
        for j in range(_n):
            v[j, i] /= s[i]

    a = []
    for i in range(_k):
        v_ = np.dot(W.T, v[:, i])
        a.append(1. / np.dot(v_.T, v_))

    a = np.array(a)

    r = vals

    return a, r, v

def solve_modal(model:Model.fem_model,k:int):
    """
    model: FEM model.
    k: number of modes to extract.
    """
    Logger.info('Solving eigen modes...')
    K_bar,M_bar,F_bar,index=model.eliminate_matrix(True)
    Dvec=model.D
    n_nodes=model.node_count
    deltas=[]
    try:
        #general eigen solution, should be optimized later!!
        A=sp.csr_matrix(K_bar)
        B=sp.csr_matrix(M_bar)

        omega2s,modes = sl.eigsh(A,k,B,which='SM')
        periods=[]

        for omega2 in omega2s:
            periods.append(2*np.pi/np.sqrt(omega2))

        #extract vibration mode
        for mode in modes.T:
            delta_bar = mode/sum(mode)
            delta=delta_bar
            #fill original displacement vector
            prev = 0
            for idx in index:
                gap=idx-prev
                if gap>0:
                    delta=np.insert(delta,prev,[0]*gap)
                prev = idx + 1               
                if idx==index[-1] and idx!=n_nodes-1:
                    delta = np.insert(delta,prev, [0]*(n_nodes*6-prev))
            delta += Dvec
            deltas.append(delta)

    except Exception as e:
        print(str(e))
        return None
    model.is_solved=True
    return periods, deltas

def SolveKM():
    """
    Loads and solves a mass and stiffness matrix from an ansys full file
    """
    import numpy as np
    from scipy.sparse import linalg
    import pyansys

    # load the mass and stiffness matrices
    full = pyansys.FullReader(pyansys.examples.fullfile)
    dofref, k, m = full.LoadKM(utri=False)

    # removing constrained nodes (this is not efficient)
    free = np.logical_not(full.const).nonzero()[0]
    k = k[free][:, free]
    m = m[free][:, free]

    # Solve
    w, v = linalg.eigsh(k, k=20, M=m, sigma=10000)

    # System natural frequencies
    f = (np.real(w))**0.5 / (2 * np.pi)

    #%% Plot result
    import vtkInterface

    # Get the 4th mode shape
    mode_shape = v[:, 3]  # x, y, z displacement for each node

    # create the full mode shape including the constrained nodes
    full_mode_shape = np.zeros(dofref.shape[0])
    full_mode_shape[np.logical_not(full.const)] = mode_shape

    # reshape and compute the normalized displacement
    disp = full_mode_shape.reshape((-1, 3))
    n = (disp * disp).sum(1)**0.5
    n /= n.max()  # normalize

    # load an archive file and create a vtk unstructured grid
    archive = pyansys.ReadArchive(pyansys.examples.hexarchivefile)
    grid = archive.ParseVTK()

    # plot the normalized displacement
#    grid.Plot(scalars=n)

    # Fancy plot the displacement
    plobj = vtkInterface.PlotClass()

    # add two meshes to the plotting class.  Meshes are copied on load
    plobj.AddMesh(grid, style='wireframe')
    plobj.AddMesh(grid, scalars=n, stitle='Normalized\nDisplacement',
                  flipscalars=True)

    # Update the coordinates by adding the mode shape to the grid
    plobj.UpdateCoordinates(grid.GetNumpyPoints() + disp / 80, render=False)
    plobj.AddText('Cantliver Beam 4th Mode Shape at {:.4f}'.format(f[3]),
                  fontsize=30)
    plobj.Plot()
    del plobj