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

File trac_10848-hermitian-matrices-v7.patch, 11.6 KB (added by rbeezer, 8 years ago)

Fixed commit message, rebased to 4.7.1.alpha4

  • sage/matrix/matrix0.pyx

    # HG changeset patch
    # User Rob Beezer <beezer@ups.edu>
    # Date 1298574394 28800
    # Node ID c810f8db0caf42000c0eef1c15707d28d15e1aa8
    # Parent  03cfa5a44bf039d2d81f1e044ee988d50e5d2411
    10848: add is_hermitian for matrices
    
    diff --git a/sage/matrix/matrix0.pyx b/sage/matrix/matrix0.pyx
    a b  
    30193019                    return False
    30203020        return True
    30213021
     3022    def is_hermitian(self):
     3023        r"""
     3024        Returns ``True`` if the matrix is equal to its conjugate-transpose.
     3025
     3026        OUTPUT:
     3027
     3028        ``True`` if the matrix is square and equal to the transpose
     3029        with every entry conjugated, and ``False`` otherwise.
     3030
     3031        Note that if conjugation has no effect on elements of the base
     3032        ring (such as for integers), then the :meth:`is_symmetric`
     3033        method is equivalent and faster.
     3034
     3035        This routine is for matrices over exact rings and so may not
     3036        work properly for matrices over ``RR`` or ``CC``.  For matrices with
     3037        approximate entries, the rings of double-precision floating-point
     3038        numbers, ``RDF`` and ``CDF``, are a better choice since the
     3039        :meth:`sage.matrix.matrix_double_dense.Matrix_double_dense.is_hermitian`
     3040        method has a tolerance parameter.  This provides control over
     3041        allowing for minor discrepancies between entries when checking
     3042        equality.
     3043
     3044        The result is cached.
     3045
     3046        EXAMPLES::
     3047
     3048            sage: A = matrix(QQbar, [[ 1 + I,  1 - 6*I, -1 - I],
     3049            ...                      [-3 - I,     -4*I,     -2],
     3050            ...                      [-1 + I, -2 - 8*I,  2 + I]])
     3051            sage: A.is_hermitian()
     3052            False
     3053            sage: B = A*A.conjugate_transpose()
     3054            sage: B.is_hermitian()
     3055            True
     3056
     3057        Sage has several fields besides the entire complex numbers
     3058        where conjugation is non-trivial. ::
     3059
     3060            sage: F.<b> = QuadraticField(-7)
     3061            sage: C = matrix(F, [[-2*b - 3,  7*b - 6, -b + 3],
     3062            ...                  [-2*b - 3, -3*b + 2,   -2*b],
     3063            ...                  [   b + 1,        0,     -2]])
     3064            sage: C.is_hermitian()
     3065            False
     3066            sage: C = C*C.conjugate_transpose()
     3067            sage: C.is_hermitian()
     3068            True
     3069
     3070        A matrix that is nearly Hermitian, but for a non-real
     3071        diagonal entry. ::
     3072
     3073            sage: A = matrix(QQbar, [[    2,   2-I, 1+4*I],
     3074            ...                      [  2+I,   3+I, 2-6*I],
     3075            ...                      [1-4*I, 2+6*I,     5]])
     3076            sage: A.is_hermitian()
     3077            False
     3078            sage: A[1,1] = 132
     3079            sage: A.is_hermitian()
     3080            True
     3081
     3082        Rectangular matrices are never Hermitian.  ::
     3083
     3084            sage: A = matrix(QQbar, 3, 4)
     3085            sage: A.is_hermitian()
     3086            False
     3087
     3088        A square, empty matrix is trivially Hermitian.  ::
     3089
     3090            sage: A = matrix(QQ, 0, 0)
     3091            sage: A.is_hermitian()
     3092            True
     3093        """
     3094        key = 'hermitian'
     3095        h = self.fetch(key)
     3096        if not h is None:
     3097            return h
     3098        if not self.is_square():
     3099            self.cache(key, False)
     3100            return False
     3101        if self._nrows == 0:
     3102            self.cache(key, True)
     3103            return True
     3104
     3105        cdef Py_ssize_t i,j
     3106        hermitian = True
     3107        for i in range(self._nrows):
     3108            for j in range(i+1):
     3109                if self.get_unsafe(i,j) != self.get_unsafe(j,i).conjugate():
     3110                    hermitian = False
     3111                    break
     3112            if not hermitian:
     3113                break
     3114        self.cache(key, hermitian)
     3115        return hermitian
     3116
    30223117    def is_skew_symmetric(self):
    30233118        """
    30243119        Returns True if this is a skew symmetric matrix.
  • sage/matrix/matrix_double_dense.pyx

    diff --git a/sage/matrix/matrix_double_dense.pyx b/sage/matrix/matrix_double_dense.pyx
    a b  
    14471447        self.cache(key, unitary)
    14481448        return unitary
    14491449
     1450    def _is_lower_triangular(self, tol):
     1451        r"""
     1452        Returns ``True`` if the entries above the diagonal are all zero.
     1453
     1454        INPUT:
     1455
     1456        - ``tol`` -  the largest value of the absolute value of the
     1457          difference between two matrix entries for which they will
     1458          still be considered equal.
     1459
     1460        OUTPUT:
     1461
     1462        Returns ``True`` if each entry above the diagonal (entries
     1463        with a row index less than the column index) is zero.
     1464
     1465        EXAMPLES::
     1466
     1467            sage: A = matrix(RDF, [[ 2.0, 0.0,  0.0],
     1468            ...                    [ 1.0, 3.0,  0.0],
     1469            ...                    [-4.0, 2.0, -1.0]])
     1470            sage: A._is_lower_triangular(1.0e-17)
     1471            True
     1472            sage: A[1,2] = 10^-13
     1473            sage: A._is_lower_triangular(1.0e-14)
     1474            False
     1475        """
     1476        global numpy
     1477        if numpy is None:
     1478            import numpy
     1479        cdef Py_ssize_t i, j
     1480        for i in range(self._nrows):
     1481            for j in range(i+1, self._ncols):
     1482                if numpy.absolute(self.get_unsafe(i,j)) > tol:
     1483                    return False
     1484        return True
     1485
     1486    def is_hermitian(self, tol = 1e-12, algorithm='orthonormal'):
     1487        r"""
     1488        Returns ``True`` if the matrix is equal to its conjugate-transpose.
     1489
     1490        INPUT:
     1491
     1492        - ``tol`` - default: ``1e-12`` - the largest value of the
     1493          absolute value of the difference between two matrix entries
     1494          for which they will still be considered equal.
     1495
     1496        - ``algorithm`` - default: 'orthonormal' - set to 'orthonormal'
     1497          for a stable procedure and set to 'naive' for a fast
     1498          procedure.
     1499
     1500        OUTPUT:
     1501
     1502        ``True`` if the matrix is square and equal to the transpose
     1503        with every entry conjugated, and ``False`` otherwise.
     1504
     1505        Note that if conjugation has no effect on elements of the base
     1506        ring (such as for integers), then the :meth:`is_symmetric`
     1507        method is equivalent and faster.
     1508
     1509        The tolerance parameter is used to allow for numerical values
     1510        to be equal if there is a slight difference due to round-off
     1511        and other imprecisions.
     1512
     1513        The result is cached, on a per-tolerance and per-algorithm basis.
     1514
     1515        ALGORITHMS::
     1516
     1517        The naive algorithm simply compares corresponing entries on either
     1518        side of the diagonal (and on the diagonal itself) to see if they are
     1519        conjugates, with equality controlled by the tolerance parameter.
     1520
     1521        The orthonormal algorithm first computes a Schur decomposition
     1522        (via the :meth:`schur` method) and checks that the result is a
     1523        diagonal matrix with real entries.
     1524
     1525        So the naive algorithm can finish quickly for a matrix that is not
     1526        Hermitian, while the orthonormal algorithm will always compute a
     1527        Schur decomposition before going through a similar check of the matrix
     1528        entry-by-entry.
     1529
     1530        EXAMPLES::
     1531
     1532            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
     1533            ...                    [-3 - I,     -4*I,     -2],
     1534            ...                    [-1 + I, -2 - 8*I,  2 + I]])
     1535            sage: A.is_hermitian(algorithm='orthonormal')
     1536            False
     1537            sage: A.is_hermitian(algorithm='naive')
     1538            False
     1539            sage: B = A*A.conjugate_transpose()
     1540            sage: B.is_hermitian(algorithm='orthonormal')
     1541            True
     1542            sage: B.is_hermitian(algorithm='naive')
     1543            True
     1544
     1545        A matrix that is nearly Hermitian, but for one non-real
     1546        diagonal entry. ::
     1547
     1548            sage: A = matrix(CDF, [[    2,   2-I, 1+4*I],
     1549            ...                    [  2+I,   3+I, 2-6*I],
     1550            ...                    [1-4*I, 2+6*I,     5]])
     1551            sage: A.is_hermitian(algorithm='orthonormal')
     1552            False
     1553            sage: A[1,1] = 132
     1554            sage: A.is_hermitian(algorithm='orthonormal')
     1555            True
     1556
     1557        We get a unitary matrix from the SVD routine and use this
     1558        numerical matrix to create a matrix that should be Hermitian
     1559        (indeed it should be the identity matrix), but with some
     1560        imprecision.  We use this to illustrate that if the tolerance
     1561        is set too small, then we can be too strict about the equality
     1562        of entries and achieve the wrong result.  ::
     1563
     1564            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
     1565            ...                    [-3 - I,     -4*I,     -2],
     1566            ...                    [-1 + I, -2 - 8*I,  2 + I]])
     1567            sage: U, _, _ = A.SVD()
     1568            sage: B=U*U.conjugate_transpose()
     1569            sage: B.is_hermitian(algorithm='naive')
     1570            True
     1571            sage: B.is_hermitian(algorithm='naive', tol=1.0e-17)
     1572            False
     1573            sage: B.is_hermitian(algorithm='naive', tol=1.0e-15)
     1574            True
     1575
     1576        A square, empty matrix is trivially Hermitian.  ::
     1577
     1578            sage: A = matrix(RDF, 0, 0)
     1579            sage: A.is_hermitian()
     1580            True
     1581
     1582        Rectangular matrices are never Hermitian, no matter which
     1583        algorithm is requested.  ::
     1584
     1585            sage: A = matrix(CDF, 3, 4)
     1586            sage: A.is_hermitian()
     1587            False
     1588
     1589        TESTS:
     1590
     1591        The tolerance must be strictly positive.  ::
     1592
     1593            sage: A = matrix(RDF, 2, range(4))
     1594            sage: A.is_hermitian(tol = -3.1)
     1595            Traceback (most recent call last):
     1596            ...
     1597            ValueError: tolerance must be positive, not -3.1
     1598
     1599        The ``algorithm`` keyword gets checked.  ::
     1600
     1601            sage: A = matrix(RDF, 2, range(4))
     1602            sage: A.is_hermitian(algorithm='junk')
     1603            Traceback (most recent call last):
     1604            ...
     1605            ValueError: algorithm must be 'naive' or 'orthonormal', not junk
     1606
     1607        AUTHOR:
     1608
     1609         - Rob Beezer (2011-03-30)
     1610        """
     1611        import sage.rings.complex_double
     1612        global numpy
     1613        tol = float(tol)
     1614        if tol <= 0:
     1615            raise ValueError('tolerance must be positive, not {0}'.format(tol))
     1616        if not algorithm in ['naive', 'orthonormal']:
     1617            raise ValueError("algorithm must be 'naive' or 'orthonormal', not {0}".format(algorithm))
     1618
     1619        key = 'hermitian_{0}_{1}'.format(algorithm, tol)
     1620        h = self.fetch(key)
     1621        if not h is None:
     1622            return h
     1623        if not self.is_square():
     1624            self.cache(key, False)
     1625            return False
     1626        if self._nrows == 0:
     1627            self.cache(key, True)
     1628            return True
     1629        if numpy is None:
     1630            import numpy
     1631        cdef Py_ssize_t i, j
     1632        cdef Matrix_double_dense T
     1633        if algorithm == 'orthonormal':
     1634            # Schur decomposition over CDF will be diagonal and real iff Hermitian
     1635            _, T = self.schur(base_ring=sage.rings.complex_double.CDF)
     1636            hermitian = T._is_lower_triangular(tol)
     1637            if hermitian:
     1638                for i in range(T._nrows):
     1639                    if numpy.absolute(numpy.imag(T.get_unsafe(i,i))) > tol:
     1640                        hermitian = False
     1641                        break
     1642        elif algorithm == 'naive':
     1643            hermitian = True
     1644            for i in range(self._nrows):
     1645                for j in range(i+1):
     1646                    if numpy.absolute(self.get_unsafe(i,j) - self.get_unsafe(j,i).conjugate()) > tol:
     1647                        hermitian = False
     1648                        break
     1649                if not hermitian:
     1650                    break
     1651        self.cache(key, hermitian)
     1652        return hermitian
     1653
    14501654    def schur(self, base_ring=None):
    14511655        r"""
    14521656        Returns the Schur decomposition of the matrix.