Faster bspline evaluation using vectors?

[schirmermischa]

Dear Bjarne,

I have looked into splinter because I need to calculate b-splines / p-splines for large data sets (astronomical images, usually 8-16 million pixels each). While the splines describe the data exactly as i had hoped for, I ran into a bottleneck evaluating them. Basically, I have sth like

for (k=0; k<1e7; ++k) { pixel[k] = bspline.eval(<pixel coordinate>); }

On a 3GHz CPU this takes about 15 seconds to evaluate. Usually, there are hundreds of images involved, and the evaluation of the spline fit by far dominates the execution time. I tried changing the source code by feeding a vector of coordinates to bspline.eval() instead of a single data point. However, I got lost in the Eigen library. It appears that eval() transforms the data to some Eigen DenseVector and then simply back, with some magic happening in between that I didn't want to unearth by digging through the Eigen library.

Do you have other ideas how the evaluation could be sped up? It would be a great help, because good C++ p-spline libraries are hard to come by.

Thank you,

Mischa Schirmer

[bgrimstad]

Hi Mischa,

Happy that you find the library useful.

I can think of two obvious ways to make your evaluations faster:

1. Evaluate in parallel (easy)
2. Implement a smart batch evaluation (hard)

The first option should be relatively easy to implement since evaluation is a read-only operation.

A naive implementation of batch evaluation - like the one I just pushed to the "batch_eval" branch - does not improve upon the iterative approach. We are already using sparse vectors and matrices and there is not much to gain from this. In fact, this batch evaluation is a tiny bit slower. It is however possible to implement a batch evaluation procedure that will be much faster for your use case. The idea is to exploit the fact that you evaluate the (monovariable) basis functions at the same points many times over, which is unnecessary. Say you evaluate at points (x1, x2) = (1, 1), (1, 2), ..., (1, n), then the basis functions for x1 are evaluated n times at 1 (which of course yields the same vector each time). To implement a batch evaluation procedure that exploits this will demand some knowledge about SPLINTER, B-splines and the Eigen library so it is a challenging task.

Hope this helps.

Have a nice weekend,
Bjarne

[schirmermischa]

Hi Bjarne, the batch evaluation indeed make sense, but I won't have the time to implement it now. My code is already parallelized (working on several images simultaneously).

[kfjahnke]

Hi collegues!
I just happened to browse your issues and saw 'vectors' mentioned. I'm the author of vspline, a FOSS C++ template library for b-spline processing, and since my library uses SIMD, I thought this might be what schirmermischa is looking for - it can in fact evaluate SIMD arguments. My library is specifically useful for processing large raster data sets, providing automatic multithreading and vectorization. I use it a lot for image processing.
Kay