Grid Interpolation layer

1

This post has been superseeded by:

I recommend you read that instead of the below which is outdated

I’ve started implementing interpolation on linearly spaced points between the a_grid points. It will be a long time while this is gradually rolled out across the toolkit, but I just wanted to briefly describe the concept here.

As user you just set
vfoptions.gridinterplayer=1
which tells toolkit to use a ‘gridded interpolation layer’, and then use
vfoptions.ngridinterp=5
to say how many points to use for the interpolation.
And then you put the same into simoptions.

Idea is that if you have, say, n_a=20 (and so a_grid is 20-by-1), then between every two consecutive points on the a_grid it will put vfoptions.ngridinterp evenly spaced points. Conceptually, the problem we want to solve is using aprime_grid=interp1(1:1:N_a,a_grid,linspace(1,N_a,N_a+(N_a-1)*vfoptions.ngridinterpt)); (N_a=prod(n_a);) We evaluate the ReturnFn exactly on this aprime_grid, but since a_grid is smaller and the value function will only be on a_grid we have to interpolate E[V] (the expected next period value function) from a_grid onto aprime_grid and this is one with linear interpolation.

Rather than do this in one hit, the codes do it in layers (to avoid the large aprime_grid causing out-of-memory errors). [Actually, I should try the one-hit as well, see how runtimes compare.]

First, the codes just do the usual thing of solving ‘pure discretization’ with the 20 points on next period assets (the same points used for current assets, so a_grid). [Or can use divide-and-conquer for this first part.] This gives us the optimal aprime index (on the a_grid) [Actually slightly different to usual, as solve for the optimal aprime conditional on the possible choices of d]

In a second ‘layer’, the codes consider the range for aprime from one a_grid point below the optimal aprime index to one point above, and puts ‘vfoptions.ngridinterp’ points between these (so we get a_grid-point-below-optimal then vfoptions.ngridinterp evenly spaced points then a_grid-point-optimal then vfoptions.ngridinterp evenly spaced points then a_grid-point-above-optimal), so consider 2*vfoptions.ngridinterp+3 points. This is done conditional on any decision variables, d. Finally find the optimal of these, using exact valuation ReturnFn, and linear interpolation of E[V].

The output for the value function is essentially the same, and the interpretation is identical. The output for the policy function now contains two indexes for the aprime value, the first is a ‘lower grid point’ and the second is from 1-to-(vfoptions.ngridinterp+2) and the two endpoints represent the lower grid point and the upper grid point (and the rest count the evenly spaced points we interpolated onto).

I’ll try write up a pseudo-code sometime soon.

This feature is likely to take at least a full year for me to roll out across VFI Toolkit, so will be a while before you can use it. But the good news is it works, and it is coming your way :smiley:

2 Likes
2

This is a great development.

In terms of performance, I suspect that the bottleneck in the code are the lines where

EVinterp=interp1(a_grid,EV,aprime_grid)

Instead of interp1, one can use the following function

 
function [aprimeIndexes,aprimeProbs] = interp_toolkit(a_grid,aprimeVals)
 
 
%aprimeVals=reshape(aprimeVals,[1,N_d*N_a]);
 
aprimeVals=aprimeVals';%,[1,N_d*N_a]);
 
 
a_griddiff=a_grid(2:end)-a_grid(1:end-1); % Distance between point and the next point
 
 
temp=a_grid-aprimeVals;
 
temp(temp>0)=1; % Equals 1 when a_grid is greater than aprimeVals
 
 
[~,aprimeIndexes]=max(temp,[],1); % Keep the dimension corresponding to aprimeVals, minimize over the a_grid dimension
 
% Note, this is going to find the 'first' grid point such that aprimeVals is smaller than or equal to that grid point
 
% This is the 'upper' grid point
 
 
% Switch to lower grid point index
 
aprimeIndexes=aprimeIndexes-1;
 
aprimeIndexes(aprimeIndexes==0)=1;
 
 
% Now, find the probabilities
 
aprime_residual=aprimeVals'-a_grid(aprimeIndexes);
 
% Probability of the 'lower' points
 
aprimeProbs=1-aprime_residual./a_griddiff(aprimeIndexes);
 
 
% Those points which tried to leave the top of the grid have probability 1 of the 'upper' point (0 of lower point)
 
offTopOfGrid=(aprimeVals>=a_grid(end));
 
aprimeProbs(offTopOfGrid)=0;
 
% Those points which tried to leave the bottom of the grid have probability 0 of the 'upper' point (1 of lower point)
 
