Why can ill-conditioned linear systems be solved precisely?

Question

According to the answer here, large condition number (for linear system solving) decreases the guaranteed number of correct digits in the floating point solution. Higher order differentiation matrices in pseudospectral methods are typically very ill-conditioned. Why is it then that they are still very accurate methods?

I understand that the low precision coming from the ill-conditioned matrices is only a guaranteed value, but still it makes me wonder why ill-conditioned matrices are accurately solved by direct methods in practice - e.g., the LCOL columns of Table 3.1 on page 11 of Wang et al., A WELL-CONDITIONED COLLOCATION METHOD USING A PSEUDOSPECTRAL INTEGRATION MATRIX, SIAM J. Sci. Comput., 36(3).

Answer 1

Added after my initial answer:

It appears to me now that the author of the referenced paper is giving condition numbers (apparently 2-norm condition numbers but possibly infinity-norm condition numbers) in the table while giving maximum absolute errors rather than norm relative errors or maximum elementwise relative errors (these are all different measures.) Note that the maximum elementwise relative error is not the same thing as the infinity-norm relative error. Furthermore, the errors in the table are relative to the exact solution to the original differential equation boundary value problem rather than the discretized linear system of equations. Thus the information provided in the paper really isn't appropriate for use with the error bound based on the condition number.

However, In my replication of the computations, I do see situations where the relative infinity norm error (or two-norm relative error) is substantially smaller than the bound set by the infinity-norm condition number (respectively 2-norm condition number.) Sometimes you just get lucky.

I used the DMSUITE MATLAB package and solve the example problem from this paper using the pseudospectral method with Chebyshev polynomials. My condition numbers and maximum absolute errors were similar to those reported in the paper.

I also saw norm relative errors that were somewhat better than one might expect based on the condition number. For example, on the example problem with $\epsilon=0.01$, using $N=1024$, I get

cond(A,2)=7.9e+8

cond(A,inf)=7.8e+8

norm(u-uexact,2)/norm(uexact,2)=3.1e-12

norm(u-uexact,inf)/norm(uexact,inf)=2.7e-12

It would appear that the solution is good to about 11-12 digits, while the condition number is on the order of 1e8.

However, The situation with element-wise errors is more interesting.

max(abs(u-uexact))=2.7e-12

That still looks good.

