| 1337 | def schur(self, base_ring=None): |
| 1338 | r""" |
| 1339 | Returns the Schur decomposition of the matrix. |
| 1340 | |
| 1341 | INPUT: |
| 1342 | |
| 1343 | - ``base_ring`` - optional, defaults to the base ring of ``self``. |
| 1344 | Use this to request the base ring of the returned matrices, which |
| 1345 | will affect the format of the results. |
| 1346 | |
| 1347 | OUTPUT: |
| 1348 | |
| 1349 | A pair of immutable matrices. The first is a unitary matrix `Q`. |
| 1350 | The second, `T`, is upper triangular when returned over the complex |
| 1351 | numbers, while it is almost upper triangular over the reals. In the |
| 1352 | latter case, there can be some `2\times 2` blocks on the diagonal |
| 1353 | which represent a pair of conjugate complex eigenvalues of ``self``. |
| 1354 | |
| 1355 | If ``self`` is the matrix `A`, then |
| 1356 | |
| 1357 | .. math:: |
| 1358 | |
| 1359 | A = QT({\overline Q})^t |
| 1360 | |
| 1361 | where the latter matrix is the conjugate-transpose of ``Q``, which |
| 1362 | is also the inverse of ``Q`` since it is unitary. |
| 1363 | |
| 1364 | Note that in the case of a normal matrix (Hermitian, symmetric, and |
| 1365 | others), the upper triangular matrix is a diagonal matrix with |
| 1366 | eigenvalues of ``self`` on the diagonal, and the unitary matrix |
| 1367 | has columns that form an orthonormal basis composed of eigenvectors |
| 1368 | of ``self``. This is known as "orthonormal diagonalization". |
| 1369 | |
| 1370 | EXAMPLES: |
| 1371 | |
| 1372 | First over the complexes. The similar matrix is always upper |
| 1373 | triangular in this case. :: |
| 1374 | |
| 1375 | sage: A = matrix(CDF, 4, 4, range(16)) + matrix(CDF, 4, 4, [x^3*I for x in range(0, 16)]) |
| 1376 | sage: Q, T = A.schur() |
| 1377 | sage: Q.round(4) |
| 1378 | [ 0.006 - 0.001*I 0.0425 + 0.1018*I -0.6446 + 0.2201*I 0.3981 + 0.6045*I] |
| 1379 | [ 0.0875 - 0.005*I 0.1984 + 0.5537*I -0.5022 + 0.1944*I -0.3414 - 0.4896*I] |
| 1380 | [ 0.3572 - 0.0159*I 0.2444 + 0.685*I 0.4499 - 0.1643*I 0.1983 + 0.2728*I] |
| 1381 | [ 0.929 - 0.0377*I -0.111 - 0.317*I -0.1212 + 0.044*I -0.0466 - 0.0627*I] |
| 1382 | sage: T.round(4) |
| 1383 | [ 30.7332 + 4648.5418*I -2342.252 + 767.2915*I -477.0197 - 1460.6287*I -532.7363 + 318.6099*I] |
| 1384 | [ 0 -0.1841 - 159.0576*I 151.3315 + 0.458*I -15.2369 - 67.0242*I] |
| 1385 | [ 0 0 -0.5239 + 11.158*I 6.2073 - 1.1208*I] |
| 1386 | [ 0 0 0 -0.0252 - 0.6422*I] |
| 1387 | sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) |
| 1388 | [1.0 0 0 0] |
| 1389 | [ 0 1.0 0 0] |
| 1390 | [ 0 0 1.0 0] |
| 1391 | [ 0 0 0 1.0] |
| 1392 | sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)]) |
| 1393 | True |
| 1394 | sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) |
| 1395 | [0 0 0 0] |
| 1396 | [0 0 0 0] |
| 1397 | [0 0 0 0] |
| 1398 | [0 0 0 0] |
| 1399 | |
| 1400 | We begin with a real matrix but ask for a decomposition over the |
| 1401 | complexes. The result will always have an upper triangular |
| 1402 | matrix for ``T``. :: |
| 1403 | |
| 1404 | sage: A = matrix(RDF, 4, 4, [x^3 for x in range(16)]) |
| 1405 | sage: Q, T = A.schur(base_ring=CDF) |
| 1406 | sage: Q.round(4) |
| 1407 | [ 0.006 0.1102 - 0.0003*I 0.6792 + 0.0468*I -0.3382 - 0.6403*I] |
| 1408 | [ 0.0877 0.5882 - 0.0015*I 0.5374 + 0.037*I 0.2787 + 0.5277*I] |
| 1409 | [ 0.3575 0.7273 - 0.0019*I -0.4781 - 0.0329*I -0.1573 - 0.2979*I] |
| 1410 | [ 0.9298 + 0.0001*I -0.3359 + 0.0009*I 0.1288 + 0.0089*I 0.0364 + 0.069*I] |
| 1411 | sage: T.round(4) |
| 1412 | [ 4648.5536 2464.7695 - 6.6819*I 1532.623 + 105.406*I -290.0435 - 549.0666*I] |
| 1413 | [ 0 -159.0726 -150.9268 - 10.7911*I 31.954 + 60.8896*I] |
| 1414 | [ 0 0 11.1621 -3.3197 - 5.3583*I] |
| 1415 | [ 0 0 0 -0.6431] |
| 1416 | sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) |
| 1417 | [1.0 0 0 0] |
| 1418 | [ 0 1.0 0 0] |
| 1419 | [ 0 0 1.0 0] |
| 1420 | [ 0 0 0 1.0] |
| 1421 | sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)]) |
| 1422 | True |
| 1423 | sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) |
| 1424 | [0 0 0 0] |
| 1425 | [0 0 0 0] |
| 1426 | [0 0 0 0] |
| 1427 | [0 0 0 0] |
| 1428 | |
| 1429 | Now totally over the reals. But with complex eigenvalues, the |
| 1430 | similar matrix may not be upper triangular. But "at worst" there |
| 1431 | may be some `2\times 2` blocks on the diagonal which represent |
| 1432 | a pair of conjugate complex eigenvalues. These blocks will then |
| 1433 | just interrupt the zeros below the main diagonal. |
| 1434 | |
| 1435 | sage: A = matrix(RDF, 4, 4, [[1, 0, -3, -1], |
| 1436 | ... [4, -16, -7, 0], |
| 1437 | ... [1, 21, 1, -2], |
| 1438 | ... [26, -1, -2, 1]]) |
| 1439 | sage: Q, T = A.schur() |
| 1440 | sage: Q.round(4) |
| 1441 | [ -0.088 0.5832 0.7724 0.2356] |
| 1442 | [ 0.1078 -0.2778 0.4746 -0.8282] |
| 1443 | [-0.2985 0.6989 -0.4072 -0.5066] |
| 1444 | [ 0.9442 0.307 -0.1109 -0.0437] |
| 1445 | sage: T.round(4) |
| 1446 | [ -0.7898 15.455 15.6639 14.4823] |
| 1447 | [ -0.3534 -0.7898 15.3 -13.7298] |
| 1448 | [ 0.0 0.0 -5.7102 16.0869] |
| 1449 | [ 0.0 0.0 -4.3681 -5.7102] |
| 1450 | sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) |
| 1451 | [1.0 0.0 0.0 0.0] |
| 1452 | [0.0 1.0 0.0 0.0] |
| 1453 | [0.0 0.0 1.0 0.0] |
| 1454 | [0.0 0.0 0.0 1.0] |
| 1455 | sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)]) |
| 1456 | False |
| 1457 | sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i-1)]) |
| 1458 | True |
| 1459 | sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) |
| 1460 | [0.0 0.0 0.0 0.0] |
| 1461 | [0.0 0.0 0.0 0.0] |
| 1462 | [0.0 0.0 0.0 0.0] |
| 1463 | [0.0 0.0 0.0 0.0] |
| 1464 | |
| 1465 | Starting with complex numbers and requesting a result over the reals |
| 1466 | will never happen. :: |
| 1467 | |
| 1468 | sage: A = matrix(CDF, 2, 2, [[2+I, -1+3*I], [5-4*I, 2-7*I]]) |
| 1469 | sage: A.schur(base_ring=RDF) |
| 1470 | Traceback (most recent call last): |
| 1471 | ... |
| 1472 | TypeError: unable to convert input matrix over CDF to a matrix over RDF |
| 1473 | |
| 1474 | If theory predicts your matrix is real, but it contains some |
| 1475 | very small imaginary parts, you can specify the cutoff for "small" |
| 1476 | imaginary parts, Then request the output as real matrices, and let |
| 1477 | the routine do the rest. :: |
| 1478 | |
| 1479 | sage: A = matrix(RDF, 2, 2, [1, 1, -1, 0]) + matrix(CDF, 2, 2, [1.0e-14*I]*4) |
| 1480 | sage: B = A.zero_at(1.0e-12) |
| 1481 | sage: B.parent() |
| 1482 | Full MatrixSpace of 2 by 2 dense matrices over Complex Double Field |
| 1483 | sage: Q, T = B.schur(RDF) |
| 1484 | sage: Q.parent() |
| 1485 | Full MatrixSpace of 2 by 2 dense matrices over Real Double Field |
| 1486 | sage: T.parent() |
| 1487 | Full MatrixSpace of 2 by 2 dense matrices over Real Double Field |
| 1488 | sage: Q.round(6) |
| 1489 | [ 0.707107 0.707107] |
| 1490 | [-0.707107 0.707107] |
| 1491 | sage: T.round(6) |
| 1492 | [ 0.5 1.5] |
| 1493 | [-0.5 0.5] |
| 1494 | sage: (Q*T*Q.conjugate().transpose()-B).zero_at(1.0e-11) |
| 1495 | [0 0] |
| 1496 | [0 0] |
| 1497 | |
| 1498 | A Hermitian matrix has real eigenvalues, so the similar matrix |
| 1499 | will be upper triangular. Furthermore, a Hermitian matrix is |
| 1500 | diagonalizable with respect to an orthonormal basis, composed |
| 1501 | of eigenvectors of the matrix. Here that basis is the set of |
| 1502 | columns of the unitary matrix. :: |
| 1503 | |
| 1504 | sage: A = matrix(CDF, [[ 52, -9*I - 8, 6*I - 187, -188*I + 2], |
| 1505 | ... [ 9*I - 8, 12, -58*I + 59, 30*I + 42], |
| 1506 | ... [-6*I - 187, 58*I + 59, 2677, 2264*I + 65], |
| 1507 | ... [ 188*I + 2, -30*I + 42, -2264*I + 65, 2080]]) |
| 1508 | sage: Q, T = A.schur() |
| 1509 | sage: T = T.zero_at(1.0e-12).change_ring(RDF) |
| 1510 | sage: T.round(6) |
| 1511 | [4680.13301 0.0 0.0 0.0] |
| 1512 | [ 0.0 102.715967 0.0 0.0] |
| 1513 | [ 0.0 0.0 35.039344 0.0] |
| 1514 | [ 0.0 0.0 0.0 3.11168] |
| 1515 | sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) |
| 1516 | [1.0 0 0 0] |
| 1517 | [ 0 1.0 0 0] |
| 1518 | [ 0 0 1.0 0] |
| 1519 | [ 0 0 0 1.0] |
| 1520 | sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) |
| 1521 | [0 0 0 0] |
| 1522 | [0 0 0 0] |
| 1523 | [0 0 0 0] |
| 1524 | [0 0 0 0] |
| 1525 | |
| 1526 | Similarly, a real symmetric matrix has only real eigenvalues, and |
| 1527 | there is an orthonormal basis composed of eigenvectors of the matrix. |
| 1528 | |
| 1529 | sage: A = matrix(RDF, [[ 1, -2, 5, -3], |
| 1530 | ... [-2, 9, 1, 5], |
| 1531 | ... [ 5, 1, 3 , 7], |
| 1532 | ... [-3, 5, 7, -8]]) |
| 1533 | sage: Q, T = A.schur() |
| 1534 | sage: Q.round(4) |
| 1535 | [-0.3027 -0.751 0.576 -0.1121] |
| 1536 | [ 0.139 -0.3892 -0.2648 0.8713] |
| 1537 | [ 0.4361 0.359 0.7599 0.3217] |
| 1538 | [ -0.836 0.3945 0.1438 0.3533] |
| 1539 | sage: T.round(4) |
| 1540 | [-13.5698 0.0 0.0 0.0] |
| 1541 | [ 0.0 -0.8508 -0.0 -0.0] |
| 1542 | [ 0.0 0.0 7.7664 0.0] |
| 1543 | [ 0.0 0.0 0.0 11.6542] |
| 1544 | sage: (Q*Q.transpose()).zero_at(1.0e-12) |
| 1545 | [1.0 0.0 0.0 0.0] |
| 1546 | [0.0 1.0 0.0 0.0] |
| 1547 | [0.0 0.0 1.0 0.0] |
| 1548 | [0.0 0.0 0.0 1.0] |
| 1549 | sage: (Q*T*Q.transpose()-A).zero_at(1.0e-11) |
| 1550 | [0.0 0.0 0.0 0.0] |
| 1551 | [0.0 0.0 0.0 0.0] |
| 1552 | [0.0 0.0 0.0 0.0] |
| 1553 | [0.0 0.0 0.0 0.0] |
| 1554 | |
| 1555 | The results are cached, both as a real factorization and also as a |
| 1556 | complex factorization. This means the returned matrices are |
| 1557 | immutable. :: |
| 1558 | |
| 1559 | sage: A = matrix(RDF, 2, 2, [[0, -1], [1, 0]]) |
| 1560 | sage: Qr, Tr = A.schur(base_ring=RDF) |
| 1561 | sage: Qc, Tc = A.schur(base_ring=CDF) |
| 1562 | sage: all([M.is_immutable() for M in [Qr, Tr, Qc, Tc]]) |
| 1563 | True |
| 1564 | sage: Tr.round(6) != Tc.round(6) |
| 1565 | True |
| 1566 | |
| 1567 | TESTS: |
| 1568 | |
| 1569 | The Schur factorization is only defined for square matrices. :: |
| 1570 | |
| 1571 | sage: A = matrix(RDF, 4, 5, range(20)) |
| 1572 | sage: A.schur() |
| 1573 | Traceback (most recent call last): |
| 1574 | ... |
| 1575 | ValueError: Schur decomposition requires a square matrix, not a 4 x 5 matrix |
| 1576 | |
| 1577 | A base ring request is checked. :: |
| 1578 | |
| 1579 | sage: A = matrix(RDF, 3, range(9)) |
| 1580 | sage: A.schur(base_ring=QQ) |
| 1581 | Traceback (most recent call last): |
| 1582 | ... |
| 1583 | ValueError: base ring of Schur decomposition matrices must be RDF or CDF, not Rational Field |
| 1584 | """ |
| 1585 | global scipy |
| 1586 | from sage.rings.real_double import RDF |
| 1587 | from sage.rings.complex_double import CDF |
| 1588 | |
| 1589 | cdef Matrix_double_dense Q, T |
| 1590 | |
| 1591 | if not self.is_square(): |
| 1592 | raise ValueError('Schur decomposition requires a square matrix, not a {0} x {1} matrix'.format(self.nrows(), self.ncols())) |
| 1593 | if base_ring == None: |
| 1594 | base_ring = self.base_ring() |
| 1595 | if not base_ring in [RDF, CDF]: |
| 1596 | raise ValueError('base ring of Schur decomposition matrices must be RDF or CDF, not {0}'.format(base_ring)) |
| 1597 | |
| 1598 | if self.base_ring() != base_ring: |
| 1599 | try: |
| 1600 | self = self.change_ring(base_ring) |
| 1601 | except TypeError: |
| 1602 | raise TypeError('unable to convert input matrix over CDF to a matrix over RDF') |
| 1603 | if base_ring == CDF: |
| 1604 | format = 'complex' |
| 1605 | else: |
| 1606 | format = 'real' |
| 1607 | |
| 1608 | schur = self.fetch('schur_factors_' + format) |
| 1609 | if not schur is None: |
| 1610 | return schur |
| 1611 | Q = self._new(self._nrows, self._nrows) |
| 1612 | T = self._new(self._nrows, self._nrows) |
| 1613 | if scipy is None: |
| 1614 | import scipy |
| 1615 | import scipy.linalg |
| 1616 | T._matrix_numpy, Q._matrix_numpy = scipy.linalg.schur(self._matrix_numpy, output=format) |
| 1617 | Q.set_immutable() |
| 1618 | T.set_immutable() |
| 1619 | # Our return order is the reverse of NumPy's |
| 1620 | schur = (Q, T) |
| 1621 | self.cache('schur_factors_' + format, schur) |
| 1622 | return schur |