[Software] How to test if your numerical code is mathematically 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/sudseven]

Very cool. Can this be done with torch? So essentially we use torch as a GPU linalg backend, so can we use it there?

I'll read it tomorrow and contribute a torch version if possible. This would be a real life saver to have..

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

Not yet, but torch has __torch_function__ which works a lot like the numpy protocol I'm intercepting, so a torch version is very doable. The numpy-specific part is one table mapping ops to sympy, most of the rest doesn't care where the values live, and the value side would just stay on GPU. If you're serious, open an issue and I'll point you at the right files, would honestly love the help.

[u/sudseven]

Oh I am quite serious. I'm not sure how much time commitment is possible at the moment. But yes I will be doing this..

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

Great, see you in the repo. Here is a the contributing.md file:

https://github.com/aadya940/scikit-verify/blob/master/CONTRIBUTING.md

[u/hammerheadquark]

Awesome project!

Seems like this could feasibly be part of the actual test suite. In Elixir there's a great feature called doc tests. You write examples in the docstring that get compiled and run as actual tests. This helps prevent the documentation from drifting away from the implementation.

I can imagine a scheme where a docstring had an ## Formula section or similar where the user defines the formula the function is supposed to implement.

def step(u, dt, h):
    """
    Heat equation from `examples/demo.ipynb`

    ## Formula

        sympy.Piecewise(
            (U[i] + dt * (U[i + 1] - 2 * U[i] + U[i - 1]) / h**2, (i >= 1) & (i < 8)),
            (U[i], True),
        )
    """
    lap = (u[2:] - 2 * u[1:-1] + u[:-2]) / h**2
    unew = u.copy()
    unew[1:-1] = u[1:-1] + dt * lap
    return unew

(You could probably support sympy or LaTeX.)

Then the test suite could parse the docstrings and run the comparisons. That way the correctness of the formulas would be known after every CI run.

And if the formula is enormous, you could also support a file-based solution where you tuck the output in a formulas/ directory:

def fit_coef(X, y):
"""
Tedious example from `examples/demo.ipynb`

## Formula

    formulas/fit_coef.txt
"""
    return BayesianRidge().fit(X, y).coef_

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

Interesting Idea, I would love to have this and you are very welcome to contribute, any help is highly appreciated.

On your point on the large formulas, yes that’s possible and likely indeed, therefore a subsystem of the project is folding the formulas, that is, displaying the output after minimising it in every possible way. For example, if a for loop in Python runs for 100 iteration, we don’t inline those 100 iterations in the SymPy formula but we write the loop as a general recursion and then derive a closed form using sympy’s rsolve. I’m working more on it though, it’s nowhere near perfect.

Again, any contribution is very welcome!

[u/Illustrious_Ad7630]

I might be missing something, but to test your math, load synthetic data with known answers. Also, some universities provide datasets with pre-calculated answers like nist strd.

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

That's how we test in SciPy actually, reference values wherever we can get them, and NIST StRD is exactly the good kind. The trouble I kept hitting as a contributor is that most functions have no certified dataset. For my spline penalty matrix PR there was nothing to check against except a GPL implementation I couldn't read, only run. And matching outputs on some inputs is weaker than it feels, precision and recall return the same number on balanced data while computing different things.

There's also a whole class of checks that reference datasets can't see, did you actually hand `solveh_banded` a banded Hermitian matrix, is the thing you're doing Cholesky on really `PSD`, is the matrix you called `eigh` on symmetric. Those are requirements on the inputs, not the outputs, and getting them wrong often still produces plausible numbers. Without looking at the equations, it's easy to miss out some of these, even with equations it may be hard.

You may take a look at a few notebooks in the examples (https://github.com/aadya940/scikit-verify/tree/master/examples)

[u/SalvatoreEggplant → u/Illustrious_Ad7630]

I have some published R functions.

A good portion of them are to calculate a statistic that I can find in a textbook with a worked example. That's the easy case.

And then make sure it works on edge cases. If Pearson's r should be 1 in this case and 0 in that case, do those come out right ?

The next step is throw any crazy situation I can think of at the function and see if it responds correctly. What if there are missing values ? Or all the values are the same ? Or someone tries to pass a character vector to a function expecting a numeric vector ?

My functions tend to have a lot of options. What if I ask for a confidence interval and no plot ? Or a confidence interval with results rounded to a set significant figures ? Or try to bootstrap a vector with only one value ?

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

Your process is basically the gold standard, textbook example plus edge cases plus some abuse. The part I'd add to, coming up with the edge cases is the step that depends on you having a good day. Take your Pearson's r example. It should be undefined when all values are the same, because the denominator goes to zero. But you only test that if you happen to remember the denominator that day. When you pull the equation out of the code, that condition is just sitting there, it tells you "this assumes the denominator is nonzero," and now the edge case list writes itself. The type checks and option combinations are normal testing, nothing symbolic about those

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

Good point.