Many recent developments in numerical algorithms concern improvements in performance and the exploitation of parallel architectures. However it's important not to lose sight of one crucial point: first and foremost, algorithms must be accurate. This begs the question, how do we know whether a routine is giving us the right answer? I'm asking this in the context of the matrix function routines I've been writing (these are found in chapter F01 of the NAG Library), but I'm sure you’ll agree that it’s an important question for all numerical software developers to consider.
First, it's important to note that the term "the right answer", is a little unfair. Given a matrix A stored in floating point, there is no guarantee that, for example, its square root A1/2 can be represented exactly in floating point arithmetic. In addition, rounding errors introduced early on can propagate through the computation so that their effects are magnified. We'd like to be able to take into account this sort of behaviour when we test our software.
For a matrix function f(A), the computed solution F can be written as F=f(A)+∆f, where ∆f is known as the forward error. In some cases it may be possible to compute the forward error, by obtaining f(A) via some other means (for example, analytically, or by using extended-precision arithmetic). Usually we are interested in the norm of ∆f , which we scale by the norm of f(A) to form a relative forward error. As an example, consider the following matrix, whose exponential is known explicitly:
The matrix A was used as a test matrix in the 2009 version of the scaling and squaring algorithm for the matrix exponential . With this algorithm, we obtain a relative forward error of 2.5x10-16. If we compute the exponentialï»¿ of A using the 2005 version of the scaling and squaring algorithm , then the relative forward error is about 3.0x10-12. Clearly the newer algorithm is giving a better result, but how do we decide whether these forward errors are small enough?
To try to answer this, it's useful to consider the backward error. We assume that the computed solution can also be written as F=f(A+∆A). Whereas the forward error tells me how far from the actual solution I am, the backward error ∆A tells me what problem I have actually solved.
A useful rule of thumb states that the forward error is approximately bounded by the product of the backward error and the condition number of the problem. Uncertainties in the input data, together with the limitations of floating point arithmetic, mean that a relative backward error in the region of unit roundoff is the best that we should expect. Thus if we can estimate the condition number, then we can get an idea of what size of forward error is acceptable. If a problem is poorly conditioned, then the forward error could be very large, even though the backward error is small and the algorithm is performing stably. At NAG we've been developing routines for estimating condition numbers of matrix functions, based on algorithms developed at the University of Manchester.
Returning to our matrix A, we find that the condition number of the exponential is very large: 1.6x1015. The product of the condition number and unit roundoff is about 1.7x10-1 so it looks like both algorithms were actually performing quite well.
The discussion above assumes that we are somehow able to compute forward errors, but this isn't always possible. One approach we been using at NAG is to test if certain matrix function identities hold. For example
are generalizations of well-known scalar identities. Of course, we are now back to our original problem: if I find that the computed residual in my identity is of the order of 10-14, is that acceptable or is it too large? We're currently working on some error analysis for such identities to answer this question, and hope to publish a paper on this soon.
So can I answer my original question; how do I know I'm getting the right answer? Well, I hope I've persuaded you that in numerical analysis, this is a thornier issue than you might have thought. One thing is clear though: testing a new library routine takes far longer than writing it in the first place!
 N.J. Higham The scaling and squaring method for the matrix exponential revisited, SIAM J. Matrix. Anal. Appl., 26 (2005), pp. 1179-1193.
 A.H. Al-Mohy and N.J. Higham A new scaling and squaring algorithm for the matrix exponential, SIAM J. Matrix. Anal. Appl., 26 (2009), pp. 970-989.