hpc.social


High Performance Computing
Practitioners
and friends /#hpc
Share: 
This is a crosspost from   Cleve’s Corner: Cleve Moler on Mathematics and Computing Cleve Moler is the author of the first MATLAB, one of the founders of MathWorks, and is currently Chief Mathematician at the company. He writes here about MATLAB, scientific computing and interesting mathematics.. See the original post here.

Singular Matrix Pencils and the QZ Algorithm

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

Ax = λBx

The algorithm computes these eigevalues without inverting either A or B. And, the QZ-algorithm can help detect and analyze exceptional situaions known as singular pencils.

Contents

Matrix pencils

If A and B are two square matrices, the linear matrix pencil is the matrix-valued function

  A - λB

A pencil is regular if there is at least one value of λ for which A - λB if not singular. The pencil is singular if both A and B are singular and, moreover, A - λB is singular for all λ. In other words,

  det(A - λB) = 0 for all λ.

Singular pencils are more insiduous than migt appear at first glance.

Example

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

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

     9     8     7
     6     5     4
     3     2     1


B =

     7     9     8
     4     6     5
     1     3     2

   syms s

   AB = A - s*B

   d = det(AB)
 
AB =
 
[9 - 7*s, 8 - 9*s, 7 - 8*s]
[6 - 4*s, 5 - 6*s, 4 - 5*s]
[  3 - s, 2 - 3*s, 1 - 2*s]
 
 
d =
 
0
 
   eig1 = eig(A,B)

   eig2 = 1./eig(B,A)

   [QAZ,QBZ,Q,Z,V,W] = qz(A,B); QAZ, QBZ
eig1 =

   -0.4071
    1.0000
    0.2439


eig2 =

   -2.0000
    1.0000
    0.3536


QAZ =

   -1.0298  -13.0363    7.7455
         0    5.6991   -4.6389
         0         0    0.0000


QBZ =

    2.4396  -11.4948    9.6394
         0    5.6991   -4.6389
         0         0    0.0000

   eig3 = eig(A',B')

   eig4 = 1./eig(B',A')

   [QATZ,QBTZ,Q,Z,V,W] = qz(A',B'); QATZ, QBTZ
eig3 =

   -0.2169
       Inf
    1.0000


eig4 =

   -0.0738
         0
    1.0000


QATZ =

   -0.0000  -15.0218    6.8390
         0    2.6729   -2.2533
         0         0    0.5922


QBTZ =

    0.0000  -15.2578    7.1280
         0         0    1.0203
         0         0    0.5922

Wilkinson example

   clear

   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

   syms s

   AB = A - s*B

   d = det(AB)
 
AB =
 
[4 - 2*s,   3 - s, 2 - 3*s, 5 - 4*s]
[6 - 3*s, 4 - 3*s, 2 - 3*s, 7 - 5*s]
[     -1,      -1, 3*s - 2, 2*s - 2]
[5 - 3*s,   3 - s, 2 - 3*s, 6 - 5*s]
 
 
d =
 
0
 
   eig1 = eig(A,B)

   eig2 = 1./eig(B,A)

   [QAZ,QBZ,Q,Z,V,W] = qz(A,B); QAZ, QBZ
eig1 =

    1.2056
    0.7055
   -1.0000
      -Inf


eig2 =

    1.5097
    0.6408
         0
   -1.0000


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

   eig3 = eig(A',B')

   eig4 = 1./eig(B',A')

   [QATZ,QBTZ,Q,Z,V,W] = qz(A',B'); QATZ, QBTZ
eig3 =

  -0.2141 + 0.2033i
  -0.2141 - 0.2033i
   0.7013 + 0.0000i
   1.4508 + 0.0000i


eig4 =

    0.3168
    0.9823
    1.2325
         0


QATZ =

   0.1281 - 0.2434i   0.2665 + 0.0169i   0.2663 + 1.4905i   0.3721 + 3.5350i
   0.0000 + 0.0000i   0.0587 + 0.1116i   5.2603 - 1.6197i  12.7878 - 4.0110i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   4.1745 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.7572 + 0.0000i


QBTZ =

   0.9052 + 0.0000i   0.6130 - 0.6141i  -0.2443 + 0.8738i   1.2233 + 2.5485i
   0.0000 + 0.0000i   0.4150 + 0.0000i   3.5658 - 1.2114i   8.0696 - 2.2671i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   6.6127 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.5220 + 0.0000i

References

C. B. Moler and G. W. Stewart, "An Algorithm for Generalized Matrix Eigenvalue Problems", SIAM J.NUMER.ANAL. 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 APPPLICATIONS, Vol. 28, 1979. Also available at Also available at wilkinson.pdf


Get the MATLAB code

Published with MATLAB® R2023a