Single precision vs. Double precision experiment

[MichaelTiemann]

I was curious to see how possible/impossible it would be to compute a FHorz value function iteration using single precision arithmetic. Answer: not that hard.

Because this was just an experiment, I cobbled together the call:

ValueFnIter_FHorz_ExpAsset_e_raw(single(n_d1),single(n_d2),single(n_a1),single(n_a2),single(n_z), single(vfoptions.n_e), single(N_j), single(d_gridvals), single(d2_gridvals), single(a1_gridvals), single(a2_grid), single(z_gridvals_J), single(vfoptions.e_gridvals_J), single(pi_z_J), single(vfoptions.pi_e_J), ReturnFn, aprimeFn, Parameters, DiscountFactorParamNames, ReturnFnParamNames, aprimeFnParamNames, vfoptions);

Then made a two fixes to the called function:

% Create a vector containing all the return function parameters (in order)
ReturnFnParamsVec=CreateVectorFromParams(Parameters, ReturnFnParamNames,N_j);
% FIX FOR SINGLE PRECISION
ReturnFnParamsVec=single(ReturnFnParamsVec); V=zeros(N_a,N_z,N_e,N_j,'single','gpuArray'); Policy=zeros(N_a,N_z,N_e,N_j,'single','gpuArray');

    % Create a vector containing all the return function parameters (in order)
    ReturnFnParamsVec=CreateVectorFromParams(Parameters, ReturnFnParamNames,jj);
    DiscountFactorParamsVec=CreateVectorFromParams(Parameters, DiscountFactorParamNames,jj);
    DiscountFactorParamsVec=prod(DiscountFactorParamsVec);

    aprimeFnParamsVec=CreateVectorFromParams(Parameters, aprimeFnParamNames,jj);
    % FIX FOR SINGLE PRECISION
    aprimeFnParamsVec=single(aprimeFnParamsVec); DiscountFactorParamsVec=single(DiscountFactorParamsVec); ReturnFnParamsVec=single(ReturnFnParamsVec);

I also needed to create single-precision versions of my ReturnFn and aprimeFn.

And voila! In this case I observed that single precision runs faster than double precision (up to 2x) and uses 2x less GPU memory. This tells me that the VFI toolkit is not itself terribly opinionated about single vs. double precision. The biggest issue in my ReturnFn and aprimeFn was the fact that assigning zero to a parameter defaults to double, so changing x=0; to x=single(0); in my own code. The biggest issue likely for the toolkit is changing its calls to zeros and ones to return single vs. double arrays, and of course the construction of parameter vectors. Parameter values can remain their feral double-precision selves.

And of course users will have to be careful to create grids of the precisions they want. Though one could allow the PType mechanism to convert grids to preferred precisions when storing them in PType structures…

Make of this what you will, but I think it wouldn’t be too difficult to support larger models / faster runtimes when one is just experimenting and super-tight tolerances are not needed for publication/replication.

[aledinola]

Yes, I think the toolkit should have an options about single vs double precision. Not sure how easy it is to implement it in the hundreds of functions that are part of the toolkit. Maybe with AI is doable. However, one needs to be careful with using single precision: it works in 99% of the cases but there are instances where it creates numerical instability in an iterative algorithm like vfi.

Low-memory option reduces memory requirements but reduces speed too, so it’s not a substitute for a single precision option.

[MichaelTiemann]

Regarding the iterative instability issue: FHorz problems are, by nature, very finite in their iteration counts. Moreover, each such iteration is actually presenting its own legible data at each iteration step. InfHorz solutions are approximated by iterative solutions that could cycle through hundreds of iterations until they are deemed converged. And such convergence may fail (which does not directly happen in the FHorz case–it’s diagnosed in other ways, see Deadbeat agents and feasibility constraints).

Since the toolkit already tells one when InfHorz fails to converge, and since it’s really not a problem in the FHorz case, the safety is there.

I also suspect that the real gains will be found in the FHorz cases, because there are so many more interesting dimensions to the Household problem than the Firm problem. But perhaps that’s just my bias/experience talking.

