Cubic Spline optimization

[kazewong]

To start with, nice package! Somehow the Jax devs don't get how often the scientific computing community uses interpolation, and this seems to fill the gap quite well.

Context:
I needed a cubic spline a while ago, so I was searching it online and there were some implementations out there. But the problem is the performance is not quite what I need, and the bigger issue is they usually have a pretty large memory footprint because at some point in their code, they usually have jnp.diag, which creates a diagonal matrix that allocates O(n^2) memory. This is problematic when my basis is 10^6 points. Then I coded a lightweight version with equinox and lineax, Here is the code that build the spline representation using lineax (https://github.com/kazewong/JaxNRSur/blob/aa24a2d5d4221f420c24d7b02391ee1c762ca9cc/jaxNRSur/Spline.py#L59). This instead only allocate O(3n) memory since it knows the system we are solving is tridiagonal.

I am looking through your implementation of the cubic spline, which seems to do similar construct around

interpax/interpax/_spline.py

Line 1048 in c5c9d4e

def _approx_df(x, f, method, axis, **kwargs):

.

Do you have a benchmark of the performance of your implementation? If you think this is relevant, I can help open a PR to put it in this code. It makes more sense that I ship my version of code out.

[f0uriest]

I haven't done any real benchmarking, most of my applications are at most a few hundred points so it hasn't been a big deal. I'd certainly welcome contributions if you want to open a PR (I've been looking for an excuse to dive into lineax).

One other point, if you're ok with C1 continuity you can use method="cubic" (or cardinal, catmull-rom, monotonic). It's only for C2 continuity (what I called cubic2) that you need to solve a linear system.