Speed up systems with interpolant

Dear all,

I have to solve a system with interpolants and I have to find the fastest way because it requires quite a lot of CPU time. I noticed that using a 2D interpolant if faster than 1D and this sounds very strange. Here an example of the code:

x1=[0:9];
y1=[0 1 2 4 6 9 13 20 40 100];

tic
for k = 5:1:15

f = @(x)(interp1(x1,y1,x));
%mat = scatteredInterpolant([x1 0 9]',[1 1 1 1 1 1 1 1 1 1 2 2]',[y1 0 100]');
%f = @(x) mat(x,1);

      F1 = @(z1,z2,z3) f(z1)+f(z2)+k*(f(z3)*f(z3))/(f(z3)+f(z3))-15;
      F2 = @(z2,z3) f(z2)+k*(f(z3)*f(z3))/(f(z3)+f(z3))-10;
      F3 = @(z2,z3) f(z2)-k*(f(z3)*f(z3))/(f(z3)+f(z3))+9;

      fun = @(z) [F1(z(1),z(2),z(3)), F2(z(2),z(3)), F3(z(2),z(3))];

      [z] = fsolve(fun, [5 5 5]);

end
toc

If now I commentate line

f = @(x)(interp1(x1,y1,x));

and uncomment the lines

mat = scatteredInterpolant([x1 0 9]',[1 1 1 1 1 1 1 1 1 1 2 2]',[y1 0 100]');
f = @(x) mat(x,1);

the system solves almost twice faster, even if the solution is the same. Even if it is faster I would not like to use this trick, because I don't see any sense to transform a 1D interpolation to a 2D just to speed up the solution.

Do you have an idea why?

Thanks and regards.

Answer 1

[Yair Altman]

scatteredInterpolant is a built-in function so we do not know how it is actually implemented internally, but it stands to reason that an internal compiled C function would run faster than an interpreted m-function (interp1.m). Also, it is quite possible that scatteredInterpolant is optimized internally for the 1D case, so mat(x,1) might well be using an optimized code branch.

If you wish to improve interpolation speed, consider using interp1q. Unlike the generic interp1 function, interp1q avoids any checks of its inputs, which are expected to be a monotonically increasing column vector (x), a column vector or matrix with length(x) rows (Y), and a column vector (xi).

An even faster implementation of the core logic was implemented by Bruno Luong using MEX, which he called nakeinterp1. Jan Simon has also implemented a MEX variant, called ScaleTime, which is many times faster than interp1 or interp1q.

Many additional performance tips can be found in my book: Accelerating MATLAB Performance.