testing derivative approximations

Question

I'm writing a library that involves some approximations of variational calculus problems. Whenever I implement routines to evaluate the derivative or Hessian of an action functional $A$, I write a test to check that these are working correctly. To do this, I evaluate the error in using the first and second derivatives of the action to compute a local quadratic approximation to the objective functional:

$E = A(u + \delta v) - \left(A(u) + \delta\cdot dA(u)\cdot v + \frac{1}{2}\delta^2(d^2A(u)v)\cdot v\right)$

and check that $E/\delta^3$ doesn't grow for $\delta = 2^{-1},\ldots,2^{-N}$ for some big enough $N$.

This works well for moderately sized $N$ (up to 16 or so), but for $N$ too large, truncation error begins to dominate and the error in the approximation gets worse.

What is a good way to test derivative approximations in the face of truncation error?

The way I see it, I have only bad options:

1. Write an initial version of the test that prints the approximation errors, find the $N$ for which truncation error takes over, then use this $N$ as the cutoff point in the final version of the unit test.
2. Have some way of automatically diagnosing where truncation error takes over and discard all the results thereafter.
3. Don't make my unit tests return failure unless it's really obvious. Graph the results and let the user decide.
4. Determine analytically where truncation error is likely to occur.

I don't like option 1 because it feels like cherry-picking my tests so that they will pass. I don't like option 2 because a genuinely failing implementation's inherent mathematical errors might be erroneously assessed as truncation error, so that automated regression testing would miss the bug. At that rate, I'd have to check the results manually, in which case I might as well go with option 3. But I don't like option 3 because I want automated testing. I think option 4 is practical for functions of a single variable -- evaluating $\sqrt{\epsilon|f(x)/f''(x)|}$ where $\epsilon$ is the machine epsilon gives you a pretty good upper bound -- but for big PDE problems where the unknown is a field with thousands of variables, not so much.

Answer 1

Well, I can tell you what I do, but it's somewhat unsatisfactory.

1. Check $dA(u)$ using $A(u)$. I prefer a fourth-order finite difference scheme, but I check $dA(u)\cdot \delta u$ against the approximation

$$ d A(u)\cdot \delta u \approx \frac{A(u - 2\epsilon \delta u)-8A(u - \epsilon \delta u)+8A(u +\epsilon \delta u)-A(u + 2\epsilon \delta u)}{12} $$

2. Once we know that $dA(u)$ is accurate, check $d^2A(u)\delta u$ versus the approximation $$ d^2 A(u)\delta u \approx \frac{dA(u - 2\epsilon \delta u)-8dA(u - \epsilon \delta u)+8dA(u +\epsilon \delta u)-dA(u + 2\epsilon \delta u)}{12} $$ This is not that different from above, except we now have a vector instead of a scalar, so the error is calculated with a norm and not an absolute value.

As you mentioned above, it's hard to know when there's a switch between truncation and round-off error. As such, I run this for several values of $\epsilon$ and look for two things

1. The minimum error in the approximation over all attempted $\epsilon$ reaches an acceptable level

2. Typically, there's a curve where the error goes high-low-high as the value of $\epsilon$ is reduced

In a unit test, I typically only use 1. When validating a function by hand, I look for 2. Technically, 1 could give an erroneous pass in case we get lucky and everything cancels. I've never seen that happen, so I don't worry about it.