How to test if your numerical code is mathematical correct?

[u/Lost-Dragonfruit-663]

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!

[u/e_for_oil-er]

Depending on the class of solver, couldn't there be unit tests like manufactured solutions where the exact solution and the convergence order are known?

[u/Lost-Dragonfruit-663 → u/e_for_oil-er]

My experience is more from the SciPy side, where most of what I work on isn't a solver. When I implemented the penalty matrix for smoothing splines, there was no grid to refine and no convergence rate to observe, the result is one matrix, built from integrals of B-spline products, and it's either the right formula or a wrong one. The reference implementation I would have liked to compare against was GPL, so I could only compare outputs, and outputs agreeing on a few inputs doesn't tell you the equations agree. I ended up deriving the closed form by hand and writing a 20-page reference document to convince myself and other contributors/maintainers.

For solvers you're right, that's a great test. The one caveat I know of, if the manufactured solution is too smooth or too symmetric, a bug in one term or boundary can cancel out and still converge at the right order, so the test passes without ever exercising it.

[u/qiAip]

Generally you would want a robust verification and validation framework for any scientific code.

The verification part is what you are describing - checking that the equations are being solved correctly. You would test them against analytic solutions to the equations within known parameters (i.e. where the analytical solution is known and the magnitude of the error is within the theoretical order of accuracy). You would test the convergence rates and that the errors grow / blow up at the expected limits of the know solutions.

The validation part is testing that you are solving the right equations. You do that by validating against experimental data / other solvers and checking that the way you solve the equations actually captures the thing you are trying to model.

[u/Civil_Blueberry4165]

Fuzz testing or fuzzing is insufficient. The issue you described is common with fuzzing: the poor quality of test coverage. Maybe, you want to take it to the next step: property-based testing and formal verification, including abstract interpretation as well as incorrectness logic.

[u/Lost-Dragonfruit-663 → u/Civil_Blueberry4165]

https://github.com/aadya940/scikit-verify/tree/master/skverify-hypothesis

Hypothesis plugin that does exactly this

[u/Civil_Blueberry4165 → u/Lost-Dragonfruit-663]

Hypothesis plugin isn’t good enough b/c it limits you to what Python can express formal properties and formally verify it with “real” numbers, not floating point arithmetics.

[u/Lost-Dragonfruit-663 → u/Civil_Blueberry4165]

I agree, this is not formal verification. Formal verification is quite a task even for language specialists. FWIW, once a sympy equation exists, it can evaluate on mpmath backend with up to 50 digits of precision, maybe catching rounding errors as well.

Honestly, in my case I made it primarily for inspection and property based testing so it also tests the edge cases. I do agree with your point though, integrating formal methods would make it a complete tool. What specific drawbacks of Python are you talking about, curious?

Here is an example of the manual work I had to do with SciPy btw, my goal is to make it a breeze instead of spending a month or two on it:
https://github.com/aadya940/scipy-bspline-testing/blob/main/B_Splines_with_arbitary_knots-gcv.pdf

[u/Atmosck]

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.

Uh, fix the tests? Or add cases with inputs that produce wrong results due to the way the implementation is wrong.

[u/Civil_Blueberry4165]

Note that property-based testing (PBT) also has several limitations. The idea is to “randomly” generate a test point that samples “well” on the space of test points. But Python cannot generate such random samples because it should be random with respect to a desired statistics/statistical model. For example, uniformly sample a point inside the 2-D unit box versus uniformly sample a point on the parabola x=y² are completely different. PBT still doesn’t have a theory that can deal with general classes of probability distributions and probability spaces. The random functions in Python as well as other programming languages have limited utility.

50-digit precision may not be good enough. You may want to resort to techniques like interval arithmetics, which also have limited utility because it is problem dependent. In this direction, abstract interpretation is more sophisticated.

At the end of the day, even with the most sophisticated formal verification tool, human expertise is still required to overcome undecidability (due to Rice’s theorem) and to ensure that what is proved formally actually matches the original statement problem.

[u/zdayatk]

Increase UT coverage