(This question might be legitimately crossposted to stackoverflow or mathoverflow or programming StackExchanges.)
Preface
I'm reading this paper on solving linear systems of equations Ax = b
using an explicit inverse (held by conventional wisdom as a big no-no): Druinsky & Toledo, "How accurate is inv(A)*b?", http://arxiv.org/abs/1201.6035
It asserts that in non-pathological cases, calculating a solution to Ax=b
via inv(A)*b
using the explicit inverse of a matrix produces solutions that are as accurate and as stable as preferred methods such as via LU factorization. It goes on to describe exactly which situations using "inv" might be bad (technically, when the right-hand side, b
, is nearly orthogonal to the subspace spanned by left-singular vectors of A
with low singular values). Even in these cases, the paper asserts, the solutions of inv(A)*b
are as accurate as preferred methods, they're just not backwards stable.
Question
My question is: what are the specific drawbacks of using an algorithm that isn't guaranteed to be backwards stable? If all I wanted was some solutions to Ax=b
, and I was guaranteed accurate results, does it matter that the solutions aren't backwards stable? Is there any specific examples of when it's ok to use an accurate but backwards-unstable algorithm?
Furthermore
My experimentation in Matlab/Octave shows me that the difference between a backwards-unstable and -stable algorithm is that slight perturbations to the solution result in bigger errors when going back through the linear system:
norm(A * (inv(A) * b) - b)
is much larger than
norm(A * (A \ b) - b)
for a "bad" b
(nearly orthogonal to the subspace of left-singular vectors with low singular values), where \
solves linear systems using Gaussian elimination (http://www.mathworks.com/help/matlab/ref/mldivide.html). The solutions themselves have the same accuracy, i.e., norm(A\b - true_x)
is about the same as norm(inv(A)*b - true_x)
.
I can imagine that in the case that
- I only cared about getting a solution to
Ax=b
and - knew I would never propagate the solution back through
A
,
could I justify using inv(A)
when I didn't want to bother checking that b
isn't "bad", i.e., without ensuring that the solution was backwards-stable. To me, this doesn't seem worth it, even without the speed and convenience advantages of using LU or Cholesky or QR decompositions.