Two Endogenous States - Divide&Conquer + Grid Interpolation Layer

robertdkirkby · November 14, 2025, 12:52am

For anyone out there using a life-cycle model with two (standard) endogenous states, such as a model with ‘assets’ and ‘house’, you can now turn on both divide-and-conquer and grid interpolation layer. Both are applied to only the first of the two endogenous states (so make sure the one that needs most grid points is the first one).

You just set,

vfoptions.divideandconquer=1;
vfoptions.gridinterplayer=1;
vfoptions.ngridinterp=15;
simoptions.gridinterplayer=vfoptions.gridinterplayer;
simoptions.ngridinterp=vfoptions.ngridinterp;

and the life-cycle models will be faster, less memory, and more accurate!

You could also play with
vfoptions.level1n=7;
which controls the number of points used in the first layer when doing divide-and-conquer. Better choices can be faster, but it has a default so setting it is optional.

Makes these models much easier and faster to solve. Works with all life-cycle models with two standard endogenous shocks.

PS. Probably works with three endogenous state (applying divide-and-conquer and grid interpolation layer to the first state), but I have not tested this so mileage may vary.

For two endogenous states I did lots of testing with all shock types etc.: GitHub - robertdkirkby/vfitoolkitTests: Various tests of VFI Toolkit features

aledinola · November 14, 2025, 9:20am

Nice! I will test this new feature on GitHub - aledinola/Chen2010housing: Rough implemention of baseline model of Chen (2010) - A Life-Cycle Analysis of Social Security with Housing

aledinola · November 14, 2025, 10:55pm

Should we also specify the number of points used in the second layer? If divide and conquer applies only to the first endogenous state, I think the second layer should use exactly n_a(2) points

I ask this question since in my repo on Chen 2010 I have this line:

vfoptions.level1n=[9; n_a(2)]

and I’m trying to remember what I wanted to do

robertdkirkby · November 14, 2025, 11:17pm

You can just set vfoptions.level1n=9 because it is only applied to the first of the two endogenous states. If you set vfoptions.level1n=[9,n_a(2)] it will anyway just do the same thing as the second entry is ignored.

Similarly, I only set vfoptions.ngridinterp=15; as it will only be applied to the first of the two endogenous states.

[I need to clean the internals up a little bit around this, but this is now the intended way to set up and any internals that don’t line up already will be brought to line up with it.]

aledinola · November 15, 2025, 12:00am

I ran again the Chen (2010) from my repo (see link above) but there is potentially a bug.

The model has n_d=0, n_a=[401,11], n_z=7, N_j=66. So the state space is (a,h,z,j) and correspondingly the size of the value function is

size(V)

ans =

   401    11     7    66

The size of Policy is

size(Policy)

ans =

     2   401    11     7    66

where Policy(1,:,:,:,:) is the policy function for next-period assets and Policy(2,:,:,:,:) is the policy function for housing. So far so good.
But when I look at the size of the agents distribution,

 size(StationaryDist)

ans =

       30877          66

which is not correct

aledinola · November 15, 2025, 12:15am

Unrelated to the previous problem: when using interpolation I got this error:

Error using gpuArray/reshape
Sparse gpuArray matrices are not supported for this function.

Error in StationaryDist_FHorz_Iteration_nProbs_raw (line 26)
    StationaryDist_jj=reshape(Gammatranspose*StationaryDist_jj,[N_a,N_z]);
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in StationaryDist_FHorz_Case1 (line 283)
            StationaryDist=StationaryDist_FHorz_Iteration_nProbs_raw(jequaloneDist,AgeWeightParamNames,Policy_aprime,PolicyProbs,2,N_a,N_z,N_j,pi_z_J,Parameters);
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Chen2010 (line 276)
StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

I think the reason is that Matlab R2024b, which I am using, does not support reshape for sparse gpuarrays. It is a bit unfortunate because I don’t want to update to R2025 since it quite buggy On a more serious note, I noticed a deterioration in performance on gpu code when using R2025b instead of R2024b. I will do some serious testing when I have time. But yeah, this is why I am reluctant to upgrade to R2025. Please do not rewrite the whole Tan improvement using sparse gpuarrays

aledinola · November 15, 2025, 12:21am

Found some inconsistency with vfoptions.divideandconquer = 1. If I set vfoptions.level1n = 9 there is an error:

Error using ValueFnIter_FHorz_DC_GI (line 119)
With two endogenous states, can only do divide-and-conquer in the first endogenous state (not in both)

Error in ValueFnIter_Case1_FHorz (line 551)
    [V,Policy]=ValueFnIter_FHorz_DC_GI(n_d, n_a, n_z, N_j, d_grid, a_grid, z_gridvals_J, pi_z_J, ReturnFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, vfoptions);
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Chen2010 (line 258)
[V,Policy]=ValueFnIter_Case1_FHorz(n_d,n_a,n_z,N_j,d_grid, a_grid, z_grid, pi_z, ReturnFn, Params, DiscountFactorParamNames, [], vfoptions);
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

