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’)
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 locateI 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 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 interp1, linear interpolation MATLAB Answers — New Questions
