f0uriest
commented
1 year ago 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.
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 withequinoxandlineax, 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 https://github.com/f0uriest/interpax/blob/c5c9d4e0203afb231f709e3537f4f935d211469c/interpax/_spline.py#L1048.
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.