I have done a bunch of runtimes for “sparse(lower grid)+sparse(upper grid)” versus “sparse([lower grid, upper grid)”
Most of the time the “sparse([lower grid, upper grid)” is faster, but not always. The average runtime gain is less than 20%, so the difference is small (in small problems the runtime gains were faster, but most of the time the runtime gain is minor). And there are plenty of runs where doing them separate was faster (I didn’t look at individual runs, but averages over 100 runs, and even there the ratio bounces all over).
So I am going to leave the codes doing “sparse([lower grid, upper grid)”. But really, it is pretty irrelevant which of the two ways you code it.
[I also did one adding 4 things, rather than just two. Much the same outcome for relative runtimes.]
That code for discretize() actually has the same problem I did, namely people off the top of the grid get bumped to the bottom.
For example,
xval=[3,9.9,10,15];
n_x=5; % number of grid points
x_grid=7+(1:1:n_x)'; % integers from 7 to 12
indx=discretize(xval,x_grid);
indx = max(indx,1); % Deal with xval<x_grid(1)
indx = min(indx,n_x-1); % Deal with xval>x_grid(end)
which is just the same thing as your above code, puts the 15 value in first index. Reason is that discretize() returns NaN for points outside the range of edges. And then indx = max(indx,1); overwrites the NaN with 1, and so indx = min(indx,n_x-1); does not end up doing anything.
Fix is easy, just replace last three lines with
indx=discretize(xval,x_grid);
indx(xval<x_grid(1))=1; % Deal with xval<x_grid(1)
indx(xval>x_grid(end))=n_x-1; % Deal with xval>x_grid(end)
Did some runtimes for discretize() versus my trick with max() [both explained above]. All runtimes are on GPU (as that is how I would always do them in practice).
For small problems, max() was about 3x faster. But for large problems discretize() is 10x faster.
I will have to think about which is going to be faster in practice, based on sizes of the things I am putting on the grid. Maybe I need an if-statement to see how big it is and then decide which to do based on that. That said, for the small problems, it is 3x of almost nothing, whereas for the big problems it is 10x of something substantial, so I might just go with discretize as it is only worse when worse is not very important.
Thanks a lot for spotting this! I have to review some of my codes then, because I was using this function in several projects
A side note: before I discovered discretize, I was using histc (it does more or less the same thing), but unfortunately histc is now deprecated by Matlab
Here is the updated version of my function
function [indx,omega] = find_loc_vec2(x_grid,xi)
%Find indx s.t. x_grid(indx)<=xi<x_grid(indx+1)
%for indx=1,..,N-1
% omega is the weight on x_grid(indx)
n_x = size(x_grid,1);
% For each 'xi', get the position of the 'x' element bounding it on the left [p x 1]
indx = discretize(xi,x_grid);
indx(xi<x_grid(1))=1; % Deal with xi<x_grid(1)
indx(xi>x_grid(end))=n_x-1; % Deal with xi>x_grid(end)
%Weight on x_grid(indx)
omega = (x_grid(indx+1)-xi)./(x_grid(indx+1)-x_grid(indx));
omega = max(min(omega,1),0);
end %end function "find_loc_vec2"
Played with runtimes for discretize() a bit more. Size of x_grid only mildly important. When xval is 1000 points or less, the max() is faster, but from about 100000 points or more on xval discretize() is faster. In between it can vary.
VFI Toolkit essentially does this operation of putting points onto grid in two places. Once ‘CreateaprimeFn’ (in various versions), which can be just a few hundred points for xi. Once ‘CreateaprimePolicy’ (in various versions) where xi is going to be millions of points. So I changed all the ‘CreateaprimePolicy’ to discretize(). Put an if-statement in the ‘CreateaprimeFn’ and choose which based on if xval has lots of points (runtimes suggest the if-statement slowdown is tiny, and avoiding being substantially wrong on which of max() and discretize() to use is worth the if-statement). [The ‘CreateaprimeFn’ are used in value function iteration, whereas ‘CreateaprimePolicy’ is in agent distribution iteration.]
That said, if you just always use discretize() the maximum potential losses in runtime are only around 10^(-3) or 10^(-4) seconds [per time you call it]. So if I wanted an easy option I would go with just using discretize() all the time. [lost runtime from using max() when discretize() is better can be hundreths of a second, even tenths of a second if xval is large]
[My exact if-statement is that numel(xval)*numel(x_grid)<1000000, then use max(), else use discretize(). Seems to always get within about 10% of the best possible runtimes.]
The comparison was [~,indx]=max(xval>x_grid); % upper grid point indx=indx-1; % lower grid point indx(indx==1); indx(xval<x_grid(1))=1; indx(xval>x_grid(end))=n_x-1;
Versus, indx=discretize(xval,x_grid); indx(xval<x_grid(1))=1; indx(xval>x_grid(end))=n_x-1;
UPDATE
I found the problem. When computing the distribution I was doing
for j=2:N_j
Distribution at age j = f(Distribution at j-1, Policy at j-1, j)
end
while it should have been
for j=2:N_j
Distribution at age j = f(Distribution at j-1, Policy at j-1, j-1)
end
The symbol f above denotes the “operations” to update the distribution.
Now the toolkit and my codes deliver exactly the same results
It’s about the law of motion for human capital, coded in the function f_HC_accum. Human capital at age j depends on human capital at age j-1, labor supply at age j-1 and age j-1. In my codes I had age j, so it was not consistent.
