# 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

## 2 CommentsOldest to Newest

**1**of 2

**2**of 2

## Recent Comments