interp1 - [Matlab R2016a] - Speed Issues

Heeeey guys,

as I already noticed in this community, I am not alone with that problem. I have tried several solutions proposed in here but it didnt help for reasons like: Old Matlab Versions (interp1 has changed significatly) or another kind of problem.

Let me first tell you what i am doing to make it easier for you to help me:

I a for-loop with approx. 30.000 iterations i use interp1 several time (5 times per iteration - see the profiler below).

Example of the call: [interp1( x, y, xq)]

Here, x is a vector (length 10.000) and y is also a vector (length 10.000) and xq is a vector (length 50.000 or more).

The problem is also that y changes in each iteration and therefore, I cannot precompute any operation (xq is also updated each iteration).

I would be very grateful if you had any suggestions to improve the computational time in this case.

Best regards, Ilya

[Walter Roberson]

Would x happen to be equally spaced? Is xq sorted?

Answer 1

[Steven Lord]

Test if creating a griddedInterpolant once before the start of the for loop, replacing the Values property of the object each time you change y, is faster.

[Ilya Klyashtornyy]

Mb i didnt get you right but i have to create a new griddedInterp. each time y is update right? (so basically in esch iteration)

[Steven Lord]

No. Create once, modify the property multiple times, use the object multiple times.

x = 1:10;
y = x;
GI = griddedInterpolant(x, y);
z = zeros(1, 5);
for k = 1:5
  GI.Values = x.^k;
  z(k) = GI(5.5);
end
disp(z)

[Ilya Klyashtornyy]

I replaced all

     int=interp1(x,y,xq);

with

     S=griddedInterpolant(x,y);
     int=S(xq);

It gave me a speed up of almost 25%

Thank you! :)

[Ilya Klyashtornyy]

Hey, thanks for this advice.

I found out that (mb only in my case) creating a new griddedInterpolant takes less time than updating the values.

 S=griddedInterpolant(x,y_new): 1.03s
 S.Values=y_new: 1.32s

Do you agree with me on that? Do you receive different results?

Answer 2

[Walter Roberson]

There are several faster interp1 in the file exchange, including

https://www.mathworks.com/matlabcentral/fileexchange/10286-fast-interpolation

https://www.mathworks.com/matlabcentral/fileexchange/43325-quicker-1d-linear-interpolation--interp1qr

[Ilya Klyashtornyy]

I have already tried the first one and it is dramatically slower (last Update 2006)

Will try the second one now.

Thank you!

[Ilya Klyashtornyy]

for the second one i get the following results:

time for interp1: Elapsed time is 0.003742 seconds. time for interp1qr: Elapsed time is 0.017551 seconds.

I did several iteration to make sure that is not an outlier.

[Walter Roberson]

Are the x equally spaced?

[Ilya Klyashtornyy]

Yes the x vector is equally spaced (it divides the region form -pi to pi in 10.000 sample)

[Walter Roberson]

Under the assumption that xq might be "exactly" -pi but is less than +pi:

xmin = -pi; xmax = +pi; span = xmax - xmin;
idx = floor( (xq - xmin) * length(x) / span + 1 );
offset = xq  - x(idx);
yq = y(idx) * (1-offset) + y(idx+1) * offset;

I would need to double-check that the rounding for the index calculation works out well enough.

The idea is that with x being equally spaced with known min and max, then algebraically you can calculate the index into the x vector through a simple scaling operation. Then you can use the index vector to look up the values to do linear interpolation.

[Ilya Klyashtornyy]

Hey Walter,

thank you very much for your time!

I ve evaluated ther method and compared it to the griddedInterpolant:

 %Time per iteration:
 griddedInterpolant: average - 0.00328s
 your suggestion: average    - 0.00672s

Answer 3

[Nikolaus Koopmann]

% Current date   = June 03, 2020
% Matlab version = 9.6.0.1072779 (R2019a)
% User name      = Nikolaus Koopmann
function [yq,p] = interp1_lin(x,y,xq)
    validateattributes(x,{'double'},{'increasing','vector'}) % slow
    %% lin. regr.
    X = [ones(length(x),1) x(:)];
    p = flipud(X\y(:));                         % see https://www.mathworks.com/help/matlab/data_analysis/linear-regression.html
    %     flin = @(x_)p(1)*x_ + p(2); % slow
    yq = p(1)*xq + p(2); 
end