# 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