but if I set vfoptions.level1n = [9,n_a(2)] then all is good. I report this problem since in your post above you wrote that I can just set vfoptions.level1n = 9, but it doesn’t seem to be the case in this code.

Divide and conquer delivers a considerable performance improvement and the run time shrinks from 14 seconds to 6 seconds on my computer, which is nice!

robertdkirkby · November 15, 2025, 12:21am

Hmm, 30877=401117, so it was calculated correctly but then left in the wrong size when outputting.

Fixed now (I hadn’t noticed as you could still input it to AllStats and it worked just fine)

robertdkirkby · November 15, 2025, 12:23am

Will fix the vfoptions.level1n setup sometime in the next week or two (seems silly to make the user declare two numbers when the second is ignored, will be easier to understand that it only applies to the 1st endo state if it is just one number).

robertdkirkby · November 15, 2025, 12:25am

Actually I decided that none of the stationary dist commands will be based on sparse gpu matrices. The speed gains are minor and you get more out of memory errors, just wasn’t worth the trade-off given how little of the runtime of an overall model this already is. I cleaned things back up so that none of them use the sparse gpu matrix but clearly I missed this particular combo.

aledinola · November 15, 2025, 12:30am

Wow, thanks for adjusting this in real time!
Overall the combination of divide and conquer and linear interpolation is good for accuracy and speed. I attach below the running times of the Chen (2010) model, only for value function iteration:

Divide and conquer reduces the run times by half and adding interpolation only slightly increases the run time. Interpolation gives more accurate results as can be seen by plotting life time profiles (not shown here)

robertdkirkby · November 15, 2025, 12:32am

Agent dist should be fixed now (use sparse matrix, not sparse gpu matrix)

robertdkirkby · November 15, 2025, 12:34am

As well as the runtime gains there are massive loosening of gpu memory bottlenecks (so you can not only solve faster, you can just solve larger models that were not previously possible).

aledinola · November 15, 2025, 12:39am

The bugs reported in this post are fixed now

Uhmm I think there are still some issues. I ran the code on R2024b and R202a. On R2024b it gives the same error as before (sparse gpuarray not supported).

Error when using R2024b

Error using gpuArray/reshape
Sparse gpuArray matrices are not supported for this function.

Error in StationaryDist_FHorz_Iteration_nProbs_raw (line 30)
    StationaryDist_jj=reshape(StationaryDist_jj*pi_z,[N_a*N_z,1]);
                      ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in StationaryDist_FHorz_Case1 (line 286)
            StationaryDist=StationaryDist_FHorz_Iteration_nProbs_raw(jequaloneDist,AgeWeightParamNames,Policy_aprime,PolicyProbs,2,N_a,N_z,N_j,pi_z_J,Parameters);
                           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Chen2010 (line 276)
StationaryDist=StationaryDist_FHorz_Case1(jequaloneDist,AgeWeightParamNames,Policy,n_d,n_a,n_z,N_j,pi_z,Params,simoptions);
               ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

Instead, on R2025a it runs OK! I have realized that there was a typo in my main script (at some point I clean simoptions and I forget to pass again

simoptions.gridinterplayer = vfoptions.gridinterplayer;
simoptions.ngridinterp = vfoptions.ngridinterp;

so the error reported below is gone

Please disregard

Error using gpuArray/reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.

Error in PolicyInd2Val_FHorz (line 178)
            Policy=reshape(Policy,[l_aprime,N_a*N_z*N_j]);
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in LifeCycleProfiles_FHorz_Case1 (line 130)
    PolicyValues=PolicyInd2Val_FHorz(Policy,n_d,n_a,n_z,N_j,d_grid,a_grid,simoptions,1);
                 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Error in Chen2010 (line 331)
AgeConditionalStats2=LifeCycleProfiles_FHorz_Case1(StationaryDist,Policy,FnsToEvaluate,Params,[],n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,simoptions);
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

aledinola · November 15, 2025, 12:41am

Just to add that I have been using these options:

vfoptions.verbose = 1;
vfoptions.lowmemory = 0;
vfoptions.divideandconquer = 1;
vfoptions.level1n = [9,n_a(2)];
vfoptions.gridinterplayer = 1;
vfoptions.ngridinterp = 15;

simoptions=struct(); % Use default options for solving for stationary distribution
simoptions.verbose = 0;
simoptions.gridinterplayer = vfoptions.gridinterplayer;
simoptions.ngridinterp = vfoptions.ngridinterp;

robertdkirkby · November 15, 2025, 12:46am

I think I just fixed the R2024a error.

The PolicyInd2Val one will have to wait until Monday as cannot do that from my laptop without GPU.

aledinola · November 15, 2025, 12:52am

Great, thanks! Now the code in R2024a runs well and the error in R2025a was due to me messing up something (I have updated the post)
So everything runs well in both versions of Matlab, with and without gridinterplayer!
Have a nice weekend!

robertdkirkby · November 15, 2025, 12:54am

Perfect! I had just fixed the PolicyInd2Vals this morning and was peeved that it still wasn’t working happy to hear it is