Hilbert Residual vs Error

A MATLAB experiment on why a tiny residual does not necessarily mean an accurate solution.

A numerical solver can produce coefficients that almost perfectly satisfy a system of equations and still be far from the true answer. This project makes that gap visible with Hilbert matrices, a classic stress test for floating-point arithmetic.

The experiment uses the simple function $f(x)=1+x$, whose exact coefficients are known, and represents it in polynomial bases of increasing size. The function and its true answer stay the same, but the resulting Hilbert system becomes extremely ill-conditioned.

I first read the results as a comparison between Gaussian elimination and Cholesky. Looking at them again, I think the more interesting question is different: how long can the equations appear to be solved while the coefficients are already becoming unreliable?

Why this setup works

For a basis $1,x,\ldots,x^{n-1}$ on $[0,1]$, the projection equations are

$$\begin{aligned} Hc &= b, \\\ H_{ij} &= \int_0^1 x^{i+j-2}\,dx = \frac{1}{i+j-1}, \\\ b_i &= \int_0^1 x^{i-1}(1+x)\,dx = \frac{1}{i}+\frac{1}{i+1}. \end{aligned}$$

Since $1+x$ is already in the basis, the exact coefficient vector is known without running either solver:

$$c_*=(1,1,0,\ldots,0)^T.$$

That is the useful part of the setup. The mathematical answer never changes. Increasing $n$ only makes the numerical representation harder, so any loss of coefficient accuracy comes from finite-precision arithmetic rather than from a difficult target function.

What I compare

For each computed vector $\hat c$, I keep three measurements separate:

Question	Measurement
Does $\hat c$ nearly satisfy the equations?	$|b-H\hat c|_2 / |b|_2$
Did the solver recover the true coefficients?	$|\hat c-c_|2 / |c|_2$
Does the polynomial still reproduce $1+x$ on $[0,1]$?	maximum absolute error on 10,001 points in $[0,1]$

The residual answers only the first question. With a badly conditioned Hilbert matrix, very different coefficient vectors can produce almost the same right-hand side. Their errors can also cancel when the polynomial is evaluated, so an inaccurate coefficient vector may still draw a convincing curve.

What happens as the system grows

The canonical script covers $n=2$ through $n=51$. It uses Gaussian elimination without pivoting and a custom Cholesky routine followed by triangular solves. An SVD supplies a double-precision estimate of $\kappa_2(H)$.

The first version used H\b as the reference for solution error. Coming back to the code, I realized the exact vector $c_$ is already available, so H\b is another numerical result rather than ground truth. The current experiment measures forward error directly against $c_$.

Around the low teens, the curves stop looking like a normal solver comparison. The Gaussian residual can remain close to machine precision even while the coefficient estimates become unreliable. In the original Cholesky run, a computed pivot eventually became nonpositive and MATLAB continued in complex arithmetic. The current implementation stops that solve and marks the affected matrix size as breakdown instead.

The condition-number plot has its own limit. At large $n$, any flattening or fluctuation in this curve is a resolution limit, not evidence that the Hilbert matrices stop getting worse. Double precision can no longer resolve the smallest singular value reliably.

The degree-50 Gaussian case is the most counterintuitive part for me. Coefficient error and function error answer different questions, so the script reports both instead of judging the result from two overlapping lines.

Repository

File	Role
run_experiment.m	Canonical experiment and figures
Gausselim.m	Gaussian elimination with optional complete pivoting
Cholesky.m	Real Cholesky factorization with an explicit breakdown check
forsub.m, backsub.m	Triangular-system solvers
labo.mlx	Original MATLAB R2024b live experiment, retained as a historical record
labo.pdf	Five-page export of the original run, kept alongside the revised experiment

Run it

The original live script records MATLAB R2024b. The canonical script uses base MATLAB functions. From the repository root, run:

run("run_experiment.m")

The script reports the first Cholesky breakdown observed in the current run. It leaves the error, residual and condition arrays in the workspace and opens the three figures generated by the experiment.

What stayed with me is that an excellent residual can still answer the wrong question. It tells me that the equations are nearly satisfied, not that I recovered the right solution.