Bohemian Matrices in the MATLAB® Gallery 2

Posted by Cleve Moler,

We will have a two-part minisymposium on "Bohemian Matrices" at ICIAM2019, the International Congress on Industrial and Applied Mathematics in Valencia, Spain, July 15-19. This is an outline of my talk.

Contents

Why "Bohemian"

The colorful name "Bohemian Matrices" is the responsibility of Rob Corless and Steven Thornton of the University of Western Ontario. Quoting the web site.

When this project was in its early stages, our focus was on random integer matrices of bounded height. We originally used the phrase "bounded height integer matrix eigenvalues" to refer to our work. This led to the acronym BHIME which isn't too far off from "Bohemian".

The MATLAB Gallery

The gallery is a collection of interesting matrices, maintained in MATLAB by Nick Higham. Many of them are Bohemian, or Bohemian wannabees. A catalog for the collection is available with

   help gallery_
 GALLERY Higham test matrices.
    [out1,out2,...] = GALLERY(matname, param1, param2, ...)
    takes matname, a string that is the name of a matrix family, and
    the family's input parameters. See the listing below for available
    matrix families. Most of the functions take an input argument
    that specifies the order of the matrix, and unless otherwise
    stated, return a single output.
 
    For additional information, type "help private/matname", where matname
    is the name of the matrix family.
 
    binomial    Binomial matrix -- multiple of involutory matrix.
    cauchy      Cauchy matrix.
    chebspec    Chebyshev spectral differentiation matrix.
    chebvand    Vandermonde-like matrix for the Chebyshev polynomials.
    chow        Chow matrix -- a singular Toeplitz lower Hessenberg matrix.
    circul      Circulant matrix.
    clement     Clement matrix -- tridiagonal with zero diagonal entries.
    compar      Comparison matrices.
    condex      Counter-examples to matrix condition number estimators.
    cycol       Matrix whose columns repeat cyclically.
    dorr        Dorr matrix -- diagonally dominant, ill-conditioned, tridiagonal.
                (One or three output arguments, sparse)
    dramadah    Matrix of ones and zeroes whose inverse has large integer entries.
    fiedler     Fiedler matrix -- symmetric.
    forsythe    Forsythe matrix -- a perturbed Jordan block.
    frank       Frank matrix -- ill-conditioned eigenvalues.
    gcdmat      GCD matrix.
    gearmat     Gear matrix.
    grcar       Grcar matrix -- a Toeplitz matrix with sensitive eigenvalues.
    hanowa      Matrix whose eigenvalues lie on a vertical line in the complex
                plane.
    house       Householder matrix. (Three output arguments)
    integerdata Array of arbitrary data from uniform distribution on
                specified range of integers
    invhess     Inverse of an upper Hessenberg matrix.
    invol       Involutory matrix.
    ipjfact     Hankel matrix with factorial elements. (Two output arguments)
    jordbloc    Jordan block matrix.
    kahan       Kahan matrix -- upper trapezoidal.
    kms         Kac-Murdock-Szego Toeplitz matrix.
    krylov      Krylov matrix.
    lauchli     Lauchli matrix -- rectangular.
    lehmer      Lehmer matrix -- symmetric positive definite.
    leslie      Leslie matrix.
    lesp        Tridiagonal matrix with real, sensitive eigenvalues.
    lotkin      Lotkin matrix.
    minij       Symmetric positive definite matrix MIN(i,j).
    moler       Moler matrix -- symmetric positive definite.
    neumann     Singular matrix from the discrete Neumann problem (sparse).
    normaldata  Array of arbitrary data from standard normal distribution
    orthog      Orthogonal and nearly orthogonal matrices.
    parter      Parter matrix -- a Toeplitz matrix with singular values near PI.
    pei         Pei matrix.
    poisson     Block tridiagonal matrix from Poisson's equation (sparse).
    prolate     Prolate matrix -- symmetric, ill-conditioned Toeplitz matrix.
    qmult       Pre-multiply matrix by random orthogonal matrix.
    randcolu    Random matrix with normalized cols and specified singular values.
    randcorr    Random correlation matrix with specified eigenvalues.
    randhess    Random, orthogonal upper Hessenberg matrix.
    randjorth   Random J-orthogonal (hyperbolic, pseudo-orthogonal) matrix.
    rando       Random matrix with elements -1, 0 or 1.
    randsvd     Random matrix with pre-assigned singular values and specified
                bandwidth.
    redheff     Matrix of 0s and 1s of Redheffer.
    riemann     Matrix associated with the Riemann hypothesis.
    ris         Ris matrix -- a symmetric Hankel matrix.
    sampling    Nonsymmetric matrix with integer, ill conditioned eigenvalues.
    smoke       Smoke matrix -- complex, with a "smoke ring" pseudospectrum.
    toeppd      Symmetric positive definite Toeplitz matrix.
    toeppen     Pentadiagonal Toeplitz matrix (sparse).
    tridiag     Tridiagonal matrix (sparse).
    triw        Upper triangular matrix discussed by Wilkinson and others.
    uniformdata Array of arbitrary data from standard uniform distribution
    wathen      Wathen matrix -- a finite element matrix (sparse, random entries).
    wilk        Various specific matrices devised/discussed by Wilkinson.
                (Two output arguments)
 
    GALLERY(3) is a badly conditioned 3-by-3 matrix.
    GALLERY(5) is an interesting eigenvalue problem.  Try to find
    its EXACT eigenvalues and eigenvectors.
 
    See also MAGIC, HILB, INVHILB, HADAMARD, PASCAL, ROSSER, VANDER, WILKINSON.

