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