[robertdkirkby]

Nice Demo!

In principle you would want to do, e.g., the stationary eqm, by using single-precision floating point for first 30 iterations, then move into double floating precision. I’ve been thinking about this concept for a few years because it is the way all the GPU hardware is going (apparently one of the big breakthroughs of DeepSeek in 2025 was how to use a lower interger precision than anyone else had pulled off).

For example, Policy should obviously just be stored as int16 or maybe int32, this would save a lot of memory.

The problem is that Matlab does not do type-setting for variables, so implementing it is really tricky.

If you can see a way to do it that isn’t just “three versions of every code” then I am very very interested, but I need to see how it can be practically done.

[When I periodically consider moving to Julia, this is the number one motive. The catch is the translation of all the code.]

PS. For anyone reading who is not very familiar with these concepts and why they matter. If you have double-point floating numbers they are using 64 bits (0s and 1s) to represent a number, while single-point floating numbers are using 32 bits. Obviously 64 is more accurate and uses more memory than 32. Also, because you do 64-times-64 when you multiply, it is roughly 4x as much runtime as 32-times-32, so speed differences should be quadratic in bits-per-number. In practice of course this is all just a rough guide to actual performance. The catch is of course that if the solution requires a level of accuracy that can be represented in 64 bits, but cannot be captured with 32, then using 32 won’t get you the right answer. [there is such a thing as 128, but it is boutique and not something you will ever use unless you have a specific reason to do so]

[MichaelTiemann]

Here’s a primer on floating point formats: IEEE 754 - Wikipedia

It has a nice chart showing that for values around zero, single-precision floating point runs out of representational precision around 1e-7 and double-precision runs out around 1e-16. For papers trying to prove things down to 1e-9, double-precision is the way to go. To test ideas where 1e-6 is our limit of concern, single-precision should be sufficient.

[MichaelTiemann]

Draft pull requests submitted for LifeCycleModel1, LifeCycleModel8 (Grid Interp), LifeCycleModel18 (slightly more complex), and VFIToolkit-matlab.

Next steps could include

• March forward through LifeCycleModel samples, fixing the toolkit as we go.
• Explain to Claude that Policy really should be a matrix of 'int32', not single-precision–and certainly not double-precision–floating point.

'int16' is way too small for many Policy functions. 'int32' is large enough because anyway MATLAB bottles out when array dimensions exceed a billion (which 'int32' easily handles).

[robertdkirkby]

Are you modifying life-cycle models or toolkit? Shouldn’t it only be the later?

[MichaelTiemann]

The ReturnFns need to work in single precision. Otherwise they would generate double-precision values that would have to be cleaned up EVERYWHERE. Note that

F=-Inf;

is a double-precision statement.

In any case, I think you’ll find that the changes to the LifeCycle codes are pretty simple. It took me 5 minutes max to do LifeCycleModel18.

[MichaelTiemann]

LifeCycleModel40 revealed something interesting. Here’s a piece of an index calculation from ValueFnIter_FHorz_DC2A_GI2A_nod_raw:

K>> [n2long,N_a2,N_a,20]

ans =

  1×4 single row vector

          43          21        6321          20

If we multiply the numbers as integers, we get 114157260. If we multiply them as single precision floats, we get 114157264 because we don’t have enough precision. Double precision handles it fine. So for large arrays (like we find in LifeCycleModel40), we cannot trust single precision floating point numbers. This is not a problem for the actual value functions–they lose precision, but not in this sort of catastrophic way.

The results of creating an int32 Policy were successful, 334MB vs 172MB of memory:

whos Policy
  Name        Size                        Bytes  Class       Attributes

  Policy      4x6321x21x81            344064672  gpuArray             

whos Policy
  Name        Size                        Bytes  Class       Attributes

  Policy      4x6321x21x81            172032336  gpuArray             

Some learnings: there are at least two ways in which grid interpolation complicate the construction and management of Policies. The first is the midpoint calculation:

