Linear interpolation with interp1 is slow: How to improve run times?

I wrote a MWE to compare the perforance of interp1 (MATLAB built-in) with my custom interp1_scal (for linear interpolation), using a simple dynamic programming problem solved with value function iteration

$$ V(k) = max_{k'} \log(k^\alpha-k')+\beta*V(k') $$

Question 1: I am a bit surprised that the custom interpolation routine I wrote is much faster than the built-in interp1 for linear interpolation (on my laptop with R2024b, it is 4 times faster, here (RUN button) the speed advantage is smaller but still significant). To find the location of the query point xi on the x grid, I am using binary search (see function locate), directly translated from an older Fortran/C routine.

Question 2: Am I using interp1 in the wrong way? would using griddedInterp or other Matlab routines significantly improve run times of this MWE?

Any help, suggestions or comments is greatly appreciated!

P.S. I am glad to provide additional info on the dynamic programming problem that this simple MWE is meant to caputure in few lines of code, if someone asks for.

Minimum Working Example here:

clear,clc,close all

%% Parameters
maxit = 50;
tol   = 1e-5;
alpha = 0.33;
beta  = 0.96;
k_ss   = (alpha * beta)^(1 / (1 - alpha));
n_k    = 300;
k_min  = 0.5 * k_ss;
k_max  = 1.5 * k_ss;
k_grid = linspace(k_min, k_max, n_k)';
% Precompute output y = k^alpha
y_grid = k_grid .^ alpha;

%% Method 1: value function iteration with Matlab built-in interpolation
tic
% Initialization
V0 = zeros(n_k, 1);
V1 = zeros(n_k, 1);

% Value Function Iteration
for iter=1:maxit
    for i = 1:n_k
        y_val     = y_grid(i);
        k_max_now = min(y_val, k_max); % upper bound for Brent minimization
        % Objective function (negative value for minimization)
        obj = @(kp) -(log(y_val-kp) + beta*interp1(k_grid,V0,kp,'linear'));
        % Optimal policy
        kp_opt = fminbnd(obj, k_min, k_max_now);
        % Updated value function
        V1(i) = -obj(kp_opt);
    end
    % Convergence check
    err = max(abs(V1 - V0));
    V0  = V1;
end %end iter
time_vfi1=toc;

%% Method 2: value function iteration with *custom* interpolation
tic
% Initialization
V0 = zeros(n_k, 1);
V2 = zeros(n_k, 1);

% Value Function Iteration
for iter=1:maxit
    for i_k = 1:n_k
        y_val     = y_grid(i_k);
        k_max_now = min(y_val, k_max);
        % Objective function (negative value for minimization)
        obj = @(kp) -( log(y_val-kp) + beta*interp1_scal(k_grid,V0,kp) );
        % Optimal policy
        kp_opt = fminbnd(obj, k_min, k_max_now);
        % Updated value function
        V2(i_k) = -obj(kp_opt);
    end
    % Convergence check
    err = max(abs(V2 - V0));
    V0  = V2;
end %end iter
time_vfi2=toc;

%% Compare vfi1 and vfi2
err_V   = max(abs(V1-V2));
disp('RESULTS OF COMPARISON')

RESULTS OF COMPARISON

fprintf('||V1-V2||         = %f \n',err_V)

||V1-V2||         = 0.000000 

fprintf('Run time of vfi1 with Matlab interp = %f \n',time_vfi1)

Run time of vfi1 with Matlab interp = 1.293248 

fprintf('Run time of vfi2 with custom interp = %f \n',time_vfi2)

Run time of vfi2 with custom interp = 0.751512 

fprintf('time vfi1 / time_vfi2 = %f \n',time_vfi1/time_vfi2)

time vfi1 / time_vfi2 = 1.720861 

%% Plot value function
figure;
plot(k_grid, V1, 'LineWidth', 1.5);
hold on
plot(k_grid, V2, 'LineWidth', 1.5);
legend('V1','V2')
xlabel('Current capital');
grid on;

% Fast linear interpolation routine
% Usage:
% yi = interp1_scal(x,y,xi)
% where x and y are column vectors with n elements, xi is a scalar 
% and yi is a scalar
% Input Arguments
% x - Sample points
%   column vector
% Y - Sample data
%   column vector
% xi - Query point
%   scalar
function yi = interp1_scal(x,y,xi)

n = size(x,1);
j = locate(x,xi);
j = max(min(j,n-1),1);

slope = (y(j+1)-y(j))/(x(j+1)-x(j));
yi = y(j)+(xi-x(j))*slope;

end %end function interp1_scal

function jl = locate(xx,x)
%function jl = locate(xx,x)
%
% x is between xx(jl) and xx(jl+1)
%
% jl = 0 and jl = n means x is out of range
%
% xx is assumed to be monotone increasing

n = length(xx);
if x<xx(1)
    jl = 0;
elseif x>xx(n)
    jl = n;
else
    jl = 1;
    ju = n;
    while (ju-jl>1)
        jm = floor((ju+jl)/2);
        if x>=xx(jm)
            jl = jm;
        else
            ju=jm;
        end
    end
end

end %end function locate

[Stephen23]

Question 1: By avoiding the input checking, error handling, and handling of different input configurations it is often possible to write code for domain-specific code which is faster than MATLAB's generic inbuilt Mfiles. This is not unexpected.

Question 2: Your use of INTERP1 seems okay. I doubt that GRIDDEDINTERPOLANT would be faster, unless you can provide the sample points once and the query points in large arrays. You might be interested in trying INTERP1Q, which essentially does the same as your custom function:

https://www.mathworks.com/help/matlab/ref/interp1q.html

Answer 1

[Paul]

Hi Alessandro,

The base version of interp1 is an ordinary m-file, so you can see all of the extra processing it needs to do to account for all possible inputs, whereas interp1_scal only works for column vector sample points with linear interpolation and no extrapolation and a scalar query point.

To that last point, consider the comparison when one wants to interpolate lots of data points and we have to write a wrapper around interp1_scal to loop over the query points.

x = (0:.01:10).';
y = sin(x);
rng(100);
xq = rand(1e4,1)*x(end);

ntrials = 40;
[texec1,texec2] = deal(nan(ntrials,1));
for ii = 1:ntrials
    tic,yi1 = interp1(x,y,xq);          texec1(ii) = toc;
    tic,yi2 = interp1_scal_vec(x,y,xq); texec2(ii) = toc;
end

figure
plot([texec1,texec2]),grid,legend('interp1','interp1_scal','Interpreter','none');

% results aren't identical, but very close
[max(abs(yi1-yi2)), norm(yi1-yi2)]  

ans = 1×2

1.0e-14 * 0.0222 0.4965

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I think timeit does not include the first call

timeit(@() interp1(x,y,xq))

ans = 4.5437e-05

timeit(@() interp1_scal_vec(x,y,xq))

ans = 0.0021

Calling with a scalar input, as in the question, does show an advantage for the specialized code, which might become more pronounced if converted to mex.

timeit(@() interp1(x,y,xq(1)))

ans = 7.5319e-06

timeit(@() interp1_scal(x,y,xq(1)))

ans = 2.2946e-06

function yi = interp1_scal_vec(x,y,xq)

yi = nan(size(xq));
for ii = 1:numel(xq)
    yi(ii) = interp1_scal(x,y,xq(ii));
end

end

function yi = interp1_scal(x,y,xi)
n = size(x,1);
j = locate(x,xi);
j = max(min(j,n-1),1);
slope = (y(j+1)-y(j))/(x(j+1)-x(j));
yi = y(j)+(xi-x(j))*slope;
end %end function interp1_scal

function jl = locate(xx,x)
%function jl = locate(xx,x)
%
% x is between xx(jl) and xx(jl+1)
%
% jl = 0 and jl = n means x is out of range
%
% xx is assumed to be monotone increasing
n = length(xx);
if x<xx(1)
    jl = 0;
elseif x>xx(n)
    jl = n;
else
    jl = 1;
    ju = n;
    while (ju-jl>1)
        jm = floor((ju+jl)/2);
        if x>=xx(jm)
            jl = jm;
        else
            ju=jm;
        end
    end
end
end %end function locate

[John D'Errico]

Good answer by Paul. I would add that you can improve interp1_scal though, since if the independent variable is known to lie on a grid, then you can eliminate the call to locate, which employs a bisection scheme. And yes, bisection is reasonably efficient, but IF the grid is KNOWN to be uniform, then there is a simple formula to decide in which interval a query point lies. All you need is the first point in the sequence, and the stride between points. And that will then be far faster than bisection, since it will require only one basic line of code. (Which is itself a form of linear interpolation.)

Another simple trick would be to use a lookup table to determine which bin any point lies in, or at least come very close. If the stride is not uniform, this will solve the problem nicely, though it will force you to precompute the lookup table.

Again, IF you are able to simplify the problem, then you can greatly decrease the time needed. But that should not be at all surprising. A general code like interp1 cannot make any assumptions, so it needs to test things. It needs to validate the data as it comes in, etc.

[Paul]

Looks like interp1 is faster with equally spaced sample points, but I'm a bit surprised it's not more faster. AFAICT, the doc doesn't say if interp1 checks for equally spaced sample points.

rng(100);
x = [0;sort(rand(1001,1));1]*10;
y = sin(x);
xq = rand(1e6,1)*x(end);

timeit(@() interp1(x,y,xq))

ans = 0.0065

x = linspace(0,10,numel(x));
timeit(@() interp1(x,y,xq))

ans = 0.0016

[Alessandro]

@John D'Errico thanks for your tips on improving interp1_scal. In my applications, the x grid is always a strictly increasing column vector, but not equally spaced. If they were equally spaced (as they are in the MWE), I could find the index jl such that x_grid(jl)<xi<x_grid(jl+1) as follows:

%jl = floor((xi-x_grid(1))/step)+1;

where step is the stride between points. However, in most cases I want to concentrate the x_grid points near the lower bound, so I would generate the grid as

x_min = 0.001; x_max = 5; n_x = 50;
x_grid = x_min + (x_max-x_min)*(linspace(0,1,n_x).^3)';
figure,plot(x_grid,x_grid,'-o'),title('x grid')

●

Is there something better than binary monotonicity for this case? Could you please expand on the idea of a lookup table? Thanks!

[Alessandro]

I have compared also interp1 with interp1q (based on @Stephen23 suggestion) using @Paul data

clear,clc,close all
x = (0:.01:10).';
y = sin(x);
rng(100);
xq = rand(1e4,1)*x(end);
timeit(@() interp1(x,y,xq))

ans = 6.6625e-05

timeit(@() interp1q(x,y,xq))

ans = 6.5462e-04

But on my laptop interp1q is not faster. Here it seems that interp1 is the winner. Maybe matlab does not maintain anymore interp1q? I read on the doc that interp1q is NOT recommended.

[John D'Errico]

@Alessandro

A lookup table is easy enough to do.

x = cumsum(randi(3,[1,10000]))

x = 1×10000

1 2 3 6 7 8 10 12 13 15 18 19 20 21 23 26 29 32 35 37 39 40 43 45 48 50 51 53 56 57

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[min(x),max(x)]

ans = 1×2

1 20109

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, suppose you will come into this with some random number in that interval? You could use bisection, but there are MANY intervals to search through. A bisection code can take roughly

log2(10000)

ans = 13.2877

tests. While that is finite, you might be better off if you precompute many possible results in advance, using a vectorized tool like discretize.

xb = x(1):1:x(end);
tic,xbindex = discretize(xb,x);toc

Elapsed time is 0.007563 seconds.

Yes, this will take some time, but as you can see, not a lot. But now for any value of xb, you already know which interval it lives in. That takes no more than an index into the vector xbindex.

If xi is a continuous variable, it may then require an extra step to refine the result, but you will be quite close.