Python linalg.eig方法代码示例

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def compare_solutions(A,B,m):
    n = A.shape[0]

    numpy.random.seed(0)

    V = rand(n,m)
    X = linalg.orth(V)

    eigs,vecs = lobpcg(A, X, B=B, tol=1e-5, maxiter=30)
    eigs.sort()

    #w,v = symeig(A,B)
    w,v = eig(A,b=B)
    w.sort()

    assert_almost_equal(w[:int(m/2)],eigs[:int(m/2)],decimal=2)

    #from pylab import plot, show, legend, xlabel, ylabel
    #plot(arange(0,len(w[:m])),w[:m],'bx',label='Results by symeig')
    #plot(arange(0,len(eigs)),eigs,'r+',label='Results by lobpcg')
    #legend()
    #xlabel(r'Eigenvalue $i$')
    #ylabel(r'$\lambda_i$')
    #show() 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def test_eigs_for_k_greater():
    # Test eigs() for k beyond limits.
    A_sparse = diags([1, -2, 1], [-1, 0, 1], shape=(4, 4))  # sparse
    A = generate_matrix(4, sparse=False)
    M_dense = np.random.random((4, 4))
    M_sparse = generate_matrix(4, sparse=True)
    M_linop = aslinearoperator(M_dense)
    eig_tuple1 = eig(A, b=M_dense)
    eig_tuple2 = eig(A, b=M_sparse)

    with suppress_warnings() as sup:
        sup.filter(RuntimeWarning)

        assert_equal(eigs(A, M=M_dense, k=3), eig_tuple1)
        assert_equal(eigs(A, M=M_dense, k=4), eig_tuple1)
        assert_equal(eigs(A, M=M_dense, k=5), eig_tuple1)
        assert_equal(eigs(A, M=M_sparse, k=5), eig_tuple2)

        # M as LinearOperator
        assert_raises(TypeError, eigs, A, M=M_linop, k=3)

        # Test 'A' for different types
        assert_raises(TypeError, eigs, aslinearoperator(A), k=3)
        assert_raises(TypeError, eigs, A_sparse, k=3) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def dlqr(a, b, q, r):
    """ Get the feedback controls from linearized system at the current time step

    for a discrete time system Ax+Bu
    find the infinite horizon optimal feedback controller
    to steer the system to the origin
    with
    u = -K*x
    """
    x = np.matrix(sLa.solve_discrete_are(a, b, q, r))

    k = np.matrix(sLa.inv(b.T * x * b + r) * (b.T * x * a))

    eigVals, eigVecs = sLa.eig(a - b * k)

    return np.asarray(k), np.asarray(x), eigVals 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def null(A,tol=1e-6):
	ee, ev = la.eig(A)
	
	#for E,V in zip(ee,ev.T):
	#	print 'Eigs:',abs(E), '\t', E#, '\t', V
	#print '\n'
	
	z = list(zip(ee,ev.T))
	zs = sorted(z, key=lambda f: abs(f[0])) # sort by absolute value of eigenvectors
	ees, evs = list(zip(*zs))
	
	#for E,V in zip(ee,ev):
	#	print abs(E), '\t', E, '::', V
	
	if abs(ees[0]<tol):
		return evs[0].T
	else:
		print('No null eigenvector found! List of eigenvalules:')
		for E,V in zip(ee,ev.T):
			print(('Eigs:',abs(E), '\t', E, '\n\t', V))
		print('\n')
		return 0 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def test_sampler(sampler, H, steps=1000):
    """Runs the sampler on the given energy
    Prints a bunch of statistics about how well it's doing
    returns t_obs, distr_obs
    """
    order = len(H) * 2
    smp = sampler(order, energies=H)

    smp.sample(steps)
    t_obs = smp.get_transition_matrix()

    print "Predicted distribution: {}  \n".format(smp.prd_distr)
    print "Observed distribution: {} \n".format(smp.get_distr())
    print "Sampling error (L1): {} \n".format(smp.sampling_err())
    print "Observed transition matrix: \n {} \n".format(t_obs)
    print "Eigenspectrum of observed transition matrix: \n"
    eigs = rectify_evecs(eig(t_obs, left=True, right=False))
    pprint_eigs(eigs)

    return t_obs, smp.get_distr() 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def rectify_evecs(eigs):
    """
    eigs: output of linalg.eig
    normalizes evecs by L1 norm, truncates small complex components,
    ensures things are positive
    """
    evecs = eigs[1].T
    l1_norm = np.abs(evecs).sum(axis=1)
    norm_evecs = evecs / l1_norm[:, np.newaxis]
    real_evals = [np.around(np.real_if_close(l), decimals=5) for l in eigs[0]]
    real_evecs = []

    for v in norm_evecs:
        real_v = np.real_if_close(v)
        if (real_v < 0).all():
            real_v *= -1
        real_evecs.append(real_v)

    # skip sorting for now: argsort is pain because numpy will typecase to complex arr
    #    desc_idx = np.argsort(real_evals)[::-1]
    #   return real_evals[desc_idx], real_evecs[desc_idx]
    return real_evals, real_evecs 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def compute_separatrices(Xss, Js, func, x_range, y_range, t=50, n_sample=500, eps=1e-6):
    ret = []
    for i, x in enumerate(Xss):
        print(x)
        J = Js[i]
        w, v = eig(J)
        I_stable = np.where(np.real(w) < 0)[0]
        print(I_stable)
        for j in I_stable:  # I_unstable
            u = np.real(v[j])
            u = u / np.linalg.norm(u)
            print("u=%f, %f" % (u[0], u[1]))

            # Parameters for building separatrix
            T = np.linspace(0, t, n_sample)
            # all_sep_a, all_sep_b = None, None
            # Build upper right branch of separatrix
            ab_upper = odeint(lambda x, _: -func(x), x + eps * u, T)
            # Build lower left branch of separatrix
            ab_lower = odeint(lambda x, _: -func(x), x - eps * u, T)

            sep = np.vstack((ab_lower[::-1], ab_upper))
            ret.append(sep)
    ret = clip_curves(ret, [x_range, y_range])
    return ret 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def pca_projection(C, L):
    """solve the problem size(C) = NxN, size(W) = NxL. max_W trace( W' C W ) : W' W = I

    Arguments
    ---------
    C: (ndarrya) The matrix of
    L: (int) The number of Eigenvalues
    Return
    ------
    W: The L largest Eigenvalues
    """

    V, U = eig(C)
    eig_idx = np.argsort(V).tolist()
    eig_idx.reverse()
    W = U.T[eig_idx[0:L]].T
    return W 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def dmd_exact(x_data, y_data):
    # decompose x
    u, s, v = lin.svd(x_data, full_matrices=False, overwrite_a=True, check_finite=False, lapack_driver='gesvd')

    # construct reduced matrix
    reduced_matrix = u.T @ y_data @ v.T @ np.diag(np.reciprocal(s))

    # find eigenvalues
    eigenvalues, eigenvectors = lin.eig(reduced_matrix, overwrite_a=True, check_finite=False)

    # sort eigenvalues
    ind = np.argsort(eigenvalues)[::-1]
    dmd_eigenvalues = eigenvalues[ind]

    # compute modes
    dmd_modes = y_data @ v.T @ np.diag(np.reciprocal(s)) @ eigenvectors[:, ind] @ np.diag(
        np.reciprocal(dmd_eigenvalues))

    return dmd_eigenvalues, dmd_modes 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def _method_2(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]
    tmp = np.dot(data, data.T)
    corr_offset = np.zeros(tmp.shape)
    corr_offset[1:] = tmp[:-1]
    corr_offset[:-1] += tmp[1:]
    if num_pcs is None:
        eivals, eivects = eig(corr_offset)
    else:
        eivals, eivects = eigs(corr_offset, num_pcs, which='LR')
    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]
    transformed_eivects = np.dot(data.T, eivects)
    for i in range(transformed_eivects.shape[1]):  # normalize the eigenvectors
        transformed_eivects[:, i] /= np.linalg.norm(transformed_eivects[:, i])
    return eivals, transformed_eivects, np.dot(data, transformed_eivects) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def gevd(x1,x2,no_pairs):
	'''Solve generalized eigenvalue decomposition
	
	Keyword arguments:
	x1 -- numpy array of size [NO_channels, NO_samples]
	x2 -- numpy array of size [NO_channels, NO_samples]
	no_pairs -- number of pairs of eigenvectors to be returned 

	Return:	numpy array of 2*No_pairs eigenvectors 
	'''
	ev,vr= linalg.eig(x1,x2,right=True) 
	evAbs = np.abs(ev)
	sort_indices = np.argsort(evAbs)
	chosen_indices = np.zeros(2*no_pairs).astype(int)
	chosen_indices[0:no_pairs] = sort_indices[0:no_pairs]
	chosen_indices[no_pairs:2*no_pairs] = sort_indices[-no_pairs:]
	
	w = vr[:,chosen_indices] # ignore nan entries 
	return w 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def symeig(mtxA, mtxB=None, eigenvectors=True, select=None):
    import scipy.linalg as sla
    if select is None:
        if np.iscomplexobj(mtxA):
            if mtxB is None:
                fun = sla.get_lapack_funcs('heev', arrays=(mtxA,))
            else:
                fun = sla.get_lapack_funcs('hegv', arrays=(mtxA,))
        else:
            if mtxB is None:
                fun = sla.get_lapack_funcs('syev', arrays=(mtxA,))
            else:
                fun = sla.get_lapack_funcs('sygv', arrays=(mtxA,))
