The benefit of LU decomposition over explicitly computing the inverse

Question

I'm going to teach a linear algebra course in the fall, and I want to motivate the topic of matrix factorizations such as the LU decomposition. A natural question one can ask is, why care about this when one already knows how to compute $A^{-1}$ and can use it to solve $Ax=b$ just by matrix-vector multiplication? The standard answer is that is that in practice one should (almost) never actually invert a matrix, because you will incur less roundoff error by backsubstituting a factorization like $LUx = b$ than by performing the multiplication $x=A^{-1}b$.

However, I recently came across the paper "How accurate is inv(A)*b?" by Druinsky and Toledo (2012) where they argue that the received wisdom is misleading, and a solution $x_{\text{inv}}=Vx$ obtained using an approximate inverse $V\approx A^{-1}$ is typically just as close to the true $x$ as a solution $x_{\text{LU}}$ obtained by backsubstitution. So I am no longer as sure about numerical accuracy as a motivation for the LU decomposition as I used to be.

I did a few numerical experiments with random ill-conditioned matrices in Matlab, based on Druinsky and Toledo's code in Sec. 6 of their paper. It seems like the forward errors $\|x_{\text{inv}}-x\|$ and $\|x_{\text{LU}}-x\|$ were indeed usually quite similar, but the backward error $\|Ax_{\text{inv}}-b\|$ could be much bigger than $\|Ax_{\text{LU}}-b\|$. I also worked out a tiny example by hand: $$\begin{align} A &= \begin{bmatrix}1.00&2.01 \\ 1.01 & 2.03\end{bmatrix} \\ &= \underbrace{\begin{bmatrix}1.00 & 2.01 \\ 0 & -1.00\times10^{-4}\end{bmatrix}}_{L}\underbrace{\begin{bmatrix}1 & 0 \\ 1.01 & 1\end{bmatrix}}_{U} \\ &= {\underbrace{\begin{bmatrix}-2.03\times10^4 & 2.01\times10^4 \\ 1.01\times10^4 & -1.00\times10^4 \end{bmatrix}}_{A^{-1}}}^{-1}, \\ b &= \begin{bmatrix}1.01 \\ 1.02\end{bmatrix}. \end{align}$$ The exact solution is $x=[-1, 1]^T$. Computing $A^{-1}b$ with three decimal digits of precision in the intermediate calculations yields the result $x_{\text{inv}} = [0, 0]^T$ due to catastrophic cancellation. The same precision in LU backsubstitution gives $x_{\text{LU}} = [1.01, 0]^T$. Both solutions clearly differ from the true solution $x$ by an error on the order of $10^0$, but $Ax_{\text{LU}}$ matches $b$ to numerical precision while $Ax_{\text{inv}}$ is totally off.

My questions are:

1. Is the explicit $2\times2$ example above representative of what happens in practice, or is it a meaningless example computed with too little precision that just coincidentally happens to match the other numerical experiments?

2. Are the numerical observations generally true? To make this precise, let $A=LU$ and define $V = U^{-1}*L^{-1}$, $x_{\text{inv}} = V*b$, and $x_{\text{LU}}= U^{-1}*(L^{-1}*b)$, where $*$ denotes multiplication with numerical roundoff error. Is it usually the case that $\|Ax_{\text{inv}}-b\|>\|Ax_{\text{LU}}-b\|$?

3. What are some examples of practical applications where forward error is more important than backward error, and vice versa?

Answer 1

Let $B * y$ denote matrix multiplication performed in finite precision arithmetics. Let $x_{LU}$ be the solution of $A x = b$ obtained by the LU method, and $x_{inv}$ be the solution computed as $x_{inv} = \text{inv}(A) * b$. Let $\text{inv}(A)$ denote the inverse of $A$ computed from the LU decomposition $A \approx \hat{L}\hat{U}$ by solving $Ax_j = e_j$. Note that $\text{inv}(A)$ is different from $\hat{U}^{-1} * \hat{L}^{-1}$ but is much simpler to analyse.

The forward error for $x_{LU}$ is $$\tag{1}\|x - x_{LU}\|_\infty \leq c_n u \||A^{-1}| |\hat{L}| |\hat{U}| |x_{LU}| \|_\infty \leq c_n u \||A^{-1}| |\hat{L}| |\hat{U}|\|_\infty \times \|x_{LU} \|_\infty$$