% Turn this into the 'midpoint'
midpoints_jj=max(min(midpoints_jj,n_a1-1),2); % avoid the top end (inner), and avoid the bottom end (outer)
% midpoint is 1-by-n_a2-by-n_a1-by-n_a2-by-n_z
a1primeindexes=(midpoints_jj+(midpoints_jj-1)*n2short)+(-n2short-1:1:1+n2short)'; % aprime points either side of midpoint
% aprime possibilities are n2long-by-n_a2-by-n_a1-by-n_a2-by-n_z

MATLAB doesn’t tolerate the concept of a negative unsigned integer, so -n2short means that n2short must be an int32. We don’t really lose much because MATLAB also won’t allow indexes on the GPU to exceed the maximum positive int32 value. So there’s no benefit to using uint32 as an index.

I was going to write about this slightly more vexing appears in the calculation of the Stationary Distribution:

elseif simoptions.gridinterplayer==1
    % (a,z,1,j)
    Policy_aprime=reshape(Policy_aprime,[N_a,N_ze,1,N_j]);
    Policy_aprime=repmat(Policy_aprime,1,1,2,1);
    PolicyProbs=ones([N_a,N_ze,2,N_j],precision,'gpuArray');
    % Policy_aprime(:,:,1,:) lower grid point for a1 is unchanged
    Policy_aprime(:,:,2,:)=Policy_aprime(:,:,2,:)+1; % add one to a1, to get upper grid point

    aprimeProbs_upper=reshape(shiftdim((Policy(end-1,:,:,:)-1)/(simoptions.ngridinterp+1),1),[N_a,N_ze,1,N_j]); % probability of upper grid point (from L2 index; end-1 because end is now L2flag)
    PolicyProbs(:,:,1,:)=PolicyProbs(:,:,1,:).*(1-aprimeProbs_upper); % lower a1
    PolicyProbs(:,:,2,:)=PolicyProbs(:,:,2,:).*aprimeProbs_upper; % upper a1

and specifically how the calculation of aprimeProbs_upper goes against our int32 Policy, but I found it was easy (with a few small adjustments) to accommodate. By casting (Policy(end-1,:,:,:)-1) to double, we can complete the calculation of our PolicyProbs and thus the Stationary Distribution. So maybe int32 Policies is another thing that can come out of this!

[aledinola]

@MichaelTiemann

Regarding the iterative instability issue: FHorz problems are, by nature, very finite in their iteration counts.

Yes, of course. I had in mind infinite horizon models.

[robertdkirkby]

Looking at the git pull request to implement single precision.

My thoughts: single vs double takes a fair bit of coding, but it is all very simple coding. My main complaint from user perspective is just that the user has to specify ‘single’ in the ReturnFn, but I guess this is easy enough to learn.

It can be worthwhile for the ValueFnIter, since that is where most runtime resides.

I don’t think it will be worth the runtime in the StationaryDist/AgentDist commands, the gains will be small for the complication involved. They are typically already such a small part of the overall runtime (with sole exception being where almost all shocks are semi-exogenous and there are lots of those).

I am undecided if the runtime is worth it for the FnsToEvaluate commands. Maybe.

Your thoughts? Do you agree?

[I can’t see much point putting it into things like discretizeAR1_FarmerToda and MarkovChainMoments. Lots of code for little impact. I’d rather focus the additional code handling on the bottlenecks where they will make a real impact.]

[MichaelTiemann]

I think the principal use case here is giving your GPU a 2x RAM upgrade when models get large (and also benefiting from a speed improvement as well). I think it’s definitely a “if it’s not broke, don’t fix it” capability. If people stick to small models, they might never see it.

The reason I put it into Stationary Dist was to see about changing the basis of Policies from double to int32. I’m seeing that it’s quite easy. Now, whether these are big memory hogs or not, I’m not quite sure. What I hope to do (when I get back) is complete as much of the merge of the changes I’ve made for my large model so that I can do proper experiments with it in terms of size and speed. The quick test I did at the start confirmed that there was a speedup and size advantage in VFI, but who knows what will break when I increase the model size…

