Example in MATLAB.

Let us consider the following $5 \times 5$ singular pencil $A - \lambda B$ with integer entries:

$$A = \left[ \begin{array}{ccccc} 22 & 34 & 31 & 31 & 17 \\ ... ... & 14 & 23 \\ 24 & 35 & 18 & 21 & 22 \\ \end{array}\right].$$

By construction the KCF of this example is ${\rm diag}( L_0, J_2(0), J_1(2), N_1, L^T_0 )$. We compute the GUPTRI form of $A - \lambda B$ by the call

$${\tt [S, T, P, Q, kstr] = guptri(A, B)};$$

i.e., we use the default values on EPSU, GAP, and ZERO. The output from this call is displayed below:

  S =
           0         0  -31.2179   69.7186 -142.6727
           0         0         0   15.8532  -42.1039
           0         0         0  -15.5639   -3.4712
           0         0         0         0   13.5979
           0         0         0         0         0

 T =
           0  -21.5942  -22.1660   44.5153 -110.9818
           0         0  -25.0581   19.3190  -43.8905
           0         0         0   -7.7820    9.2041
           0         0         0         0         0
           0         0         0         0         0

  P =
      0.2697    0.4520    0.8072    0.0226   -0.2660
      0.4045    0.5101   -0.4985    0.5400   -0.1900
      0.6742   -0.3818    0.1705    0.2136    0.5701
      0.1348    0.5620   -0.1771   -0.6188    0.5017
      0.5394   -0.2718   -0.1986   -0.5285   -0.5625

  Q =
     -0.4044    0.1510    0.4636    0.7253   -0.2697
      0.6985   -0.4097    0.1200    0.1975   -0.5394
      0.1470    0.6896   -0.5633    0.1481   -0.4045
     -0.4044   -0.0310    0.1859   -0.5886   -0.6742
     -0.4044   -0.5769   -0.6471    0.2582   -0.1348

  kstr =

       2     1     0    -1     2     0    -1     1    -1
       1     1     0    -1     1     0    -1     1    -1

The first part of kstr (all the columns to the left of the first -1) shall be interpreted as follows:

• The number of $L_{k-1}$ blocks is $n_k - r_k \equiv {\tt kstr}(1,k+1) - {\tt kstr}(2,k+1)$.
• The number of $J_k(0)$ blocks is $r_k - n_{k+1} \equiv {\tt kstr}(2,k) - {\tt kstr}(1,k+1)$.

The second part of kstr contains the corresponding information about $L^T_i$ and $N_i$ blocks and the third contains information about the size of the regular part. For this example we see that the computed GUPTRI form has the expected Kronecker structure, including one $L_0$ block ( ${\tt kstr}(1,1) - {\tt kstr}(2,1) = 1$), one $J_2(0)$ block ( ${\tt kstr}(2,2) - {\tt kstr}(1,3) = 1$), one $L_0^T$ block ( ${\tt kstr}(1,5) - {\tt kstr}(2,5) = 1$), one $N_1$ block ( ${\tt kstr}(2,5) - {\tt kstr}(1,6) = 1$), and finally a $1 \times 1$ regular part corresponding to the finite nonzero eigenvalue 2 ( ${\tt kstr}(1,8) = {\tt kstr}(2,8) = 1$). The $L_0$ block corresponds to $A$ and $B$ having a common column null space (both $S$ and $T$ have a first column with zero entries). Similarly, the $L^T_0$ block corresponds to $A$ and $B$ having a common row null space (both $S$ and $T$ have a last row with zero entries). Since $L_0$ and $L^T_0$ are the only singular blocks, $A_r - \lambda B_r$ is of size $0 \times 1$ and $A_l - \lambda B_l$ is of size $1 \times 0$. The remaining Kronecker structure in the GUPTRI form is contained in $S(\mbox{1:4, 2:5})$ and $T(\mbox{1:4, 2:5})$, i.e., rows 1 to 4 and columns 2 to 5 of $S$ and $T$: $A_z - \lambda B_z = S(\mbox{1:2, 2:3}) - \lambda T(\mbox{1:2, 2:3})$, $A_f - \lambda B_f = S(3,4) - \lambda T(3,4)$, and finally $A_i - \lambda B_i = S(4,5) - \lambda T(4,5)$.

The computed transformation matrices $P$ and $Q$ are orthogonal to machine precision accuracy. Moreover,

$$\delta = \Vert(A - A', B - B') \Vert _F = 2.8184e-13$$

is an upper bound on the distance to the closest $\{A + \delta A, B + \delta B\}$ with the KCF of $\{A', B'\}$, where $A' = PSQ^{\ast}$ and $B' = PTQ^{\ast}$. For this example, $\delta$ is of size $O( {\eps}\Vert(A, B) \Vert _F )$, where ${\eps}$ is the machine precision.

By setting ${\tt ZERO} = 0$ and making the call (with the default values on EPSU and GAP)

$${\tt [S, T, P, Q, kstr] = guptri(A, B, EPSU, GAP, ZERO)},$$

we perform a true equivalence transformation of $A - \lambda B$:

S =

  -2.9227e-15            0   3.1218e+01   6.9719e+01  -1.4267e+02
   2.0749e-15  -3.0175e-15   1.7377e-14   1.5853e+01  -4.2104e+01
  -5.4031e-15  -1.7747e-15   8.8818e-16  -1.5564e+01  -3.4712e+00
   1.6694e-15  -4.4157e-15   2.2204e-15            0   1.3598e+01
  -8.3206e-16   3.1071e-16   1.9984e-15  -2.2204e-16            0

T =

  -6.9809e-31   2.1594e+01   2.2166e+01   4.4515e+01  -1.1098e+02
   3.4328e-31            0   2.5058e+01   1.9319e+01  -4.3890e+01
   1.7764e-15            0            0  -7.7820e+00   9.2041e+00
  -5.5511e-16            0            0            0  -7.5221e-15
            0            0            0            0   3.1465e-15

The tiny nonzero elements in the lower triangular parts of $S$ and $T$ correspond to the singular values that are interpreted as numerically zero in the rank determination process of the GUPTRI algorithm. The results with the default value on ZERO correspond to the regularized pencil that has the computed $S$ and $T$ as its exact GUPTRI form with Kronecker structure as reported in kstr.

We end the discussion by exposing our seemingly harmless $5 \times 5$ example to the MATLAB function eig. The call [V,D] = eig(A,B) is supposed to compute a diagonal eigenvalue matrix $D$ and a full matrix $V$ whose columns are the corresponding eigenvectors so that $A V = B V D$. MATLAB computes a $D$ with the diagonal entries

        -1.8351e+16  
         2.0000e+00 
         7.2695e-01 - 4.1359e-25i             
        -6.2535e-16 + 2.4399e-08i
        -5.9077e-16 - 2.4399e-08i

and an eigenvector matrix $V$ with condition number $\Vert V \Vert \cdot \Vert V^{-1} \Vert = 8.8817e+08$. The large condition number and the residual $\Vert A V - B V D \Vert = 15.595$ signal that this is not a harmless example and further investigations are necessary. We conclude that it is only by software tools like guptri that we can get a more complete understanding of such ill-posed problems. And they do exist in real applications! Examples include controllability and observability issues in linear systems theory (see [447,120]).