From bug-octave-request at bevo dot che dot wisc dot edu Mon Mar 5 09:09:21 2001 Subject: RE: numerical artifacts by defective input to matrix functions From: "Lippert, Ross A." To: "'Rolf Fabian'" Cc: "'bug-octave UWISC'" Date: Mon, 5 Mar 2001 09:09:10 -0600 The Schur decomposition has more uses that you might think. >A Schur decomposition is only appropriate for the >LIMITED subclass of 'normal' matrices, which are known to >be 'unitarily diagonalizable' ( with all distinct -or- multiple >eigenvalues). Normal matrices A have the property that they I did not say that the schur decomposition should be used to perform the diagonalization. I do not recommend performing diagonaization at all. funm can be evaluated using the complex schur decomposition. The Schur decomposition, as I understand it, can be applied to any square matrix and produces a decomposition A = U * R * U', where U is unitary and R is upper triangular. > A * A' - A' * A == ZERO (**) does imply that the schur decomposition does diagonalize the matrix bc it is equivalent to R *R' - R' *R == ZERO which can only happen when R is diagonal (a stronger condition that upper triangularity). Once again, diagonalization is not a step needed to do funm. **** math theory I guess I can state the basic idea for the algorithm since it is a mathematical theory and not protected by anyone's copyright. A = U*R*U', you do complex schur (using, perhaps, the trick you suggest) I think all parties will agree that R is a upper-tri matrix with the eigenvalues of A on the diagonal. now, to apply f spectrally to A, one recognizes that f(A) = U * f(R) * U' (proof: spectral application of a function is invariant under all similarity transformations. With a little more proof it can be shown that f(R) must also be upper-tri.) What is now needed is just f(R). But f(R) can be computed without diagonalization. You just need to use the following facts: 1) diag(f(R)) = f(diag(R)) 2) f(R)*R - R*f(R) = 0 **** why is this better? 1) the Schur decomposition exists and is stable for all matrices, while the eigen decomposition is not. 2) it is possible to do some extra massaging of the f(R) algorithm so it handles multiple eigenvalues pretty well-- at the very least, it is possible to get good bounds on f(R) so you can know whether or not to trust the f(R) you compute. -r -----Original Message----- From: Rolf Fabian [mailto:fabian at TU-Cottbus dot De] Sent: Monday, March 05, 2001 4:44 AM To: Lippert, Ross A. Cc: 'bug-octave UWISC' Subject: RE: numerical artifacts by defective input to matrix functions Ross A. Lippert wrote |Just because a matrix has multiple eigenvalues does not mean |that it does not have a valid square root. |The same is true for many of the other functions. I've already pointed this out to Paul Kienzle. |The eigen decomposition looks good in texts for function application, |but it is really the schur decomposition which is numerically stable. |Use that and you will be better off. A Schur decomposition is only appropriate for the LIMITED subclass of 'normal' matrices, which are known to be 'unitarily diagonalizable' ( with all distinct -or- multiple eigenvalues). Normal matrices A have the property that they >> commutate with their own (complex conjugated) transposed << A * A' - A' * A == ZERO (**) That's the major decision criterium to check whether a schur decomposition is appropriate. But we've to be careful using Octave '[U,S]=schur(A)' function. Even if A is 'normal', S is ONLY garantueed to be the diagonal matrix of eigenvalues if A is COMPLEX in addition! Otherwise, A is REAL Octave's '[U,S]=schur(A)' call performs a '2x2 BLOCK REAL SCHUR' decomposition (see Octave documentation), which generally does NOT return S as the diagonal matrix of eigenvalues! A workaround for normal real A is to FORCE a complex schur call with an i*A argument [ U,IS ] = schur ( i*A ) and subsequently transform the result IS by S = -i*IS This trick guarantees S to be diagonal if A is normal and real. If the input to [U,S]=schur(A) is NOT normal, you'll get S as a triangular matrix, which is of limited help concerning the calculation of general square matrix functions. Hence, usage of schur decomposition is ONLY appropriate for normal matrices, which pocess an orthogonal system of eigenvectors. The normality check (**) should be the first step in ANY matrix function like 'logm, sqrtm, polyvalm, etc. If it returns TRUE, it is guaranteed, that A can be presented by A == U * S * U' A matrix function F(A) then evaluates to F( A ) = U * diag( F(diag(S)) ) * U' where U and S are results of the Schur decomposition, as long as F(diag(S)) exists. Besides normal matrices, the schur decomposition is NOT APPROPRIATE in the general case. Rolf ------------------------------------------------------------- Octave is freely available under the terms of the GNU GPL. Octave's home on the web: http://www.octave.org How to fund new projects: http://www.octave.org/funding.html Subscription information: http://www.octave.org/archive.html -------------------------------------------------------------