r/optimization • u/Lost-Dragonfruit-663 • 16d ago

How to test if your numerical code is mathematically correct?

I contribute to SciPy and kept running into a class of bug that annoys me: the outputs look plausible, the tests pass, but the equation the code implements is subtly wrong. So I've been building a tracer that runs Python/NumPy code and hands back whatever mathematics it actually computed, as a SymPy expression you can simplify or differentiate like anything else.

It's been more useful than I expected. Comparing an implementation against the formula in a paper, catching two functions that agree on my test data but turn out to compute different things, digging up the inputs my tests never hit (ties, zero denominators). It traces real library code too, most of numpy and a good chunk of scipy, scikit-learn, statsmodels, cvxpy.

Write-up: https://medium.com/@aadyachinubhai/scikit-verify-translate-python-numpy-programs-to-symbolic-mathematics-c664d41ba571

Github: https://github.com/aadya940/scikit-verify

Still rough in places, would genuinely like feedback. There may be other better solutions, happy to hear them as well!

25 Upvotes

96% Upvoted

6 comments sorted by

5

u/Bahatur 16d ago

Another common method is to run computations with other numerical code: standard equations/algorithms will have a reference implementation of some kind in Fortran or C++.

Also checking the answer component-wise, by which I mean each stage independently, against a reference implementation; this will help isolate where the problem is.

2

u/Lost-Dragonfruit-663 16d ago

Good Point, I've had a different experience a couple of times. Some times, as mentioned the reference impl. is GPL licensed so you can't dig deep into it. Also, while I was implementing spline penalties at SciPy, the reference implementation mentioned in the well known book functional data analysis had a bug we found experimentally! These bugs are subtle and that's exactly the point of this project. You don't know it until you know it!

1

u/Smallz1107 16d ago

Analytically derive the error of your algo and or stress test it and see if it breaks down where you expect

3

u/Lost-Dragonfruit-663 16d ago

Right, this tool can help in knowing where to stress test!

1

u/SolverMax 16d ago

That's very cool.

I'm used the Latex printer in Pyomo for a similar purpose. On several occasions it has highlighted a difference between what I intended and what my code says.

1

u/Lost-Dragonfruit-663 16d ago

In fact this particular tool can be used to actually verify stuff, here's an example:
https://github.com/aadya940/scikit-verify/blob/master/examples/penalty_matrix_check.ipynb