Singular Matrix Pencils and the QZ Algorithm, Update - Cleve Moler on Mathematics and Computing

[poster]

(The January 5 posting was premature and incomplete.)

This year, 2023, is the 50-th anniversary of the QZ algorithm for the generalized matrix eigenvalue problem,

Ax = λBxThe algorithm avoids inverting either A or B. And, importantly, the QZ algorithm can be used to detect and analyze exceptional instances of the problem known as singular pencils. These pencils do not occur in the standard eigenvalue problem where B is the identity matrix.

Contents

• Singular pencils
• 3-by-3 example
• Wilkinson example
• References

Singular pencils

A matrix pencil is singular if both A and B are singular and, moreover, A - λ B is singular for all λ. Equivalently,

det(A - λB) = 0 for all λ.Singular pencils are more insidious than might appear at first glance. In some sense, the spectrum is the entire complex plane.

There are no singular pencils with the standard eigenvalue problem

Ax = λxwhere B is the identity matrix and certainly nonsingular.

If A and B are both real symmetric or complex Hermitian, but neither is positive definite, the pencil may or may not be singular. Symmetric problems occur frequently and there are separate algorithms and perturbation and convergence theories.

The QZ algorithm does not compute the eigenvalues λ until the very last step. It stably reduces A and B to triangular form with diagonals α and β. The eigenvalues are then the ratios

λ = α./βIsolated zeros in α or β yield zero or infinite eigenvalues with no special difficulties. Singular pencils occur if and only if zeros in α and β occur with the same index and (with exact arithmetic) lead to λ = 0/0 = NaN.

With a singular pencil, some, perhaps all, of the diagonal values in α and β are highly sensitive to perturbation, and all of the eigenvalues computed from the ratios are suspect. Theoretically the ill-conditioned eigenvalues can be determined from the Kronecker Canonical Form, a generalization of the notorious Jordan Canonical Form. But, like the Jordan Form, the Kronecker Form cannot provide a stable numerical algorithm.

3-by-3 example

A = [9 8 7; 6 5 4; 3 2 1] B = [1 3 2; 4 6 5; 7 9 8]A = 9 8 7 6 5 4 3 2 1B = 1 3 2 4 6 5 7 9 8Let's verify that this is a singular pencil. Use the Symbolic Math Toolbox to introduce a free variable.

syms lambda AB = A - lambda*B AB = [ 9 - lambda, 8 - 3*lambda, 7 - 2*lambda][6 - 4*lambda, 5 - 6*lambda, 4 - 5*lambda][3 - 7*lambda, 2 - 9*lambda, 1 - 8*lambda] Without further computation, we can see that the second row is the average of the first and third rows for all lambda and consequently that the determinant must be identically zero.

With exact arithmetic, each of these statements would produce the same eigenvalues. After the introduction of some roundoff error, two of the lambdas are indeterminant, but lambda = -1 is present in all four results. Is lambda = -1 a stable eigenvalue?

lambda_AB = eig(A,B) lambda_BA = 1./eig(B,A) lambda_ATBT = eig(A',B') lambda_BTAT = 1./eig(B',A')lambda_AB = 1.8984 -1.0000 -0.0807lambda_BA = -1.0000 0.5837 -0.9274lambda_ATBT = 0.0829 Inf -1.0000lambda_BTAT = -0.9661 0 -1.0000The triangular matrices for lambda_AB.

[QAZ,QBZ] = qz(A,B); QAZ QBZQAZ = 1.6131 10.2664 -11.0905 0 -4.2969 5.9613 0 0 -0.0000QBZ = 0.7898 6.8901 -13.5242 0 4.2969 -5.9613 0 0 0.0000Careful examination of the diagonals reveals that alfa(2)/beta(2) is producing the -1, while alfa(3)/beta(3) is roundoff over roundoff.

format long e alfa = diag(QAZ) beta = diag(QBZ) format shortalfa = 1.613087771308989e+00 -4.296911800112353e+00 -1.965207685813115e-15beta = 7.898460671891234e-01 4.296911800112357e+00 1.359052275299816e-15Wilkinson example

Jim Wilkinson published a survey paper about QZ and Kronecker products in 1979. One of his examples is

A = [4 3 2 5; 6 4 2 7; -1 -1 -2 -2; 5 3 2 6] B = [2 1 3 4; 3 3 3 5; 0 0 -3 -2; 3 1 3 5]A = 4 3 2 5 6 4 2 7 -1 -1 -2 -2 5 3 2 6B = 2 1 3 4 3 3 3 5 0 0 -3 -2 3 1 3 5Use the Symbolic Math Toolbox to verify that this is a singular pencil.

syms lambda AB = A - lambda*B AB = [4 - 2*lambda, 3 - lambda, 2 - 3*lambda, 5 - 4*lambda][6 - 3*lambda, 4 - 3*lambda, 2 - 3*lambda, 7 - 5*lambda][ -1, -1, 3*lambda - 2, 2*lambda - 2][5 - 3*lambda, 3 - lambda, 2 - 3*lambda, 6 - 5*lambda] The determinant is identically zero for all lambda.

d = det(AB) d = 0 With exact arithmetic, each of these statements would produce the same eigenvalues, but in practice each set is different. None of the eigenvalues is stable.

lambda_AB = eig(A,B) lambda_BA = 1./eig(B,A) lambda_ATBT = eig(A',B') lambda_BTAT = 1./eig(B',A')lambda_AB = 1.2056 0.7055 -1.0000 -Inflambda_BA = 1.5097 0.6408 0 -1.0000lambda_ATBT = -0.2141 + 0.2033i -0.2141 - 0.2033i 0.7013 + 0.0000i 1.4508 + 0.0000ilambda_BTAT = 0.3168 0.9823 1.2325 0The triangular matrices for lambda_AB are

[QAZ,QBZ] = qz(A,B); QAZ QBZQAZ = 0.7437 4.1769 -12.7279 -5.5000 0 0.0000 5.2328 2.1602 0 0 0.7857 0.0123 0 0 0 -0.2887QBZ = 0.5005 6.6143 -8.4853 -2.5000 0 0.0000 3.2668 2.0105 0 0 1.1525 -0.7904 0 0 0 0.2887Examine the diagonals more carefully. alfa(2)/beta(2) is the only roundoff over roundoff, but all four eigenvalues are unstable.

format long e alfa = diag(QAZ) beta = diag(QBZ) format shortalfa = 7.437114999643711e-01 1.216947725307920e-14 7.857314232211017e-01 -2.886751345948121e-01beta = 5.005405248737872e-01 1.021080292327182e-13 1.152509249099882e+00 2.886751345948153e-01References

C. B. Moler and G. W. Stewart, "An Algorithm for Generalized Matrix Eigenvalue Problems", SIAM J. Numerical Analysis, Vol.10, No.2, April 1973. Also available at cbm_gws.pdf

J. H. Wilkinson, "Kronecker's Canonical Form and the QZ Algorithm", Linear Algebra and its Applications, Vol. 28, 1979. Also available at wilkinson.pdf

Get the MATLAB code (requires JavaScript)

Published with MATLAB® R2023a