The main reason for touching discretizeAR1_FarmerToda and MarkovChainMoments was to get clean z_grid, pi_z, and the like. I really wanted to minimize the extra interventions that users would have to do (beyond the custom ReturnFn/aprimeFn/etc). So it’s not for performance, but invisibility.

Let’s take a harder look at it once I’ve got a large model to throw at it. As you can see from the many files lightly touched, Claude should be good at adding it more generally, if and when.

[aledinola]

I haven’t had time to look into the code in the pull request, but it sounds like a nice addition to the toolkit: decrease memory requirements and make the code faster.
I would look a bit into accuracy, though, in the context of a “real” model, not just simple/toy examples. For example, would it be possible to run these two examples on Robert’s github, based on two published papers and check if the results look the same:

1. Huggett and Ventura (2000) https://github.com/robertdkirkby/HuggettVentura2000
2. Chen (2010) GitHub - robertdkirkby/Chen2010housing: Rough implemention of baseline model of Chen (2010) - A Life-Cycle Analysis of Social Security with Housing · GitHub

In the interest of time, one could skip general equilibrium computation.

[robertdkirkby]

I both am and am not worried about accuracy.

If it does hurt accuracy, it can still be very useful for solving general eqm and transition paths and calibration/estimation, because you just use the fast-but-low-accuracy to get near the optimum, and then switch back to the full accuracy to get the final solution.

So while it is important to know if it hurts accuracy and how much, it can still be a feature worth implementing even though it is not as accurate. In this sense, I am happy to implement first and worry about when and where it can be used second. But you are absolutely right that looking at when and where it can be used has to be done on ‘production’ models.

It will definitely always be a feature that is optional and is off by default.

[robertdkirkby]

@MichaelTiemann

The main reason for touching discretizeAR1_FarmerToda and MarkovChainMoments was to get clean z_grid, pi_z, and the like. I really wanted to minimize the extra interventions that users would have to do (beyond the custom ReturnFn/aprimeFn/etc). So it’s not for performance, but invisibility.

I had figured that the best way to make it ‘invisible’ and minimize the extra interventions required of user was that z_grid and pi_z would have to be converted to single precision inside the ValueFnIter command. It was in this sense that I didn’t seem much point implementing it to things like the discretization commands. We are going to have to check and convert to single anyway.

[MichaelTiemann]

Would you rather see an implementation that supports the construction of single precision/int32 policies (by keeping going when it sees them), or introducing some vfoptions/simoptions/heteroagent/transpathoptions parameters that tells the toolkit “Make it so”?

The one thing we cannot do (unfortunately) is query type within the arrayfun context, so there we really do require the user to code special functions that don’t break the format mixing rules.

BTW, I’m running a baseline of HV2000 and I’m struck by how low the %GPU utilization is in the baseline condition. Model4 takes 24GB of GPU ram. Since a lot of time is spent iterating through the StationaryDist calculations, I’ll be interested to see whether int32 Policies shine compared with the default double policies.

[robertdkirkby]

The HV2000 codes for model 4 have vfoptions.lowmemory=1 turned on as otherwise it errors on my desktop. I suspect this is the cause of the " I’m struck by how low the %GPU utilization is in the baseline condition", which is to say it is wasting utilization to avoid going out-of-memory. I’m just guessing here.

[robertdkirkby]

@MichaelTiemann

introducing some vfoptions/simoptions/heteroagent/transpathoptions parameters that tells the toolkit “Make it so”?

I was thinking along these lines, a vfoptions.singlefloat=0 by default, can can set to 1. Seems to fit with rest of the toolkit approach.

[MichaelTiemann]

Easily done in the next round of iterations. I’m continuing to learn where promotion rules (single vs double, float vs int) work either in my favor or against. I’ll be stripping out dozens of superfluous promotions and adding a few needed ones to hopefully get things nice and tight.

[MichaelTiemann]

Here’s a look at all the time spent in StationaryDist calculations compared with VFI during the GE calculations.

The horizonal scale of the windows is 4 minutes wide, so 6 VFIs per 4 minutes is 45 sec for each GE estimation. And of those 45 seconds, looks like more than half that time is non-VFI time.