Cleve Moler is the author of the first MATLAB, one of the founders of MathWorks, and is currently Chief Mathematician at the company. He is the author of two books about MATLAB that are available online. He writes here about MATLAB, scientific computing and interesting mathematics.

Contents

David Hilbert

Around the turn of the 20th century, David Hilbert was the world's most famous mathematician. He introduced the matrix that now bears his name in a paper in 1895. The elements of the matrix, which are reciprocals of consecutive positive integers, are constant along the antidiagonals.

Least squares

Hilbert was interested in this matrix because it comes up in the least squares approximation of a continuous function on the unit interval by polynomials, expressed in the conventional basis as linear combinations of monomials.

$$ p(x) = \sum_{j=1}^n c_j x^{j-1} $$

The coefficient matrix for the normal equations has elements

$$ \int_0^1 x^{i+j-2} dx \ = \ \frac{1}{i+j-1} $$

Properties

A Hilbert matrix has many useful properties.

Symmetric.

Positive definite.

Hankel, $a_{i,j}$ is a function of $i+j$.

Cauchy, $a_{i,j} = 1/(x_i - y_j)$.

Nearly singular.

Each column of a Hilbert matrix is nearly a multiple of the other columns. So the columns are nearly linearly dependent and the matrix is close to, but not exactly, singular.

hilb

MATLAB has always had functions hilb and invhilb that compute the Hilbert matrix and its inverse. The body of hilb is now only two lines.

J = 1:n;
H = 1./(J'+J-1);

We often cite this as a good example of singleton expansion. A column vector is added to a row vector to produce a matrix, then a scalar is subtracted from that matrix, and finally the reciprocals of the elements produce the result.

Inverse

It is possible to express the elements of the inverse of the Hilbert matrix in terms of binomial coefficients. For reasons that I've now forgotten, I always use $T$ for $H^{-1}$.

A checkerboard sign pattern with large integers in the inverse cancels the smooth surface of the Hilbert matrix itself.

T12 = invhilb(12);
spy(T12 > 0)

A log scale mitigates the jaggedness.

surf(sign(T12).*log(abs(T12)))
view(60,60)

Rookie experiment

Now using MATLAB, I am going to repeat the experiment that I did on the Burroughs 205 when I was still a rookie. I had just written my first program that used Gaussian elimination to invert matrices. I proceeded to test it by inverting Hilbert matrices and comparing the results with the exact inverses. (Today's code uses this utility function that picks out the largest element in a matrix.)

Wow! The error is $10^8$. That's a pretty big error. Can I trust my matrix inversion code? What went wrong?

But I knew the elements of the inverse are huge. We should be looking at relative error.

relerr = maxabs(E)/maxabs(T)

relerr =
1.4439e-04

OK. The relative error is $10^{-4}$. That still seems like a lot.

I knew that the Hilbert matrix is nearly singular. That's why I was using it. John Todd was one of the first people to write about condition numbers. An error estimate that involves nearness to singularity and the floating point accuracy would be expressed today by

esterr = cond(H)*eps

esterr =
0.0036

That was about all I understood at the time. The roundoff error in the inversion process is magnified by the condition number of the matrix. And, the error I observe is less than the estimate that this simple analysis provides. So my inversion code passes this test.

Deeper explanation

I met Jim Wilkinson a few years later and came to realize that there is more to the story. I'm not actually inverting the Hilbert matrix. There are roundoff errors involved in computing H even before it is passed to the inversion routine.

Today, the Symbolic Math Toolbox helps provide a deeper explanation. The 'f' flag on the sym constructor says to convert double precision floating point arguments exactly to their rational representation. Here's how it works in this situation. To save space, I'll print just the first column.

The elements of H that are not exact binary fractions become ratios of large integers. The denominators are powers of two; in this case $2^{54}$ and $2^{55}$. The numerators are these denominators divided by $3$, $5$, etc. and then rounded to the nearest integer. The elements of F are as close to the exact elements of a Hilbert matrix as we can get in binary floating point.

Let's invert F, using the exact rational arithmetic provided by the Symbolic Toolbox. (I couldn't do this in 1960.)

S = inv(F);

We now have three inverse Hilbert matrices, X, S, and T.

X is the approximate inverse computed with floating point arithmetic by the routine I was testing years ago, or by MATLAB inv function today.

S is the exact inverse of the floating point matrix that was actually passed to the inversion routine.

T is the exact Hilbert inverse, obtained from the binomial coefficient formula.

The error in the computed inverse comes from two sources -- generating the matrix in the first place and then computing the inverse. The first of these is actually larger than the second, although the two are comparable.

Note

6 CommentsOldest to Newest

Dear Professor Moler: I would like to recall an additional useful property of the Hilbert matrix, namely its total positivity. Interestingly, this property is also related to the work of John Todd and Olga Taussky, since the reference given by Nick Higham in his “Test Matrix Toolbox for MATLAB” for the total positivity of Cauchy and Hilbert matrices was

But maybe I didn’t make this point strongly enough — we can’t use these high accuracy algorithms, for two reasons. First, we only start with a floating point approximation to what would be a totally positive matrix. The roundoff errors made in generating the matrix in the first place have a bigger effect on the inverse than those generated during the inversion process. Second, to take advantage of totally positivity, it is necessary to have the representation as a product of exact, rational, bidiagonal matrices.

You made the remark that on the Burroughs 205, you did better on inverting the Hilber matrix than the condition number seemed to imply. cond(H) gives an upper estimate of the error, but this upper estimate can be way off. Suppose I have two matrices: the identity matrix and the diagonal matrix with elements 1,2,4,8,16…2^n; their condition number is 1 and 2^n respectively. I don’t see any reason why in the second case I would lose log10(2^n) digits.
Can you comment on when the condition number gives a tight estimate and whether there is a better estimator?

Good question, Michele. The answer deserves a proper explanation, so I’ll make it the subject of a blog post. I’ve already got something cooking for next week, so it will be a couple of weeks before I get to this topic. Thanks for asking. — Cleve