(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