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

sage/matrix/matrix0.pyx

@@ -2896 +2896 @@
                     return False
         return True
 
+    def is_hermitian(self):
+        r"""
+        Returns ``True`` if the matrix is equal to its conjugate-transpose.
+
+        OUTPUT:
+
+        ``True`` if the matrix is square and equal to the tranpose
+        with every entry conjugated, and ``False`` otherwise.
+
+        Note that if conjugation has no effect on elements of the base
+        ring (such as for integers), then the :meth:`is_symmetric`
+        method is equivalent and faster.
+
+        For numerical matrices a specialized routine available
+        over ``RDF`` and ``CDF`` is a good choice.
+
+        EXAMPLES::
+
+            sage: A = matrix(QQbar, [[ 1 + I,  1 - 6*I, -1 - I],
+            ...                      [-3 - I,     -4*I,     -2],
+            ...                      [-1 + I, -2 - 8*I,  2 + I]])
+            sage: A.is_hermitian()
+            False
+            sage: B = A*A.conjugate_transpose()
+            sage: B.is_hermitian()
+            True
+
+        Sage has several fields besides the entire complex numbers
+        where conjugation is non-trivial. ::
+
+            sage: F.<b> = QuadraticField(-7)
+            sage: C = matrix(F, [[-2*b - 3,  7*b - 6, -b + 3],
+            ...                  [-2*b - 3, -3*b + 2,   -2*b],
+            ...                  [   b + 1,        0,     -2]])
+            sage: C.is_hermitian()
+            False
+            sage: C = C*C.conjugate_transpose()
+            sage: C.is_hermitian()
+            True
+
+        A matrix that is nearly Hermitian, but for a non-real
+        diagonal entry. ::
+
+            sage: A = matrix(QQbar, [[    2,   2-I, 1+4*I],
+            ...                      [  2+I,   3+I, 2-6*I],
+            ...                      [1-4*I, 2+6*I,     5]])
+            sage: A.is_hermitian()
+            False
+            sage: A[1,1] = 132
+            sage: A.is_hermitian()
+            True
+
+        Rectangular matrices are never Hermitian.  ::
+
+            sage: A = matrix(QQbar, 3, 4)
+            sage: A.is_hermitian()
+            False
+        """
+        if not self.is_square():
+            return False
+
+        cdef Py_ssize_t i,j
+
+        for i from 0 <= i < self._nrows:
+            for j from 0 <= j <= i:
+                if self.get_unsafe(i,j) != self.get_unsafe(j,i).conjugate():
+                    return False
+        return True
+
     def is_skew_symmetric(self):
         """
         Returns True if this is a skew symmetric matrix.

sage/matrix/matrix_double_dense.pyx

@@ -1327 +1327 @@
             return False
         b = True
         for i from 0 < i < self._nrows:
-            for j from 0 <= j < i:                
+            for j from 0 <= j < i:
                 if math.fabs(self.get_unsafe(i,j) - self.get_unsafe(j,i)) > tol:
                     b = False
                     break
         self.cache(key, b)
         return b
 
+    def is_hermitian(self, tol = 1e-12):
+        r"""
+        Return ``True`` if the matrix is equal to its conjugate-transpose.
+
+        INPUT:
+
+        - ``tol`` - default: ``1e-12`` - the largest value of the
+          absolute value of the difference between two matrix entries
+          for which they will still be considered equal.
+
+        OUTPUT:
+
+        ``True`` if the matrix is square and equal to the tranpose
+        with every entry conjugated, and ``False`` otherwise.
+
+        Note that if conjugation has no effect on elements of the base
+        ring (such as for integers), then the :meth:`is_symmetric`
+        method is equivalent and faster.
+
+        The tolerance parameter is used to allow for numerical values
+        to be equal if there is a slight difference due to round-off
+        and other imprecisions.
+
+        The result is cached, on a per-tolerance basis.
+
+        EXAMPLES::
+
+            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
+            ...                    [-3 - I,     -4*I,     -2],
+            ...                    [-1 + I, -2 - 8*I,  2 + I]])
+            sage: A.is_hermitian()
+            False
+            sage: B = A*A.conjugate_transpose()
+            sage: B.is_hermitian()
+            True
+
+        We get a unitary matrix from the SVD routine and use this
+        numerical matrix to create a matrix that should be Hermitian
+        (indeed it should be the identity matrix), but with some
+        imprecision.  We use this to illustrate that if the tolerance
+        is set too small, then we can be too strict about the equality
+        of entries and achieve the wrong result.  ::
+
+            sage: A = matrix(CDF, [[ 1 + I,  1 - 6*I, -1 - I],
+            ...                    [-3 - I,     -4*I,     -2],
+            ...                    [-1 + I, -2 - 8*I,  2 + I]])
+            sage: U, _, _ = A.SVD()
+            sage: B=U*U.conjugate_transpose()
+            sage: B.is_hermitian()
+            True
+            sage: B.is_hermitian(tol=5.0e-17)
+            True
+            sage: B.is_hermitian(tol=3.0e-17)
+            False
+
+        A matrix that is nearly Hermitian, but for a non-real
+        diagonal entry. ::
+
+            sage: A = matrix(CDF, [[    2,   2-I, 1+4*I],
+            ...                    [  2+I,   3+I, 2-6*I],
+            ...                    [1-4*I, 2+6*I,     5]])
+            sage: A.is_hermitian()
+            False
+            sage: A[1,1] = 132
+            sage: A.is_hermitian()
+            True
+
+        Rectangular matrices are never Hermitian.  ::
+
+            sage: A = matrix(CDF, 3, 4)
+            sage: A.is_hermitian()
+            False
+        """
+        global numpy
+        tol = float(tol)
+        key = 'hermitian_%s'%tol
+        b = self.fetch(key)
+        if not b is None:
+            return b
+        if not self.is_square():
+            self.cache(key, False)
+            return False
+        if numpy is None:
+            import numpy
+        cdef Py_ssize_t i, j
+        hermitian = True
+        for i from 0 < i < self._nrows:
+            for j from 0 <= j <= i:
+                if numpy.absolute(self.get_unsafe(i,j) - self.get_unsafe(j,i).conjugate()) > tol:
+                    hermitian = False
+                    break
+        self.cache(key, hermitian)
+        return hermitian
 
     def cholesky(self):
         r"""