Efficient replacement for ppval

Question

I have a loop in which I use ppval to evaluate a set of values from a piecewise polynomial spline. The interpolation is easily the most time consuming part of the loop and I am looking for a way improve the function's efficiency.

More specifically, I'm using a finite difference scheme to calculate transient temperature distributions in friction welds. To do this I need to recalculate the material properties (as a function of temperature and position) at each time step. The rate limiting factor is the interpolation of these values. I could use an alternate finite difference scheme (less restrictive in the time domain) but would rather stick with what I have if at all possible.

I've included a MWE below:

x=0:.1:10;
y=sin(x);
pp=spline(x,y);
tic
for n=1:10000
    x_int=10*rand(1000,1);
    y_int=ppval(pp,x_int);
end
toc

plot(x,y,x_int,y_int,'*') % plot for sanity of data

Elapsed time is 1.265442 seconds.

Edit - I should probably mention that I would be more than happy with a simple linear interpolation between values but the interp1 function is slower than ppval

x=0:.1:10;
y=sin(x);
tic
for n=1:10000
    x_int=10*rand(1000,1);
    y_int=interp1(x,y,x_int,'linear');
end
toc

plot(x,y,x_int,y_int,'*') % plot for sanity of data

Elapsed time is 1.957256 seconds.

Answer 1

This is slow, because you're running into the single most annoying limitation of JIT. It's the cause of many many many oh so many questions in the MATLAB tag here on SO:

MATLAB's JIT accelerator cannot accelerate loops that call non-builtin functions.

Both ppval and interp1 are not built in (check with type ppval or edit interp1). Their implementation is not particularly slow, they just aren't fast when placed in a loop.

Now I have the impression it's getting better in more recent versions of MATLAB, but there are still quite massive differences between "inlined" and "non-inlined" loops. Why their JIT doesn't automate this task by simply recursing into non-builtins, I really have no idea.

Anyway, to fix this, you should copy-paste the essence of what happens in ppval into the loop body:

% Example data
x = 0:.1:10;
y = sin(x);
pp = spline(x,y);


% Your original version
tic
for n = 1:10000
    x_int = 10*rand(1000,1);
    y_int = ppval(pp, x_int);
end
toc


% "inlined" version

tic

br = pp.breaks.';
cf = pp.coefs;

for n = 1:10000

    x_int = 10*rand(1000,1);

    [~, inds] = histc(x_int, [-inf; br(2:end-1); +inf]); 

    x_shf = x_int - br(inds);    
    zero  = ones(size(x_shf));
    one   = x_shf;
    two   = one .* x_shf;
    three = two .* x_shf;

    y_int = sum( [three two one zero] .* cf(inds,:), 2);
end
toc

Profiler:

Results on my crappy machine:

Elapsed time is 2.764317 seconds.  % ppval
Elapsed time is 1.695324 seconds.  % "inlined" version

The difference is actually less than what I expected, but I think that's mostly due to the sum() -- for this ppval case, I usually only need to evaluate a single site per iteration, which you can do without histc (but with simple vectorized code) and matrix/vector multiplication x*y (BLAS) instead of sum(x.*y) (fast, but not BLAS-fast).

Oh well, a ~60% reduction is not bad :)

Answer 2

It is a bit surprising that interp1 is slower than ppval, but having a quick look at its source code, it seems that it has to check for many special cases and has to loop over all the points since it it cannot be sure if the step-size is constant.

I didn't check the timing, but I guess you can speed up the linear interpolation by a lot if you can guarantee that steps in x of your table are constant, and that the values to be interpolated are stricktly within the given range, so that you do not have to do any checking. In that case, linear interpolation can be converted to a simple lookup problem like so:

%data to be interpolated, on grid with constant step
x = 0:0.5:10;
y = sin(x);

x_int = 0:0.1:9.9;

%make sure it is interpolation, not extrapolation
assert(all(x(1) <= x_int & x_int < x(end)));

% compute mapping, this can be precomputed for constant grid
slope = (length(x) - 1) / (x(end) - x(1));
offset = 1 - slope*x(1); 

%map x_int to interval 1..lenght(i)
xmapped = offset + slope * x_int;
ind = floor(xmapped);
frac = xmapped - ind;
%interpolate by taking weighted sum of neighbouring points
y_int = y(ind) .* (1 - frac) + y(ind+1) .* frac;

% make plot to check correctness
plot(x, y, 'o-', x_int, y_int, '.')