The error of a linear solve follows the condition number

Every computed answer here is compared with an exact one. The Hilbert systems are scaled to have integer entries, so the system stored in the computer is the exact system and its solution is exactly all ones; the Python reference confirms it by solving in rational arithmetic with fractions. The exact condition numbers come from mpmath at 80 digits, and numpy's linalg.cond agrees with them to $3\times 10^{-6}$ through $n=10$. The reference, posts/refs/the-error-follows-the-condition-number.py, also repeats every solve in Python in the same order of operations (output beside it).

Conditioning

The condition number $\kappa (A)={\sigma }_{\max }/{\sigma }_{\min }$, the ratio of the largest to the smallest singular value, bounds how much a relative change in $b$ or $A$ can change the solution of $Ax=b$: by up to $\kappa$ times as much. Rounding a number to double precision changes it, relative to its size, by at most half of $\epsilon =2^{-52}\approx 2.2\times 10^{-16}$ (machine epsilon, the gap between 1 and the next double), so no method can promise better than about $\kappa \epsilon$.

A method is backward stable if its answer is the exact solution of a nearby system, one whose data differ from $A$ and $b$ by a few $\epsilon$ relative. Its residual $b-A\hat{x}$ is then a few $\epsilon$ relative to $\Vert A\Vert \Vert \hat{x}\Vert$, and its error is about $\kappa \epsilon$. Gaussian elimination with partial pivoting (GEPP) behaves this way in practice. So a small residual says the method did its job, and it says nothing about the error when $\kappa$ is large.

The Hilbert matrix $H_{ij}=1/(i+j+1)$ is the standard ill-conditioned example; $\kappa$ grows about 32 times with each row added. Multiplying by $L=\mathrm{lcm}(1,\dots ,2n-1)$ makes every entry an integer below $2^{53}$, stored exactly, with the same $\kappa$. The right side $b=H(1,\dots ,1)$ is a sum of integers, also exact.

The condition number has to be computed too, and it runs into the same limit. eig_sym finds the smallest eigenvalue to within about $\epsilon$ times the largest, so $\kappa$ comes out with a relative error below $\kappa \epsilon$: $1.1\times 10^{-15}$ at $n=3$ and 4.4% at $n=12$, where numpy's cond is off by 1.6%. The tests use the exact values.

In least squares the method sets the condition number

A least-squares fit asks for the $c$ that minimizes $\Vert Vc-y\Vert$ for a tall matrix $V$. Here $V$ is the Vandermonde matrix of a polynomial fit, $V_{ik}=t_i^k$, at $m=40$ points $t_i=i/39$, and the data are $y=V(1,\dots ,1)$, so the coefficients sought are all ones.

The normal equations $V^{T}Vc=V^{T}y$ are the textbook route: a square system, solved here with the same elimination. But $\kappa (V^{T}V)=\kappa (V{)}^2$, so the error is about ${\kappa }^2\epsilon$. Householder QR multiplies $V$ by reflections $I-2v{v}^T/v^{T}v$, each of which zeroes one column below the diagonal, until $V$ is upper triangular, $R$. Reflections are orthogonal and change no lengths, so the problem becomes $Rc=Q^{T}y$ with $\kappa (R)=\kappa (V)$, and the error is about $\kappa \epsilon$.

The data are rounded when they are formed, which moves the exact least-squares answer of the stored problem slightly away from all ones. The Python reference computes that answer in rational arithmetic: it is at most $2.2\times 10^{-8}$ from ones (at degree 13), a small part of the errors below. The singular values for $\kappa (V)$ come from one-sided Jacobi (lib/svd).