6
$\begingroup$

MATLAB has the function qz() such that [AA, BB, Q, Z, V, W] = qz(A, B) produces matrices V and W whose columns are generalized eigenvectors. This function gives the right and left eigenvectors. I didn't find a similar function in Mathematica. Does it exist?

Thanks.

$\endgroup$
  • 3
    $\begingroup$ I'm pretty sure the built-ins Eigensystem, Eigenvalues, or Eigenvectors can do this. See the documentation. $\endgroup$ – march Nov 29 '17 at 16:39
  • 2
    $\begingroup$ This page has maping of some matlab functions to Mathematica, including qz. But Matlab qz returns number of values. So may be there is some functionality not there. But f you have specific example you want to see it in Matlab, you can provide one. $\endgroup$ – Nasser Nov 29 '17 at 17:01
16
$\begingroup$

As already hinted, Mathematica is perfectly capable of computing the eigensystem of a pencil. However, none of the commenters gave the function that is the closest to the functionality of MATLAB's qz(), and that is Mathematica's SchurDecomposition[].

Briefly put, when given a pencil (as a list of two square matrices of the same dimensions), SchurDecomposition[] will return the generalized Schur decomposition (apparently termed the "QZ decomposition" in some circles). More concretely, given a pencil $(\mathbf A,\mathbf B)$, SchurDecomposition[] will return a similar (i.e., same eigenvalues) pencil $(\mathbf T,\mathbf U)$, where $\mathbf T$ and $\mathbf U$ are quasitriangular matrices, and are related to $\mathbf A$ and $\mathbf B$ as $\mathbf Q\mathbf T\mathbf Z^\top=\mathbf A$, and $\mathbf Q\mathbf U\mathbf Z^\top=\mathbf B$, where $\mathbf Q$ and $\mathbf Z$ are orthogonal matrices. (For complex matrices, replace transpose with conjugate transpose, quasitriangular with triangular, and orthogonal with unitary.)

(If you are reading Golub and Van Loan, note that my convention is slightly different from theirs.)

If the resulting pencil is triangular, the eigenvalues of the pencil can be obtained as the ratio of the diagonal elements of $\mathbf T$ and $\mathbf U$. Zero diagonal elements in $\mathbf U$ correspond to infinite eigenvalues of the pencil.

Let's look first at what MATLAB does:

A = [3 2 1; 2 3 1; 5 1 2];

B = [1 2 3; 4 5 6; 7 8 9];

[at, bt, q, z] = qz(A, B)

at =

   -2.1394    6.0702    2.8152
         0    0.8982    1.8737
         0         0   -2.0817


bt =

    1.9287   14.3282    2.2505
         0    8.2195    1.8326
         0         0         0


q =

   -0.6272   -0.5758   -0.5244
    0.6633   -0.0420   -0.7472
    0.4082   -0.8165    0.4082


z =

    0.4729   -0.4012   -0.7845
    0.3630   -0.7226    0.5883
   -0.8029   -0.5630   -0.1961

Here is the Mathematica equivalent:

amat = N[{{3, 2, 1}, {2, 3, 1}, {5, 1, 2}}];
bmat = N[{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}];

MatrixForm /@ ({q, at, z, bt} = SchurDecomposition[{amat, bmat}])

$$\left\{\begin{pmatrix} -0.627226 & 0.663265 & 0.408248 \\ -0.575823 & -0.0419709 & -0.816497 \\ -0.524419 & -0.747207 & 0.408248 \\ \end{pmatrix}, \begin{pmatrix} -2.13941 & 6.07017 & -2.81516 \\ 0. & 0.898164 & -1.87371 \\ 0. & 0. & 2.08167 \\ \end{pmatrix}, \begin{pmatrix} 0.47294 & -0.401178 & 0.784465 \\ 0.362966 & -0.722566 & -0.588348 \\ -0.802859 & -0.562988 & 0.196116 \\ \end{pmatrix}, \begin{pmatrix} 1.92866 & 14.3282 & -2.25052 \\ 0. & 8.21954 & -1.83255 \\ 0. & 0. & 0. \\ \end{pmatrix}\right\}$$

Note the similarities and differences in the output. For MATLAB, we have the relations q.' * at * z.' == A and q.' * bt * z.' == B, while for Mathematica, we have q.at.Transpose[z] == amat and q.bt.Transpose[z] == bmat.


Additionally, MATLAB's qz() provides a form where the (right and left) eigenvectors of the given pencil are also produced. In Mathematica, Eigenvectors[] (and thus Eigensystem[] as well) can work on a pencil. Using the same matrices as above in both MATLAB and Mathematica:

[at, bt, q, z, vl, wl] = qz(A, B)

...

vl =

    0.5891    0.2404   -0.5000
    0.4521   -0.0128    1.0000
   -1.0000   -1.0000   -0.5000


wl =

   -1.0000    1.0000    0.5000
    0.2915   -0.7044   -1.0000
    0.1749   -0.4226    0.5000

where I had omitted the output of the other four matrices. In Mathematica, the equivalent for vl is

MatrixForm[vecs = Eigenvectors[{amat, bmat}]]

$$\begin{pmatrix} 0.408248 & -0.816497 & 0.408248 \\ 0.47294 & 0.362966 & -0.802859 \\ 0.233732 & -0.0124306 & -0.972221 \\ \end{pmatrix}$$

It looks a little different from what MATLAB produced, but we can verify that certain properties are satisfied. In MATLAB:

adi = diag(diag(at));
bdi = diag(diag(bt));

norm(A * vl * bdi - B * vl * adi)

ans =

   1.3591e-14

and in Mathematica:

amat.Transpose[vecs[[{2, 3, 1}]]].DiagonalMatrix[Diagonal[bt]] - 
bmat.Transpose[vecs[[{2, 3, 1}]]].DiagonalMatrix[Diagonal[at]]
   7.65362*10^-15

(Note that I needed to permute the eigenvectors to agree with the order in the generalized Schur decomposition.)

As for wl, that is the set of left eigenvectors for the pencil. Unfortunately, in Mathematica, there is no convenient way to get the left eigenvectors, so you have to do a second computation:

MatrixForm[lvecs = Eigenvectors[Transpose /@ {amat, bmat}]]

$$\begin{pmatrix} 0.408248 & -0.816497 & 0.408248 \\ -0.946782 & 0.276006 & 0.165603 \\ -0.772724 & 0.544287 & 0.326572 \\ \end{pmatrix}$$

Check in MATLAB:

norm(bdi * wl.' * A - adi * wl.' * B)

ans =

   1.0668e-14

and in Mathematica:

Norm[DiagonalMatrix[Diagonal[bt]].lvecs[[{2, 3, 1}]].amat - 
     DiagonalMatrix[Diagonal[at]].lvecs[[{2, 3, 1}]].bmat]
   5.36595*10^-15
| improve this answer | |
$\endgroup$
  • 2
    $\begingroup$ Perhaps a naive question, but where does the term pencil come from here? $\endgroup$ – b3m2a1 Mar 21 '18 at 7:06
  • 3
    $\begingroup$ "Pencil" is actually a more general term, but most conventionally applied to a pair of matrices. It comes from projective geometry; see this MO thread. $\endgroup$ – J. M.'s discontentment Mar 21 '18 at 7:11
  • $\begingroup$ @J.M. Thanks for the sophisticated answer. $\endgroup$ – qahtah Mar 22 '18 at 16:49

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.