scipy.interpolate.PPoly becomes numerically unstable for high order

[aarchiba]

The Bulirsch-Stoer integrator can easily produce dense-output polynomials of order higher than 20. This is known to be a problem if you represent them in the power basis, as scipy.interpolate.PPoly does. scipy.interpolate.KroghInterpolator uses an algorithm that is considerablye more stable, numerically. So forcing solvers to output a PPoly "spline" is a numerical problem. Perhaps it would be easier if the solver simply provided a callable for each step, which OdeSolution could then forward calls to?

[drhagen]

The docs clearly need some more work because this is exactly how it works right now. Each solver has a spline method. Perhaps it should be renamed to interpolator.

[aarchiba]

What is the spline object supposed to return? I saw that method but assumed it was required to return a PPolu object for OdeSolution to use.

[drhagen]

It is supposed to return a callable (t) -> y. PPoly is one type of callable. interp1 is another.

[aarchiba]

Then perhaps spline is a misleading name?

More generally, is it helpful to have a function that takes a collection of steps and returns a single callable? If we're using OdeSolution anyway to wrap the callable objects from multiple steps, what is the advantage to using interp1d and PPoly objects (as opposed to one callable per step)? Or just providing an evaluation method in the integrator that takes a state object as information?

[drhagen]

I have already renamed it locally to interpolator.

Both of the alternatives you mention could be done. If there is any kind of setup cost for the interpolator, both of those will be a lot slower. The list of interpolators will require thousands of small setups. The evaluation method would require setting up the interpolator each time points are requested. But I don't know if the setup cost is meaningful, so those remain valid alternatives.

[aarchiba]

If there's per-segment setup cost for the integrator, do that when state is requested. The only thing the interpolator interface does differently is handle interpolators that need information from many state objects. I suspect that no ODE solver will use an object with such global structure, because most are implemented like our solver objects: they can interpolate within the last state, and keep handy all the data to do this. So perhaps simply allow 'state' or 'interpolator' to return an interpolator good for the current step when requested?

[pv]

The advantage of using a generic spline object is in that you have access to all derivatives, integration, etc. with an API that can be "standard" --- which you don't have with a "black-box" representation. This is not where scipy.interpolate is currently, but it is the direction where I want it to go. (C.f. Mathematica's NDSolve returns InterpolatingFunction which I've found very useful in several cases.)

In addition to Krogh (which has a sparse feature set, no integration etc), there's also bernstein basis (don't know about its stability, though). I don't know exactly the representation used in Bulirsch-Stoer, but perhaps it can be converted in a stable way to a representation that either exists in scipy.interpolate --- or, maybe a new representation could be added?

[aarchiba]

The problem is, depending on the ODE solver, you get (potentially) different representations of polynomials, and converting between polynomial representations is numerically troublesome. Collocation methods, for example, might naturally yield polynomials specified in the Lagrange basis, or specified by one value and derivatives but not values at many points (which Krogh can't handle). Plus if we want to use FORTRAN or C solvers from elsewhere - the SUNDIALS CVODE springs to mind - the most natural approach is to use their own dense output routines. So maybe the solution is, like scipy.stats, to specify the interface of a "function" object and maybe provide generic vectorization, numerical integration, differentiation, and root-finding routines that can be overloaded when something better is available? Obviously scipy.interpolate would produce objects that follow the standard.

[aarchiba]

Bulirsch-Stoer naturally produces Krogh-type interpolation - a polynomial specified by value and derivative at both ends, and value and (lots of) derivatives at the middle. So KroghInterpolator works well, constructing the output values from the input data in a fairly clean and stable way. I don't know that there is a clean way to get the Bernstein basis out of this, and with orders that can be in the 20s, sticking to a stable representation is kind of important.

[pv]

Different representations are fine IMO. I would however prefer well-defined and documented representations that can be converted to each other (modulo numerical stability) and inspected directly over structures with opaque internal state, even if provided with a common API. The point would be to get out objects that you can work with potentially in ways for which there is no method in the API.

I don't know if fallbacks to "generic" routines are very good --- numerical root finding, integration, and differentiation has performance and accuracy that is not necessarily so easy to predict, and you e.g. have difficulty in deciding whether all roots were found. Given that the users can apply generic routines here when necessary, it's maybe OK to limit to operations that can be done "exactly" in some sense.

Re Bernstein: http://docs.scipy.org/doc/scipy-0.18.0/reference/generated/scipy.interpolate.BPoly.from_derivatives.html#scipy.interpolate.BPoly.from_derivatives --- whether this is stable at high orders, I don't know.

[aarchiba]

Unfortunately that function only allows derivatives and values to be specified at breakpoints, not in the middle of an interval.

Opaque representations are indeed unpleasant, but I think it's important that we be able to work with well-established and robust ODE solver codes. I'm content using my own Bulirsch-Stoer integrator for my particular problem, but for general use I think it's hard to write a solver that is as bulletproof as something that's been out there being hammered on by users for a couple of decades. And that includes the dense output routines, even if it is feasible to delve into the guts and construct a nominally equivalent polynomial. Though that may be necessary to allow storage of the relevant state.

[pv]

Yes, I agree figuring out the internal representation used by a given ODE solver, and potentially having to reimplement the evaluation routines can be lots of error-prone extra work (even B-splines are not fully trivial), whose usefulness in the end depends on the particular use case. This is not necessarily in conflict in being transparent about the representation and in trying to stick to a common API --- by which I don't necessarily mean a class hierarchy or using a common "OdeInterpolant" high-level interface object, but more about trying to be consistent with the operating philosophy, method names and arguments, behavior vs. vectorization etc.

[aarchiba]

Just so. I think scipy.interpolate is the only place we construct function objects, and unfortunately there isn't a very consistent interface. Part of that is my fault: for KroghInterpolator there's an easy way calculate derivatives up to degree k, so I implemented a derivative function, but for BarycentricInterpolator there isn't a nice way to calculate derivatives. In neither case is there a good (that is, without losing accuracy, maybe catastrophically) way to produce new polynomials representing derivatives or antiderivatives. Those methods are more reasonable for low-order polynomials, where numerical issues aren't as severe. The fact that some methods are missing for some objects doesn't prevent standardizing what they should look like when they exist.

There's also the issue of single-polynomial objects versus piecewise-polynomial objects; for ODEs I think it makes sense to allow infinite domains by computing the interpolating polynomials lazily, but the container object should still support the same interface.