offBottomOfGrid=(aprimeVals<=a_grid(1));
 
aprimeProbs(offBottomOfGrid)=1;
 
 
 
end %end function

In a two step way:

[ind,prob] = interp_toolkit(a_grid,aprime_grid);
EV_interp = prob.*EV(ind) + (1-prob).*EV(ind+1)

There are two advantages in this approach:

(1) Avoid the overhead of calling the Matlab function interp1, which has a lot of internal checks that make it slow
(2) By splitting the interpolation in two steps, one can often precompute ind and prob

3

I’ll try this out next week, set up some runtime tests on gpu and see how it goes (and if a minor modification does what I am after, which is that I have EV over a_rough-by-z and I interpolate onto a_fine-by-z).

PS. Should be easier to improve a little further as I know that the ‘additional’ interpolation points in a_fine are evenly spaced between the points in a_rough, so I already essentially know what the probabilities will be and so can calculate them in a trivial manner.

1 Like
4

That’s nice: if you know that the points are equally spaced, you can save some time in the calculation of ind and prob.

The function that I shared in the previous post works for an arbitrarily spaced grid. I had actually copied it some time ago from the toolkit.

For equally spaced points, I typically use this other function:

	
% table lookup: find position of X in grid gridX;
 
% gridX must be EQUIDISTANT !!
 
% output:
 
%     iPos: index of grid point left to X
 
%     weightHi: (X-gridX(iPos))./ (gridX(iPos+1) - gridX(iPos))
 
function [iPos weightHi] = lookup_equidist(gridX,X)
 
n = length(gridX);
 
stepx = gridX(2)-gridX(1);
 
% position of largest grid point that is smaller than X:
 
iPos = max(1,floor( (X-gridX(1))/stepx ) + 1);
 
iPos = min(iPos,n-1); %iPos must not be n
 
% probability to go to high end of interval:
 
weightHi = (X-gridX(iPos))./stepx;
 
% test whether X are between gridX(iPos) and gridX(iPos+1):
 
% this tests whether X must be within grid
 
assert(all(weightHi>=-1e-10 & weightHi<=1+1e-10));
 
end
5

I did a bunch of runtime tests. Idea is that I have a matrix EV, and I want to interpolate the first dimension (as EV is over (a,z), and I just want to interpolate the a). interp1() is fastest, when I precompute the index and probability the runtimes are slower.

This was true on cpu (which is really slow) and on gpu (which is what I actually care about). I tried adding a third dimension to EV (still just interpolate the first dimension), and it remained the case that interp1() was fastest.

You can see all the tests here:

PS. I wouldn’t be surprised if this was untrue in Matlab a decade ago, and then interp1() got improved since.

1 Like
6

This is really interesting, thanks for setting up these tests.

I’m surprised about this finding for the cpu. I will go over your tests later, since this might give some information useful for my codes on the cpu (for other projects, unrelated to the toolkit).

7

I didn’t try much with the cpu (just one or two combos of grid sizes), mostly focused on the gpu. So may or may not be more generally the case for cpu.

PS. In principle the way interp1() is used here would allow you to do Akima spline interpolation, but apparently that is not yet supported for gpu.

1 Like
8

I have redone the test for cpu comparing interp1, interp1q, griddedInterpolant and interpn. It turns out that griddedInterpolant and interpn are faster than interp1. interp1q is really much slower.
However my test is a bit different from yours. I will do the GPU and report the run times.

disp('interp1 -- interp1q -- griddedInterpolant -- interpn')
mean(timer_interp)
ans =

    0.0104    0.0224    0.0076    0.0089

clear,clc,close all

S = 100;
w = 1;
r = 0.04;
sigma = 2;
f_util = @(c) c.^(1-sigma)/(1-sigma);
n_a = 5000;
n_aprime = 30000;
n_z = 99;

a_grid = linspace(0,30,n_a)';    % (a,1)
z_grid = linspace(0.5,1.5,n_z);  % (1,z)
EV     = f_util(w*z_grid+r*a_grid);
aprime_grid = linspace(0,30,n_aprime)';

%figure
%plot(a_grid,EV(:,1),a_grid,EV(:,n_z))

timer_interp=zeros(S,4);

%% interp1
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp1=interp1(a_grid,EV,aprime_grid,'linear');
    timer_interp(ii,1)=toc;
end

if ~isequal(size(EV_interp1),[n_aprime,n_z])
    error('EV_interp1 has wrong size!')