##         print fun
        if mtxB is None:
            out = fun(mtxA)
        else:
            out = fun(mtxA, mtxB)
            out = out[1], out[0], out[2]
##         print w
##         print v
##         print info
##         from symeig import symeig
##         print symeig( mtxA, mtxB )
    else:
        out = sla.eig(mtxA, mtxB, right=eigenvectors)
        w = out[0]
        ii = np.argsort(w)
        w = w[slice(*select)]
        if eigenvectors:
            v = out[1][:,ii]
            v = v[:,slice(*select)]
            out = w, v, 0
        else:
            out = w, 0

    return out[:-1] 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def getFiedlerVectorDense(self):
        """
         Returns the (dense)Fiedler vector of the higher-order network. The Fiedler 
         vector can be used for a spectral bisectioning of the network.             
        """
    
        # NOTE: The Laplacian is transposed for the sparse case to get the left
        # NOTE: eigenvalue.
        L = self.getLaplacianMatrix()
        # convert to dense matrix and transpose again to have the untransposed
        # laplacian again.
        w, v = _la.eig(L.todense().transpose(), right=False, left=True)

        return v[:,_np.argsort(_np.absolute(w))][:,1] 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def compare_solutions(A,B,m):
    n = A.shape[0]

    np.random.seed(0)

    V = rand(n,m)
    X = linalg.orth(V)

    eigs,vecs = lobpcg(A, X, B=B, tol=1e-5, maxiter=30)
    eigs.sort()

    w,v = eig(A,b=B)
    w.sort()

    assert_almost_equal(w[:int(m/2)],eigs[:int(m/2)],decimal=2) 

# 需要导入模块: from scipy import linalg [as 别名]
# 或者: from scipy.linalg import eig [as 别名]
def pagerank_weighted_scipy(graph, damping=0.85):
    adjacency_matrix = build_adjacency_matrix(graph)
    probability_matrix = build_probability_matrix(graph)

    # Suppress deprecation warnings from numpy.
    # See https://github.com/summanlp/textrank/issues/57
    import warnings
    with warnings.catch_warnings():
        from numpy import VisibleDeprecationWarning
        warnings.filterwarnings("ignore", category=VisibleDeprecationWarning)
        warnings.filterwarnings("ignore", category=PendingDeprecationWarning)
        pagerank_matrix = damping * adjacency_matrix.todense() + (1 - damping) * probability_matrix

    vals, vecs = eig(pagerank_matrix, left=True, right=False)
    return process_results(graph, vecs)