QZ factorization in Mathematica
- To: mathgroup@smc.vnet.net
- Subject: [mg11982] QZ factorization in Mathematica
- From: elisha@dot.net.au (Luci Ellis)
- Date: Fri, 17 Apr 1998 03:40:19 -0400
- Organization: Elisha, Ink!
I am trying to track down a Mathematica version of the QZ algorithm. This algorithm was first published by Moler and Stewart (1973), "An algorithm for generalised matrix eigenvalue problems" SIAM Journal on Numerical Analysis 10. pp 241-256, and apparently it also appears in Coleman and van Loan (1988) Handbook for matrix computation. According to the paper I am reading, the factorization works like this \!\(\[CapitalLambda]\ = \ Q\ \(\[CapitalGamma]\_0\) Z\) \!\(\[CapitalOmega]\ = \ Q\ \(\[CapitalGamma]\_1\) Z\) ( Lambda = Q. Gamma0. Z and Omega = Q. Gamma1. Z ) where Gamma1 and Gamma0 are known square numerical matrices. CapitalLambda and CapitalOmega are upper triangular. Transpose[Q].Q == IdentityMatrix[n] and also Transpose[Z].Z == IdentityMatrix[n] As far as I can work out from the numerical example in the paper that I'm reading this from, Z is the SchurDecomposition of Inverse[Gamma0].Gamma1 (the paper is H.M. Amman and D.A. Kendrick (1998) "Computing the steady state of linear quadratic optimization models with rational expectations", Economics Letters 58 pp 185-191, if anyone's interested) But I can't work out what Q is and therefore I can't work out what CapitalLambda and CapitalOmega are, which renders the rest of the paper pretty useless. The only other pieces of information I have about the problem is that the ratios of the elements on the lead diagonals of CapitalOmega/CapitalLambda equal Eigenvalues[Inverse[Gamma0].Gamma1]. I tried using the fact that various elements of CapitalLambda and CapitalOmega are zero (because they are upper triangular) to Solve for the elements of Q, but I couldn't get it to work and it must be a bit ugly to write a general version of such a procedure. Has anyone implemented this factorization in Mathematica? I know it exists in another program, but I don't have it and would rather not give money to a company that has broadcast its intention to abandon my platform of choice. A complex-matrix version of the algorithm was published in "Transactions in Mathematical Software" as Algorithm 535: it is available on the web at http://www.netlib.org/toms/535 and probably elsewhere. Unfortunately that version is in Fortran, and although I can vaguely work out what's going on, I don't trust my understanding of either the mathematics or the Fortran language to translate it adequately. Please email me directly as well as posting any reply. Thanks in advance, Luci Ellis elisha@remove-to-reply.dot.net.au