Ticket #10848: trac_10848-hermitian-matrices-v5.patch

File trac_10848-hermitian-matrices-v5.patch, 11.0 KB (added by rbeezer, 9 years ago)
  • sage/matrix/matrix0.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1298574394 28800
    # Node ID 9e8fe224a52936a26d940d3012fa945e45a541a1
    # Parent  f884aba42d5ca2cce5f0ba196e4cf755a1424e79
    10848: add is_hermitian for matrices
    * * *
    [mq]: hermitian-two-speed
    
    diff -r f884aba42d5c -r 9e8fe224a529 sage/matrix/matrix0.pyx
    a b  
    28962896                    return False
    28972897        return True
    28982898
     2899    def is_hermitian(self):
     2900        r"""
     2901        Returns ``True`` if the matrix is equal to its conjugate-transpose.
     2902
     2903        OUTPUT:
     2904
     2905        ``True`` if the matrix is square and equal to the transpose
     2906        with every entry conjugated, and ``False`` otherwise.
     2907
     2908        Note that if conjugation has no effect on elements of the base
     2909        ring (such as for integers), then the :meth:`is_symmetric`
     2910        method is equivalent and faster.
     2911
     2912        This routine is for matrices over exact rings and so may not
     2913        work properly for matrices over ``RR`` or ``CC``.  For matrices with
     2914        approximate entries, the rings of double-precision floating-point
     2915        numbers, ``RDF`` and ``CDF``, are a better choice since the
     2916        :meth:`sage.matrix.matrix_double_dense.Matrix_double_dense.is_hermitian`
     2917        method has a tolerance parameter.  This provides control over
     2918        allowing for minor discrepancies between entries when checking
     2919        equality.
     2920
     2921        EXAMPLES::
     2922
     2923            sage: A = matrix(QQbar, [[ 1 + I,  1 - 6*I, -1 - I],
     2924            ...                      [-3 - I,     -4*I,     -2],
     2925            ...                      [-1 + I, -2 - 8*I,  2 + I]])
     2926            sage: A.is_hermitian()
     2927            False
     2928            sage: B = A*A.conjugate_transpose()
     2929            sage: B.is_hermitian()
     2930            True
     2931
     2932        Sage has several fields besides the entire complex numbers
     2933        where conjugation is non-trivial. ::
     2934
     2935            sage: F.<b> = QuadraticField(-7)
     2936            sage: C = matrix(F, [[-2*b - 3,  7*b - 6, -b + 3],
     2937            ...                  [-2*b - 3, -3*b + 2,   -2*b],
     2938            ...                  [   b + 1,        0,     -2]])
     2939            sage: C.is_hermitian()
     2940            False
     2941            sage: C = C*C.conjugate_transpose()
     2942            sage: C.is_hermitian()
     2943            True
     2944
     2945        A matrix that is nearly Hermitian, but for a non-real
     2946        diagonal entry. ::
     2947
     2948            sage: A = matrix(QQbar, [[    2,   2-I, 1+4*I],
     2949            ...                      [  2+I,   3+I, 2-6*I],
     2950            ...                      [1-4*I, 2+6*I,     5]])
     2951            sage: A.is_hermitian()
     2952            False
     2953            sage: A[1,1] = 132
     2954            sage: A.is_hermitian()
     2955            True
     2956
     2957        Rectangular matrices are never Hermitian.  ::
     2958
     2959            sage: A = matrix(QQbar, 3, 4)
     2960            sage: A.is_hermitian()
     2961            False
     2962        """
     2963        if not self.is_square():
     2964            return False
     2965
     2966        cdef Py_ssize_t i,j
     2967
     2968        for i from 0 <= i < self._nrows:
     2969            for j from 0 <= j <= i:
     2970                if self.get_unsafe(i,j) != self.get_unsafe(j,i).conjugate():
     2971                    return False
     2972        return True
     2973
    28992974    def is_skew_symmetric(self):
    29002975        """
    29012976        Returns True if this is a skew symmetric matrix.
  • sage/matrix/matrix_double_dense.pyx

    diff -r f884aba42d5c -r 9e8fe224a529 sage/matrix/matrix_double_dense.pyx
    a b  
    16581658        self.cache(key, b)
    16591659        return b
    16601660
     1661    def _is_lower_triangular(self, tol):
     1662        r"""
     1663        Returns ``True`` if the entries above the diagonal are all zero.
     1664
     1665        INPUT:
     1666
     1667        - ``tol`` -  the largest value of the absolute value of the
     1668          difference between two matrix entries for which they will
     1669          still be considered equal.
     1670
     1671        OUTPUT:
     1672
     1673        Returns ``True`` if each entry above the diagonal (row index
     1674        is less than the column index) is zero.
     1675
     1676        EXAMPLES::
     1677
     1678            sage: A = matrix(RDF, [[ 2.0, 0.0,  0.0],
     1679            ...                    [ 1.0, 3.0,  0.0],
     1680            ...                    [-4.0, 2.0, -1.0]])
     1681            sage: A._is_lower_triangular(1.0e-17)
     1682            True
     1683            sage: A[1,2] = 10^-13
     1684            sage: A._is_lower_triangular(1.0e-14)
     1685            False
     1686        """
     1687        global numpy
     1688        if numpy is None:
     1689            import numpy
     1690        cdef Py_ssize_t i, j
     1691        for i from 0 < i < self._nrows:
     1692            for j from i+1 <= j < self._ncols:
     1693                if numpy.absolute(self.get_unsafe(i,j)) > tol:
     1694                    return False
     1695        return True
     1696
     1697    def is_hermitian(self, tol = 1e-12, algorithm='orthonormal'):
     1698        r"""
     1699        Returns ``True`` if the matrix is equal to its conjugate-transpose.
     1700
     1701        INPUT:
     1702
     1703        - ``tol`` - default: ``1e-12`` - the largest value of the
     1704          absolute value of the difference between two matrix entries
     1705          for which they will still be considered equal.
     1706
     1707        - ``algorithm`` - default: 'orthonormal' - set to 'orthonormal'
     1708          for a stable procedure and set to 'naive' for a fast
     1709          procedure.
     1710
     1711        OUTPUT:
     1712
     1713        ``True`` if the matrix is square and equal to the transpose
     1714        with every entry conjugated, and ``False`` otherwise.
     1715
     1716        Note that if conjugation has no effect on elements of the base
     1717        ring (such as for integers), then the :meth:`is_symmetric`
     1718        method is equivalent and faster.
     1719
     1720        The tolerance parameter is used to allow for numerical values
     1721        to be equal if there is a slight difference due to round-off
     1722        and other imprecisions.
     1723
     1724        The result is cached, on a per-tolerance and per-algorithm basis.
     1725
     1726        ALGORITHMS::
     1727
     1728        The naive algorithm simply compares corresponing entries on either
     1729        side of the diagonal (and on the diagonal itself) to see if they are
     1730        conjugates, with equality controlled by the tolerance parameter.
     1731
     1732        The orthonormal algorithm first computes a Schur decomposition
     1733        (via the :meth:`schur` method) and checks that the result is a
     1734        diagonal matrix with real entries.
     1735
     1736        So the naive algorithm can finish quickly for a matrix that is not
     1737        Hermitian, while the orthonormal algorithm will always compute a
     1738        Schur decomposition before going through a similar check of the matrix
     1739        entry-by-entry.
     1740
     1741        EXAMPLES::
     1742
     1743            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
     1744            ...                    [-3 - I,     -4*I,     -2],
     1745            ...                    [-1 + I, -2 - 8*I,  2 + I]])
     1746            sage: A.is_hermitian(algorithm='orthonormal')
     1747            False
     1748            sage: A.is_hermitian(algorithm='naive')
     1749            False
     1750            sage: B = A*A.conjugate_transpose()
     1751            sage: B.is_hermitian(algorithm='orthonormal')
     1752            True
     1753            sage: B.is_hermitian(algorithm='naive')
     1754            True
     1755
     1756        A matrix that is nearly Hermitian, but for one non-real
     1757        diagonal entry. ::
     1758
     1759            sage: A = matrix(CDF, [[    2,   2-I, 1+4*I],
     1760            ...                    [  2+I,   3+I, 2-6*I],
     1761            ...                    [1-4*I, 2+6*I,     5]])
     1762            sage: A.is_hermitian(algorithm='orthonormal')
     1763            False
     1764            sage: A[1,1] = 132
     1765            sage: A.is_hermitian(algorithm='orthonormal')
     1766            True
     1767
     1768        We get a unitary matrix from the SVD routine and use this
     1769        numerical matrix to create a matrix that should be Hermitian
     1770        (indeed it should be the identity matrix), but with some
     1771        imprecision.  We use this to illustrate that if the tolerance
     1772        is set too small, then we can be too strict about the equality
     1773        of entries and achieve the wrong result.  ::
     1774
     1775            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
     1776            ...                    [-3 - I,     -4*I,     -2],
     1777            ...                    [-1 + I, -2 - 8*I,  2 + I]])
     1778            sage: U, _, _ = A.SVD()
     1779            sage: B=U*U.conjugate_transpose()
     1780            sage: B.is_hermitian(algorithm='naive')
     1781            True
     1782            sage: B.is_hermitian(algorithm='naive', tol=1.0e-17)
     1783            False
     1784            sage: B.is_hermitian(algorithm='naive', tol=1.0e-15)
     1785            True
     1786
     1787        Rectangular matrices are never Hermitian, no matter which
     1788        algorithm is requested.  ::
     1789
     1790            sage: A = matrix(CDF, 3, 4)
     1791            sage: A.is_hermitian()
     1792            False
     1793
     1794        TESTS:
     1795
     1796        The tolerance must be strictly positive.  ::
     1797
     1798            sage: A = matrix(RDF, 2, range(4))
     1799            sage: A.is_hermitian(tol = -3.1)
     1800            Traceback (most recent call last):
     1801            ...
     1802            ValueError: tolerance must be positive, not -3.1
     1803
     1804        The ``algorithm`` keyword gets checked.  ::
     1805
     1806            sage: A = matrix(RDF, 2, range(4))
     1807            sage: A.is_hermitian(algorithm='junk')
     1808            Traceback (most recent call last):
     1809            ...
     1810            ValueError: algorithm must be 'naive' or 'orthonormal', not junk
     1811
     1812        AUTHOR:
     1813
     1814         - Rob Beezer (2011-03-30)
     1815        """
     1816        import sage.rings.complex_double
     1817        global numpy
     1818        tol = float(tol)
     1819        if tol <= 0:
     1820            raise ValueError('tolerance must be positive, not {0}'.format(tol))
     1821        if not algorithm in ['naive', 'orthonormal']:
     1822            raise ValueError("algorithm must be 'naive' or 'orthonormal', not {0}".format(algorithm))
     1823
     1824        key = 'hermitian_{0}_{1}'.format(algorithm, tol)
     1825        b = self.fetch(key)
     1826        if not b is None:
     1827            return b
     1828        if not self.is_square():
     1829            self.cache(key, False)
     1830            return False
     1831        if numpy is None:
     1832            import numpy
     1833        cdef Py_ssize_t i, j
     1834        cdef Matrix_double_dense T
     1835        if algorithm == 'orthonormal':
     1836            # Schur decomposition over CDF will be diagonal and real iff Hermitian
     1837            _, T = self.schur(base_ring=sage.rings.complex_double.CDF)
     1838            hermitian = T._is_lower_triangular(tol)
     1839            if hermitian:
     1840                for i from 0 < i < T._nrows:
     1841                    if numpy.absolute(numpy.imag(T.get_unsafe(i,i))) > tol:
     1842                        hermitian = False
     1843                        break
     1844        elif algorithm == 'naive':
     1845            hermitian = True
     1846            for i from 0 < i < self._nrows:
     1847                for j from 0 <= j <= i:
     1848                    if numpy.absolute(self.get_unsafe(i,j) - self.get_unsafe(j,i).conjugate()) > tol:
     1849                        hermitian = False
     1850                        break
     1851                if not hermitian:
     1852                    break
     1853        self.cache(key, hermitian)
     1854        return hermitian
     1855
    16611856    def schur(self, base_ring=None):
    16621857        r"""
    16631858        Returns the Schur decomposition of the matrix.