| 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 | |