Ticket #10795: trac_10795-QR-decomposition-double-dense-v2.patch

File trac_10795-QR-decomposition-double-dense-v2.patch, 8.9 KB (added by rbeezer, 9 years ago)
  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1297983491 28800
    # Node ID 3fe937195454ce8dacdd2ebbf4453a8ab5d7f756
    # Parent  cc97b3f58cd2ecf9ab29ed1eb8b079126680ef47
    10795: fix trivial case, and upgrade docs for double dense matrix QR decomposition
    * * *
    10795: numerical noise fixes
    
    diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
    a b  
    22152215        return USV
    22162216
    22172217    def QR(self):
    2218         """
    2219         Return the Q,R factorization of a real matrix A.
    2220 
    2221         The computed decomposition is cached and returned on subsequent calls.
     2218        r"""
     2219        Returns a factorization into a unitary matrix and an
     2220        upper-triangular matrix.
    22222221
    22232222        INPUT:
    22242223
    2225         - self -- a real matrix A
     2224        Any matrix over ``RDF`` or ``CDF``.
    22262225
    22272226        OUTPUT:
    22282227
    2229         - Q, R -- immutable matrices such that A = Q*R such that the columns of Q are
    2230           orthogonal (i.e., $Q^t Q = I$), and R is upper triangular.
    2231 
    2232         EXAMPLES::
     2228        `Q`, `R` -- a pair of matrices such that `A = QR` such that
     2229        the columns of Q are orthogonal (i.e., `Q^\ast Q = I`),
     2230        and R is upper triangular, where `Q^\ast` is the
     2231        conjugate-transpose in the complex case and the transpose
     2232        in the real case. If `A` is an `m\times n` matrix and
     2233        ``full=True`` then this method returns a pair of matrices:
     2234        `Q` is an `m\times m` unitary matrix (meaning its inverse
     2235        is its conjugate-transpose) and `R` is an `m\times n`
     2236        upper-triangular matrix.  For a matrix of full rank this
     2237        factorization is unique up to the sign of the diagonal entries.
     2238
     2239        The resulting decomposition is cached.
     2240
     2241        ALGORITHM:
    22332242       
    2234             sage: m = matrix(RDF,3,range(0, 12)); m
    2235             [ 0.0  1.0  2.0  3.0]
    2236             [ 4.0  5.0  6.0  7.0]
    2237             [ 8.0  9.0 10.0 11.0]
    2238             sage: Q,R = m.QR()
    2239             sage: Q*R
    2240             [ 0.0  1.0  2.0  3.0]
    2241             [ 4.0  5.0  6.0  7.0]
    2242             [ 8.0  9.0 10.0 11.0]
    2243 
    2244         Note that Q is an orthogonal matrix::
    2245 
    2246             sage: (Q*Q.transpose()).zero_at(1e-10)
     2243        Calls "linalg.qr" from SciPy, which is in turn an
     2244        interface to LAPACK routines.
     2245
     2246        EXAMPLES:
     2247
     2248        Over the reals, the inverse of `Q` is its transpose,
     2249        since including a conjugate has no effect.  In the real
     2250        case, we say `Q` is orthogonal. ::
     2251
     2252            sage: A = matrix(RDF, [[-2, 0, -4, -1, -1],
     2253            ...                    [-2, 1, -6, -3, -1],
     2254            ...                    [1, 1, 7, 4, 5],
     2255            ...                    [3, 0, 8, 3, 3],
     2256            ...                    [-1, 1, -6, -6, 5]])
     2257            sage: Q, R = A.QR()
     2258            sage: Q.round(6).zero_at(10^-6)
     2259            [-0.458831  0.126051 -0.381212 -0.394574  -0.68744]
     2260            [-0.458831  -0.47269  0.051983  0.717294 -0.220963]
     2261            [ 0.229416 -0.661766 -0.661923 -0.180872  0.196411]
     2262            [ 0.688247 -0.189076  0.204468   0.09663 -0.662889]
     2263            [-0.229416 -0.535715  0.609939 -0.536422  0.024551]
     2264            sage: R.round(6).zero_at(10^-6)
     2265            [ 4.358899 -0.458831 13.076697  6.194225  2.982405]
     2266            [      0.0 -1.670172 -0.598741   1.29202 -6.207997]
     2267            [      0.0       0.0 -5.444402 -5.468661  0.682716]
     2268            [      0.0       0.0       0.0  1.027626   -3.6193]
     2269            [      0.0       0.0       0.0       0.0  0.024551]
     2270            sage: (Q*Q.transpose()).zero_at(10^-14)
     2271            [1.0 0.0 0.0 0.0 0.0]
     2272            [0.0 1.0 0.0 0.0 0.0]
     2273            [0.0 0.0 1.0 0.0 0.0]
     2274            [0.0 0.0 0.0 1.0 0.0]
     2275            [0.0 0.0 0.0 0.0 1.0]
     2276            sage: (Q*R - A).zero_at(10^-14)
     2277            [0.0 0.0 0.0 0.0 0.0]
     2278            [0.0 0.0 0.0 0.0 0.0]
     2279            [0.0 0.0 0.0 0.0 0.0]
     2280            [0.0 0.0 0.0 0.0 0.0]
     2281            [0.0 0.0 0.0 0.0 0.0]
     2282
     2283        Now over the complex numbers, demonstrating that the SciPy libraries
     2284        are (properly) using the Hermitian inner product, so that `Q` is
     2285        a unitary matrix (its inverse is the conjugate-transpose).  ::
     2286
     2287            sage: A = matrix(CDF, [[-8, 4*I + 1, -I + 2, 2*I + 1],
     2288            ...                    [1, -2*I - 1, -I + 3, -I + 1],
     2289            ...                    [I + 7, 2*I + 1, -2*I + 7, -I + 1],
     2290            ...                    [I + 2, 0, I + 12, -1]])
     2291            sage: Q, R = A.QR()
     2292            sage: Q.round(6).zero_at(10^-6)
     2293            [             -0.730297  0.207057 + 0.538347*I -0.246305 + 0.076446*I   0.238162 - 0.10366*I]
     2294            [              0.091287 -0.207057 - 0.377878*I -0.378656 + 0.195222*I  0.701244 - 0.364371*I]
     2295            [  0.63901 + 0.091287*I  0.170822 + 0.667758*I  0.034115 - 0.040902*I  0.314017 - 0.082519*I]
     2296            [ 0.182574 + 0.091287*I  -0.036235 + 0.07247*I -0.863228 - 0.063228*I -0.449969 - 0.011612*I]
     2297            sage: R.round(6).zero_at(10^-6)
     2298            [             10.954451            -1.917029*I   5.385938 - 2.19089*I  -0.273861 - 2.19089*I]
     2299            [                   0.0               4.829596 -0.869638 - 5.864879*I  0.993872 - 0.305409*I]
     2300            [                   0.0                    0.0             -12.001608  0.270953 - 0.442063*I]
     2301            [                   0.0                    0.0                    0.0               1.942964]
     2302            sage: (Q.conjugate().transpose()*Q).zero_at(10^-15)
     2303            [1.0 0.0 0.0 0.0]
     2304            [0.0 1.0 0.0 0.0]
     2305            [0.0 0.0 1.0 0.0]
     2306            [0.0 0.0 0.0 1.0]
     2307            sage: (Q*R - A).zero_at(10^-14)
     2308            [0.0 0.0 0.0 0.0]
     2309            [0.0 0.0 0.0 0.0]
     2310            [0.0 0.0 0.0 0.0]
     2311            [0.0 0.0 0.0 0.0]
     2312
     2313        An example of a rectangular matrix that is also rank-deficient.
     2314        If you run this example yourself, you may see a very small, nonzero
     2315        entries in the third row, in the third column, even though the exact
     2316        version of the matrix has rank 2.  The final two columns of ``Q`` span the left kernel of ``A`` and different platforms will compute different
     2317        bases so we do not exhibit the actual matrix.  ::
     2318
     2319            sage: Arat = matrix(QQ, [[2, -3, 3],
     2320            ...                      [-1, 1, -1],
     2321            ...                      [-1, 3, -3],
     2322            ...                      [-5, 1, -1]])
     2323            sage: Arat.rank()
     2324            2
     2325            sage: A = Arat.change_ring(CDF)
     2326            sage: Q, R = A.QR()
     2327            sage: R.round(6).zero_at(10^-6)
     2328            [-5.567764   2.69408  -2.69408]
     2329            [      0.0 -3.569585  3.569585]
     2330            [      0.0       0.0       0.0]
     2331            [      0.0       0.0       0.0]
     2332            sage: (Q.conjugate_transpose()*Q).zero_at(10^-14)
     2333            [1.0 0.0 0.0 0.0]
     2334            [0.0 1.0 0.0 0.0]
     2335            [0.0 0.0 1.0 0.0]
     2336            [0.0 0.0 0.0 1.0]
     2337
     2338        Results are cached, meaning they are immutable matrices.
     2339        Make a copy if you need to manipulate a result. ::
     2340
     2341            sage: A = random_matrix(CDF, 2, 2)
     2342            sage: Q, R = A.QR()
     2343            sage: Q.is_mutable()
     2344            False
     2345            sage: R.is_mutable()
     2346            False
     2347            sage: Q[0,0] = 0
     2348            Traceback (most recent call last):
     2349            ...
     2350            ValueError: matrix is immutable; please change a copy instead (i.e., use copy(M) to change a copy of M).
     2351            sage: Qcopy = copy(Q)
     2352            sage: Qcopy[0,0] = 679
     2353            sage: Qcopy[0,0]
     2354            679.0
     2355
     2356        TESTS:
     2357
     2358        Trivial cases return trivial results of the correct size,
     2359        and we check `Q` itself in one case, verifying a fix for
     2360        :trac:`10795`.  ::
     2361
     2362            sage: A = zero_matrix(RDF, 0, 10)
     2363            sage: Q, R = A.QR()
     2364            sage: Q.nrows(), Q.ncols()
     2365            (0, 0)
     2366            sage: R.nrows(), R.ncols()
     2367            (0, 10)
     2368            sage: A = zero_matrix(RDF, 3, 0)
     2369            sage: Q, R = A.QR()
     2370            sage: Q.nrows(), Q.ncols()
     2371            (3, 3)
     2372            sage: R.nrows(), R.ncols()
     2373            (3, 0)
     2374            sage: Q
    22472375            [1.0 0.0 0.0]
    22482376            [0.0 1.0 0.0]
    22492377            [0.0 0.0 1.0]
    2250 
    2251         The result is immutable::
    2252        
    2253             sage: Q[0,0] = 0
    2254             Traceback (most recent call last):
    2255                 ...
    2256             ValueError: matrix is immutable; please change a copy instead (i.e., use copy(M) to change a copy of M).
    2257             sage: R.is_immutable()
    2258             True
    2259            
    22602378        """
    22612379        global scipy
    22622380        cdef Matrix_double_dense Q,R
    22632381
    22642382        if self._nrows == 0 or self._ncols == 0:
    2265             return self.new_matrix(self._nrows, self._nrows), self.new_matrix()
     2383            return self.new_matrix(self._nrows, self._nrows, entries=1), self.new_matrix()
    22662384
    22672385        QR = self.fetch('QR_factors')
    22682386        if QR is None: