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

File trac_10848-hermitian-matrices-v6.patch, 11.8 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 1e9dd27c0c50524af699c534538f9977df03c5f8
    # Parent  55cd048d7689e648c728df31a2369a23cac34350
    10848: add is_hermitian for matrices
    * * *
    [mq]: hermitian-two-speed
    
    diff -r 55cd048d7689 -r 1e9dd27c0c50 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        The result is cached.
     2922
     2923        EXAMPLES::
     2924
     2925            sage: A = matrix(QQbar, [[ 1 + I,  1 - 6*I, -1 - I],
     2926            ...                      [-3 - I,     -4*I,     -2],
     2927            ...                      [-1 + I, -2 - 8*I,  2 + I]])
     2928            sage: A.is_hermitian()
     2929            False
     2930            sage: B = A*A.conjugate_transpose()
     2931            sage: B.is_hermitian()
     2932            True
     2933
     2934        Sage has several fields besides the entire complex numbers
     2935        where conjugation is non-trivial. ::
     2936
     2937            sage: F.<b> = QuadraticField(-7)
     2938            sage: C = matrix(F, [[-2*b - 3,  7*b - 6, -b + 3],
     2939            ...                  [-2*b - 3, -3*b + 2,   -2*b],
     2940            ...                  [   b + 1,        0,     -2]])
     2941            sage: C.is_hermitian()
     2942            False
     2943            sage: C = C*C.conjugate_transpose()
     2944            sage: C.is_hermitian()
     2945            True
     2946
     2947        A matrix that is nearly Hermitian, but for a non-real
     2948        diagonal entry. ::
     2949
     2950            sage: A = matrix(QQbar, [[    2,   2-I, 1+4*I],
     2951            ...                      [  2+I,   3+I, 2-6*I],
     2952            ...                      [1-4*I, 2+6*I,     5]])
     2953            sage: A.is_hermitian()
     2954            False
     2955            sage: A[1,1] = 132
     2956            sage: A.is_hermitian()
     2957            True
     2958
     2959        Rectangular matrices are never Hermitian.  ::
     2960
     2961            sage: A = matrix(QQbar, 3, 4)
     2962            sage: A.is_hermitian()
     2963            False
     2964
     2965        A square, empty matrix is trivially Hermitian.  ::
     2966
     2967            sage: A = matrix(QQ, 0, 0)
     2968            sage: A.is_hermitian()
     2969            True
     2970        """
     2971        key = 'hermitian'
     2972        h = self.fetch(key)
     2973        if not h is None:
     2974            return h
     2975        if not self.is_square():
     2976            self.cache(key, False)
     2977            return False
     2978        if self._nrows == 0:
     2979            self.cache(key, True)
     2980            return True
     2981
     2982        cdef Py_ssize_t i,j
     2983        hermitian = True
     2984        for i in range(self._nrows):
     2985            for j in range(i+1):
     2986                if self.get_unsafe(i,j) != self.get_unsafe(j,i).conjugate():
     2987                    hermitian = False
     2988                    break
     2989            if not hermitian:
     2990                break
     2991        self.cache(key, hermitian)
     2992        return hermitian
     2993
    28992994    def is_skew_symmetric(self):
    29002995        """
    29012996        Returns True if this is a skew symmetric matrix.
  • sage/matrix/matrix_double_dense.pyx

    diff -r 55cd048d7689 -r 1e9dd27c0c50 sage/matrix/matrix_double_dense.pyx
    a b  
    13341334        self.cache(key, b)
    13351335        return b
    13361336
     1337    def _is_lower_triangular(self, tol):
     1338        r"""
     1339        Returns ``True`` if the entries above the diagonal are all zero.
     1340
     1341        INPUT:
     1342
     1343        - ``tol`` -  the largest value of the absolute value of the
     1344          difference between two matrix entries for which they will
     1345          still be considered equal.
     1346
     1347        OUTPUT:
     1348
     1349        Returns ``True`` if each entry above the diagonal (entries
     1350        with a row index less than the column index) is zero.
     1351
     1352        EXAMPLES::
     1353
     1354            sage: A = matrix(RDF, [[ 2.0, 0.0,  0.0],
     1355            ...                    [ 1.0, 3.0,  0.0],
     1356            ...                    [-4.0, 2.0, -1.0]])
     1357            sage: A._is_lower_triangular(1.0e-17)
     1358            True
     1359            sage: A[1,2] = 10^-13
     1360            sage: A._is_lower_triangular(1.0e-14)
     1361            False
     1362        """
     1363        global numpy
     1364        if numpy is None:
     1365            import numpy
     1366        cdef Py_ssize_t i, j
     1367        for i in range(self._nrows):
     1368            for j in range(i+1, self._ncols):
     1369                if numpy.absolute(self.get_unsafe(i,j)) > tol:
     1370                    return False
     1371        return True
     1372
     1373    def is_hermitian(self, tol = 1e-12, algorithm='orthonormal'):
     1374        r"""
     1375        Returns ``True`` if the matrix is equal to its conjugate-transpose.
     1376
     1377        INPUT:
     1378
     1379        - ``tol`` - default: ``1e-12`` - the largest value of the
     1380          absolute value of the difference between two matrix entries
     1381          for which they will still be considered equal.
     1382
     1383        - ``algorithm`` - default: 'orthonormal' - set to 'orthonormal'
     1384          for a stable procedure and set to 'naive' for a fast
     1385          procedure.
     1386
     1387        OUTPUT:
     1388
     1389        ``True`` if the matrix is square and equal to the transpose
     1390        with every entry conjugated, and ``False`` otherwise.
     1391
     1392        Note that if conjugation has no effect on elements of the base
     1393        ring (such as for integers), then the :meth:`is_symmetric`
     1394        method is equivalent and faster.
     1395
     1396        The tolerance parameter is used to allow for numerical values
     1397        to be equal if there is a slight difference due to round-off
     1398        and other imprecisions.
     1399
     1400        The result is cached, on a per-tolerance and per-algorithm basis.
     1401
     1402        ALGORITHMS::
     1403
     1404        The naive algorithm simply compares corresponing entries on either
     1405        side of the diagonal (and on the diagonal itself) to see if they are
     1406        conjugates, with equality controlled by the tolerance parameter.
     1407
     1408        The orthonormal algorithm first computes a Schur decomposition
     1409        (via the :meth:`schur` method) and checks that the result is a
     1410        diagonal matrix with real entries.
     1411
     1412        So the naive algorithm can finish quickly for a matrix that is not
     1413        Hermitian, while the orthonormal algorithm will always compute a
     1414        Schur decomposition before going through a similar check of the matrix
     1415        entry-by-entry.
     1416
     1417        EXAMPLES::
     1418
     1419            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
     1420            ...                    [-3 - I,     -4*I,     -2],
     1421            ...                    [-1 + I, -2 - 8*I,  2 + I]])
     1422            sage: A.is_hermitian(algorithm='orthonormal')
     1423            False
     1424            sage: A.is_hermitian(algorithm='naive')
     1425            False
     1426            sage: B = A*A.conjugate_transpose()
     1427            sage: B.is_hermitian(algorithm='orthonormal')
     1428            True
     1429            sage: B.is_hermitian(algorithm='naive')
     1430            True
     1431
     1432        A matrix that is nearly Hermitian, but for one non-real
     1433        diagonal entry. ::
     1434
     1435            sage: A = matrix(CDF, [[    2,   2-I, 1+4*I],
     1436            ...                    [  2+I,   3+I, 2-6*I],
     1437            ...                    [1-4*I, 2+6*I,     5]])
     1438            sage: A.is_hermitian(algorithm='orthonormal')
     1439            False
     1440            sage: A[1,1] = 132
     1441            sage: A.is_hermitian(algorithm='orthonormal')
     1442            True
     1443
     1444        We get a unitary matrix from the SVD routine and use this
     1445        numerical matrix to create a matrix that should be Hermitian
     1446        (indeed it should be the identity matrix), but with some
     1447        imprecision.  We use this to illustrate that if the tolerance
     1448        is set too small, then we can be too strict about the equality
     1449        of entries and achieve the wrong result.  ::
     1450
     1451            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
     1452            ...                    [-3 - I,     -4*I,     -2],
     1453            ...                    [-1 + I, -2 - 8*I,  2 + I]])
     1454            sage: U, _, _ = A.SVD()
     1455            sage: B=U*U.conjugate_transpose()
     1456            sage: B.is_hermitian(algorithm='naive')
     1457            True
     1458            sage: B.is_hermitian(algorithm='naive', tol=1.0e-17)
     1459            False
     1460            sage: B.is_hermitian(algorithm='naive', tol=1.0e-15)
     1461            True
     1462
     1463        A square, empty matrix is trivially Hermitian.  ::
     1464
     1465            sage: A = matrix(RDF, 0, 0)
     1466            sage: A.is_hermitian()
     1467            True
     1468
     1469        Rectangular matrices are never Hermitian, no matter which
     1470        algorithm is requested.  ::
     1471
     1472            sage: A = matrix(CDF, 3, 4)
     1473            sage: A.is_hermitian()
     1474            False
     1475
     1476        TESTS:
     1477
     1478        The tolerance must be strictly positive.  ::
     1479
     1480            sage: A = matrix(RDF, 2, range(4))
     1481            sage: A.is_hermitian(tol = -3.1)
     1482            Traceback (most recent call last):
     1483            ...
     1484            ValueError: tolerance must be positive, not -3.1
     1485
     1486        The ``algorithm`` keyword gets checked.  ::
     1487
     1488            sage: A = matrix(RDF, 2, range(4))
     1489            sage: A.is_hermitian(algorithm='junk')
     1490            Traceback (most recent call last):
     1491            ...
     1492            ValueError: algorithm must be 'naive' or 'orthonormal', not junk
     1493
     1494        AUTHOR:
     1495
     1496         - Rob Beezer (2011-03-30)
     1497        """
     1498        import sage.rings.complex_double
     1499        global numpy
     1500        tol = float(tol)
     1501        if tol <= 0:
     1502            raise ValueError('tolerance must be positive, not {0}'.format(tol))
     1503        if not algorithm in ['naive', 'orthonormal']:
     1504            raise ValueError("algorithm must be 'naive' or 'orthonormal', not {0}".format(algorithm))
     1505
     1506        key = 'hermitian_{0}_{1}'.format(algorithm, tol)
     1507        h = self.fetch(key)
     1508        if not h is None:
     1509            return h
     1510        if not self.is_square():
     1511            self.cache(key, False)
     1512            return False
     1513        if self._nrows == 0:
     1514            self.cache(key, True)
     1515            return True
     1516        if numpy is None:
     1517            import numpy
     1518        cdef Py_ssize_t i, j
     1519        cdef Matrix_double_dense T
     1520        if algorithm == 'orthonormal':
     1521            # Schur decomposition over CDF will be diagonal and real iff Hermitian
     1522            _, T = self.schur(base_ring=sage.rings.complex_double.CDF)
     1523            hermitian = T._is_lower_triangular(tol)
     1524            if hermitian:
     1525                for i in range(T._nrows):
     1526                    if numpy.absolute(numpy.imag(T.get_unsafe(i,i))) > tol:
     1527                        hermitian = False
     1528                        break
     1529        elif algorithm == 'naive':
     1530            hermitian = True
     1531            for i in range(self._nrows):
     1532                for j in range(i+1):
     1533                    if numpy.absolute(self.get_unsafe(i,j) - self.get_unsafe(j,i).conjugate()) > tol:
     1534                        hermitian = False
     1535                        break
     1536                if not hermitian:
     1537                    break
     1538        self.cache(key, hermitian)
     1539        return hermitian
     1540
    13371541    def schur(self, base_ring=None):
    13381542        r"""
    13391543        Returns the Schur decomposition of the matrix.