Python Util.indSvd方法代码示例

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def addCols2(U, s, V, B):
        """
        Find the SVD of a matrix [A, B] where  A = U diag(s) V.T. Uses the SVD 
        decomposition to find an orthogonal basis on B. 
        
        :param U: The left singular vectors of A  
        
        :param s: The singular values of A 
        
        :param V: The right singular vectors of A 
        
        :param B: The matrix to append to A 
        
        """
        if U.shape[0] != B.shape[0]:
            raise ValueError("U must have same number of rows as B")
        if s.shape[0] != U.shape[1]:
            raise ValueError("Number of cols of U must be the same size as s")
        if s.shape[0] != V.shape[1]:
            raise ValueError("Number of cols of V must be the same size as s")

        m, k = U.shape
        r = B.shape[1]
        n = V.shape[0]

        C = numpy.dot(numpy.eye(m) - numpy.dot(U, U.T), B)
        Ubar, sBar, Vbar = numpy.linalg.svd(C, full_matrices=False)
        inds = numpy.flipud(numpy.argsort(sBar))[0:k]
        Ubar, sBar, Vbar = Util.indSvd(Ubar, sBar, Vbar, inds)

        rPrime = Ubar.shape[1]

        D = numpy.r_[numpy.diag(s), numpy.zeros((rPrime, k))]
        E = numpy.r_[numpy.dot(U.T, B), numpy.diag(sBar).dot(Vbar.T)]
        D = numpy.c_[D, E]

        Uhat, sHat, Vhat = numpy.linalg.svd(D, full_matrices=False)
        inds = numpy.flipud(numpy.argsort(sHat))[0:k]
        Uhat, sHat, Vhat = Util.indSvd(Uhat, sHat, Vhat, inds)

        #The best rank k approximation of [A, B]
        Utilde = numpy.dot(numpy.c_[U, Ubar], Uhat)
        sTilde = sHat

        G1 = numpy.r_[V, numpy.zeros((r, k))]
        G2 = numpy.r_[numpy.zeros((n ,r)), numpy.eye(r)]
        Vtilde = numpy.dot(numpy.c_[G1, G2], Vhat)

        return Utilde, sTilde, Vtilde

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def testSvdSoft(self): 
        A = scipy.sparse.rand(10, 10, 0.2)
        A = A.tocsc()
        
        lmbda = 0.2
        U, s, V = SparseUtils.svdSoft(A, lmbda)
        ATilde = U.dot(numpy.diag(s)).dot(V.T)     
        
        #Now compute the same matrix using numpy
        A = A.todense() 
        
        U2, s2, V2 = numpy.linalg.svd(A)
        inds = numpy.flipud(numpy.argsort(s2))
        inds = inds[s2[inds] > lmbda]
        U2, s2, V2 = Util.indSvd(U2, s2, V2, inds) 
        
        s2 = s2 - lmbda 
        s2 = numpy.clip(s, 0, numpy.max(s2)) 

        ATilde2 = U2.dot(numpy.diag(s2)).dot(V2.T)
        
        nptst.assert_array_almost_equal(s, s)
        nptst.assert_array_almost_equal(ATilde, ATilde2)
        
        #Now run svdSoft with a numpy array 
        U3, s3, V3 = SparseUtils.svdSoft(A, lmbda)
        ATilde3 = U.dot(numpy.diag(s)).dot(V.T)  
        
        nptst.assert_array_almost_equal(s, s3)
        nptst.assert_array_almost_equal(ATilde3, ATilde2)

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def _addSparseRSVD(U, s, V, X, k=10, kX=None, kRand=None, q=None):
        """
        Perform a randomised SVD of the matrix X + U diag(s) V.T. We use th
        """
        if kX==None:
            kX=k
        if kRand==None:
            kRand=k
        if q==None:
            q=1

        m, n = X.shape
        Us = U*s

        kX = numpy.min([m, n, kX])
        UX, sX, VX = SparseUtils.svdPropack(X, kX)
        omega = numpy.c_[V, VX, numpy.random.randn(n, kRand)]
        
        def rMultA(x):
            return Us.dot(V.T.dot(x)) + X.dot(x)
        def rMultAT(x):
            return V.dot(Us.T.dot(x)) + X.T.dot(x)
        
        Y = rMultA(omega)
        for i in range(q): 
            Y = rMultAT(Y)
            Y = rMultA(Y)
        
        Q, R = numpy.linalg.qr(Y)
        B = rMultAT(Q).T   
        U, s, VT = numpy.linalg.svd(B, full_matrices=False)
        U, s, V = Util.indSvd(U, s, VT, numpy.flipud(numpy.argsort(s))[:k])
        U = Q.dot(U)
        
        return U, s, V 

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def svdSoft(X, lmbda, kmax=None):
        """
        Find the partial SVD of the sparse or dense matrix X, for which singular
        values are >= lmbda. Soft threshold the resulting singular values
        so that s <- max(s - lambda, 0)
        """
        if scipy.sparse.issparse(X):
            k = min(X.shape[0], X.shape[1])
            L = scipy.sparse.linalg.aslinearoperator(X)

            U, s, V = SparseUtils.svdPropack(L, k, kmax=kmax)
            V = V.T
        else:
            U, s, V = numpy.linalg.svd(X)

        inds = numpy.flipud(numpy.argsort(s))
        inds = inds[s[inds] >= lmbda]
        U, s, V = Util.indSvd(U, s, V, inds)

        #Soft threshold
        if s.shape[0] != 0:
            s = s - lmbda
            s = numpy.clip(s, 0, numpy.max(s))

        return U, s, V

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
 def testSvd(self): 
     n = 100 
     m = 80
     A = scipy.sparse.rand(m, n, 0.1)
     
     ks = [10, 20, 30, 40] 
     q = 2 
     
     lastError = numpy.linalg.norm(A.todense())        
     
     for k in ks: 
         U, s, V = RandomisedSVD.svd(A, k, q)
         
         nptst.assert_array_almost_equal(U.T.dot(U), numpy.eye(k))
         nptst.assert_array_almost_equal(V.T.dot(V), numpy.eye(k))
         A2 = (U*s).dot(V.T)
         
         error = numpy.linalg.norm(A - A2)
         self.assertTrue(error <= lastError)
         lastError = error 
         
         #Compare versus exact svd 
         U, s, V = numpy.linalg.svd(numpy.array(A.todense()))
         inds = numpy.flipud(numpy.argsort(s))[0:k*2]
         U, s, V = Util.indSvd(U, s, V, inds)
         
         Ak = (U*s).dot(V.T)
         
         error2 = numpy.linalg.norm(A - Ak)
         self.assertTrue(error2 <= error)

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def eigenAdd(omega, Q, Y, k):
        """
        Perform an eigen update of the form A*A + Y*Y in which Y is a low-rank matrix
        and A^*A = Q Omega Q*. We use the rank-k approximation of A:  Q_k Omega_k Q_k^*
        and then approximate [A^*A_k Y^*Y]_k.
        """
        #logging.debug("< eigenAdd >")
        Parameter.checkInt(k, 0, omega.shape[0])
        #if not numpy.isrealobj(omega) or not numpy.isrealobj(Q):
        #    raise ValueError("Eigenvalues and eigenvectors must be real")
        if omega.ndim != 1:
            raise ValueError("omega must be 1-d array")
        if omega.shape[0] != Q.shape[1]:
            raise ValueError("Must have same number of eigenvalues and eigenvectors")

        if __debug__:
            Parameter.checkOrthogonal(Q, tol=EigenUpdater.tol, softCheck=True, arrayInfo="input Q in eigenAdd()")

        #Taking the abs of the eigenvalues is correct
        inds = numpy.flipud(numpy.argsort(numpy.abs(omega)))

        omega, Q = Util.indEig(omega, Q, inds[numpy.abs(omega)>EigenUpdater.tol])
        Omega = numpy.diag(omega)

        YY = Y.conj().T.dot(Y)
        QQ = Q.dot(Q.conj().T)
        Ybar = Y - Y.dot(QQ)

        Pbar, sigmaBar, Qbar = numpy.linalg.svd(Ybar, full_matrices=False)
        inds = numpy.flipud(numpy.argsort(numpy.abs(sigmaBar)))
        inds = inds[numpy.abs(sigmaBar)>EigenUpdater.tol]
        Pbar, sigmaBar, Qbar = Util.indSvd(Pbar, sigmaBar, Qbar, inds)
        
        SigmaBar = numpy.diag(sigmaBar)
        Qbar = Ybar.T.dot(Pbar)
        Qbar = Qbar.dot(numpy.diag(numpy.diag(Qbar.T.dot(Qbar))**-0.5))

        r = sigmaBar.shape[0]

        YQ = Y.dot(Q)
        Zeros = numpy.zeros((r, omega.shape[0]))
        D = numpy.c_[Q, Qbar]

        YYQQ = YY.dot(QQ)
        Z = D.conj().T.dot(YYQQ + YYQQ.conj().T).dot(D)
        F = numpy.c_[numpy.r_[Omega - YQ.conj().T.dot(YQ), Zeros], numpy.r_[Zeros.T, SigmaBar.conj().dot(SigmaBar)]]
        F = F + Z 

        pi, H = scipy.linalg.eigh(F)
        inds = numpy.flipud(numpy.argsort(numpy.abs(pi)))

        H = H[:, inds[0:k]]
        pi = pi[inds[0:k]]

        V = D.dot(H)
        #logging.debug("</ eigenAdd >")
        return pi, V

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def addRows(U, s, V, B, k=None): 
        """
        Find the SVD of a matrix [A ; B] where  A = U diag(s) V.T. Uses the QR 
        decomposition to find an orthogonal basis on B. 
        
        :param U: The left singular vectors of A  
        
        :param s: The singular values of A 
        
        :param V: The right singular vectors of A 
        
        :param B: The matrix to append to A 
        """
        if V.shape[0] != B.shape[1]:
            raise ValueError("U must have same number of rows as B cols")
        if s.shape[0] != U.shape[1]:
            raise ValueError("Number of cols of U must be the same size as s")
        if s.shape[0] != V.shape[1]:
            raise ValueError("Number of cols of V must be the same size as s")
    
        if k == None: 
            k = U.shape[1]
        m, p = U.shape
        r = B.shape[0]
        
        C = B.T - V.dot(V.T).dot(B.T)
        Q, R = numpy.linalg.qr(C)

        rPrime = Util.rank(C)
        Q = Q[:, 0:rPrime]
        R = R[0:rPrime, :]

        D = numpy.c_[numpy.diag(s), numpy.zeros((p, rPrime))]
        E = numpy.c_[B.dot(V), R.T]
        D = numpy.r_[D, E]
        
        G1 = numpy.c_[U, numpy.zeros((m, r))]
        G2 = numpy.c_[numpy.zeros((r, p)), numpy.eye(r)]
        G = numpy.r_[G1, G2]
        
        H = numpy.c_[V, Q]
        
        nptst.assert_array_almost_equal(G.T.dot(G), numpy.eye(G.shape[1])) 
        nptst.assert_array_almost_equal(H.T.dot(H), numpy.eye(H.shape[1])) 
        nptst.assert_array_almost_equal(G.dot(D).dot(H.T), numpy.r_[(U*s).dot(V.T), B])

        Uhat, sHat, Vhat = numpy.linalg.svd(D, full_matrices=False)
        inds = numpy.flipud(numpy.argsort(sHat))[0:k]
        Uhat, sHat, Vhat = Util.indSvd(Uhat, sHat, Vhat, inds)

        #The best rank k approximation of [A ; B]
        Utilde = G.dot(Uhat)
        Stilde = sHat
        Vtilde = H.dot(Vhat)

        return Utilde, Stilde, Vtilde

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
 def testAddRows(self): 
     
     #Test case when k = rank 
     Utilde, Stilde, Vtilde = SVDUpdate.addRows(self.U, self.s, self.V, self.C)
     
     nptst.assert_array_almost_equal(Utilde.T.dot(Utilde), numpy.eye(Utilde.shape[1]))
     nptst.assert_array_almost_equal(Vtilde.T.dot(Vtilde), numpy.eye(Vtilde.shape[1]))
     
     self.assertEquals(Stilde.shape[0], self.k)
     
     #Check we get the original solution with full SVD 
     U, s, V = numpy.linalg.svd(self.A)
     inds = numpy.flipud(numpy.argsort(s))
     U, s, V = Util.indSvd(U, s, V, inds)
     
     Utilde, Stilde, Vtilde = SVDUpdate.addRows(U, s, V, self.C)
     D = numpy.r_[self.A, self.C]
     
     nptst.assert_array_almost_equal(D, (Utilde*Stilde).dot(Vtilde.T), 4)
     
     #Check solution for partial rank SVD 
     k = 20 
     U, s, V = numpy.linalg.svd(self.A)
     inds = numpy.flipud(numpy.argsort(s))[0:k]
     U, s, V = Util.indSvd(U, s, V, inds)
     
     Utilde, Stilde, Vtilde = SVDUpdate.addRows(U, s, V, self.C)
     D = numpy.r_[(U*s).dot(V.T), self.C]
     U, s, V = numpy.linalg.svd(D)
     inds = numpy.flipud(numpy.argsort(s))[0:k]
     U, s, V = Util.indSvd(U, s, V, inds)
     
     nptst.assert_array_almost_equal((U*s).dot(V.T), (Utilde*Stilde).dot(Vtilde.T), 4)
     
     #Test if same as add cols 
     U, s, V = numpy.linalg.svd(self.A)
     inds = numpy.flipud(numpy.argsort(s))[0:k]
     U, s, V = Util.indSvd(U, s, V, inds)
     Utilde, sTilde, Vtilde = SVDUpdate.addRows(U, s, V, self.C)
     Vtilde2, sTilde2, Utilde2 = SVDUpdate.addCols(V, s, U, self.C.T)
     
     nptst.assert_array_almost_equal((Utilde*sTilde).dot(Vtilde.T),  (Utilde2*sTilde2).dot(Vtilde2.T))

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def svdArpack(X, k, kmax=None):
        """
        Perform the SVD of a sparse matrix X using ARPACK for the largest k
        singular values. Note that the input matrix should be of float dtype.

        :param X: The input matrix as scipy.sparse.csc_matrix or a LinearOperator

        :param k: The number of singular vectors/values for None for all

        :param kmax: The maximal number of iterations / maximal dimension of Krylov subspace.
        """
        if k==None:
            k = min(X.shape[0], X.shape[1])
        if kmax==None:
            kmax = SparseUtils.kmaxMultiplier*k

        if scipy.sparse.isspmatrix(X):
            L = scipy.sparse.linalg.aslinearoperator(X)
        else:
            L = X

        m, n = L.shape

        def matvec_AH_A(x):
            Ax = L.matvec(x)
            return L.rmatvec(Ax)

        AH_A = LinearOperator(matvec=matvec_AH_A, shape=(n, n), dtype=L.dtype)

        eigvals, eigvec = scipy.sparse.linalg.eigsh(AH_A, k=k, ncv=kmax)
        s2 = scipy.sqrt(eigvals)
        V2 = eigvec
        U2 = L.matmat(V2)/s2

        inds = numpy.flipud(numpy.argsort(s2))
        U2, s2, V2 = Util.indSvd(U2, s2, V2.T, inds)

        return U2, s2, V2

