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:

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?

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\|$?

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