An interface to Matlab in terms of a MEX-file is available for the routine GUPTRI. In the following we give a brief description of the arguments and output of the guptri function and demonstrate its usage on a small numerical example.

SYNOPSIS

guptri(A,B)

guptri(A,B,EPSU)

guptri(A,B,EPSU,GAP)

guptri(A,B,EPSU,GAP,ZERO)

[ S, T ] = guptri(...)

[ S, T, kstr ] = guptri(...)

[ S, T, P, Q ] = guptri(...)

[ S, T, P, Q, kstr ] = guptri(...)

The most general call [S, T, P, Q, kstr] = guptri(A, B, EPSU, GAP, ZERO) stores the transformation matrices in P and Q and P^HAQ and P^HBQ in S and T, respectively. The computed Kronecker structure is revealed by kstr, as described in the example dicussed later.

The guptri function reduces an m × n pencil A - zB to generalized upper triangular form S - zT = P^H(A - zB)Q:

S=                          T=
  A_r   *    *    *    *       B_r   *    *    *    *  
   0   A_z   *    *    *        0   B_z   *    *    *
   0    0   A_f   *    *        0    0   B_f   *    *
   0    0    0   A_i   *        0    0    0   B_i   *
   0    0    0    0   A_l       0    0    0    0   B_l

where the diagonal blocks of S and T in staircase forms describe the Kronecker structure of A - zB.

Besides (A, B), the user can optionally provide three input parameters (EPSU, GAP and ZERO) that control the computation of the GUPTRI form. EPSU (relative uncertainty in data) and GAP (should be at least 1 and nominally 1000) are used used to make rank decisions in order to determine the Kronecker structure of an input pencil. The default value is GAP = 1000, but other values may be necessary for certain examples. ZERO is used as a logical variable which when set to true (ZERO = nonzero value) forces guptri to z ero out small singular values during the reduction process so that the returned pencil really has the computed GUPTRI form. Otherwise, the returned pencil is a true equivalence transformation of the input pencil and all zero blocks in the GUPTRI form will contain small entries (normally of size EPSU at most).

Example

Let us consider the following 5 × 5 singular pencil A - zB with integer entries:

A =                      B = 
    22  34  31   31  17      13  26  25  17  24
    45  45  42   19  29      31  46  40  26  37 
    39  47  49   26  34      26  40  19  25  25 
    27  31  26   21  15      16  25  27  14  23 
    38  44  44   24  30      24  35  18  21  22 

By construction the KCF of this example is diag(L₀, J₂(0), J₁(2), N₁, L^T₀). We compute the GUPTRI form of A - zB by the call [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 (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 <=> kstr(1,k+1) - kstr(2,k+1).
• The number of J_k(0) blocks is r_k - n_k+1 <=> kstr(2,k) - 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₀ block (kstr(1,1) - kstr(2,1) = 1), one J₂(0) block (kstr(2,2) - kstr(1,3) = 1), one L^T(0) block (kstr(1,5) - kstr(2,5) = 1), one N₁ block (kstr(2,5) - kstr(1,6) = 1), and finally a 1 × 1 regular part corresponding to the finite nonzero eigenvalue 2 (kstr(1,8) = kstr(2,8) = 1). The L₀ block corresponds to that A and B have a commmon column nullspace (both S and T have a first column with zero entries). Similarly, the L^T₀ block corresponds to that A and B have a common row nullspace (both S and T have a last row with zero entries). Since L₀ and L^T₀ are the only singular blocks, A_r - zB_r is of size 0 × 1 and A_l - zB_l is of size 1 × 0. The remaining Kronecker structure in the GUPTRI form is contained in S(1:4, 2:5) and T(1:4, 2:5), i.e., rows 1 to 4 and columns 2 to 5 of S and T: A_z - zB_z = S(1:2, 2:3) - zT(1:2, 2:3), A_f - zB_f = S(3,4) - zT(3,4), and finally A_i - zB_i = S(4,5) - zT(4,5).

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

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^H and B' = PTQ^H. For this example, delta is of size O(max(||A||_F, ||B||_F)*eps), where eps is the relative machine accuracy.

By setting ZERO = 0 and making the call (with the default values on EPSU and GAP)

we perform a true equivalence transformation of A - zB:

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 presented previous 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 × 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 ||V|| * ||V^-1|| = 8.8817e+08. The large condition number and a large residual ||A V - B V D|| 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 problem. And they do exist in real applications! Examples include controllability and observability issues in linear systems theory.