# Why can ill-conditioned linear systems be solved precisely?

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).

• My intuition is that solvability/accuracy of an Ax=b system is tied to the forcing vector b, not just the matrix A. Perhaps if b doesn't "probe" or "excite" the ill conditioned modes of A, then accurate solution remains possible. As a limiting example, A can be exactly singular (infinite condition number), yet Ax=b may still possess a solution, that can be computed accurately, if the forcing data b resides in the range of A. I admit this is pretty hand-wavy, which is why I only comment instead of answer. – rchilton1980 Sep 29 '17 at 17:37
• @rchilton1980 "yet Ax=b may still possess a solution" But that solution is not unique. And the examples I'm referring to possess a unique solution. – Zoltán Csáti Sep 29 '17 at 17:57
• That's a fair counterpoint - perhaps an artifact of picking infinite condition number (an exactly zero eigenvalue). However I think you can replace that zero eigenvalue with machine epsilon and my point still stands. (That is, the sytem has very large condition number, the system is nonsingular with a a unique solution, which we can compute very accurately provided b has no component along that tiny eigenpair). – rchilton1980 Sep 29 '17 at 18:11
• To be more specific, my thought experiment here is something like A = diag([1 1 1 1 1 eps]), b = [b1 b2 b3 b4 b5 0]. It is contrived, but I think it is suffcient to justify the original claim: "sometimes ill conditioned A's can be solved accurately for particular choices of b" – rchilton1980 Sep 29 '17 at 18:16
• Just give another example from Moler's blog blogs.mathworks.com/cleve/2015/02/16/… – percusse Oct 3 '17 at 11:17

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.

• Wouldn't this argument imply that it is possible (even if unlikely) to reach the worst-case bound of $\kappa_\infty(A)$ for the matrix in the example? AFAIU, based on my answer and based on ?getrs's documentation this shouldn't be possible. – Kirill Sep 29 '17 at 20:06
• @BrianBorchers Could you please elaborate why "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$." holds? From the example below it is logical, but I would need some formulas. Let the SVD of $A$ be $A = U \Sigma V^T$. In the first case, $A = \Delta b \sigma_1 v_1^T + \sum_{i=2}^N u_i \sigma_i v_i^T$. How to proceed? – Zoltán Csáti Sep 29 '17 at 21:16
• I haven't discussed round-off errors in the $A$ matrix, but the general effect is similar- unless you get really unlucky in the round-off errors, you typically do somewhat better than the pessimistic worst-case bound. – Brian Borchers Sep 29 '17 at 23:23
• (-1) The discussion of component-wise relative errors in the output is seriously misleading. – Kirill Oct 3 '17 at 0:16

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

• The condition numbers in the paper referred to by the OP are two-norm condition numbers. If you go back to reference [17] by ElBarbary you'll see that in the earlier paper these were two-norm condition numbers. Also, I setup the examples from this paper using DMsuite, and got nearly exactly the same 2-norm condition numbers as reported in the paper. – Brian Borchers Sep 29 '17 at 22:54
• The infinity norm condition norm numbers for these examples that I got using dmsuite and Chebyshev interpolation were similar in magnitude to the two-norm condition numbers. I don't think that the difference between 2-norm in infinity-norm condition numbers is that important for this particular example. – Brian Borchers Sep 29 '17 at 22:58
• I believe that the errors reported in the paper are absolute rather than relative errors (it doesn't make too much difference except for $\epsilon=0.01$, where the solution dips down close to 0. – Brian Borchers Sep 29 '17 at 23:19
• For $\epsilon=0.01$ and $N=1024$, the relative errors for the parts of the solution that are near 0 are huge, but the absolute errors are small. I agree that the paper was very vague about what condition number was used and about what the "errors" were exactly (relative or absolute errors.) – Brian Borchers Sep 29 '17 at 23:24
• @BrianBorchers I'm not sure what you mean: this isn't the difference between 2-norm and infty-norm condition numbers, but rather normwise- and component-wise condition numbers (component-wise relative perturbations in the input, not component-wise relative errors in the output as in your answer). – Kirill Sep 30 '17 at 17:42