How to deal with complexity in numerical code, for example, when dealing with large Jacobian matrices?

Question

I am solving a non-linear system of coupled equations, and have calculated the Jacobian of the discretised system. The result is really complicated, below are (only!) the first 3 columns of a $3\times 9$ matrix,

(The complexity arises, in part, because the numerical scheme requires exponential fitting for stability.)

I have quite a general question about implementation of numerical codes using Jacobians.

I can go ahead and implement this matrix in code. But my intuition is telling me to expect a few days (maybe weeks!) of tedious debugging due to the sheer complexity and the unavoidability of introducing errors. How does one cope with complexity such as this in numerical code, it seems unavoidable?! Do you use automatic code-generation from symbolic packages (then tweak the code by hand)?

First I plan to debug the analytical Jacobian with a finite difference approximation, should I be aware of any pitfalls? How do you deal with similar problems in your code?

Update

I am coding this in Python and used sympy to generate the Jacobian. Maybe I can used the code-generation feature?

Answer 1

One word: Modularity.

There are a lot of repeated expressions in your Jacobian that could be written as their own function. There's no reason for you to write the same operation more than once, and that will make debugging easier; if you only write it once there's only one place for an error (in theory).

Modular code will also make testing easier; You can write tests for each component of your Jacobian as opposed to trying to test the whole matrix. For example, if you write your function am() in a modular fashion, you can easily write sanity tests for it, check if you're differentiating it properly, etc.

Another suggestion would be to look at automatic differentiation libraries for assembling the Jacobian. There's no guarantee they're error-free, but there will probably be less debugging/fewer errors than writing your own. Here are a few you may want to look at it:

• Sacado (Sandia Labs)
• ADIC (Argonne)

Sorry, just saw that you're using python. ScientificPython has support for AD.

• ScientificPython

Answer 2

Let me weigh in here with a few words of caution, prefaced with a story. Long ago, I worked with a fellow when I was just starting out. He had an optimization problem to solve, with a rather messy objective. His solution was to generate the analytical derivatives for an optimization.

The problem that I saw was these derivatives were nasty. Generated using Macsyma, converted to fortran code, they were each dozens of continuation statements long. In fact, the Fortran compiler got upset at that, as it exceeded the maximum number of continuation statements. While we found a flag that allowed us to get around that problem, there were other issues.

• In long expressions, as are commonly generated by CA systems, there is a risk of massive subtractive cancellation. Compute lots of big numbers, only to find they all cancel each other out to yield a small number.

• Often analytically generated derivatives are actually more costly to evaluate than are numerically generated derivatives using finite differences. A gradient for n variables may take more than n times the cost of evaluating your objective function. (You may be able to save some time because many of the terms can be re-used across the various derivatives, but that will also force you to do careful hand coding, instead of using computer generated expressions. And any time you code nasty mathematical expressions, the probability of an error is not trivial. Make sure you verify these derivatives for accuracy.)

The point of my story is these CA generated expressions have issues of their own. The funny thing is my colleague was actually proud of the complexity of the problem, that he was clearly solving a really difficult problem because the algebra was so nasty. What I don't think he considered was if that algebra was actually computing the correct thing, was it doing so accurately, and was it doing so efficiently.

Had I been the senior person at the time on this project, I would have read him the riot act. His pride caused him to use a solution that was probably unnecessarily complex, without even checking that a finite difference based gradient was adequate. I'll bet we had spent perhaps a man-week of time to get this optimization running. At the very least, I would have counseled him to carefully test the gradient produced. Was it accurate? How accurate was it, compared to finite difference derivatives? In fact, there are tools around today that will also return an estimate of the error in their derivative prediction. This is certainly true for the adaptive differentiation code, (derivest) I've written in MATLAB.

Test the code. Verify the derivatives.

But before you do ANY of this, consider if other, better optimization schemes are an option. For example, if you are doing exponential fitting, then there is a very good chance you can use a partitioned nonlinear least squares (sometimes called separable least squares. I think that was the term used by Seber and Wild in their book.) The idea is to break the set of parameters into intrinsically linear and intrinsically nonlinear sets. Use an optimization that works on only the nonlinear parameters. Given those parameters are "known", then the intrinsically linear parameters can be estimated using simple linear least squares. This scheme will reduce the parameter space in the optimization. It makes the problem more robust, since you don't need to find starting values for the linear parameters. It reduces the dimensionality of your search space, so making the problem run more quickly. Again I have supplied a tool for this purpose, but only in MATLAB.

If you do use the analytical derivatives, code them to reuse terms. This can be a serious time savings, and may actually reduce the bugs, saving your own time. But then check those numbers!

Answer 3

There are several strategies to consider:

1. Find the derivatives in symbolic form using a CAS, then export code for computing the derivatives.

2. Use an automatic differentiation (AD) tool to produce code that computes the derivatives from code to compute the functions.

3. Use finite difference approximations to approximate the Jacobian.

Automatic differentiation could produce more efficient code for computing the entire Jacobian then using symbolic computation to produce a formula for each entry in the matrix. Finite differences are a good way to double check your derivatives.

Answer 4

Here is an example of where we have used automatic differentiation using Sacado in one code: http://www.dealii.org/developer/doxygen/deal.II/step_33.html

Answer 5

In addition to BrianBorcher's excellent suggestions, another possible approach for real-valued functions is to use the complex-step derivative approximation (see this article (paywalled) and this article). In some cases, this approach yields more accurate numerical derivatives at the cost of changing the values of variables in your function from real to complex. The second article lists some cases where the complex step function approximation might break down.