# 需要导入模块: from sandbox.util.Util import Util [as 别名]
# 或者: from sandbox.util.Util.Util import indSvd [as 别名]
    def eigenAdd2(omega, Q, Y1, Y2, k, debug= False):
        """
        Compute an approximation of the eigendecomposition A^*A + Y1Y2^* +Y2Y1^*
        in which Y1, Y2 are low rank matrices, Y1^*Y2=0 and A^*A = Q Omega Q*. We 
        use the rank-k approximation of A^*A: Q_k Omega_k Q_k^* and then find
        [A^*A_k + Y1Y2^* + Y2Y1^*]. If debug=False then pi, V are returned which 
        respectively correspond to all the eigenvalues/eigenvectors of 
        [A^*A_k + Y1Y2^* + Y2Y1^*]. 
        """
        #logging.debug("< eigenAdd2 >")
        Parameter.checkInt(k, 0, float('inf'))
        Parameter.checkClass(omega, numpy.ndarray)
        Parameter.checkClass(Q, numpy.ndarray)
        Parameter.checkClass(Y1, numpy.ndarray)
        Parameter.checkClass(Y2, numpy.ndarray)
        if not numpy.isrealobj(omega) or not numpy.isrealobj(Q):
            logging.warn("Eigenvalues or eigenvectors are not real")
        if not numpy.isrealobj(Y1) or not numpy.isrealobj(Y2):
            logging.warn("Y1 or Y2 are not real")
        if omega.ndim != 1:
            raise ValueError("omega must be 1-d array")
        if omega.shape[0] != Q.shape[1]:
            raise ValueError("Must have same number of eigenvalues and eigenvectors")
        if Q.shape[0] != Y1.shape[0]:
            raise ValueError("Q must have the same number of rows as Y1 rows")
        if Q.shape[0] != Y2.shape[0]:
            raise ValueError("Q must have the same number of rows as Y2 rows")
        if Y1.shape[1] != Y2.shape[1]:
            raise ValueError("Y1 must have the same number of columns as Y2 columns")

        if __debug__:
            Parameter.checkArray(omega, softCheck=True, arrayInfo="omega as input in eigenAdd2()")
            Parameter.checkArray(Q, softCheck=True, arrayInfo="Q as input in eigenAdd2()")
            Parameter.checkOrthogonal(Q, tol=EigenUpdater.tol, softCheck=True, arrayInfo="Q as input in eigenAdd2()")
            Parameter.checkArray(Y1, softCheck=True, arrayInfo="Y1 as input in eigenAdd2()")
            Parameter.checkArray(Y2, softCheck=True, arrayInfo="Y2 as input in eigenAdd2()")
            


        #Get first k eigenvectors/values of A^*A
        omega, Q = Util.indEig(omega, Q, numpy.flipud(numpy.argsort(omega))[0:k])

        QY1 = Q.conj().T.dot(Y1)
        Y1bar = Y1 - Q.dot(QY1)

        P1bar, sigma1Bar, Q1bar = Util.safeSvd(Y1bar)
        inds = numpy.arange(sigma1Bar.shape[0])[numpy.abs(sigma1Bar)>EigenUpdater.tol]
        P1bar, sigma1Bar, Q1bar = Util.indSvd(P1bar, sigma1Bar, Q1bar, inds)
        # checks on SVD decomposition of Y1bar
        if __debug__:
            Parameter.checkArray(QY1, softCheck=True, arrayInfo="QY1 in eigenAdd2()")
            Parameter.checkArray(Y1bar, softCheck=True, arrayInfo="Y1bar in eigenAdd2()")
            Parameter.checkArray(P1bar, softCheck=True, arrayInfo="P1bar in eigenAdd2()")
            if not Parameter.checkOrthogonal(P1bar, tol=EigenUpdater.tol, softCheck=True, arrayInfo="P1bar in eigenAdd2()", investigate=True):
                print ("corresponding sigma: ", sigma1Bar)
            Parameter.checkArray(sigma1Bar, softCheck=True, arrayInfo="sigma1Bar in eigenAdd2()")
            Parameter.checkArray(Q1bar, softCheck=True, arrayInfo="Q1bar in eigenAdd2()")
            if not Parameter.checkOrthogonal(Q1bar, tol=EigenUpdater.tol, softCheck=True, arrayInfo="Q1bar in eigenAdd2()"):
                print ("corresponding sigma: ", sigma1Bar)

        del Y1bar

        P1barY2 = P1bar.conj().T.dot(Y2)
        QY2 = Q.conj().T.dot(Y2)
        Y2bar = Y2 - Q.dot(QY2) - P1bar.dot(P1barY2)
        
        P2bar, sigma2Bar, Q2bar = Util.safeSvd(Y2bar)
        inds = numpy.arange(sigma2Bar.shape[0])[numpy.abs(sigma2Bar)>EigenUpdater.tol]
        P2bar, sigma2Bar, Q2bar = Util.indSvd(P2bar, sigma2Bar, Q2bar, inds)
        # checks on SVD decomposition of Y1bar
        if __debug__:
            Parameter.checkArray(P1barY2, softCheck=True, arrayInfo="P1barY2 in eigenAdd2()")
            Parameter.checkArray(QY2, softCheck=True, arrayInfo="QY2 in eigenAdd2()")
            Parameter.checkArray(Y2bar, softCheck=True, arrayInfo="Y2bar in eigenAdd2()")
            Parameter.checkArray(P2bar, softCheck=True, arrayInfo="P2bar in eigenAdd2()")
            Parameter.checkOrthogonal(P2bar, tol=EigenUpdater.tol, softCheck=True, arrayInfo="P2bar in eigenAdd2()")
            Parameter.checkArray(sigma2Bar, softCheck=True, arrayInfo="sigma2Bar in eigenAdd2()")
            Parameter.checkArray(Q2bar, softCheck=True, arrayInfo="Q2bar in eigenAdd2()")
            Parameter.checkOrthogonal(Q2bar, tol=EigenUpdater.tol, softCheck=True, arrayInfo="Q2bar in eigenAdd2()")

        del Y2bar 

        r = omega.shape[0]
        p = Y1.shape[1]
        p1 = sigma1Bar.shape[0]
        p2 = sigma2Bar.shape[0]

        D = numpy.c_[Q, P1bar, P2bar]
        del P1bar
        del P2bar 
        # rem: A*s = A.dot(diag(s)) ; A*s[:,new] = diag(s).dot(A)
        DStarY1 = numpy.r_[QY1, sigma1Bar[:,numpy.newaxis] * Q1bar.conj().T, numpy.zeros((p2, p))]
        DStarY2 = numpy.r_[QY2, P1barY2, sigma2Bar[:,numpy.newaxis] * Q2bar.conj().T]
        DStarY1Y2StarD = DStarY1.dot(DStarY2.conj().T)

        del DStarY1
        del DStarY2
        
        r = omega.shape[0]
        F = numpy.zeros((r+p1+p2, r+p1+p2))
#.........这里部分代码省略.........