Ticket #10794: trac_10794-QR-decomposition-exact.patch

File trac_10794-QR-decomposition-exact.patch, 16.0 KB (added by Rob Beezer, 12 years ago)

  • sage/matrix/matrix2.pyx 

    # HG changeset patch
# User Rob Beezer <beezer@ups.edu>
# Date 1297978335 28800
# Node ID cd4917ecc814ce7aa1a53aea09d684297fa86118
# Parent  3aa276649760e86f6a86eee05ccaf059ade6a31b
10794: QR matrix decomposition for exact rings

diff -r 3aa276649760 -r cd4917ecc814 sage/matrix/matrix2.pyx
    a b  
    58315831       
    58325832        return A
    58335833
     5834    def QR(self, full=True):
     5835        r"""
     5836        Returns a factorization of ``self`` as a unitary matrix
     5837        and an upper-triangular matrix.
     5838
     5839        INPUT:
     5840
     5841        - ``full`` - default: ``True`` - if ``True`` then the
     5842          returned matrices have dimensions as described below.
     5843          If ``False`` the ``R`` matrix has no zero rows and the
     5844          columns of ``Q`` are a basis for the column space of
     5845          ``self``.
     5846
     5847        OUTPUT:
     5848
     5849        If ``self`` is an `m\times n` matrix and ``full=True`` then this
     5850        method returns a pair of matrices: `Q` is an `m\times m` unitary
     5851        matrix (meaning its inverse is its conjugate-transpose) and `R`
     5852        is an `m\times n` upper-triangular matrix with non-negative entries
     5853        on the diagonal.  For a matrix of full rank this factorization is
     5854        unique (due to the restriction to positive entries on the diagonal).
     5855
     5856        If ``full=False`` then `Q` has `m` rows and the columns form an
     5857        orthonormal basis for the column space of ``self``.  So, in particular,
     5858        the conjugate-transpose of `Q` times `Q` will be an identity matrix.
     5859        The matrix `R` will still be upper-triangular but will also have full
     5860        rank, in particular it will lack the zero rows present in a full
     5861        factorization of a rank-deficient matrix.
     5862
     5863        The results obtained when ``full=True`` are cached,
     5864        hence `Q` and `R` are immutable matrices in this case.
     5865
     5866        .. note::
     5867
     5868            This is an exact computation, so limited to exact rings.
     5869            Also the base ring needs to have a fraction field implemented
     5870            in Sage and this field must contain square roots.  One example
     5871            is the field of algebraic numbers, ``QQbar``, as used in the
     5872            examples below.  If you need numerical results, convert the
     5873            base ring to the field of complex double numbers, ``CDF``,
     5874            which will use a faster routine that is careful about
     5875            numerical subtleties.
     5876
     5877        ALGORITHM:
     5878
     5879        "Modified Gram-Schmidt," Algorithm 8.1 of [TREFETHEN-BAU]_.
     5880
     5881        EXAMPLES:
     5882
     5883        For a nonsingular matrix, the QR decomposition is unique. ::
     5884
     5885            sage: A = matrix(QQbar, [[-2, 0, -4, -1, -1],
     5886            ...                      [-2, 1, -6, -3, -1],
     5887            ...                      [1, 1, 7, 4, 5],
     5888            ...                      [3, 0, 8, 3, 3],
     5889            ...                      [-1, 1, -6, -6, 5]])
     5890            sage: Q, R = A.QR()
     5891            sage: Q
     5892            [ -0.4588314677411235?  -0.1260506983326509?   0.3812120831224489?  -0.3945737113384181?  -0.6874400625963308?]
     5893            [ -0.4588314677411235?   0.4726901187474409? -0.05198346588033394?   0.7172941251646595?  -0.2209628772631064?]
     5894            [  0.2294157338705618?   0.6617661662464172?   0.6619227988762521?  -0.1808720937375480?   0.1964114464560946?]
     5895            [  0.6882472016116853?   0.1890760474989764?  -0.2044682991293135?   0.0966302966543065?  -0.6628886317893191?]
     5896            [ -0.2294157338705618?   0.5357154679137663?  -0.6099393329959182?  -0.5364220314271115?  0.02455143080701182?]
     5897            sage: R
     5898            [  4.358898943540674? -0.4588314677411235?   13.07669683062202?   6.194224814505168?   2.982404540317303?]
     5899            [                   0   1.670171752907625?  0.5987408170800917?  -1.292019657909672?   6.207996892883057?]
     5900            [                   0                    0   5.444401659866974?   5.468660610611130? -0.6827161852283857?]
     5901            [                   0                    0                    0   1.027626039419836?  -3.619300149686620?]
     5902            [                   0                    0                    0                    0 0.02455143080701182?]
     5903            sage: Q.conjugate_transpose()*Q
     5904            [1.000000000000000?            0.?e-37            0.?e-36            0.?e-34            0.?e-31]
     5905            [           0.?e-37 1.000000000000000?            0.?e-36            0.?e-34            0.?e-31]
     5906            [           0.?e-36            0.?e-36 1.000000000000000?            0.?e-34            0.?e-31]
     5907            [           0.?e-34            0.?e-34            0.?e-34 1.000000000000000?            0.?e-31]
     5908            [           0.?e-31            0.?e-31            0.?e-31            0.?e-31 1.000000000000000?]
     5909            sage: Q*R == A
     5910            True
     5911
     5912
     5913        An example with complex numbers in ``QQbar``, the field of algebraic
     5914        numbers. ::
     5915
     5916            sage: A = matrix(QQbar, [[-8, 4*I + 1, -I + 2, 2*I + 1],
     5917            ...                      [1, -2*I - 1, -I + 3, -I + 1],
     5918            ...                      [I + 7, 2*I + 1, -2*I + 7, -I + 1],
     5919            ...                      [I + 2, 0, I + 12, -1]])
     5920            sage: Q, R = A.QR()
     5921            sage: Q
     5922            [                          -0.7302967433402215?    0.2070566455055649? + 0.5383472783144687?*I    0.2463049809998642? - 0.0764456358723292?*I   0.2381617683194332? - 0.10365960327796942?*I]
     5923            [                           0.0912870929175277?   -0.2070566455055649? - 0.3778783780476559?*I    0.3786559533863033? - 0.1952221495524667?*I    0.7012444502144683? - 0.3643711650986595?*I]
     5924            [   0.6390096504226938? + 0.0912870929175277?*I    0.1708217325420910? + 0.6677576817554466?*I -0.03411475806452072? + 0.04090198741767143?*I    0.3140171085506764? - 0.0825191718705412?*I]
     5925            [   0.1825741858350554? + 0.0912870929175277?*I  -0.03623491296347385? + 0.0724698259269477?*I   0.8632284069415110? + 0.06322839976356195?*I  -0.4499694867611521? - 0.01161191812089175?*I]
     5926            sage: R
     5927            [                          10.95445115010333?               0.?e-37 - 1.917028951268082?*I    5.385938482134133? - 2.190890230020665?*I  -0.2738612787525831? - 2.190890230020665?*I]
     5928            [                                           0                           4.829596256417300?   -0.869637911123373? - 5.864879483945125?*I   0.993871898426712? - 0.3054085521207082?*I]
     5929            [                                           0                                            0                           12.00160760935814? -0.2709533402297274? + 0.4420629644486323?*I]
     5930            [                                           0                                            0                                            0               1.942963944258992? + 0.?e-33*I]
     5931            sage: Q.conjugate_transpose()*Q
     5932            [1.000000000000000? + 0.?e-37*I            0.?e-36 + 0.?e-36*I            0.?e-33 + 0.?e-33*I            0.?e-33 + 0.?e-33*I]
     5933            [           0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-36*I            0.?e-33 + 0.?e-33*I            0.?e-33 + 0.?e-33*I]
     5934            [           0.?e-33 + 0.?e-33*I            0.?e-33 + 0.?e-33*I 1.000000000000000? + 0.?e-34*I            0.?e-33 + 0.?e-33*I]
     5935            [           0.?e-33 + 0.?e-33*I            0.?e-33 + 0.?e-33*I            0.?e-33 + 0.?e-33*I 1.000000000000000? + 0.?e-33*I]
     5936            sage: Q*R - A
     5937            [                  0 0.?e-36 + 0.?e-36*I 0.?e-32 + 0.?e-32*I 0.?e-33 + 0.?e-33*I]
     5938            [                  0 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-33 + 0.?e-33*I]
     5939            [0.?e-36 + 0.?e-36*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-33 + 0.?e-33*I]
     5940            [0.?e-37 + 0.?e-37*I 0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-33 + 0.?e-33*I]
     5941
     5942        A rank-deficient rectangular matrix, with both values of the ``full`` keyword.  ::
     5943
     5944            sage: A = matrix(QQbar, [[2, -3, 3],
     5945            ...                      [-1, 1, -1],
     5946            ...                      [-1, 3, -3],
     5947            ...                      [-5, 1, -1]])
     5948            sage: Q, R = A.QR()
     5949            sage: Q
     5950            [  0.3592106040535498?  -0.5693261797050169?   0.7239227659930268?   0.1509015305256380?]
     5951            [ -0.1796053020267749?   0.1445907757980996?                     0   0.9730546968377341?]
     5952            [ -0.1796053020267749?   0.7048800320157352?    0.672213996993525?  -0.1378927778941174?]
     5953            [ -0.8980265101338745?  -0.3976246334447737?   0.1551263069985058? -0.10667177157846818?]
     5954            sage: R
     5955            [ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
     5956            [                  0  3.569584777515583? -3.569584777515583?]
     5957            [                  0                   0                   0]
     5958            [                  0                   0                   0]
     5959            sage: Q.conjugate_transpose()*Q
     5960            [                 1            0.?e-37            0.?e-37            0.?e-38]
     5961            [           0.?e-37                  1            0.?e-37            0.?e-38]
     5962            [           0.?e-37            0.?e-37                  1            0.?e-38]
     5963            [           0.?e-38            0.?e-37            0.?e-38 1.000000000000000?]
     5964
     5965            sage: Q, R = A.QR(full=False)
     5966            sage: Q
     5967            [ 0.3592106040535498? -0.5693261797050169?]
     5968            [-0.1796053020267749?  0.1445907757980996?]
     5969            [-0.1796053020267749?  0.7048800320157352?]
     5970            [-0.8980265101338745? -0.3976246334447737?]
     5971            sage: R
     5972            [ 5.567764362830022? -2.694079530401624?  2.694079530401624?]
     5973            [                  0  3.569584777515583? -3.569584777515583?]
     5974            sage: Q.conjugate_transpose()*Q
     5975            [      1 0.?e-37]
     5976            [0.?e-37       1]
     5977
     5978        Another rank-deficient rectangular matrix, with complex entries,
     5979        as a reduced decomposition. ::
     5980
     5981            sage: A = matrix(QQbar, [[-3*I - 3, I - 3, -12*I + 1, -2],
     5982            ...                      [-I - 1, -2, 5*I - 1, -I - 2],
     5983            ...                      [-4*I - 4, I - 5, -7*I, -I - 4]])
     5984            sage: Q, R = A.QR(full=False)
     5985            sage: Q
     5986            [ -0.4160251471689219? - 0.4160251471689219?*I   0.5370861555295747? + 0.1790287185098583?*I]
     5987            [ -0.1386750490563073? - 0.1386750490563073?*I  -0.7519206177414046? - 0.2506402059138015?*I]
     5988            [ -0.5547001962252291? - 0.5547001962252291?*I -0.2148344622118299? - 0.07161148740394329?*I]
     5989            sage: R
     5990            [                        7.211102550927979?  3.328201177351375? - 5.269651864139676?*I   7.904477796209515? + 8.45917799243475?*I  4.021576422632911? - 2.634825932069838?*I]
     5991            [                                         0                         1.074172311059150?  -1.611258466588724? - 9.13046464400277?*I 1.611258466588724? + 0.5370861555295747?*I]
     5992            sage: Q.conjugate_transpose()*Q
     5993            [1.000000000000000? + 0.?e-37*I            0.?e-36 + 0.?e-36*I]
     5994            [           0.?e-36 + 0.?e-36*I 1.000000000000000? + 0.?e-36*I]
     5995            sage: Q*R-A
     5996            [0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-35*I]
     5997            [0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-34*I 0.?e-35 + 0.?e-35*I]
     5998            [0.?e-37 + 0.?e-37*I 0.?e-36 + 0.?e-36*I 0.?e-35 + 0.?e-35*I 0.?e-35 + 0.?e-36*I]
     5999
     6000        Results of full decompositions are cached and thus returned
     6001        immutable.  ::
     6002
     6003            sage: A = random_matrix(QQbar, 2, 2)
     6004            sage: Q, R = A.QR()
     6005            sage: Q.is_mutable()
     6006            False
     6007            sage: R.is_mutable()
     6008            False
     6009
     6010        Trivial cases return trivial results of the correct size,
     6011        and we check `Q` itself in one case.  ::
     6012
     6013            sage: A = zero_matrix(QQbar, 0, 10)
     6014            sage: Q, R = A.QR()
     6015            sage: Q.nrows(), Q.ncols()
     6016            (0, 0)
     6017            sage: R.nrows(), R.ncols()
     6018            (0, 10)
     6019            sage: A = zero_matrix(QQbar, 3, 0)
     6020            sage: Q, R = A.QR()
     6021            sage: Q.nrows(), Q.ncols()
     6022            (3, 3)
     6023            sage: R.nrows(), R.ncols()
     6024            (3, 0)
     6025            sage: Q
     6026            [1 0 0]
     6027            [0 1 0]
     6028            [0 0 1]
     6029
     6030        TESTS:
     6031
     6032        Inexact rings are caught and ``CDF`` suggested.  ::
     6033
     6034            sage: A = matrix(RealField(100), 2, range(4))
     6035            sage: A.QR()
     6036            Traceback (most recent call last):
     6037            ...
     6038            NotImplementedError: QR decomposition is implemented over exact rings, try CDF for numerical results, not Real Field with 100 bits of precision
     6039
     6040        Without a fraction field, we cannot hope to run the algorithm. ::
     6041
     6042            sage: A = matrix(Integers(6), 2, range(4))
     6043            sage: A.QR()
     6044            Traceback (most recent call last):
     6045            ...
     6046            ValueError: QR decomposition needs a fraction field of Ring of integers modulo 6
     6047
     6048        The biggest obstacle is making unit vectors, thus requiring square
     6049        roots, though some small cases pass through.  ::
     6050
     6051            sage: A = matrix(ZZ, 3, range(9))
     6052            sage: A.QR()
     6053            Traceback (most recent call last):
     6054            ...
     6055            TypeError: QR decomposition unable to compute square roots in Rational Field
     6056
     6057            sage: A = matrix(ZZ, 2, range(4))
     6058            sage: Q, R = A.QR()
     6059            sage: Q
     6060            [0 1]
     6061            [1 0]
     6062            sage: R
     6063            [2 3]
     6064            [0 1]
     6065
     6066        REFERENCES:
     6067
     6068        .. [TREFETHEN-BAU] Trefethen, Lloyd N., Bau, David, III
     6069           "Numerical Linear Algebra"
     6070           SIAM, Philadelphia, 1997.
     6071
     6072        AUTHOR:
     6073
     6074        - Rob Beezer (2011-02-17)
     6075        """
     6076        from sage.modules.free_module_element import zero_vector
     6077        from sage.matrix.constructor import zero_matrix, matrix
     6078        from sage.functions.other import sqrt
     6079
     6080        if full:
     6081            QR = self.fetch('QR_factors')
     6082            if QR is not None:
     6083                return QR
     6084        R = self.base_ring()
     6085        if not R.is_exact():
     6086            raise NotImplementedError('QR decomposition is implemented over exact rings, try CDF for numerical results, not %s' % R)
     6087        try:
     6088            F = R.fraction_field()
     6089        except:
     6090            raise ValueError("QR decomposition needs a fraction field of %s" % R)
     6091        m = self.nrows()
     6092        n = self.ncols()
     6093
     6094        R = zero_matrix(F, m, n)
     6095        V = self.columns(copy=True)
     6096        Q = []
     6097        row = 0  # row of R being filled
     6098        for i in range(n):
     6099            v = V[i]
     6100            hip = v.hermitian_inner_product(v)
     6101            if hip != 0:
     6102                try:
     6103                    scale = sqrt(hip)
     6104                    q = (1/scale)*v
     6105                    Q.append(q)
     6106                    R[row,i] = scale
     6107                    for j in range(i+1, n):
     6108                        R[row,j] = q.hermitian_inner_product(V[j])
     6109                        V[j] = V[j] - R[row,j]*q
     6110                    row = row + 1
     6111                except TypeError:
     6112                    raise TypeError('QR decomposition unable to compute square roots in %s' % F)
     6113        # complete to full orthonormal basis, or reduce to truncated R
     6114        if full:
     6115            Qt = matrix(Q) # as rows here
     6116            if Qt.nrows() == 0:
     6117                Qt = zero_matrix(F, 0, m)
     6118            orthogonal = Qt.right_kernel().basis_matrix().transpose()
     6119            Qperp, _ = orthogonal.QR(full=False)
     6120            Q = Q + Qperp.columns()
     6121        else:
     6122            R = R[0:len(Q), 0:n]
     6123        Q = matrix(Q).transpose()
     6124        # Adjust rows of Q if empty
     6125        if Q.ncols() == 0:
     6126            Q = zero_matrix(F, m, 0)
     6127        QR = (Q, R)
     6128        if full:
     6129            Q.set_immutable()
     6130            R.set_immutable()
     6131            self.cache('QR_factors', QR)
     6132        return QR
     6133
    58346134    def gram_schmidt(self):
    58356135        r"""
    58366136        Return the matrix G whose rows are obtained from the rows of self