Too many derivatives

For a model with a large number of variables, the number of columns of 3rd order derivatives used at the end of dynamic.json and static.json is larger than an Int32 and appears as a negative number.

This raises a more general question: is there in the code a protection against computing an unreasonable number of (higher order0 derivatives when asking by mistake or lack of knowledge too high an order of approximation for a large model?

assigned to @sebastien

There is already a check in the Dynare++ code, see https://git.dynare.org/Dynare/dynare/-/blob/master/dynare++/tl/cc/tl_static.cc#L53

I guess we could implement a similar check at the preprocessor level. What do you think?

Yes, that should do the job.

@sebastien I may not understand the check in Dynare++, but I have the impression that the test does not depend on the nature of the model at hand. Does it take into account the sparsity, w.r.t to the nonlinearity, of the model? If we have a very large model with very few nonlinear equations, we should be able to solve models at higher orders, is this true here?

There is nothing in Dynare++ that separates the linear from the nonlinear part of the model. Ferhat implemented something like that in the block-decomposition / bytecode stuff, but it is unrelated.

I did not mean that we should separate the linear and nonlinear parts. My point is that this test seems, but I may be wrong, overly restrictive. If most of the higher derivatives are nil, why should the preprocessor throw an error on the sole basis of the model's dimension (I am not 100% sure that it is what is done here)?

Because indices in the matrix will overflow, whether the matrix is dense or sparse.

OK. I did not understand how these objects are stored. Why do we only consider Int32 integers for the indices?

Historically BLAS, LAPACK and SuiteSparse were only providing an interface with 32-bit integers for indices. So we use 32-bit indices everywhere in our code. Nowadays, there exist versions of those libraries with 64-bit integers for indices, so it would be possible to upgrade all our indices. This is a project in itself, because you basically have to go through all our C++ and Fortran codebase (preprocessor, MEX files) to make the change at all relevant places. It may even be necessary to update some MATLAB code, if the Int32 type is used.

Some more thoughts on this issue:

• If
  n
  is the number of dynamic variables, and
  k
  is the derivation order, then the overflow happens when
  n^k ⩾ 2^{31}
  (we’re using signed integers, so the maximum value is
  2^{31}-1
  ). With
  k=3
  , the maximum
  n
  is as low as
  1290
  .
• In the JSON format, there is no integer type; all numbers are double precision floats. So the maximum integer that can be represented exactly is slightly above
  2^{53}
  , and we are hit by this limit since we currently pass column indices as numbers. This limit could be lifted by passing column indices as strings, but the consumer code would of course have to be adapted.
• There is a similar issue in the MATLAB/Octave dynamic file (both the .m and the use_dll versions), since in MATLAB/Octave the default numeric type is a double precision float. This could be worked around by using the int64 type instead.
• Julia by default handles big integers, so the above problem does not affect it. However, since the preprocessor internally uses 32-bit integers to compute the column indices before printing them to the dynamic .jl file, the indices are wrong. This should be easy to fix in the preprocessor.
• The bytecode is not really affected, since only first-order derivatives can be stored in this format.
• Even if we manage to raise the limit to 64-bit in the dynamic file in all formats, this won’t be the end of the story. The Dynare++ codebase also needs to be adapted, which is a big task. Then, we need to make sure that the lower-level libraries (BLAS, LAPACK, SuiteSparse) accept 64-bit indices. In particular, Octave in Debian currently links against 32-bit versions of those; for Windows, the default Octave package also uses 32-bit indices, but there is a 64-bit alternative.
• The block decomposition could in theory help when the constraint is binding. But the current implementation of block decomposition only computes the first order derivatives, so in practice it does not help.

I think that for the time being I’m simply going to verify that

n^k < 2^{31}

, and exit with an error in that case (as the Dynare++ codebase already does). If at some point we feel that this constraint is too restrictive, we could try to raise the limit using 64-bit integers, taking into account the above points.

Thank you @sebastien. I agree with the above. Additional thoughts:

• We should check how long it takes to compute the 3rd order approximation of a model with 1290 variables
• A model of this size would need to be about entirely linear for the non-zero 3rd order derivatives to be stored in the RAM of a regular computer.

@MichelJuillard We had this with a NK-HANK model of 900 equations. That one was not sparse and took about 45 minutes to preprocess and resulted in static and dynamic-files of about 45MB. The main issue was the preprocessor trying to optimize. So the bottleneck may not be the higher order stuff.

Interesting @JohannesPfeifer

closed with commit 5e8b478c

mentioned in commit 1683baee

mentioned in issue dynare#1859