The matrices in the See also line were available before gallery was introduced. I think of them as part of the gallery.

gallery5

My favorite matrix in the gallery is

   G = gallery(5)
G =

          -9          11         -21          63        -252
          70         -69         141        -421        1684
        -575         575       -1149        3451      -13801
        3891       -3891        7782      -23345       93365
        1024       -1024        2048       -6144       24572

Let's compute its eigenvalues.

   format long e
   eig(G)
ans =

     -3.472940132398842e-02 + 2.579009841174434e-02i
     -3.472940132398842e-02 - 2.579009841174434e-02i
      1.377760760018629e-02 + 4.011025813393478e-02i
      1.377760760018629e-02 - 4.011025813393478e-02i
      4.190358744843689e-02 + 0.000000000000000e+00i

How accurate are these? What are the exact eigenvalues?

I will give readers of this post, and attendees at my talk, some time to contemplate these questions by delaying the answers until the end.

triw

The triw matrix demonstrates that Gaussian elimination cannot reliably detect near singularity.

   dbtype 1:2 private/triw
1     function T = triw(n, alpha, k, classname)
2     %TRIW   Upper triangular matrix discussed by Wilkinson and others.
   n = 12;
   T = gallery('triw',n)
T =

     1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1
     0     1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1
     0     0     1    -1    -1    -1    -1    -1    -1    -1    -1    -1
     0     0     0     1    -1    -1    -1    -1    -1    -1    -1    -1
     0     0     0     0     1    -1    -1    -1    -1    -1    -1    -1
     0     0     0     0     0     1    -1    -1    -1    -1    -1    -1
     0     0     0     0     0     0     1    -1    -1    -1    -1    -1
     0     0     0     0     0     0     0     1    -1    -1    -1    -1
     0     0     0     0     0     0     0     0     1    -1    -1    -1
     0     0     0     0     0     0     0     0     0     1    -1    -1
     0     0     0     0     0     0     0     0     0     0     1    -1
     0     0     0     0     0     0     0     0     0     0     0     1

T is already upper triangular. It is the U of its own LU-decomposition. Since there are no small pivots on the diagonal, we might be tempted to conclude that T is well conditioned. However, let

   x = 2.^(-(1:n-1))';
   format rat
   x(n) = x(n-1)
