How do I know I'm getting the right answer?

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 [2]. 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 [1], 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-1so 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
exp(log(A))=A,
sin2A+cos2A=I,

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!

[1] N.J. Higham The scaling and squaring method for the matrix exponential revisited, SIAM J. Matrix. Anal. Appl., 26 (2005), pp. 1179-1193.
[2] 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.

1. Maybe you need more testing.

2. The forward error of 2.5e-16 was computed making use of the known result. For this high precision, can you still evaluate the algorithm forward and backward error using the rules for IEEE arithmetics? Or are you forced to compute the errors of specific matrices? Also you talk about size, which I take it means order of magnitude. Again, when you reach 1e-16, can one still assume that the coefficient (for the order of magnitude) is of order 1?

3. To get an accurate value for the forward error in general (as
opposed to for specific matrices with known exponential) you will
need to use higher precision computation, or somehow simulate it.
If the original computation was done in IEEE single precision,
then using double for the calculation of the forward error would
usually be fine. If the original was done in double precision,
you'd either need IEEE extended precision (which is hard to get
hold of reliably, since most computing systems don't give direct
access to it), or some sort of quad precision hardware or emulation.

1. Thank you.

Your say that you can check the algorithm's result on your test matrix A by using numerics in higher precision. But the claim is that you have developed an algorithm with better precision, i.e. on a _generic_ input matrix.

Is there an analytical way of testing the algorithm precision, (even to a precision closed to machine precision) using the machine arithmetic model, such as the method described in Higham's Accuracy and Stability book (chapter 2)?

Do people in your field carry out this sort of analysis? I ask because I have no experience in numerical analysis.

2. Hi Michele,
For things like matrix multiplication, the model described in the Accuracy and Stability book enables us to bound the forward error analytically, as you've described, by looking at how rounding errors propagate through the computation.

In some cases, the same ideas can be applied to matrix functions. For example the 'Schur' algorithm for computing the square root of a matrix can have its error bounded in this way. Unfortunately, for more complicated algorithms, such as those used to compute matrix exponentials and logarithms, full rounding error analyses giving useful error bounds are not yet available. This is why numerical testing is so important.

So to answer your final question, people in the field certainly do this sort of analysis, but it's tricky for really sophisticated algorithms.

3. Thank you, Edvin.