Linear interpolation with interp1 is slow: How to improve run times?
Show older comments
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')
fprintf('||V1-V2|| = %f \n',err_V)
fprintf('Run time of vfi1 with Matlab interp = %f \n',time_vfi1)
fprintf('Run time of vfi2 with custom interp = %f \n',time_vfi2)
fprintf('time vfi1 / time_vfi2 = %f \n',time_vfi1/time_vfi2)
%% 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
1 Comment
Stephen23
on 27 Dec 2025
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:
Accepted Answer
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)]
timeit(@() interp1(x,y,xq))
timeit(@() interp1_scal_vec(x,y,xq))
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)))
timeit(@() interp1_scal(x,y,xq(1)))
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
5 Comments
John D'Errico
on 27 Dec 2025
Edited: John D'Errico
on 27 Dec 2025
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.
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))
x = linspace(0,10,numel(x));
timeit(@() interp1(x,y,xq))
@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!
A lookup table is easy enough to do.
x = cumsum(randi(3,[1,10000]))
[min(x),max(x)]
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)
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
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.
More Answers (0)
Categories
Find more on Interpolation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
