# Ticket #11027: trac_11027-schur-decomposition.patch

File trac_11027-schur-decomposition.patch, 13.2 KB (added by Rob Beezer, 12 years ago)
• ## sage/matrix/matrix_double_dense.pyx

# HG changeset patch
# User Rob Beezer <beezer@ups.edu>
# Date 1301065587 25200
# Node ID bd05f5341465c21a937a026a98584faffa449bbe
# Parent  6a679959b54b7975af1841df4dd77d28b0faef28
11027: Schur matrix decomposition

diff -r 6a679959b54b -r bd05f5341465 sage/matrix/matrix_double_dense.pyx
 a self.cache(key, b) return b def schur(self, base_ring=None): r""" Returns the Schur decomposition of the matrix. INPUT: - base_ring - optional, defaults to the base ring of self. Use this to request the base ring of the returned matrices, which will affect the format of the results. OUTPUT: A pair of immutable matrices.  The first is a unitary matrix Q. The second, T, is upper triangular when returned over the complex numbers, while it is almost upper triangular over the reals.  In the latter case, there can be some 2\times 2 blocks on the diagonal which represent a pair of conjugate complex eigenvalues of self. If self is the matrix A, then .. math:: A = QT({\overline Q})^t where the latter matrix is the conjugate-transpose of Q, which is also the inverse of Q since it is unitary. Note that in the case of a normal matrix (Hermitian, symmetric, and others), the upper triangular matrix is  a diagonal matrix with eigenvalues of self on the diagonal, and the unitary matrix has columns that form an orthonormal basis composed of eigenvectors of self.  This is known as "orthonormal diagonalization". EXAMPLES: First over the complexes.  The similar matrix is always upper triangular in this case.  :: sage: A = matrix(CDF, 4, 4, range(16)) + matrix(CDF, 4, 4, [x^3*I for x in range(0, 16)]) sage: Q, T = A.schur() sage: Q.round(4) [   0.006 - 0.001*I  0.0425 + 0.1018*I -0.6446 + 0.2201*I  0.3981 + 0.6045*I] [  0.0875 - 0.005*I  0.1984 + 0.5537*I -0.5022 + 0.1944*I -0.3414 - 0.4896*I] [ 0.3572 - 0.0159*I   0.2444 + 0.685*I  0.4499 - 0.1643*I  0.1983 + 0.2728*I] [  0.929 - 0.0377*I   -0.111 - 0.317*I  -0.1212 + 0.044*I -0.0466 - 0.0627*I] sage: T.round(4) [  30.7332 + 4648.5418*I  -2342.252 + 767.2915*I -477.0197 - 1460.6287*I  -532.7363 + 318.6099*I] [                      0    -0.1841 - 159.0576*I      151.3315 + 0.458*I    -15.2369 - 67.0242*I] [                      0                       0      -0.5239 + 11.158*I       6.2073 - 1.1208*I] [                      0                       0                       0      -0.0252 - 0.6422*I] sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) [1.0   0   0   0] [  0 1.0   0   0] [  0   0 1.0   0] [  0   0   0 1.0] sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)]) True sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) [0 0 0 0] [0 0 0 0] [0 0 0 0] [0 0 0 0] We begin with a real matrix but ask for a decomposition over the complexes.  The result will always have an upper triangular matrix for T. :: sage: A = matrix(RDF, 4, 4, [x^3 for x in range(16)]) sage: Q, T = A.schur(base_ring=CDF) sage: Q.round(4) [             0.006  0.1102 - 0.0003*I  0.6792 + 0.0468*I -0.3382 - 0.6403*I] [            0.0877  0.5882 - 0.0015*I   0.5374 + 0.037*I  0.2787 + 0.5277*I] [            0.3575  0.7273 - 0.0019*I -0.4781 - 0.0329*I -0.1573 - 0.2979*I] [ 0.9298 + 0.0001*I -0.3359 + 0.0009*I  0.1288 + 0.0089*I   0.0364 + 0.069*I] sage: T.round(4) [             4648.5536   2464.7695 - 6.6819*I   1532.623 + 105.406*I -290.0435 - 549.0666*I] [                     0              -159.0726  -150.9268 - 10.7911*I     31.954 + 60.8896*I] [                     0                      0                11.1621     -3.3197 - 5.3583*I] [                     0                      0                      0                -0.6431] sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) [1.0   0   0   0] [  0 1.0   0   0] [  0   0 1.0   0] [  0   0   0 1.0] sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)]) True sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) [0 0 0 0] [0 0 0 0] [0 0 0 0] [0 0 0 0] Now totally over the reals.  But with complex eigenvalues, the similar matrix may not be upper triangular. But "at worst" there may be some 2\times 2 blocks on the diagonal which represent a pair of conjugate complex eigenvalues. These blocks will then just interrupt the zeros below the main diagonal. sage: A = matrix(RDF, 4, 4, [[1, 0, -3, -1], ...                          [4, -16, -7, 0], ...                          [1, 21, 1, -2], ...                          [26, -1, -2, 1]]) sage: Q, T = A.schur() sage: Q.round(4) [ -0.088  0.5832  0.7724  0.2356] [ 0.1078 -0.2778  0.4746 -0.8282] [-0.2985  0.6989 -0.4072 -0.5066] [ 0.9442   0.307 -0.1109 -0.0437] sage: T.round(4) [ -0.7898   15.455  15.6639  14.4823] [ -0.3534  -0.7898     15.3 -13.7298] [     0.0      0.0  -5.7102  16.0869] [     0.0      0.0  -4.3681  -5.7102] sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) [1.0 0.0 0.0 0.0] [0.0 1.0 0.0 0.0] [0.0 0.0 1.0 0.0] [0.0 0.0 0.0 1.0] sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i)]) False sage: all([T.zero_at(1.0e-12)[i,j] == 0 for i in range(4) for j in range(i-1)]) True sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) [0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0] Starting with complex numbers and requesting a result over the reals will never happen.  :: sage: A = matrix(CDF, 2, 2, [[2+I, -1+3*I], [5-4*I, 2-7*I]]) sage: A.schur(base_ring=RDF) Traceback (most recent call last): ... TypeError: unable to convert input matrix over CDF to a matrix over RDF If theory predicts your matrix is real, but it contains some very small imaginary parts, you can specify the cutoff for "small" imaginary parts, Then request the output as real matrices, and let the routine do the rest. :: sage: A = matrix(RDF, 2, 2, [1, 1, -1, 0]) + matrix(CDF, 2, 2, [1.0e-14*I]*4) sage: B = A.zero_at(1.0e-12) sage: B.parent() Full MatrixSpace of 2 by 2 dense matrices over Complex Double Field sage: Q, T = B.schur(RDF) sage: Q.parent() Full MatrixSpace of 2 by 2 dense matrices over Real Double Field sage: T.parent() Full MatrixSpace of 2 by 2 dense matrices over Real Double Field sage: Q.round(6) [ 0.707107  0.707107] [-0.707107  0.707107] sage: T.round(6) [ 0.5  1.5] [-0.5  0.5] sage: (Q*T*Q.conjugate().transpose()-B).zero_at(1.0e-11) [0 0] [0 0] A Hermitian matrix has real eigenvalues, so the similar matrix will be upper triangular.  Furthermore, a Hermitian matrix is diagonalizable with respect to an orthonormal basis, composed of eigenvectors of the matrix.  Here that basis is the set of columns of the unitary matrix.  :: sage: A = matrix(CDF, [[        52,   -9*I - 8,    6*I - 187,  -188*I + 2], ...                    [   9*I - 8,         12,   -58*I + 59,   30*I + 42], ...                    [-6*I - 187,  58*I + 59,         2677, 2264*I + 65], ...                    [ 188*I + 2, -30*I + 42, -2264*I + 65,       2080]]) sage: Q, T = A.schur() sage: T = T.zero_at(1.0e-12).change_ring(RDF) sage: T.round(6) [4680.13301        0.0        0.0        0.0] [       0.0 102.715967        0.0        0.0] [       0.0        0.0  35.039344        0.0] [       0.0        0.0        0.0    3.11168] sage: (Q*Q.conjugate().transpose()).zero_at(1.0e-12) [1.0   0   0   0] [  0 1.0   0   0] [  0   0 1.0   0] [  0   0   0 1.0] sage: (Q*T*Q.conjugate().transpose()-A).zero_at(1.0e-11) [0 0 0 0] [0 0 0 0] [0 0 0 0] [0 0 0 0] Similarly, a real symmetric matrix has only real eigenvalues, and there is an orthonormal basis composed of eigenvectors of the matrix. sage: A = matrix(RDF, [[ 1, -2, 5, -3], ...                    [-2,  9, 1,  5], ...                    [ 5,  1, 3 , 7], ...                    [-3,  5, 7, -8]]) sage: Q, T = A.schur() sage: Q.round(4) [-0.3027  -0.751   0.576 -0.1121] [  0.139 -0.3892 -0.2648  0.8713] [ 0.4361   0.359  0.7599  0.3217] [ -0.836  0.3945  0.1438  0.3533] sage: T.round(4) [-13.5698      0.0      0.0      0.0] [     0.0  -0.8508     -0.0     -0.0] [     0.0      0.0   7.7664      0.0] [     0.0      0.0      0.0  11.6542] sage: (Q*Q.transpose()).zero_at(1.0e-12) [1.0 0.0 0.0 0.0] [0.0 1.0 0.0 0.0] [0.0 0.0 1.0 0.0] [0.0 0.0 0.0 1.0] sage: (Q*T*Q.transpose()-A).zero_at(1.0e-11) [0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0] [0.0 0.0 0.0 0.0] The results are cached, both as a real factorization and also as a complex factorization.  This means the returned matrices are immutable.  :: sage: A = matrix(RDF, 2, 2, [[0, -1], [1, 0]]) sage: Qr, Tr = A.schur(base_ring=RDF) sage: Qc, Tc = A.schur(base_ring=CDF) sage: all([M.is_immutable() for M in [Qr, Tr, Qc, Tc]]) True sage: Tr.round(6) != Tc.round(6) True TESTS: The Schur factorization is only defined for square matrices.  :: sage: A = matrix(RDF, 4, 5, range(20)) sage: A.schur() Traceback (most recent call last): ... ValueError: Schur decomposition requires a square matrix, not a 4 x 5 matrix A base ring request is checked. :: sage: A = matrix(RDF, 3, range(9)) sage: A.schur(base_ring=QQ) Traceback (most recent call last): ... ValueError: base ring of Schur decomposition matrices must be RDF or CDF, not Rational Field """ global scipy from sage.rings.real_double import RDF from sage.rings.complex_double import CDF cdef Matrix_double_dense Q, T if not self.is_square(): raise ValueError('Schur decomposition requires a square matrix, not a {0} x {1} matrix'.format(self.nrows(), self.ncols())) if base_ring == None: base_ring = self.base_ring() if not base_ring in [RDF, CDF]: raise ValueError('base ring of Schur decomposition matrices must be RDF or CDF, not {0}'.format(base_ring)) if self.base_ring() != base_ring: try: self = self.change_ring(base_ring) except TypeError: raise TypeError('unable to convert input matrix over CDF to a matrix over RDF') if base_ring == CDF: format = 'complex' else: format = 'real' schur = self.fetch('schur_factors_' + format) if not schur is None: return schur Q = self._new(self._nrows, self._nrows) T = self._new(self._nrows, self._nrows) if scipy is None: import scipy import scipy.linalg T._matrix_numpy, Q._matrix_numpy = scipy.linalg.schur(self._matrix_numpy, output=format) Q.set_immutable() T.set_immutable() # Our return order is the reverse of NumPy's schur = (Q, T) self.cache('schur_factors_' + format, schur) return schur def cholesky(self): r"""