end

%% interp1q
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp2=interp1q(a_grid,EV,aprime_grid);
    timer_interp(ii,2)=toc;
end

if ~isequal(size(EV_interp2),[n_aprime,n_z])
    error('EV_interp2 has wrong size!')
end

%% griddedInterpolant
base_interp = griddedInterpolant(a_grid,EV,'linear','linear');
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp3=base_interp(aprime_grid);
    timer_interp(ii,3)=toc;
end

if ~isequal(size(EV_interp3),[n_aprime,n_z])
    error('EV_interp3 has wrong size!')
end

%% interpn
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp4=interpn(a_grid,EV,aprime_grid,'linear');
    timer_interp(ii,4)=toc;
end

if ~isequal(size(EV_interp4),[n_aprime,n_z])
    error('EV_interp3 has wrong size!')
end

err_1_2 = max(abs(EV_interp2-EV_interp1),[],"all");
err_1_3 = max(abs(EV_interp3-EV_interp1),[],"all");
err_1_4 = max(abs(EV_interp4-EV_interp1),[],"all");
if max([err_1_4,err_1_3,err_1_2])>1e-10
    error('Results are not consistent!')
end

disp('interp1 -- interp1q -- griddedInterpolant -- interpn')
mean(timer_interp)
1 Like
9

Results on GPU

I dropped interpn since it does not seem to work with gpu arrays. It seems that griddedInterpolant is the fastest, also on GPU. Surprisingly, interp1q becomes faster than interp1 (but slower than griddedInterpolant).
Now the question is why these results for the GPU are different from yours :slight_smile: It would be because (1) different computers (the gpu on my laptop is not so good), (2) different parameters (grid sizes, dimension of arrays, etc.) I tried only EV with two dimensions, a and z. My test is very simple:

EV_interp = interp(a_grid,EV,aprime_grid)

where a_grid is [n_a,1], EV is [n_a,n_z], aprime_grid is [n_aprime,1] where n_aprime>>n_a and hence EV_interp is [n_aprime,n_z].

interp1 -- interp1q -- griddedInterpolant

ans =

                9.6662e-05                3.7056e-05                1.2957e-05

Test code

clear,clc,close all

S = 100;
w = 1;
r = 0.04;
sigma = 2;
f_util = @(c) c.^(1-sigma)/(1-sigma);
n_a = 5000;
n_aprime = 30000;
n_z = 99;

