I am having issues understanding why different languages are producing different answers for the Schur(QZ) decomposition. I am working on writing some old stuff from Matlab into Julia and Python and am running into some issues with differences between Matlab's
schurfact, and Scipy's
qz (Julia and Scipy agree and Matlab is different). After some digging, it seems the root of the issue is that Matlab uses the LAPACK
zgges by default (and not even sure if you can use
dgges, but didn't look very hard) whereas both Julia and Scipy default to LAPACK
dgges and can reach
zgges by requesting things to be complex. I can get the answer that I am "looking" for, but am somewhat unsatisfied with simply "arriving" at the answer.
Is there a reason that Matlab sticks to
zgges? For real numbers, why do
dgges produce different answers(the answers differ by signs on only a few numbers that aren't necessarily related by row or column)?
Here is an example of what I am talking about (In all three languages for completeness and to facilitate looking at it).
A = reshape(1:9, 3, 3) B = [ 0.47143516 -1.19097569 1.43270697 -0.3126519 -0.72058873 0.88716294 0.85958841 -0.6365235 0.01569637] # This is the answer without using complex numbers Freal = schurfact(A, B) AAreal, BBreal = Freal[:S], Freal[:T] Qreal, Zreal = Freal[:Q], Freal[:Z] # This is the complex answer turned into reals # and matches Matlab's answer Fcomplex = schurfact(A+0im, B+0im) AAcomplex, BBcomplex = real(Fcomplex[:S]), real(Fcomplex[:T]) Qcomplex, Zcomplex = real(Fcomplex[:Q]), real(Fcomplex[:Z])
import numpy as np import scipy.linalg as la A = np.arange(1, 10).reshape(3, 3).T B = np.array([[ 0.47143516, -1.19097569, 1.43270697], [-0.3126519 , -0.72058873, 0.88716294], [ 0.85958841, -0.6365235 , 0.01569637]]) # Answer without using complex AAreal, BBreal, Qreal, Zreal = la.qz(A, B) # Answer using complex and # matches Matlab's answer Fcomplex = la.qz(A, B, output="complex") AAcomplex, BBcomplex, Qcomplex, Zcomplex = map(np.real, Fcomplex)
A = reshape(1:9, 3, 3) B = [ 0.47143516 -1.19097569 1.43270697 -0.3126519 -0.72058873 0.88716294 0.85958841 -0.6365235 0.01569637] [AAm, BBm, Qm, Zm] = qz(A, B)
Sorry for being a little long, but I'm interested in why this is happening.