Skip to content

How the solver works

pysnape's solver, Nested, answers two nested questions:

  • Inner: if the coefficients had these values, what smooth signal (a B-spline) best fits the measurements while respecting the equation? The equation enters as a weighted penalty, and this is solved to convergence.
  • Outer: which coefficients give the signal closest to the measurements? They are judged by the data alone; a Levenberg–Marquardt search moves them, using exact derivatives of the signal with respect to the coefficients.

Statisticians know this structure as profiled estimation (Ramsay, Hooker, Campbell and Cao, 2007). fit sets it up for you; this page explains the choices it makes and how to change them.

The basis

The fitted signal is a B-spline, set by a number of knots and an order per coordinate.

  • Order: the highest derivative the equation takes in that coordinate, plus 3. One order above the minimum improves accuracy markedly (about 30× on the KdV and Kuramoto–Sivashinsky tests).
  • Knots: by default about one per three grid points (knots="auto").

With the equation enforced, the basis must be fine enough to follow the solution's fastest feature: under-resolved bases bias the coefficients even without noise. Set it yourself when the automatic choice is wrong:

import pysnape as ps
from pysnape.data.generators import ode

data = ode.vanderpol(t_end=20.0, n=1001, mu=2.0)
eq = "x_tt + a*x_t + b*x**2*x_t + c*x = 0"
start = {"a": 5.0, "b": 5.0, "c": 5.0}

one = ps.fit(eq, data, unknowns=start, knots=200, order=6, multistart=False)
best = ps.fit(eq, data, unknowns=start, knots=[100, 200, 300], multistart=False)
print(best.summary())  # lists the bases tried and keeps the one that fits the data best

knots also takes one count per coordinate, {"x": 100, "t": 60}. A chaotic Kuramoto–Sivashinsky field, for instance, needed about 0.2 knots per point in x and 0.4 in t, which no single fraction gives. When the fitted signal misses the data by more than 25 %, fit warns that the basis is probably too coarse.

The equation weight

How strongly to enforce the equation is not known in advance. Nested raises the weight in stages, 0.1 to 10⁷ in units where 1 weighs data and equation equally, each stage starting from the last, and keeps the stage whose neighbours agree best: the estimate that is stable in the weight.

It stops early once a stable answer exists and the next stages stop being informative: when the data misfit has doubled two stages later (the basis can no longer follow the data under the equation) or a stage makes no progress. Before an answer has settled nothing stops it, since early stages may move the coefficients far. summary() says when and why it stopped. weights= sets the stages yourself.

Multi-start

By default each stage also tries a grid of starting points scaled to the data (0.01, 0.1 and 1 times the size each coefficient needs to balance the equation), so a poor start does not trap the search. This is what makes very noisy problems work, and it costs time: 3 to 30 times a single search. With many coefficients the grid is sampled, 256 points at most.

Turn it off with multistart=False when the data is clean or the start is good.

Speed

  • verbose=True prints each stage as it finishes, with the estimates and the elapsed time.
  • multistart=False and fewer weights are the largest savings.
  • Models linear in the fields (no term multiplies two fields: waves, mass–spring chains, the heat equation) are set up once and run faster than nonlinear ones.
  • On a large grid, fit a window of it, or thin it, to choose the basis, then fit the whole.
  • With several fits in parallel (jobs= in ensemble and bootstrap), set OMP_NUM_THREADS=1 so the sparse solvers do not compete for cores.
  • pip install "pysnape[sparse]" adds CHOLMOD, about 14× faster than the default sparse solver on large 2-D problems.

Measured on one laptop: Burgers on its 256 × 101 PDE-FIND grid, two minutes with the defaults; the 2-D wave on the paper's 50 × 50 × 100 grid, 26 minutes with multi-start off; a five-mass chain from noisy accelerations, 87 minutes with the defaults, of which most is multi-start.

Limits

  • Real-valued fields only, for now.
  • The basis must be able to represent the solution: see above.
  • The coefficients must be constants and enter linearly (Writing equations).