python类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 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)