x =

       1/2     
       1/4     
       1/8     
       1/16    
       1/32    
       1/64    
       1/128   
       1/256   
       1/512   
       1/1024  
       1/2048  
       1/2048  

This vector is nearly a null vector for T.

   Tx = T*x
Tx =

       0       
       0       
       0       
       0       
       0       
       0       
       0       
       0       
       0       
       0       
       0       
       1/2048  

The 1-norm condition number of T is at least as large as

   norm(x,1)/norm(Tx,1)
ans =

    2048       

That is growing like 2^n.

moler

The moler matrix demonstrates that Cholesky decomposition of a symmetric, positive definite matrix cannot reliably detect near singularity either. The matrix can either be obtained from the gallery,

   M = gallery('moler',n);

or generated from T,

   format short
   M = T'*T
M =

     1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1
    -1     2     0     0     0     0     0     0     0     0     0     0
    -1     0     3     1     1     1     1     1     1     1     1     1
    -1     0     1     4     2     2     2     2     2     2     2     2
    -1     0     1     2     5     3     3     3     3     3     3     3
    -1     0     1     2     3     6     4     4     4     4     4     4
    -1     0     1     2     3     4     7     5     5     5     5     5
    -1     0     1     2     3     4     5     8     6     6     6     6
    -1     0     1     2     3     4     5     6     9     7     7     7
    -1     0     1     2     3     4     5     6     7    10     8     8
    -1     0     1     2     3     4     5     6     7     8    11     9
    -1     0     1     2     3     4     5     6     7     8     9    12

I was surprised when Nick named the matrix after me. The comments in the code indicate he found the matrix in a 30-year old book by a friend of mine, John Nash.

   dbtype 9:14 private/moler
9     %   Nash [1] attributes the ALPHA = -1 matrix to Moler.
10    %
11    %   Reference:
12    %   [1] J. C. Nash, Compact Numerical Methods for Computers: Linear
13    %   Algebra and Function Minimisation, second edition, Adam Hilger,
14    %   Bristol, 1990 (Appendix 1).

John doesn't remember when he heard about the example from me and I don't remember where or where I first came across it.

We know that T should be the Cholesky decomposition of M, but let's make sure.

   assert(isequal(chol(M),T))

We have already seen that, despite the fact that it has no small elements, T is badly conditioned. M must be at least as badly conditioned as T. In fact, it is far worse. A good null vector is provided by

   [~,v] = condest(M);

Let's look at v alongside our x. They are nearly the same.

   format long
   [v, x]
ans =

   0.500000178813913   0.500000000000000
   0.250000089406957   0.250000000000000
   0.125000223517391   0.125000000000000
   0.062500469386522   0.062500000000000
   0.031250949948913   0.031250000000000
   0.015626905485761   0.015625000000000
   0.007816313765488   0.007812500000000
   0.003913878927961   0.003906250000000
   0.001968383554413   0.001953125000000
   0.001007079958072   0.000976562500000
   0.000549316340766   0.000488281250000
   0.000366210893844   0.000488281250000

The 1-norm condition of M is at least as large as

   norm(v,1)/norm(M*v,1)
ans =

     2.796202999840670e+06

It turns out that is growing faster than 4^n.

kahan

But wait, you say, use column pivoting. Indeed, QR with column pivoting will detect the ill-condition of T or M. But the kahan matrix in the Gallery is a version of T, rescaled so that the columns all have nearly the same norm. Even column pivoting is fooled.

   format short
   K = gallery('kahan',7)
K =

    1.0000   -0.3624   -0.3624   -0.3624   -0.3624   -0.3624   -0.3624
         0    0.9320   -0.3377   -0.3377   -0.3377   -0.3377   -0.3377
         0         0    0.8687   -0.3148   -0.3148   -0.3148   -0.3148
         0         0         0    0.8097   -0.2934   -0.2934   -0.2934
         0         0         0         0    0.7546   -0.2734   -0.2734
         0         0         0         0         0    0.7033   -0.2549
         0         0         0         0         0         0    0.6555

