Julia Issue Note #5705¶

Linear algebra test failures for linear solve with triangular matrices¶

Recent Julia issues like #5705 have adjusted the bounds on linear algebra tests. The purpose of this notebook is to provide more extensive numerical demonstrations of why the current linear algebra tests are fundamentally broken and need to be systematically replaced by more theoretically justified error bounds, as explained in #5605.

The current testing strategy¶

Here is a snippet of code representing a typical linear algebra test in Julia:

Mathematically, the test being performed corresponds to asserting that the inequality

$$\left\vert \left\Vert A x\right\Vert - \left\Vert b \right\Vert \right\vert \le C \epsilon$$

is valid, where $\epsilon$ is machine epsilon and $C = 20000$ is an arbitrary magic constant.

Most of the time, this test will appear to be well-behaved. The deviation between $Ax$ and $b$ appears to be small in norm, and it looks like the magic number 20000 would obviously cover any contingency. But how robust is this observation?

Distribution of magic numbers¶

The statistics show that every now and then, a nearly singular matrix is sampled, causing the magic number to blow up.

Where do the old and new bounds $15000\epsilon$ and $20000\epsilon$ lie on the distribution?

Despite the large change in the absolute magnitude of the magic constant, we see that this only shifts the threshold slightly on this histogram of $log(c)$, thanks to the extremely long tail in the magic number $c$. What is the probability that a given random matrix will fail the old and new cutoffs?

The change of the magic number in PR #5707 decreases the probability of failure by an essentially negligible amount.

Plotting the failure rate vs the magic constant $log(c)$ shows that in the grand scheme of things, the change in the magic constant does very little, and it would take an enormous shift in $c$ to make the failure rate negligible, so much so that the test would be practically useless.

Theoretical error bounds on linear solvers¶

Let's have a look at which method is being dispatched by the linear solver.

For a triangular matrix $A$, this either dispatches to the LAPACK xTRTRS subroutine, which in turn calls the substitution solver xTRSM, or the naïve substitution solver Base.LinAlg.naivesub!, depending on the element type.

What does theory tell us about the error bounds on substitution? You can derive an error bound using first-order matrix perturbation theory, or if you're lazy, consult a reference like Higham's Accuracy and Stability of Numerical Algorithms, Ch. 8 (pdf).

Error bound results usually come in two flavors. First, we have forward error bounds, which bound the norm of the difference between the computed solution $\hat{x}$ and the exact solution $x$.

For the triangular system, we have the forward error bound (Higham, Theorem 8.5):

$$\frac{\left\Vert x - \hat{x} \right\Vert_\infty}{\left\Vert x \right\Vert_\infty} \le \frac{\textrm{cond}(A,x) \gamma_n}{1 - \textrm{cond}(T) \gamma_n}$$

with constant

$$\gamma_n = \frac{n\epsilon}{1-n\epsilon}$$

and the Skeel condition numbers

$$\textrm{cond}(A,x) = \frac{\left\Vert \left\vert A \right\vert \left\vert A^{-1} \right\vert \left\vert x \right\vert \right\Vert_\infty}{\left\Vert x \right\Vert_\infty}$$

and

$$\textrm{cond}(A) = \left\Vert \left\vert A \right\vert \left\vert A^{-1} \right\vert \right\Vert_\infty$$

where $|A|$ means that the modulus is taken componentwise for every element of $A$.

Second, we also have backward error bounds, which assert that the computed solution $\hat{x} = A\backslash b$ is an exact solution to a slightly perturbed problem $(A + \delta A) \hat{x} = b + \delta b$. Sometimes results are also known if only $A$ or only $b$ is assumed to be perturbed.

For the triangular system, we have the backward error bound (Higham, Theorem 8.3)

$$\left\vert \delta A \right\vert \le \gamma_n \left\vert A \right\vert$$

where only the matrix $A$ is perturbed.

Forward error bound¶

Ordinarily, this error bound is useless in practice since the exact solution $x$ is not known. However, thanks to Julia's BigFloats, we can approximate the exact solution with that returned by big(A) \ big(b).