Ticket #11027: trac_11027-schur-decomposition.patch

File trac_11027-schur-decomposition.patch, 13.2 KB (added by rbeezer, 10 years ago)
  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1301065587 25200
    # Node ID bd05f5341465c21a937a026a98584faffa449bbe
    # Parent  6a679959b54b7975af1841df4dd77d28b0faef28
    11027: Schur matrix decomposition
    
    diff -r 6a679959b54b -r bd05f5341465 sage/matrix/matrix_double_dense.pyx
    a b  
    13341334        self.cache(key, b)
    13351335        return b
    13361336
     1337    def schur(self, base_ring=None):
     1338        r"""
     1339        Returns the Schur decomposition of the matrix.
     1340
     1341        INPUT:
     1342
     1343        - ``base_ring`` - optional, defaults to the base ring of ``self``.
     1344          Use this to request the base ring of the returned matrices, which
     1345          will affect the format of the results.
     1346
     1347        OUTPUT:
     1348
     1349        A pair of immutable matrices.  The first is a unitary matrix `Q`.
     1350        The second, `T`, is upper triangular when returned over the complex
     1351        numbers, while it is almost upper triangular over the reals.  In the
     1352        latter case, there can be some `2\times 2` blocks on the diagonal
     1353        which represent a pair of conjugate complex eigenvalues of ``self``.
     1354
     1355        If ``self`` is the matrix `A`, then
     1356
     1357        .. math::
     1358
     1359            A = QT({\overline Q})^t
     1360
     1361        where the latter matrix is the conjugate-transpose of ``Q``, which
     1362        is also the inverse of ``Q`` since it is unitary.
     1363
     1364        Note that in the case of a normal matrix (Hermitian, symmetric, and
     1365        others), the upper triangular matrix is  a diagonal matrix with
     1366        eigenvalues of ``self`` on the diagonal, and the unitary matrix
     1367        has columns that form an orthonormal basis composed of eigenvectors
     1368        of ``self``.  This is known as "orthonormal diagonalization".
     1369
     1370        EXAMPLES:
     1371
     1372        First over the complexes.  The similar matrix is always upper
     1373        triangular in this case.  ::
     1374
     1375            sage: A = matrix(CDF, 4, 4, range(16)) + matrix(CDF, 4, 4, [x^3*I for x in range(0, 16)])
     1376            sage: Q, T = A.schur()
     1377            sage: Q.round(4)
     1378            [   0.006 - 0.001*I  0.0425 + 0.1018*I -0.6446 + 0.2201*I  0.3981 + 0.6045*I]
     1379            [  0.0875 - 0.005*I  0.1984 + 0.5537*I -0.5022 + 0.1944*I -0.3414 - 0.4896*I]
     1380            [ 0.3572 - 0.0159*I   0.2444 + 0.685*I  0.4499 - 0.1643*I  0.1983 + 0.2728*I]
     1381            [  0.929 - 0.0377*I   -0.111 - 0.317*I  -0.1212 + 0.044*I -0.0466 - 0.0627*I]
     1382            sage: T.round(4)
     1383            [  30.7332 + 4648.5418*I  -2342.252 + 767.2915*I -477.0197 - 1460.6287*I  -532.7363 + 318.6099*I]
     1384            [                      0    -0.1841 - 159.0576*I      151.3315 + 0.458*I    -15.2369 - 67.0242*I]
     1385            [                      0                       0      -0.5239 + 11.158*I       6.2073 - 1.1208*I]
     1386            [                      0                       0                       0      -0.0252 - 0.6422*I]
     1387            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1388            [1.0   0   0   0]
     1389            [  0 1.0   0   0]
     1390            [  0   0 1.0   0]
     1391            [  0   0   0 1.0]
     1392            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1393            True
     1394            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1395            [0 0 0 0]
     1396            [0 0 0 0]
     1397            [0 0 0 0]
     1398            [0 0 0 0]
     1399
     1400        We begin with a real matrix but ask for a decomposition over the
     1401        complexes.  The result will always have an upper triangular
     1402        matrix for ``T``. ::
     1403       
     1404            sage: A = matrix(RDF, 4, 4, [x^3 for x in range(16)])
     1405            sage: Q, T = A.schur(base_ring=CDF)
     1406            sage: Q.round(4)
     1407            [             0.006  0.1102 - 0.0003*I  0.6792 + 0.0468*I -0.3382 - 0.6403*I]
     1408            [            0.0877  0.5882 - 0.0015*I   0.5374 + 0.037*I  0.2787 + 0.5277*I]
     1409            [            0.3575  0.7273 - 0.0019*I -0.4781 - 0.0329*I -0.1573 - 0.2979*I]
     1410            [ 0.9298 + 0.0001*I -0.3359 + 0.0009*I  0.1288 + 0.0089*I   0.0364 + 0.069*I]
     1411            sage: T.round(4)
     1412            [             4648.5536   2464.7695 - 6.6819*I   1532.623 + 105.406*I -290.0435 - 549.0666*I]
     1413            [                     0              -159.0726  -150.9268 - 10.7911*I     31.954 + 60.8896*I]
     1414            [                     0                      0                11.1621     -3.3197 - 5.3583*I]
     1415            [                     0                      0                      0                -0.6431]
     1416            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1417            [1.0   0   0   0]
     1418            [  0 1.0   0   0]
     1419            [  0   0 1.0   0]
     1420            [  0   0   0 1.0]
     1421            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1422            True
     1423            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1424            [0 0 0 0]
     1425            [0 0 0 0]
     1426            [0 0 0 0]
     1427            [0 0 0 0]
     1428
     1429        Now totally over the reals.  But with complex eigenvalues, the
     1430        similar matrix may not be upper triangular. But "at worst" there
     1431        may be some `2\times 2` blocks on the diagonal which represent
     1432        a pair of conjugate complex eigenvalues. These blocks will then
     1433        just interrupt the zeros below the main diagonal.
     1434
     1435            sage: A = matrix(RDF, 4, 4, [[1, 0, -3, -1],
     1436            ...                          [4, -16, -7, 0],
     1437            ...                          [1, 21, 1, -2],
     1438            ...                          [26, -1, -2, 1]])
     1439            sage: Q, T = A.schur()
     1440            sage: Q.round(4)
     1441            [ -0.088  0.5832  0.7724  0.2356]
     1442            [ 0.1078 -0.2778  0.4746 -0.8282]
     1443            [-0.2985  0.6989 -0.4072 -0.5066]
     1444            [ 0.9442   0.307 -0.1109 -0.0437]
     1445            sage: T.round(4)
     1446            [ -0.7898   15.455  15.6639  14.4823]
     1447            [ -0.3534  -0.7898     15.3 -13.7298]
     1448            [     0.0      0.0  -5.7102  16.0869]
     1449            [     0.0      0.0  -4.3681  -5.7102]
     1450            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1451            [1.0 0.0 0.0 0.0]
     1452            [0.0 1.0 0.0 0.0]
     1453            [0.0 0.0 1.0 0.0]
     1454            [0.0 0.0 0.0 1.0]
     1455            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1456            False
     1457            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i-1)])
     1458            True
     1459            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1460            [0.0 0.0 0.0 0.0]
     1461            [0.0 0.0 0.0 0.0]
     1462            [0.0 0.0 0.0 0.0]
     1463            [0.0 0.0 0.0 0.0]
     1464
     1465        Starting with complex numbers and requesting a result over the reals
     1466        will never happen.  ::
     1467
     1468            sage: A = matrix(CDF, 2, 2, [[2+I, -1+3*I], [5-4*I, 2-7*I]])
     1469            sage: A.schur(base_ring=RDF)
     1470            Traceback (most recent call last):
     1471            ...
     1472            TypeError: unable to convert input matrix over CDF to a matrix over RDF
     1473
     1474        If theory predicts your matrix is real, but it contains some
     1475        very small imaginary parts, you can specify the cutoff for "small"
     1476        imaginary parts, Then request the output as real matrices, and let
     1477        the routine do the rest. ::
     1478
     1479            sage: A = matrix(RDF, 2, 2, [1, 1, -1, 0]) + matrix(CDF, 2, 2, [1.0e-14*I]*4)
     1480            sage: B = A.zero_at(1.0e-12)
     1481            sage: B.parent()
     1482            Full MatrixSpace of 2 by 2 dense matrices over Complex Double Field
     1483            sage: Q, T = B.schur(RDF)
     1484            sage: Q.parent()
     1485            Full MatrixSpace of 2 by 2 dense matrices over Real Double Field
     1486            sage: T.parent()
     1487            Full MatrixSpace of 2 by 2 dense matrices over Real Double Field
     1488            sage: Q.round(6)
     1489            [ 0.707107  0.707107]
     1490            [-0.707107  0.707107]
     1491            sage: T.round(6)
     1492            [ 0.5  1.5]
     1493            [-0.5  0.5]
     1494            sage: (Q*T*Q.conjugate().transpose()-B).zero_at(1.0e-11)
     1495            [0 0]
     1496            [0 0]
     1497
     1498        A Hermitian matrix has real eigenvalues, so the similar matrix
     1499        will be upper triangular.  Furthermore, a Hermitian matrix is
     1500        diagonalizable with respect to an orthonormal basis, composed
     1501        of eigenvectors of the matrix.  Here that basis is the set of
     1502        columns of the unitary matrix.  ::
     1503
     1504            sage: A = matrix(CDF, [[        52,   -9*I - 8,    6*I - 187,  -188*I + 2],
     1505            ...                    [   9*I - 8,         12,   -58*I + 59,   30*I + 42],
     1506            ...                    [-6*I - 187,  58*I + 59,         2677, 2264*I + 65],
     1507            ...                    [ 188*I + 2, -30*I + 42, -2264*I + 65,       2080]])
     1508            sage: Q, T = A.schur()
     1509            sage: T = T.zero_at(1.0e-12).change_ring(RDF)
     1510            sage: T.round(6)
     1511            [4680.13301        0.0        0.0        0.0]
     1512            [       0.0 102.715967        0.0        0.0]
     1513            [       0.0        0.0  35.039344        0.0]
     1514            [       0.0        0.0        0.0    3.11168]
     1515            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1516            [1.0   0   0   0]
     1517            [  0 1.0   0   0]
     1518            [  0   0 1.0   0]
     1519            [  0   0   0 1.0]
     1520            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1521            [0 0 0 0]
     1522            [0 0 0 0]
     1523            [0 0 0 0]
     1524            [0 0 0 0]
     1525
     1526        Similarly, a real symmetric matrix has only real eigenvalues, and
     1527        there is an orthonormal basis composed of eigenvectors of the matrix.
     1528
     1529            sage: A = matrix(RDF, [[ 1, -2, 5, -3],
     1530            ...                    [-2,  9, 1,  5],
     1531            ...                    [ 5,  1, 3 , 7],
     1532            ...                    [-3,  5, 7, -8]])
     1533            sage: Q, T = A.schur()
     1534            sage: Q.round(4)
     1535            [-0.3027  -0.751   0.576 -0.1121]
     1536            [  0.139 -0.3892 -0.2648  0.8713]
     1537            [ 0.4361   0.359  0.7599  0.3217]
     1538            [ -0.836  0.3945  0.1438  0.3533]
     1539            sage: T.round(4)
     1540            [-13.5698      0.0      0.0      0.0]
     1541            [     0.0  -0.8508     -0.0     -0.0]
     1542            [     0.0      0.0   7.7664      0.0]
     1543            [     0.0      0.0      0.0  11.6542]
     1544            sage: (Q*Q.transpose()).zero_at(1.0e-12)
     1545            [1.0 0.0 0.0 0.0]
     1546            [0.0 1.0 0.0 0.0]
     1547            [0.0 0.0 1.0 0.0]
     1548            [0.0 0.0 0.0 1.0]
     1549            sage: (Q*T*Q.transpose()-A).zero_at(1.0e-11)
     1550            [0.0 0.0 0.0 0.0]
     1551            [0.0 0.0 0.0 0.0]
     1552            [0.0 0.0 0.0 0.0]
     1553            [0.0 0.0 0.0 0.0]
     1554
     1555        The results are cached, both as a real factorization and also as a
     1556        complex factorization.  This means the returned matrices are
     1557        immutable.  ::
     1558
     1559            sage: A = matrix(RDF, 2, 2, [[0, -1], [1, 0]])
     1560            sage: Qr, Tr = A.schur(base_ring=RDF)
     1561            sage: Qc, Tc = A.schur(base_ring=CDF)
     1562            sage: all([M.is_immutable() for M in [Qr, Tr, Qc, Tc]])
     1563            True
     1564            sage: Tr.round(6) != Tc.round(6)
     1565            True
     1566
     1567        TESTS:
     1568
     1569        The Schur factorization is only defined for square matrices.  ::
     1570
     1571            sage: A = matrix(RDF, 4, 5, range(20))
     1572            sage: A.schur()
     1573            Traceback (most recent call last):
     1574            ...
     1575            ValueError: Schur decomposition requires a square matrix, not a 4 x 5 matrix
     1576
     1577        A base ring request is checked. ::
     1578
     1579            sage: A = matrix(RDF, 3, range(9))
     1580            sage: A.schur(base_ring=QQ)
     1581            Traceback (most recent call last):
     1582            ...
     1583            ValueError: base ring of Schur decomposition matrices must be RDF or CDF, not Rational Field
     1584        """
     1585        global scipy
     1586        from sage.rings.real_double import RDF
     1587        from sage.rings.complex_double import CDF
     1588
     1589        cdef Matrix_double_dense Q, T
     1590
     1591        if not self.is_square():
     1592            raise ValueError('Schur decomposition requires a square matrix, not a {0} x {1} matrix'.format(self.nrows(), self.ncols()))
     1593        if base_ring == None:
     1594            base_ring = self.base_ring()
     1595        if not base_ring in [RDF, CDF]:
     1596            raise ValueError('base ring of Schur decomposition matrices must be RDF or CDF, not {0}'.format(base_ring))
     1597
     1598        if self.base_ring() != base_ring:
     1599            try:
     1600                self = self.change_ring(base_ring)
     1601            except TypeError:
     1602                raise TypeError('unable to convert input matrix over CDF to a matrix over RDF')
     1603        if base_ring == CDF:
     1604            format = 'complex'
     1605        else:
     1606            format = 'real'
     1607
     1608        schur = self.fetch('schur_factors_' + format)
     1609        if not schur is None:
     1610            return schur
     1611        Q = self._new(self._nrows, self._nrows)
     1612        T = self._new(self._nrows, self._nrows)
     1613        if scipy is None:
     1614            import scipy
     1615        import scipy.linalg
     1616        T._matrix_numpy, Q._matrix_numpy = scipy.linalg.schur(self._matrix_numpy, output=format)
     1617        Q.set_immutable()
     1618        T.set_immutable()
     1619        # Our return order is the reverse of NumPy's
     1620        schur = (Q, T)
     1621        self.cache('schur_factors_' + format, schur)
     1622        return schur
    13371623
    13381624    def cholesky(self):
    13391625        r"""