rosser

I am very fond of the Rosser matrix.

   R = rosser
R =

   611   196  -192   407    -8   -52   -49    29
   196   899   113  -192   -71   -43    -8   -44
  -192   113   899   196    61    49     8    52
   407  -192   196   611     8    44    59   -23
    -8   -71    61     8   411  -599   208   208
   -52   -43    49    44  -599   411   208   208
   -49    -8     8    59   208   208    99  -911
    29   -44    52   -23   208   208  -911    99

Before the success of the QR algorithm, the Rosser matrix was a challenge for matrix eigenvalue routines. Today we can use it to show off the Symbolic Math toolbox.

   R = sym(R);
   p = charpoly(R,'x')
   f = factor(p)
   e = solve(p)
 
p =
 
x^8 - 4040*x^7 + 5080000*x^6 + 82518000*x^5 - 5327676250000*x^4 + 4287904631000000*x^3 - 1082852512000000000*x^2 + 106131000000000000*x
 
 
f =
 
[ x, x - 1020, x^2 - 1040500, x^2 - 1020*x + 100, x - 1000, x - 1000]
 
 
e =
 
                  0
               1000
               1000
               1020
 510 - 100*26^(1/2)
 100*26^(1/2) + 510
    -10*10405^(1/2)
     10*10405^(1/2)
 

parter

Several years ago, I was totally surprised when I happened to compute the SVD of a nonsymmetric version of the Hilbert matrix.

   format long
   n = 20;
   [I,J] = ndgrid(1:n);
   P = 1./(I - J + 1/2);
   svd(P)
ans =

   3.141592653589794
   3.141592653589794
   3.141592653589794
   3.141592653589793
   3.141592653589793
   3.141592653589792
   3.141592653589791
   3.141592653589776
   3.141592653589108
   3.141592653567518
   3.141592652975738
   3.141592639179606
   3.141592364762737
   3.141587705293030
   3.141520331851230
   3.140696048746875
   3.132281782811693
   3.063054039559126
   2.646259806143500
   1.170504602453951

Where did all those pi 's come from? Seymour Parter was able to explain the connection between this matrix and a classic theorem of Szego about the coefficients in the Fourier series of a square wave. The gallery has this matrix as gallery('parter',n).

pascal

I must include the Pascal triangle.

   help pascal_
 PASCAL Pascal matrix.
    P = PASCAL(N) is the Pascal matrix of order N. P a symmetric positive
    definite matrix with integer entries, made up from Pascal's triangle.
    Its inverse has integer entries. PASCAL(N).^r is symmetric positive
    semidefinite for all nonnegative r.
 
    P = PASCAL(N,1) is the lower triangular Cholesky factor (up to signs of
    columns) of the Pascal matrix. P is involutory, that is, it is its own
    inverse.
 
    P = PASCAL(N,2) is a rotated version of PASCAL(N,1), with sign flipped
    if N is even. P is a cube root of the identity matrix.
 
    Reference:
    N. J. Higham, Accuracy and Stability of Numerical Algorithms,
    Second edition, Society for Industrial and Applied Mathematics,
    Philadelphia, 2002, Sec. 28.4.

   P = pascal(12,2)
P =

    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1    -1
    11    10     9     8     7     6     5     4     3     2     1     0
   -55   -45   -36   -28   -21   -15   -10    -6    -3    -1     0     0
   165   120    84    56    35    20    10     4     1     0     0     0
  -330  -210  -126   -70   -35   -15    -5    -1     0     0     0     0
   462   252   126    56    21     6     1     0     0     0     0     0
  -462  -210   -84   -28    -7    -1     0     0     0     0     0     0
   330   120    36     8     1     0     0     0     0     0     0     0
  -165   -45    -9    -1     0     0     0     0     0     0     0     0
    55    10     1     0     0     0     0     0     0     0     0     0
   -11    -1     0     0     0     0     0     0     0     0     0     0
     1     0     0     0     0     0     0     0     0     0     0     0

