Singular Matrix Pencils and the QZ Algorithm, Update
(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 = λBx
The 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
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 = λx
where 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 1
B =
1 3 2
4 6 5
7 9 8
Let'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.0807
lambda_BA =
-1.0000
0.5837
-0.9274
lambda_ATBT =
0.0829
Inf
-1.0000
lambda_BTAT =
-0.9661
0
-1.0000
The triangular matrices for lambda_AB.
[QAZ,QBZ] = qz(A,B); QAZ QBZ
QAZ =
1.6131 10.2664 -11.0905
0 -4.2969 5.9613
0 0 -0.0000
QBZ =
0.7898 6.8901 -13.5242
0 4.2969 -5.9613
0 0 0.0000
Careful 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 short
alfa =
1.613087771308989e+00
-4.296911800112353e+00
-1.965207685813115e-15
beta =
7.898460671891234e-01
4.296911800112357e+00
1.359052275299816e-15
Wilkinson 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 6
B =
2 1 3 4
3 3 3 5
0 0 -3 -2
3 1 3 5
Use 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
-Inf
lambda_BA =
1.5097
0.6408
0
-1.0000
lambda_ATBT =
-0.2141 + 0.2033i
-0.2141 - 0.2033i
0.7013 + 0.0000i
1.4508 + 0.0000i
lambda_BTAT =
0.3168
0.9823
1.2325
0
The triangular matrices for lambda_AB are
[QAZ,QBZ] = qz(A,B); QAZ QBZ
QAZ =
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.2887
QBZ =
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.2887
Examine 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 short
alfa =
7.437114999643711e-01
1.216947725307920e-14
7.857314232211017e-01
-2.886751345948121e-01
beta =
5.005405248737872e-01
1.021080292327182e-13
1.152509249099882e+00
2.886751345948153e-01
References
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
Published with MATLAB® R2023a