Linear Algebra Review: Hilbert Matrix

The Hilbert matrix, \({\bf H}\) is defined such that its elements are:

\[H_{ij} = (i + j -1)^{-1}\]

for \(i = 1, \ldots, N\) and \(j = 1, \ldots, N\).

Warning

This matrix is known to have a large condition number.

Here we will explore the influence of the large condition number by solving progressively larger linear systems.

Define a vector:

\[{\bf x}^{(N)} = (0, 1, \ldots, N-1)^T\]

and call the \(N\times N\) Hilbert matrix \({\bf H}^{(N)}\).

Define the righthand side of a linear system simply as:

\[{\bf b}^{(N)} = {\bf H}^{(N)} {\bf x}^{(N)}\]

Using our Gaussian elimination routine or a linear algebra library, solve the system

\[{\bf H}^{(N)} \tilde{\bf x} = {\bf b}^{(N)}\]

for \(N = 2, \ldots, 15\).

Define the error in your solution as

\[\epsilon = \max |\tilde{\bf x} - {\bf x}^{(N)} |\]

For what value of \(N\) does this error become \(\mathcal{O}(1)\)?—i.e. even the first digit in your computed solution is wrong.

We decided to implement this in C++. We’ll reuse the Array class and the other matrix routines we implemented previously in C++. They are:

• array.H : the Array class

• matmul.H : a maxtrix-vector product function for doing \({\bf b} = {\bf A}{\bf x}\)

• gauss.H : contains the Gaussian elimination solver for solving \({\bf A}{\bf x} = {\bf b}\)

Here’s the code we came up with:

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>

#include "array.H"
#include "gauss.H"
#include "matmul.H"

int main() {

    // loop over N
    for (int n = 2; n < 100; ++n) {

        // create Hilbert matrix of size NxN
        Array H(n, n);

        for (int irow = 0; irow < n; ++irow) {
            for (int jcol = 0; jcol < n; ++jcol) {
                H(irow, jcol) = 1.0 / static_cast<double>(irow + jcol + 1);
            }
        }

        // create vector x
        std::vector x(n, 0.0);
        for (int i = 0; i < n; ++i) {
            x[i] = static_cast<double>(i);
        }

        // create b = H x
        auto b = matmul(H, x);

        // solve H xtilde = b
        auto x_tilde = gauss_elim(H, b, true);

        // compute error | x - xtilde |

        double error{-1.0};
        for (int i = 0; i < n; ++i) {
            error = std::max(error, std::abs(x[i] - x_tilde[i]));
        }

        //std::cout << std::setprecision(10);

        //for (int i = 0; i < n; ++i) {
        //    std::cout << i << " " << x[i] << " " << x_tilde[i] << std::endl;
        //}

        // output
        std::cout << "for n = " << n << " " << "error = " << error << std::endl;

        // if error ~ 1 then exit
        if (error > 1.0) {
            break;
        }
    }

}

And our output:

for n = 2 error = 0
for n = 3 error = 2.22045e-16
for n = 4 error = 6.62581e-13
for n = 5 error = 1.11449e-11
for n = 6 error = 7.78773e-10
for n = 7 error = 6.92262e-09
for n = 8 error = 1.06671e-06
for n = 9 error = 7.72834e-06
for n = 10 error = 0.00195039
for n = 11 error = 0.0329536
for n = 12 error = 0.806736
for n = 13 error = 5.67127

We see that once we hit \(N = 12\) we are making nearly a 100% error.