where $c_n$ denotes a constant of order $n$ and $u$ is the machine precision. Similar forward error bounds for $x_{inv}$ is $$\tag{2} \|x - x_{inv}\|_{\infty} \leq c_n u \||A^{-1}| |\hat{L}| |\hat{U}| |\text{inv}(A)| |b|\|_{\infty} \leq c_n u \||A^{-1}| |\hat{L}| |\hat{U}|\|_\infty \times \|\text{inv}(A)| |b|\|_{\infty}$$

Comparing (1) and (2) and noting, that for sufficiently well conditioned matrix $A$, $$\tag{3} \|x_{LU}\|_\infty\approx \|x\|_\infty = \|A^{-1}b\|_\infty \leq \||A^{-1}||b|\|_\infty\approx \||\text{inv}(A)| |b|\|_\infty$$ we see, that as long as $\|x\|_\infty\approx \|A^{-1}||b\|_\infty$, then $x_{LU}$ and $x_{inv}$ have similar forward error. Notice, that this condition holds for general $b$.

Therefore forward error bounds (1) and (2) always prefer $x_{LU}$ over $x_{inv}$ but usually the forward error for $x_{inv}$ is not much worse.

Error bound (1) can be found in Higham's "Accuracy and Stability of Numerical Algorithms". Error bound (2) can be computed from forward error bound for $\text{inv}(A)$: $$|\text{inv}(A) - A^{-1}| \leq c_n u |A^{-1}| |\hat{L}| |\hat{U}| |\text{inv}(A)|$$ and from observation, that error introduces by matrix multiplication $\text{inv}(A) * b$ is always smaller than error associated with computing $\text{inv}(A)$.

Now suppose, that $\bar{x}_{inv} = A^{-1} *b$, i.e. then inverse $A^{-1}$ is computed exactly and the only source of error is matrix multiplication. Then residuals $b - A\bar{x}_{inv}$ satisfy $$\tag{4} |b - A \bar{x}_{inv}| \leq c_n u |A| |A^{-1}| |b|$$ for $x_{LU}$ we have $$\tag{5}|b - A x_{LU}| \leq c_n u |\hat{L}| |\hat{U}| |x_{LU}|$$ Since it is usually true, that $\||\hat{L}| |\hat{U}|\|_\infty \approx \|A\|_\infty$, therefore $\bar{x}_{inv}$ and $x_{LU}$ can have comparable residual norms only when $\|x_{LU}\|_\infty \approx \||A^{-1}| |b|\|_\infty$. On the other hand, from (3) we have (for sufficiently well conditioned $A$), $\|x_{LU}\|_\infty \leq \||A^{-1}| |b|\|_\infty$, therefore residuals for $x_{LU}$ are usually smaller than residuals for $x_{inv}$ even for well conditioned matrix $A$. If $A$ has large condition number, then $\text{inv}(A)$ is inaccurate, we have additional source of error in (4), therefore we can expect much larger residuals for $x_{inv}$.

Therefore:

1. Presented $2\times 2$ example is generic, i.e. usually $x_{inv}$ and $x_{LU}$ have similar forward error (even for ill conditioned matrix $A$) unless $|x|_\infty \ll \||A^{-1}| |b|\|_\infty$.

2. Residuals for $x_{inv}$ is usually bigger than for $x_{LU}$ unless $A$ is well conditioned and $|x|_\infty \approx \||A^{-1}| |b|\|_\infty$.

I am not aware of any practical applications where minimizing the backward error is the ultimate goal. However the backward error and the forward error are closely related according to inequality $$\text{forward_error} \leq \text{condition_number} \times \text{backward_error}$$

Since condition number does not depend on algorithm used to solve given problem, one may hope, that choosing an algorithm will smaller backward error will lead to lower forward error as well.

Answer 2

This is more a comment than an answer to the questions, but too long for the comment format

I answered here to a question about the use of inverse. In the answer I use the same article How accurate is inv(A)*b?.

I start with the consideration that use the inverse is not always bad, see also the comment in the question linked, but. From the article you obtain the two estimates report in the answer, and the better is the backward-stable. You can improve the stime for the inverse case but under some condition. So in a case with a general matrix, and no other info, it is even better to use solve the linear system, for example with LU. In my opinion numerical accuracy is still a good reason.

More than this there are other aspect to consider in the preference of the LU respect the inverse. For example the performance, ok use the matrix-vector multiplication but I must calculate the inverse. Or the consideration for the sparse case, where is not possible use an explicit inverse. (see comments and other answer in the question linked).