[unitaryhack] Improve the calculation of multidimensional hermite polynomials

[nquesada]

This issue has been tagged for a bounty during unitaryHACK

Currently all calculations of multidimensional hermite polynomials are carried in double precision even though the C++ functions to do these calculations are templated to allow for any datatype.

A useful addition to The Walrus would be to allow the Python modules to have the C++ functions do calculations in long double precision. Note that even if the input and output of the C++ functions is ultimately casted to double it is still useful to internally do the calculations in long double given that numerical errors can accumulate very fast. Similar issues for the calculations of hafnians have been solved by doing its internal calculation in long double.

Two alternatives come to mind to achieve this:

1. Make a wrapper function that takes as input the matrix R and vector y as double precision arrays and then converts them to long double and pass them to the templated version of hermite_multidimensional_cpp. An array in long double would be returned and then the wrapper function would cast this array into double precision and return it to Python using the ArrayWrapper for double datatypes in libwalrus.pyx.

2. A second option that would avoid any type conversion is to write new array wrappers that allow to pass transparently long double arrays into numpy np.longdouble.

An important caveat for testing is that in MS Visual C++ long double is an alias for double so no gain in precision can be obtained by using this compiler.

[nquesada]

A second improvement would be to find a way to parallelize the loop in line

thewalrus/include/hermite_multidimensional.hpp

Line 114 in bc53689

for (ullint jj = 0; jj < Hdim-1; jj++) {

.

This might require a nontrivial amount of algorithmic thinking since this loop implements a recursion relation.

[nquesada]

Finally, it would also be interesting to think about porting this part of the code to numba.

[TripleRD]

Hi @nquesada, I didn't get the second improvement you proposed. Can you elaborate on it? Thanks!

[nquesada]

I just meant, to write an implementation of Hermite multidimensional in numba.

[TripleRD]

A second improvement would be to find a way to parallelize the loop in line

thewalrus/include/hermite_multidimensional.hpp

Line 114 in bc53689

for (ullint jj = 0; jj < Hdim-1; jj++) {

.
This might require a nontrivial amount of algorithmic thinking since this loop implements a recursion relation.

I was referring to this comment, can you elaborate on it? Thanks!

[co9olguy]

Hi @TripleR47, the idea here is to use numba to replace that line (a serial for loop) with a parallelized version

[nquesada]

To add to @co9olguy : the challenge is that the for loop in line 114 is very long (dimensions Hdim) and implements a recursion relation, where the next element is calculated based on the values of previous ones. My comment is just to highlight that breaking a parallelzing a recursive loop might be non-trivial!

[e-eight]

Would love to take a stab at this. Okay if I work on this?

[co9olguy]

Hi @e-eight, yes this issue is still open. Feel free to tackle it if interested :)

[e-eight]

Alright. Will give it a shot.