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