max(abs((u-uexact)./uexact)=6.1e+9

Wow- there's a very large relative error in at least one component of the solution.

What happened? The exact solution of this equation has components that are tiny (e.g. 1.9e-22) while the approximate solution bottoms out at a much larger value of 9e-14. This is hidden by the norm relative error measurement (whether it be the 2-norm or infinity-norm) and only becomes visible when you look at the elementwise relative errors and take the maximum.

My original answer below explains why you can get a norm relative error in the solution that is less than the bound given by the condition number.

---

As you've noted in the question, the condition number, $\kappa(A)$, of a non-singular matrix gives a worst-case relative error bound for the solution to a perturbed system of equations. That is, if we solve $A(x+\Delta x)=b+\Delta b$ exactly and solve $Ax=b$ exactly, then

$\frac{\| \Delta x \|}{\| x \|} \leq \kappa(A) \frac{\| \Delta b \|}{\| b \|}$

Condition numbers can be computed with respect to a variety of norms, but the two-norm condition number is often used, and that's the condition number used in the paper that you refer to.

The worst case error occurs when $\Delta b$ is a left singular vector of $A$ corresponding to the smallest singular value of $A$. The best case occurs when $\Delta b$ is a left singular vector of $A$ corresponding to the largest singular value of $A$. When $\Delta b$ is random, then you have to look at the projections of $\Delta b$ onto all of the left singular vectors of $A$ and the corresponding singular values. Depending on the spectrum of $A$, things may go very badly or very well.

Consider two matrices $A$, both with 2-norm condition number $1.0 \times 10^{10}$. The first matrix has the singular values $1$, $1 \times 10^{-10}$, $\ldots$, $1 \times 10^{-10}$. The second matrix has singular values $1$, $1$, $\ldots$, $1$, $1 \times 10^{-10}$.

In the first case, a random perturbation is unlikely to be in the direction of the first left singular vector, and more likely to be close to one of the singular vectors with singular value $1 \times 10^{-10}$. Thus the relative change in the solution is likely to be very large. In the second case, almost any perturbation will be close in direction to a singular vector with singular value $1$, and the relative change in the solution will be small.

P.S. (added later after I got back from yoga class...)

The formula for the solution to $A\Delta x = \Delta b$ is

$\Delta x = V \Sigma^{-1} U^{T} \Delta b=\sum_{i=1}^{n} \frac{U_{i}^{T}\Delta b}{\sigma_{i}} V_{i}$

By the Pythagorean theorem,

$\| \Delta x \|_{2}^{2}= \sum_{i=1}^{n} \left( \frac{U_{i}^{T}\Delta b}{\sigma_{i}} \right)^{2} $

If we keep $\| \Delta b \|_{2}=1$, then this sum is maximized when $\Delta b=U_{n}$ and minimized when $\Delta b=U_{1}$.

In the situation considered here, $\Delta b$ is the result of random round-off errors, so the $U_{i}^{T}\Delta b$ values should all be of roughly the same magnitude. The terms with smaller values of $\sigma_{i}$ will contribute a lot to the error, while terms with larger values of $\sigma_{i}$ won't contribute much. Depending on the spectrum, this could easily be much smaller than the worst case bound.

Answer 2

tl;dr They reported a condition number, not necessarily the right condition number for the matrix, because there is a difference.

This is specific to the matrix and the right hand side vector. If you look at the documentation for *getrs, it says the forward error bound is $$ \frac{\|x-x_0\|_\infty}{\|x\|_\infty} \lesssim \mathrm{cond}(A,x)u \leq \mathrm{cond}(A)u. $$ Here $\mathrm{cond}(A,x)$ is not quite the usual condition number $\kappa_\infty(A)$, but rather $$ \mathrm{cond}(A,x) = \frac{\||A^{-1}||A||x|\|_\infty}{\|x\|_\infty},\\ \mathrm{cond}(A) = \||A^{-1}||A|\|. $$ (Here inside the norm these are component-wise absolute values.) See, for example, Iterative refinement for linear systems and LAPACK by Higham, or Higham's Accuracy and Stability of Numerical Algorithms (7.2).

For your example, I took a pseudospectral differential operator for a similar problem with $n=128$, and there is in fact a big difference between $\||A^{-1}||A|\|$ and $\kappa_\infty(A)$, I computed $7\times 10^3$ and $2.6\times 10^7$, which is enough to explain the observation that this happens for all right hand sides, because the orders of magnitudes roughly match what is seen in Table 3.1 (3-4 orders better errors). This doesn't work when I try the same for just a random ill-conditioned matrix, so it has to be a property of $A$.

An explicit example for which the two condition numbers don't match, which I took from Higham (7.17, p.124), due to Kahan is $$ \begin{pmatrix}2&-1&1\\-1&\epsilon&\epsilon\\1&\epsilon&\epsilon\end{pmatrix}, \qquad \begin{pmatrix}2+2\epsilon\\-\epsilon\\\epsilon\end{pmatrix}. $$ Another example I found is just the plain Vandermonde matrix on [1:10] with random $b$. I went through MatrixDepot.jl and some other ill-conditioned matrices also produce this type of result, like triw and moler.

Essentially, what's going on is that when you analyze the stability of solving linear systems with respect to perturbations, you first have to specify which perturbations you are considering. When solving linear systems with LAPACK, this error bound considers component-wise perturbations in $A$, but no perturbation in $b$. So this is different from the usual $\kappa(A) = \|A^{-1}\|\|A\|$, which considers normwise perturbations in both $A$ and $b$.

Consider (as a counterexample) also what would happen if you don't make the distinction. We know that using iterative refinement with double precision (see link above) we can get the best possible forward relative error of $O(u)$ for those matrices with $\kappa(A)\ll 1/u$. So if we consider the idea that linear systems can't be solved to accuracy better than $\kappa(A)u$, how would refining solutions possibly work?

P.S. It matters that ?getrs says the computed solution is the true solution of (A + E)x = b with a perturbation $E$ in $A$, but no perturbation in $b$. Things would be different if perturbations were allowed in $b$.

Edit To show this working more directly, in code, that this is not a fluke or a matter of luck, but rather the (unusual) consequence of two condition numbers being very different for some specific matrices, i.e., $$ \mathrm{cond}(A,x) \approx \mathrm{cond}(A) \ll \kappa(A). $$

function main2(m=128)
    A = matrixdepot("chebspec", m)^2
    A[1,:] = A[end,:] = 0
    A[1,1] = A[end,end] = 1
    best, worst = Inf, -Inf
    for k=1:2^5
        b = randn(m)
        x = A \ b
        x_exact = Float64.(big.(A) \ big.(b))
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        best, worst = min(best, err), max(worst, err)
    end
    @printf "Best relative error:       %.3e\n" best
    @printf "Worst relative error:      %.3e\n" worst
    @printf "Predicted error κ(A)*ε:    %.3e\n" cond(A, Inf)*eps()
    @printf "Predicted error cond(A)*ε: %.3e\n" norm(abs.(inv(A))*abs.(A), Inf)*eps()
end

julia> main2()
Best relative error:       2.156e-14
Worst relative error:      2.414e-12
Predicted error κ(A)*ε:    8.780e-09
Predicted error cond(A)*ε: 2.482e-12

Edit 2 Here is another example of the same phenomenon where the different conditions numbers unexpectedly differ by a lot. This time, $$ \mathrm{cond}(A, x) \ll \mathrm{cond}(A) \approx \kappa(A). $$ Here $A$ is the 10×10 Vandermonde matrix on $1:10$, and when $x$ is chosen randomly, $\mathrm{cond}(A,x)$ is noticably smaller than $\kappa(A)$, and the worst case $x$ is given by $x_i = i^a$ for some $a$.

function main4(m=10)
    A = matrixdepot("vand", m)
    lu = lufact(A)
    lu_big = lufact(big.(A))
    AA = abs.(inv(A))*abs.(A)
    for k=1:12
        # b = randn(m) # good case
        b = (1:m).^(k-1) # worst case
        x, x_exact = lu \ b, lu_big \ big.(b)
        err = norm(x - x_exact, Inf) / norm(x_exact, Inf)
        predicted = norm(AA*abs.(x), Inf)/norm(x, Inf)*eps()
        @printf "relative error[%2d]    = %.3e (predicted cond(A,x)*ε = %.3e)\n" k err predicted
    end
    @printf "predicted κ(A)*ε      = %.3e\n" cond(A)*eps()
    @printf "predicted cond(A)*ε   = %.3e\n" norm(AA, Inf)*eps()
end

Average case (almost 9 orders of magnitude better error):

julia> T.main4()
relative error[1]     = 6.690e-11 (predicted cond(A,x)*ε = 2.213e-10)
relative error[2]     = 6.202e-11 (predicted cond(A,x)*ε = 2.081e-10)
relative error[3]     = 2.975e-11 (predicted cond(A,x)*ε = 1.113e-10)
relative error[4]     = 1.245e-11 (predicted cond(A,x)*ε = 6.126e-11)
relative error[5]     = 4.820e-12 (predicted cond(A,x)*ε = 3.489e-11)
relative error[6]     = 1.537e-12 (predicted cond(A,x)*ε = 1.729e-11)
relative error[7]     = 4.885e-13 (predicted cond(A,x)*ε = 8.696e-12)
relative error[8]     = 1.565e-13 (predicted cond(A,x)*ε = 4.446e-12)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Worst case ($a=1,\ldots,12$):

julia> T.main4()
relative error[ 1]    = 0.000e+00 (predicted cond(A,x)*ε = 6.608e-13)
relative error[ 2]    = 1.265e-13 (predicted cond(A,x)*ε = 3.382e-12)
relative error[ 3]    = 5.647e-13 (predicted cond(A,x)*ε = 1.887e-11)
relative error[ 4]    = 8.895e-74 (predicted cond(A,x)*ε = 1.127e-10)
relative error[ 5]    = 4.199e-10 (predicted cond(A,x)*ε = 7.111e-10)
relative error[ 6]    = 7.815e-10 (predicted cond(A,x)*ε = 4.703e-09)
relative error[ 7]    = 8.358e-09 (predicted cond(A,x)*ε = 3.239e-08)
relative error[ 8]    = 1.174e-07 (predicted cond(A,x)*ε = 2.310e-07)
relative error[ 9]    = 3.083e-06 (predicted cond(A,x)*ε = 1.700e-06)
relative error[10]    = 1.287e-05 (predicted cond(A,x)*ε = 1.286e-05)
relative error[11]    = 3.760e-10 (predicted cond(A,x)*ε = 1.580e-09)
relative error[12]    = 3.903e-10 (predicted cond(A,x)*ε = 1.406e-09)
predicted κ(A)*ε      = 4.677e-04
predicted cond(A)*ε   = 1.483e-05

Edit 3 Another example is the Forsythe matrix, which is a perturbed Jordan block of any size of the form $$A = \begin{pmatrix} 0&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \\ \epsilon&0&0&0 \end{pmatrix}. $$ This has $\|A\|=1$, $\|A^{-1}\|=\epsilon^{-1}$, so $\kappa_\infty(A) = \epsilon^{-1}$, but $|A^{-1}| = A^{-1} = |A|^{-1}$, so $\mathrm{cond}(A) = 1$. And as can be verified by hand, solving systems of linear equations like $Ax=b$ with pivoting is extremely accurate, despite the potentially unbounded $\kappa_\infty(A)$. So this matrix too will yield unexpectedly precise solutions.

Edit 4 Kahan matrices are also like this, with $\mathrm{cond}(A)\ll \kappa(A)$:

A = matrixdepot("kahan", 48)
κ, c = cond(A, Inf), norm(abs.(inv(A))*abs.(A), Inf)
@printf "κ=%.3e c=%.3e ratio=%g\n" κ c (c/κ)

κ=8.504e+08 c=4.099e+06 ratio=0.00482027