a_grid = gpuArray(linspace(0,30,n_a)');    % (a,1)
z_grid = gpuArray(linspace(0.5,1.5,n_z));  % (1,z)
EV     = f_util(w*z_grid+r*a_grid);
aprime_grid = gpuArray(linspace(0,30,n_aprime)');

%figure
%plot(a_grid,EV(:,1),a_grid,EV(:,n_z))

timer_interp=zeros(S,3);

%% interp1
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp1=interp1(a_grid,EV,aprime_grid);
    timer_interp(ii,1)=toc;
end

if ~isequal(size(EV_interp1),[n_aprime,n_z])
    error('EV_interp1 has wrong size!')
end

%% interp1q
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp2=interp1q(a_grid,EV,aprime_grid);
    timer_interp(ii,2)=toc;
end

if ~isequal(size(EV_interp2),[n_aprime,n_z])
    error('EV_interp2 has wrong size!')
end

%% griddedInterpolant
base_interp = griddedInterpolant(a_grid,EV,'linear','linear');
for ii=1:S
    tic;
    % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
    EV_interp3=base_interp(aprime_grid);
    timer_interp(ii,3)=toc;
end

if ~isequal(size(EV_interp3),[n_aprime,n_z])
    error('EV_interp3 has wrong size!')
end

%% interpn
% for ii=1:S
%     tic;
%     % a_grid(n_a,1), EV(n_a,n_z), aprime_grid(n_aprime,1) ---> EV_interp1(n_aprime,n_z)
%     EV_interp4=interpn(a_grid,EV,aprime_grid,'linear');
%     timer_interp(ii,4)=toc;
% end

% if ~isequal(size(EV_interp4),[n_aprime,n_z])
%     error('EV_interp3 has wrong size!')
% end

err_1_2 = max(abs(EV_interp2-EV_interp1),[],"all");
err_1_3 = max(abs(EV_interp3-EV_interp1),[],"all");
%err_1_4 = max(abs(EV_interp4-EV_interp1),[],"all");
if max([err_1_3,err_1_2])>1e-10
    error('Results are not consistent!')
end

disp('interp1 -- interp1q -- griddedInterpolant')
mean(timer_interp)
10

I tried griddedInterpolant(). This has two steps, first step ‘creates’ and second step ‘evaluates’. The ‘evaluates’ step of griddedInterpolant() is faster than interp1(), but interp1() is faster than doing both steps of griddedInterpolant() together.

Updated my codes to include this case.

For how the toolkit uses interpolation, it is the ‘both steps’ that is relevant, and so interp1() is faster for the purposes of toolkit. [I do interpolation repeatedly, but each time it is for a different expected value fn, e.g., a different age in FHorz. For any specific expected value fn, there is just a single interpolation operation being performed.]

At least at present, interpn() and interp1q() don’t do GPUs (and later doesn’t do interpolating a 2D matrix in the first dimension). So I didn’t try those out.

1 Like
11

Thanks for the clarification. Not sure that interp1q() does not do GPU, though. In my test code I call interp1q() passing gpu arrays as inputs and it works fine.

12

Comment on grid interpolation in models with no shocks: when using grid interpolation, because some policy choices are ‘between’ grid points, when you put them back onto the grid (using linear interpolation to put them onto both the two nearest points, with weights based on where between those points you are) as part of the agent distribution you ‘create’ variance across agents where in principle there is none. This means that the model, instead of having, e.g., 0 variance of assets conditional on age (if there is 0 in the initial dist) will have something like a variance of assets conditional on age of 10^(-3) or 10^(-6). Just something to be aware of when using grid interpolation in models without shocks.

1 Like
13

Grid interpolation layer now works with standard one endogenous state FHorz models. So with/without d,z,e,semiz. You can use grid interpolation layer in combination with divide-and-conquer, or just on its own.

1 Like
14

Great news! Quick question: for a finite horizon model with two d variables, one a variable, two semiz variables, one exogenous z and N_j periods we have

size(Policy)= [4,n_a,n_semiz(1),n_semiz(2),n_z,N_j]

?
The interpretation is

Policy_1 = \arg \max_{d_1} F_j(d_1,d_2,a',a,sz_1,sz_2,z)+\beta EV_{j+1}(a',sz_1',sz_2',z')
Policy_2 = \arg \max_{d_2} F_j(d_1,d_2,a',a,sz_1,sz_2,z)+\beta EV_{j+1}(a',sz_1',sz_2',z')
Policy_3 = \arg \max_{a'} F_j(d_1,d_2,a',a,sz_1,sz_2,z)+\beta EV_{j+1}(a',sz_1',sz_2',z')

and \arg \max_{a'} is the left grid point for the optimal choice of a'. But what is Policy_4? There should also be a “weight” on the left grid point, no? So that if a'^* is the optimal choice, then weight*agrid(Policy_4)+(1-weight)*agrid((Policy_4+1)) = a'^*.
Maybe my interpretation is not correct since you are restricting a’ to be on a finer grid.

15

When using grid interpolation layer, where Policy would normally contain the index for ‘aprime’, now it will contain both the ‘lower grid point index for aprime’ and an ‘index on the interpolation layer’.

So the ‘lower grid point index for aprime’ will be a number from 1 to n_a-1 (-1 as obviously the max grid point is not able to be the ‘lower’ grid point). While the ‘index on the interpolation layer’ will be a number from 1 (the lower grid point) to ngridinterp+1 (the upper grid point); there are vfoptions.ngridinterp between the two grid points.

In your case, Policy(3,:) is the lower grid point, and Policy(4,:) is the interpolation-layer index. Note that Policy(4,:) can be switched to a ‘probability of the lower grid point’ simply by calculating Prob_lowergridpoint=1-(Policy(4,:)-1)/ngridinterp. [The ‘upper grid point’ is just ‘1+the lower grid point’, and the ‘probability of the upper grid point’ is just ‘1-prob of the lower grid point’.]

Toolkit keeps the ‘interpolation-layer index’ rather than just the probability, as the index is what gets created originally and you still want it later for evaluating model stats, or switching to policy values. So makes more sense to keep it and only change to probability inside the agent dist commands (so far at least this is the only place you want the probabilities/weights).

[Internally I create the ‘fine’ grid, and then just index this. Rather than play around with rough grid and weighted sums.]

1 Like
16

It would be great to have a pseudo-code available, indeed.

I was planning to write my own code for interpolation in infinite horizon models and the pseudo code would help

17

psuedocode in pdf linked at

1 Like