我们从Python开源项目中，提取了以下13个代码示例，用于说明如何使用scipy.linalg.qr()。
def clean_basis(basis,traj_edges,delay,orthogonalize=True): # Normalize the eigenvectors N,k = np.shape(basis) t_0_indices, t_lag_indices = start_stop_indices(traj_edges,delay) evec_norm = np.linalg.norm(basis,axis=0) basis *= np.sqrt(N)/evec_norm # Calculate orthogonal coefficients if orthogonalize: basis_t_0 = basis[t_0_indices] Q,R = spl.qr(basis_t_0) R_sub = R[:k,:k] basis = np.dot(basis,R_sub) return basis
def _unitary_haar(dim, randstate=None): """Returns a sample from the Haar measure of the unitary group of given dimension. :param int dim: Dimension :param randn: Function to create real N(0,1) distributed random variables. It should take the shape of the output as numpy.random.randn does (default: numpy.random.randn) """ z = _zrandn((dim, dim), randstate) / np.sqrt(2.0) q, r = qr(z) d = np.diagonal(r) ph = d / np.abs(d) return q * ph
def __init__(self, n): self.n = n self.m = n/5 M = np.random.rand(self.m, self.m) M = linalg.qr(M)[0] A = linalg.block_diag(M,M,M,M,M) # A = linalg.block_diag(M,M,M,M,M,M,M,M,M,M, M,M,M,M,M,M,M,M,M,M) self.A = A
def __init__(self): M = np.random.rand(5,5) M = linalg.qr(M)[0] A = linalg.block_diag(M,M,M,M) A = linalg.block_diag(M,M,M,M,M,M,M,M,M,M, M,M,M,M,M,M,M,M,M,M) self.A = A
def rq(A): from scipy.linalg import qr Q,R = qr(flipud(A).T) R = flipud(R.T) Q = Q.T return R[:,::-1],Q[::-1,:]
def principal_angles(X, Y, n): ''' compute the principal angles between two subspaces return: np.array of principal angles, orthogonal matrics U and V ''' QX, RX = la.qr(X) QY, RY = la.qr(Y) if X.shape[1] >= Y.shape[1]: C = QX.conjugate().T.dot(QY) M, cos, Nh = la.svd(C) U = QX.dot(M) V = QY.dot(Nh.conjugate().T) angles = np.arccos(cos) return angles, U, V else: C = QY.conjugate().T.dot(QX) M, cos, Nh = la.svd(C) U = QX.dot(M) V = QY.dot(Nh.conjugate().T) angles = np.arccos(cos) return angles, U, V # Similarity between subspaces by Yamaguchi's definition
def _qr_economic_old(A, **kwargs): """ Compat function for the QR-decomposition in economic mode Scipy 0.9 changed the keyword econ=True to mode='economic' """ with warnings.catch_warnings(record=True): return linalg.qr(A, econ=True, **kwargs)
def _qr_economic_new(A, **kwargs): return linalg.qr(A, mode='economic', **kwargs)
def _overlap_projector(data_int, data_res, corr): """Calculate projector for removal of subspace intersection in tSSS""" # corr necessary to deal with noise when finding identical signal # directions in the subspace. See the end of the Results section in [2]_ # Note that the procedure here is an updated version of [2]_ (and used in # Elekta's tSSS) that uses residuals instead of internal/external spaces # directly. This provides more degrees of freedom when analyzing for # intersections between internal and external spaces. # Normalize data, then compute orth to get temporal bases. Matrices # must have shape (n_samps x effective_rank) when passed into svd # computation # we use np.linalg.norm instead of sp.linalg.norm here: ~2x faster! n = np.linalg.norm(data_int) Q_int = linalg.qr(_orth_overwrite((data_int / n).T), overwrite_a=True, mode='economic', **check_disable)[0].T n = np.linalg.norm(data_res) Q_res = linalg.qr(_orth_overwrite((data_res / n).T), overwrite_a=True, mode='economic', **check_disable)[0] assert data_int.shape[1] > 0 C_mat = np.dot(Q_int, Q_res) del Q_int # Compute angles between subspace and which bases to keep S_intersect, Vh_intersect = linalg.svd(C_mat, overwrite_a=True, full_matrices=False, **check_disable)[1:] del C_mat intersect_mask = (S_intersect >= corr) del S_intersect # Compute projection operator as (I-LL_T) Eq. 12 in [2]_ # V_principal should be shape (n_time_pts x n_retained_inds) Vh_intersect = Vh_intersect[intersect_mask].T V_principal = np.dot(Q_res, Vh_intersect) return V_principal
def groupInverse(M): """ Computes the group inverse of stochastic matrix using the algorithm given by Golub and Meyer in: G. H. Golub and C. D. Meyer, Jr, SIAM J. Alg. Disc. Meth. 7, 273- 281 (1986) Parameters ---------- M : ndarray A square matrix with index 1. Returns ------- grpInvM : ndarray The group inverse of M. """ L=np.shape(M)[1] q,r=qr(M) piDist=q[:,L-1] piDist=(1/np.sum(piDist))*piDist specProjector=np.identity(L)-np.outer(np.ones(L),piDist) u=r[0:(L-1),0:(L-1)]#remember 0:(L-1) actually means 0 to L-2! uInv= inv(u)#REPLACE W. lapack, invert triangular matrix ROUTINE uInv = np.real(uInv) grpInvM=np.zeros((L,L)) grpInvM[0:(L-1),0:(L-1)]=uInv grpInvM=np.dot(specProjector,np.dot(grpInvM,np.dot(q.transpose(),specProjector))) return grpInvM
def solve(A,b,rtol=10**-8): ''' Solve the matrix equation A*x=b by QR decomposition. Parameters ---------- A : 2d ndarray The coefficient matrix. b : 1d ndarray The ordinate values. rtol : np.float64 The relative tolerance of the solution. Returns ------- 1d ndarray The solution. Raises ------ LinAlgError When no solution exists. ''' assert A.ndim==2 nrow,ncol=A.shape if nrow>=ncol: result=np.zeros(ncol,dtype=np.find_common_type([],[A.dtype,b.dtype])) q,r=sl.qr(A,mode='economic',check_finite=False) temp=q.T.dot(b) for i,ri in enumerate(r[::-1]): result[-1-i]=(temp[-1-i]-ri[ncol-i:].dot(result[ncol-i:]))/ri[-1-i] else: temp=np.zeros(nrow,dtype=np.find_common_type([],[A.dtype,b.dtype])) q,r=sl.qr(dagger(A),mode='economic',check_finite=False) for i,ri in enumerate(dagger(r)): temp[i]=(b[i]-ri[:i].dot(temp[:i]))/ri[i] result=q.dot(temp) if not np.allclose(A.dot(result),b,rtol=rtol): raise sl.LinAlgError('solve error: no solution.') return result
def stationary_distrib(F,residtol = 1.E-10,max_iter=100): """ Calculates the eigenvector of the matrix F with eigenvalue 1 (if it exists). Parameters ---------- F : ndarray A matrix known to have a single left eigenvector with eigenvalue 1. residtol : float or scalar To improve the accuracy of the computation, the algorithm will "polish" the final result using several iterations of the power method, z^T F = z^T. Residtol gives the tolerance for the associated relative residual to determine convergence. maxiter : int Maximum number of iterations to use the power method to reduce the residual. In practice, should never be reached. Returns ------- z : ndarray The eigenvector of the matrix F with eigenvalue 1. For a Markov chain stationary distribution, this is the stationary distribution. Normalization is chosen s.t. entries sum to one. """ L = len(F) # Number of States M = np.eye(L)-F q,r=qr(M) z=q[:,-1] # Stationary dist. is last column of QR fact z/=np.sum(z) # Normalize Trajectory # Polish solution using power method. for itr in xrange(max_iter): znew = np.dot(z,F) maxresid = np.max(np.abs(znew - z)/z) # Convergence Criterion if maxresid < residtol: break else: z = znew return z/np.sum(z) # Return normalized (by convention)
