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

File trac_11027-schur-decomposition-v3.patch, 12.9 KB (added by jdemeyer, 10 years ago)

Patch rebased to sage-4.7

  • sage/matrix/matrix_double_dense.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1301065587 25200
    # Node ID 24d52b1990009edf1cd157336c5a3c6dfa18059e
    # Parent  611e83c50dc9bb56a04a52ae752b118a36fd0213
    11027: Schur matrix decomposition
    
    diff -r 611e83c50dc9 -r 24d52b199000 sage/matrix/matrix_double_dense.pyx
    a b  
    14471447        self.cache(key, unitary)
    14481448        return unitary
    14491449
     1450    def schur(self, base_ring=None):
     1451        r"""
     1452        Returns the Schur decomposition of the matrix.
     1453
     1454        INPUT:
     1455
     1456        - ``base_ring`` - optional, defaults to the base ring of ``self``.
     1457          Use this to request the base ring of the returned matrices, which
     1458          will affect the format of the results.
     1459
     1460        OUTPUT:
     1461
     1462        A pair of immutable matrices.  The first is a unitary matrix `Q`.
     1463        The second, `T`, is upper-triangular when returned over the complex
     1464        numbers, while it is almost upper-triangular over the reals.  In the
     1465        latter case, there can be some `2\times 2` blocks on the diagonal
     1466        which represent a pair of conjugate complex eigenvalues of ``self``.
     1467
     1468        If ``self`` is the matrix `A`, then
     1469
     1470        .. math::
     1471
     1472            A = QT({\overline Q})^t
     1473
     1474        where the latter matrix is the conjugate-transpose of ``Q``, which
     1475        is also the inverse of ``Q``, since ``Q`` is unitary.
     1476
     1477        Note that in the case of a normal matrix (Hermitian, symmetric, and
     1478        others), the upper-triangular matrix is  a diagonal matrix with
     1479        eigenvalues of ``self`` on the diagonal, and the unitary matrix
     1480        has columns that form an orthonormal basis composed of eigenvectors
     1481        of ``self``.  This is known as "orthonormal diagonalization".
     1482
     1483        .. warning::
     1484
     1485            The Schur decomposition is not unique, as there may be numerous
     1486            choices for the vectors of the orthonormal basis, and consequently
     1487            different possibilities for the upper-triangular matrix.  However,
     1488            the diagonal of the upper-triangular matrix will always contain the
     1489            eigenvalues of the matrix (in the complex version), or `2\times 2`
     1490            block matrices in the real version representing pairs of conjugate
     1491            complex eigenvalues.
     1492
     1493            In particular, results may vary across systems and processors.
     1494
     1495        EXAMPLES:
     1496
     1497        First over the complexes.  The similar matrix is always
     1498        upper-triangular in this case.  ::
     1499
     1500            sage: A = matrix(CDF, 4, 4, range(16)) + matrix(CDF, 4, 4, [x^3*I for x in range(0, 16)])
     1501            sage: Q, T = A.schur()
     1502            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1503            [1.0   0   0   0]
     1504            [  0 1.0   0   0]
     1505            [  0   0 1.0   0]
     1506            [  0   0   0 1.0]
     1507            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1508            True
     1509            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1510            [0 0 0 0]
     1511            [0 0 0 0]
     1512            [0 0 0 0]
     1513            [0 0 0 0]
     1514            sage: eigenvalues = [T[i,i] for i in range(4)]; eigenvalues
     1515            [30.733... + 4648.541...*I, -0.184... - 159.057...*I, -0.523... + 11.158...*I, -0.025... - 0.642...*I]
     1516            sage: A.eigenvalues()
     1517            [30.733... + 4648.541...*I, -0.184... - 159.057...*I, -0.523... + 11.158...*I, -0.025... - 0.642...*I]
     1518            sage: abs(A.norm()-T.norm()) < 1e-10
     1519            True
     1520
     1521        We begin with a real matrix but ask for a decomposition over the
     1522        complexes.  The result will yield an upper-triangular matrix over
     1523        the complex numbers for ``T``. ::
     1524
     1525            sage: A = matrix(RDF, 4, 4, [x^3 for x in range(16)])
     1526            sage: Q, T = A.schur(base_ring=CDF)
     1527            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1528            [1.0   0   0   0]
     1529            [  0 1.0   0   0]
     1530            [  0   0 1.0   0]
     1531            [  0   0   0 1.0]
     1532            sage: T.parent()
     1533            Full MatrixSpace of 4 by 4 dense matrices over Complex Double Field
     1534            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1535            True
     1536            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1537            [0 0 0 0]
     1538            [0 0 0 0]
     1539            [0 0 0 0]
     1540            [0 0 0 0]
     1541
     1542        Now totally over the reals.  But with complex eigenvalues, the
     1543        similar matrix may not be upper-triangular. But "at worst" there
     1544        may be some `2\times 2` blocks on the diagonal which represent
     1545        a pair of conjugate complex eigenvalues. These blocks will then
     1546        just interrupt the zeros below the main diagonal.  This example
     1547        has a pair of these of the blocks. ::
     1548
     1549            sage: A = matrix(RDF, 4, 4, [[1, 0, -3, -1],
     1550            ...                          [4, -16, -7, 0],
     1551            ...                          [1, 21, 1, -2],
     1552            ...                          [26, -1, -2, 1]])
     1553            sage: Q, T = A.schur()
     1554            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1555            [1.0 0.0 0.0 0.0]
     1556            [0.0 1.0 0.0 0.0]
     1557            [0.0 0.0 1.0 0.0]
     1558            [0.0 0.0 0.0 1.0]
     1559            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)])
     1560            False
     1561            sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i-1)])
     1562            True
     1563            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1564            [0.0 0.0 0.0 0.0]
     1565            [0.0 0.0 0.0 0.0]
     1566            [0.0 0.0 0.0 0.0]
     1567            [0.0 0.0 0.0 0.0]
     1568            sage: eigenvalues = sorted(T[0:2,0:2].eigenvalues() + T[2:4,2:4].eigenvalues())
     1569            sage: eigenvalues.reverse(); eigenvalues
     1570            [-0.789... + 2.336...*I, -0.789... - 2.336...*I, -5.710... + 8.382...*I, -5.710... - 8.382...*I]
     1571            sage: A.eigenvalues()
     1572            [-0.789... + 2.336...*I, -0.789... - 2.336...*I, -5.710... + 8.382...*I, -5.710... - 8.382...*I]
     1573            sage: abs(A.norm()-T.norm()) < 1e-12
     1574            True
     1575
     1576        Starting with complex numbers and requesting a result over the reals
     1577        will never happen.  ::
     1578
     1579            sage: A = matrix(CDF, 2, 2, [[2+I, -1+3*I], [5-4*I, 2-7*I]])
     1580            sage: A.schur(base_ring=RDF)
     1581            Traceback (most recent call last):
     1582            ...
     1583            TypeError: unable to convert input matrix over CDF to a matrix over RDF
     1584
     1585        If theory predicts your matrix is real, but it contains some
     1586        very small imaginary parts, you can specify the cutoff for "small"
     1587        imaginary parts, then request the output as real matrices, and let
     1588        the routine do the rest. ::
     1589
     1590            sage: A = matrix(RDF, 2, 2, [1, 1, -1, 0]) + matrix(CDF, 2, 2, [1.0e-14*I]*4)
     1591            sage: B = A.zero_at(1.0e-12)
     1592            sage: B.parent()
     1593            Full MatrixSpace of 2 by 2 dense matrices over Complex Double Field
     1594            sage: Q, T = B.schur(RDF)
     1595            sage: Q.parent()
     1596            Full MatrixSpace of 2 by 2 dense matrices over Real Double Field
     1597            sage: T.parent()
     1598            Full MatrixSpace of 2 by 2 dense matrices over Real Double Field
     1599            sage: Q.round(6)
     1600            [ 0.707107  0.707107]
     1601            [-0.707107  0.707107]
     1602            sage: T.round(6)
     1603            [ 0.5  1.5]
     1604            [-0.5  0.5]
     1605            sage: (Q*T*Q.conjugate().transpose()-B).zero_at(1.0e-11)
     1606            [0 0]
     1607            [0 0]
     1608
     1609        A Hermitian matrix has real eigenvalues, so the similar matrix
     1610        will be upper-triangular.  Furthermore, a Hermitian matrix is
     1611        diagonalizable with respect to an orthonormal basis, composed
     1612        of eigenvectors of the matrix.  Here that basis is the set of
     1613        columns of the unitary matrix.  ::
     1614
     1615            sage: A = matrix(CDF, [[        52,   -9*I - 8,    6*I - 187,  -188*I + 2],
     1616            ...                    [   9*I - 8,         12,   -58*I + 59,   30*I + 42],
     1617            ...                    [-6*I - 187,  58*I + 59,         2677, 2264*I + 65],
     1618            ...                    [ 188*I + 2, -30*I + 42, -2264*I + 65,       2080]])
     1619            sage: Q, T = A.schur()
     1620            sage: T = T.zero_at(1.0e-12).change_ring(RDF)
     1621            sage: T.round(6)
     1622            [4680.13301        0.0        0.0        0.0]
     1623            [       0.0 102.715967        0.0        0.0]
     1624            [       0.0        0.0  35.039344        0.0]
     1625            [       0.0        0.0        0.0    3.11168]
     1626            sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12)
     1627            [1.0   0   0   0]
     1628            [  0 1.0   0   0]
     1629            [  0   0 1.0   0]
     1630            [  0   0   0 1.0]
     1631            sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11)
     1632            [0 0 0 0]
     1633            [0 0 0 0]
     1634            [0 0 0 0]
     1635            [0 0 0 0]
     1636
     1637        Similarly, a real symmetric matrix has only real eigenvalues,
     1638        and there is an orthonormal basis composed of eigenvectors of
     1639        the matrix.  ::
     1640
     1641            sage: A = matrix(RDF, [[ 1, -2, 5, -3],
     1642            ...                    [-2,  9, 1,  5],
     1643            ...                    [ 5,  1, 3 , 7],
     1644            ...                    [-3,  5, 7, -8]])
     1645            sage: Q, T = A.schur()
     1646            sage: Q.round(4)
     1647            [-0.3027  -0.751   0.576 -0.1121]
     1648            [  0.139 -0.3892 -0.2648  0.8713]
     1649            [ 0.4361   0.359  0.7599  0.3217]
     1650            [ -0.836  0.3945  0.1438  0.3533]
     1651            sage: T.round(4)
     1652            [-13.5698      0.0      0.0      0.0]
     1653            [     0.0  -0.8508     -0.0     -0.0]
     1654            [     0.0      0.0   7.7664      0.0]
     1655            [     0.0      0.0      0.0  11.6542]
     1656            sage: (Q*Q.transpose()).zero_at(1.0e-12)
     1657            [1.0 0.0 0.0 0.0]
     1658            [0.0 1.0 0.0 0.0]
     1659            [0.0 0.0 1.0 0.0]
     1660            [0.0 0.0 0.0 1.0]
     1661            sage: (Q*T*Q.transpose()-A).zero_at(1.0e-11)
     1662            [0.0 0.0 0.0 0.0]
     1663            [0.0 0.0 0.0 0.0]
     1664            [0.0 0.0 0.0 0.0]
     1665            [0.0 0.0 0.0 0.0]
     1666
     1667        The results are cached, both as a real factorization and also as a
     1668        complex factorization.  This means the returned matrices are
     1669        immutable.  ::
     1670
     1671            sage: A = matrix(RDF, 2, 2, [[0, -1], [1, 0]])
     1672            sage: Qr, Tr = A.schur(base_ring=RDF)
     1673            sage: Qc, Tc = A.schur(base_ring=CDF)
     1674            sage: all([M.is_immutable() for M in [Qr, Tr, Qc, Tc]])
     1675            True
     1676            sage: Tr.round(6) != Tc.round(6)
     1677            True
     1678
     1679        TESTS:
     1680
     1681        The Schur factorization is only defined for square matrices.  ::
     1682
     1683            sage: A = matrix(RDF, 4, 5, range(20))
     1684            sage: A.schur()
     1685            Traceback (most recent call last):
     1686            ...
     1687            ValueError: Schur decomposition requires a square matrix, not a 4 x 5 matrix
     1688
     1689        A base ring request is checked. ::
     1690
     1691            sage: A = matrix(RDF, 3, range(9))
     1692            sage: A.schur(base_ring=QQ)
     1693            Traceback (most recent call last):
     1694            ...
     1695            ValueError: base ring of Schur decomposition matrices must be RDF or CDF, not Rational Field
     1696
     1697        AUTHOR:
     1698
     1699            - Rob Beezer (2011-03-31)
     1700        """
     1701        global scipy
     1702        from sage.rings.real_double import RDF
     1703        from sage.rings.complex_double import CDF
     1704
     1705        cdef Matrix_double_dense Q, T
     1706
     1707        if not self.is_square():
     1708            raise ValueError('Schur decomposition requires a square matrix, not a {0} x {1} matrix'.format(self.nrows(), self.ncols()))
     1709        if base_ring == None:
     1710            base_ring = self.base_ring()
     1711        if not base_ring in [RDF, CDF]:
     1712            raise ValueError('base ring of Schur decomposition matrices must be RDF or CDF, not {0}'.format(base_ring))
     1713
     1714        if self.base_ring() != base_ring:
     1715            try:
     1716                self = self.change_ring(base_ring)
     1717            except TypeError:
     1718                raise TypeError('unable to convert input matrix over CDF to a matrix over RDF')
     1719        if base_ring == CDF:
     1720            format = 'complex'
     1721        else:
     1722            format = 'real'
     1723
     1724        schur = self.fetch('schur_factors_' + format)
     1725        if not schur is None:
     1726            return schur
     1727        if scipy is None:
     1728            import scipy
     1729        import scipy.linalg
     1730        Q = self._new(self._nrows, self._nrows)
     1731        T = self._new(self._nrows, self._nrows)
     1732        T._matrix_numpy, Q._matrix_numpy = scipy.linalg.schur(self._matrix_numpy, output=format)
     1733        Q.set_immutable()
     1734        T.set_immutable()
     1735        # Our return order is the reverse of NumPy's
     1736        schur = (Q, T)
     1737        self.cache('schur_factors_' + format, schur)
     1738        return schur
     1739
    14501740    def cholesky(self):
    14511741        r"""
    14521742        Return the cholesky factorization of this matrix.  The input