Ticket #11027: trac_11027-schur-decomposition-v2.patch

File trac_11027-schur-decomposition-v2.patch, 12.8 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 55cd048d7689e648c728df31a2369a23cac34350
    # Parent  6a679959b54b7975af1841df4dd77d28b0faef28
    11027: Schur matrix decomposition
    
    diff -r 6a679959b54b -r 55cd048d7689 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 ``Q`` 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        .. warning::
     1371
     1372            The Schur decomposition is not unique, as there may be numerous
     1373            choices for the vectors of the orthonormal basis, and consequently
     1374            different possibilities for the upper-triangular matrix.  However,
     1375            the diagonal of the upper-triangular matrix will always contain the
     1376            eigenvalues of the matrix (in the complex version), or `2\times 2`
     1377            block matrices in the real version representing pairs of conjugate
     1378            complex eigenvalues.
     1379
     1380            In particular, results may vary across systems and processors.
     1381
     1382        EXAMPLES:
     1383
     1384        First over the complexes.  The similar matrix is always
     1385        upper-triangular in this case.  ::
     1386
     1387            sage: A = matrix(CDF, 4, 4, range(16)) + matrix(CDF, 4, 4, [x^3*I for x in range(0, 16)])
     1388            sage: Q, T = A.schur()
     1389            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1390            [1.0   0   0   0]
     1391            [  0 1.0   0   0]
     1392            [  0   0 1.0   0]
     1393            [  0   0   0 1.0]
     1394            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1395            True
     1396            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1397            [0 0 0 0]
     1398            [0 0 0 0]
     1399            [0 0 0 0]
     1400            [0 0 0 0]
     1401            sage: eigenvalues = [T[i,i] for i in range(4)]; eigenvalues
     1402            [30.733... + 4648.541...*I, -0.184... - 159.057...*I, -0.523... + 11.158...*I, -0.025... - 0.642...*I]
     1403            sage: A.eigenvalues()
     1404            [30.733... + 4648.541...*I, -0.184... - 159.057...*I, -0.523... + 11.158...*I, -0.025... - 0.642...*I]
     1405            sage: abs(A.norm()-T.norm()) < 1e-10
     1406            True
     1407
     1408        We begin with a real matrix but ask for a decomposition over the
     1409        complexes.  The result will yield an upper-triangular matrix over
     1410        the complex numbers for ``T``. ::
     1411
     1412            sage: A = matrix(RDF, 4, 4, [x^3 for x in range(16)])
     1413            sage: Q, T = A.schur(base_ring=CDF)
     1414            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1415            [1.0   0   0   0]
     1416            [  0 1.0   0   0]
     1417            [  0   0 1.0   0]
     1418            [  0   0   0 1.0]
     1419            sage: T.parent()
     1420            Full MatrixSpace of 4 by 4 dense matrices over Complex Double Field
     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.  This example
     1434        has a pair of these of the blocks. ::
     1435
     1436            sage: A = matrix(RDF, 4, 4, [[1, 0, -3, -1],
     1437            ...                          [4, -16, -7, 0],
     1438            ...                          [1, 21, 1, -2],
     1439            ...                          [26, -1, -2, 1]])
     1440            sage: Q, T = A.schur()
     1441            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1442            [1.0 0.0 0.0 0.0]
     1443            [0.0 1.0 0.0 0.0]
     1444            [0.0 0.0 1.0 0.0]
     1445            [0.0 0.0 0.0 1.0]
     1446            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1447            False
     1448            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i-1)])
     1449            True
     1450            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1451            [0.0 0.0 0.0 0.0]
     1452            [0.0 0.0 0.0 0.0]
     1453            [0.0 0.0 0.0 0.0]
     1454            [0.0 0.0 0.0 0.0]
     1455            sage: eigenvalues = sorted(T[0:2,0:2].eigenvalues() + T[2:4,2:4].eigenvalues())
     1456            sage: eigenvalues.reverse(); eigenvalues
     1457            [-0.789... + 2.336...*I, -0.789... - 2.336...*I, -5.710... + 8.382...*I, -5.710... - 8.382...*I]
     1458            sage: A.eigenvalues()
     1459            [-0.789... + 2.336...*I, -0.789... - 2.336...*I, -5.710... + 8.382...*I, -5.710... - 8.382...*I]
     1460            sage: abs(A.norm()-T.norm()) < 1e-12
     1461            True
     1462
     1463        Starting with complex numbers and requesting a result over the reals
     1464        will never happen.  ::
     1465
     1466            sage: A = matrix(CDF, 2, 2, [[2+I, -1+3*I], [5-4*I, 2-7*I]])
     1467            sage: A.schur(base_ring=RDF)
     1468            Traceback (most recent call last):
     1469            ...
     1470            TypeError: unable to convert input matrix over CDF to a matrix over RDF
     1471
     1472        If theory predicts your matrix is real, but it contains some
     1473        very small imaginary parts, you can specify the cutoff for "small"
     1474        imaginary parts, then request the output as real matrices, and let
     1475        the routine do the rest. ::
     1476
     1477            sage: A = matrix(RDF, 2, 2, [1, 1, -1, 0]) + matrix(CDF, 2, 2, [1.0e-14*I]*4)
     1478            sage: B = A.zero_at(1.0e-12)
     1479            sage: B.parent()
     1480            Full MatrixSpace of 2 by 2 dense matrices over Complex Double Field
     1481            sage: Q, T = B.schur(RDF)
     1482            sage: Q.parent()
     1483            Full MatrixSpace of 2 by 2 dense matrices over Real Double Field
     1484            sage: T.parent()
     1485            Full MatrixSpace of 2 by 2 dense matrices over Real Double Field
     1486            sage: Q.round(6)
     1487            [ 0.707107  0.707107]
     1488            [-0.707107  0.707107]
     1489            sage: T.round(6)
     1490            [ 0.5  1.5]
     1491            [-0.5  0.5]
     1492            sage: (Q*T*Q.conjugate().transpose()-B).zero_at(1.0e-11)
     1493            [0 0]
     1494            [0 0]
     1495
     1496        A Hermitian matrix has real eigenvalues, so the similar matrix
     1497        will be upper-triangular.  Furthermore, a Hermitian matrix is
     1498        diagonalizable with respect to an orthonormal basis, composed
     1499        of eigenvectors of the matrix.  Here that basis is the set of
     1500        columns of the unitary matrix.  ::
     1501
     1502            sage: A = matrix(CDF, [[        52,   -9*I - 8,    6*I - 187,  -188*I + 2],
     1503            ...                    [   9*I - 8,         12,   -58*I + 59,   30*I + 42],
     1504            ...                    [-6*I - 187,  58*I + 59,         2677, 2264*I + 65],
     1505            ...                    [ 188*I + 2, -30*I + 42, -2264*I + 65,       2080]])
     1506            sage: Q, T = A.schur()
     1507            sage: T = T.zero_at(1.0e-12).change_ring(RDF)
     1508            sage: T.round(6)
     1509            [4680.13301        0.0        0.0        0.0]
     1510            [       0.0 102.715967        0.0        0.0]
     1511            [       0.0        0.0  35.039344        0.0]
     1512            [       0.0        0.0        0.0    3.11168]
     1513            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1514            [1.0   0   0   0]
     1515            [  0 1.0   0   0]
     1516            [  0   0 1.0   0]
     1517            [  0   0   0 1.0]
     1518            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1519            [0 0 0 0]
     1520            [0 0 0 0]
     1521            [0 0 0 0]
     1522            [0 0 0 0]
     1523
     1524        Similarly, a real symmetric matrix has only real eigenvalues,
     1525        and there is an orthonormal basis composed of eigenvectors of
     1526        the matrix.  ::
     1527
     1528            sage: A = matrix(RDF, [[ 1, -2, 5, -3],
     1529            ...                    [-2,  9, 1,  5],
     1530            ...                    [ 5,  1, 3 , 7],
     1531            ...                    [-3,  5, 7, -8]])
     1532            sage: Q, T = A.schur()
     1533            sage: Q.round(4)
     1534            [-0.3027  -0.751   0.576 -0.1121]
     1535            [  0.139 -0.3892 -0.2648  0.8713]
     1536            [ 0.4361   0.359  0.7599  0.3217]
     1537            [ -0.836  0.3945  0.1438  0.3533]
     1538            sage: T.round(4)
     1539            [-13.5698      0.0      0.0      0.0]
     1540            [     0.0  -0.8508     -0.0     -0.0]
     1541            [     0.0      0.0   7.7664      0.0]
     1542            [     0.0      0.0      0.0  11.6542]
     1543            sage: (Q*Q.transpose()).zero_at(1.0e-12)
     1544            [1.0 0.0 0.0 0.0]
     1545            [0.0 1.0 0.0 0.0]
     1546            [0.0 0.0 1.0 0.0]
     1547            [0.0 0.0 0.0 1.0]
     1548            sage: (Q*T*Q.transpose()-A).zero_at(1.0e-11)
     1549            [0.0 0.0 0.0 0.0]
     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
     1554        The results are cached, both as a real factorization and also as a
     1555        complex factorization.  This means the returned matrices are
     1556        immutable.  ::
     1557
     1558            sage: A = matrix(RDF, 2, 2, [[0, -1], [1, 0]])
     1559            sage: Qr, Tr = A.schur(base_ring=RDF)
     1560            sage: Qc, Tc = A.schur(base_ring=CDF)
     1561            sage: all([M.is_immutable() for M in [Qr, Tr, Qc, Tc]])
     1562            True
     1563            sage: Tr.round(6) != Tc.round(6)
     1564            True
     1565
     1566        TESTS:
     1567
     1568        The Schur factorization is only defined for square matrices.  ::
     1569
     1570            sage: A = matrix(RDF, 4, 5, range(20))
     1571            sage: A.schur()
     1572            Traceback (most recent call last):
     1573            ...
     1574            ValueError: Schur decomposition requires a square matrix, not a 4 x 5 matrix
     1575
     1576        A base ring request is checked. ::
     1577
     1578            sage: A = matrix(RDF, 3, range(9))
     1579            sage: A.schur(base_ring=QQ)
     1580            Traceback (most recent call last):
     1581            ...
     1582            ValueError: base ring of Schur decomposition matrices must be RDF or CDF, not Rational Field
     1583
     1584        AUTHOR:
     1585
     1586            - Rob Beezer (2011-03-31)
     1587        """
     1588        global scipy
     1589        from sage.rings.real_double import RDF
     1590        from sage.rings.complex_double import CDF
     1591
     1592        cdef Matrix_double_dense Q, T
     1593
     1594        if not self.is_square():
     1595            raise ValueError('Schur decomposition requires a square matrix, not a {0} x {1} matrix'.format(self.nrows(), self.ncols()))
     1596        if base_ring == None:
     1597            base_ring = self.base_ring()
     1598        if not base_ring in [RDF, CDF]:
     1599            raise ValueError('base ring of Schur decomposition matrices must be RDF or CDF, not {0}'.format(base_ring))
     1600
     1601        if self.base_ring() != base_ring:
     1602            try:
     1603                self = self.change_ring(base_ring)
     1604            except TypeError:
     1605                raise TypeError('unable to convert input matrix over CDF to a matrix over RDF')
     1606        if base_ring == CDF:
     1607            format = 'complex'
     1608        else:
     1609            format = 'real'
     1610
     1611        schur = self.fetch('schur_factors_' + format)
     1612        if not schur is None:
     1613            return schur
     1614        if scipy is None:
     1615            import scipy
     1616        import scipy.linalg
     1617        Q = self._new(self._nrows, self._nrows)
     1618        T = self._new(self._nrows, self._nrows)
     1619        T._matrix_numpy, Q._matrix_numpy = scipy.linalg.schur(self._matrix_numpy, output=format)
     1620        Q.set_immutable()
     1621        T.set_immutable()
     1622        # Our return order is the reverse of NumPy's
     1623        schur = (Q, T)
     1624        self.cache('schur_factors_' + format, schur)
     1625        return schur
    13371626
    13381627    def cholesky(self):
    13391628        r"""