4
$\begingroup$

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 qz, Julia'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 zgges and 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).

Julia

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])

Python

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)

Matlab

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.

$\endgroup$
  • $\begingroup$ FYI, qz(A,B,'real') is how you get dgges in Matlab. $\endgroup$ – Federico Poloni Sep 12 '14 at 6:26
4
$\begingroup$

The MATLAB syntax qz(A,B,'real') is consistent with schur(A,'real'), so we might as well ask why the default is complex in the Schur form.

Two reasons spring to mind.

  1. Backward compatibility. Probably there was a time when only the complex Schur form was implemented in Matlab (possibly from the pre-LAPACK times), and the default is retaining that behavior, not to break existing code. Apparently backward compatibility is a big deal in Matlab, even when it generates abominations (cfr. the syntax of logspace for b=pi).

  2. Least surprise. When people call schur, they expect the triangular factor to be, well, triangular. Those who want a real Schur form know that there are two variants, at least.

$\endgroup$
  • $\begingroup$ Thank you. I saw that Matlab allowed for qz(A, B, 'real') last night and produced the same answer as the other languages. This provides some insight into perhaps why it has a different default. $\endgroup$ – cc7768 Sep 12 '14 at 13:54
  • 1
    $\begingroup$ +1 you are right about the backward compatibility thing. Before LAPACK, MATLAB used C-translated routine from EISPACK for QZ which always assumed complex arrays: mathworks.com/company/newsletters/articles/… $\endgroup$ – Amro Nov 12 '14 at 6:54
1
$\begingroup$

Frequently, this sort of "disagreement" occurs between languages, libraries, versions of libraries, etc. because different conventions are being enforced, and the requested decompositions are equivalent (but not identical). I encountered this sort of issue in implementing an incremental singular value decomposition. The result I obtained with the algorithm I implemented differed from the standard singular value decomposition in the signs of the singular vectors, because the decomposition is unique up to the signs of those vectors.

You're best off asking a MATLAB-specific forum or mailing list why MATLAB chooses zgges.

$\endgroup$
  • $\begingroup$ Thanks. I think the key insight is that the decompositions are equivalent but not identical. The only thing that differed was the signs of some of the numbers, and it just "magically" turns out that the signs on Q, Z, AA, and BB all differ in exactly the right places for them to cancel out. Added those details in case others run into a similar issue in the future. Accepted and would upvote, but first question on scicomp. $\endgroup$ – cc7768 Sep 11 '14 at 23:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.