3
$\begingroup$

(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

  1. I only cared about getting a solution to Ax=b and
  2. 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.

$\endgroup$
1
  • 1
    $\begingroup$ You might ask the moderators to migrate your question to scicomp.stackexchange.com. I think where you care most about backward stability is when solving stiff ordinary differential equations. This query turned up some stuff that might (or might not) be relevant to you: scicomp.stackexchange.com/…. $\endgroup$ Commented Jun 10, 2013 at 14:49

1 Answer 1

1
$\begingroup$

Sivan Toledo responded to my query with the following:

In general, the significance and the broad-context meaning of mathematical results is a matter of opinion and debate, so perhaps not every researcher would agree to our explanation below.

In most of the literature, backward stability is a technical and/or analytic tool that makes it relatively easy to prove that an algorithm produces accurate results (or to bound the inaccuracy). In these contexts, the user of the numerical method worries about accuracy, not about backward stability, and the numerical analyst uses backward stability just as a mechanism to prove the bound that the user cares about.

(The analysis in our paper proves accuracy without relying directly on backward stability, but this is atypical.)

When an algorithm is backward stable, the user gets not only a bound on the magnitude of the error, but also an estimate of its direction; it lies in the direction of right singular vectors that correspond to small singular values (for the case of solving linear equations). In many cases this is completely meaningless to the user.

Backward stability may be meaningful to the end user if the parameters of the problem (say the coefficients of a linear equation) are only known approximately, say because they are measured or estimated using some limited-accuracy physical instrument. In some of these cases it may be good for the user to know that the output she/he got is an exact solution to a nearby problem.

I'd still like to better understand theoretically what the drawbacks of using inv(A)*b are when b is "bad". That is, if the solution using inv is as accurate as with a "better" algorithm, what really matters. (This is hinted at computationally in my question: the projection matrix A * inv(A' * A) * A' is at least unreliable.)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.