With this sign pattern and arrangement of binomial coefficients, P is a cube root of the identity.

   I = P*P*P
I =

     1     0     0     0     0     0     0     0     0     0     0     0
     0     1     0     0     0     0     0     0     0     0     0     0
     0     0     1     0     0     0     0     0     0     0     0     0
     0     0     0     1     0     0     0     0     0     0     0     0
     0     0     0     0     1     0     0     0     0     0     0     0
     0     0     0     0     0     1     0     0     0     0     0     0
     0     0     0     0     0     0     1     0     0     0     0     0
     0     0     0     0     0     0     0     1     0     0     0     0
     0     0     0     0     0     0     0     0     1     0     0     0
     0     0     0     0     0     0     0     0     0     1     0     0
     0     0     0     0     0     0     0     0     0     0     1     0
     0     0     0     0     0     0     0     0     0     0     0     1

grcar

Joe Grcar will give us an interesting graphic.

   dbtype 7:10 private/grcar
7     %  References:
8     %    [1] J. F. Grcar, Operator coefficient methods for linear equations,
9     %    Report SAND89-8691, Sandia National Laboratories, Albuquerque,
10    %    New Mexico, 1989 (Appendix 2).
   n = 12;
   Z = gallery('grcar',n)
Z =

     1     1     1     1     0     0     0     0     0     0     0     0
    -1     1     1     1     1     0     0     0     0     0     0     0
     0    -1     1     1     1     1     0     0     0     0     0     0
     0     0    -1     1     1     1     1     0     0     0     0     0
     0     0     0    -1     1     1     1     1     0     0     0     0
     0     0     0     0    -1     1     1     1     1     0     0     0
     0     0     0     0     0    -1     1     1     1     1     0     0
     0     0     0     0     0     0    -1     1     1     1     1     0
     0     0     0     0     0     0     0    -1     1     1     1     1
     0     0     0     0     0     0     0     0    -1     1     1     1
     0     0     0     0     0     0     0     0     0    -1     1     1
     0     0     0     0     0     0     0     0     0     0    -1     1

   n = 100;
   Z = gallery('grcar',n);
   e = eig(Z);
   plot(e,'o')

gallery5

Did you compute the exact eigenvalues of gallery(5) yet? Remember, the matrix is

   G
G =

          -9          11         -21          63        -252
          70         -69         141        -421        1684
        -575         575       -1149        3451      -13801
        3891       -3891        7782      -23345       93365
        1024       -1024        2048       -6144       24572

MATLAB says the eigenvalues are

   e = eig(G)
e =

 -0.034729401323988 + 0.025790098411744i
 -0.034729401323988 - 0.025790098411744i
  0.013777607600186 + 0.040110258133935i
  0.013777607600186 - 0.040110258133935i
  0.041903587448437 + 0.000000000000000i

But, look at the powers of G. Since G has integer entries (it is Bohemian), its entries can be computed without roundoff. The fourth power is

   G^4
ans =

           0           0           0           0           0
         -84         168        -420        1344       -5355
         568       -1136        2840       -9088       36210
       -3892        7784      -19460       62272     -248115
       -1024        2048       -5120       16384      -65280

And the fifth power is

   G^5
ans =

     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0
     0     0     0     0     0

All zero.

The Cayley-Hamilton theorem tells us the characteristic polynomial must be

x^5

The eigenvalue, the only root of this polynomial, is zero with algebraic multiplicity 5.

MATLAB is effectively solving

x^5 = roundoff

Hence the computed eigenvalues come from

x = (roundoff)^(1/5)

It's hard to recognize them as perturbations of zero.


Get the MATLAB code

Published with MATLAB® R2018b

37 views (last 30 days)  | |

Comments

To leave a comment, please click here to sign in to your MathWorks Account or create a new one.