Matrix_exp for complex numbers

Hi all,

I am excited to see that Stan is beginning to provide support for complex numbers. I am interested in using Stan for a quantum model, but this requires additional features. I was wondering whether Stan will support `matrix_exp` and user defined functions with complex number support in the future. If so, is there an approximate timeline? Thank you.

Complex vectors/matrices are a pretty high priority item for the language team. I think the hope is to have them in the next version.

User defined functions can already use complex numbers as arguments and return types

3 Likes

Thank you for the update! Iâ€™m looking forward to the next release.

I would not be so sure about supporting the matrix exponentialâ€¦itâ€™s certainly a great thing. So maybe search for complex number support for different operations on GitHub. Or is the implementation for complex numbers so general such that the matrix_exp function as currently implemented will â€śjust workâ€ť? @Bob_Carpenter ? I would think that it needs adjustments to be able to deal with complex numbers, right?

It appears that Stan math libraryâ€™s `matrix_exp` function is based on a modified copy of eigenâ€™s matrix exponential: math/MatrixExponential.h at master Â· stan-dev/math Â· GitHub

The underlying Eigen matrix exponential works for complex numbers:

tiny demo with debian stable `clang-11` & `libeigen3-dev`

``````#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>

#include <complex>
#include <iostream>

using namespace std::literals;

int main()
{
Eigen::MatrixXcd a(2,2);
a << 1.0 + 2.0i, 2.0 + 1.0i, 3.0 - 1.0i, 4.0 - 2.0i;

std::cout << "The matrix A is:\n" << a << "\n\n";
std::cout << "The matrix exponential of A is:\n" << a.exp() << "\n\n";
}
``````
``````\$ clang++-11 -I /usr/include/eigen3/ main.cpp -o main

\$ ./main
The matrix A is:
(1,2)  (2,1)
(3,-1) (4,-2)

The matrix exponential of A is:
(39.1276,7.06055) (60.9904,-9.96502)
(51.0254,-70.9555) (41.6007,-131.104)
``````

Still a decision for Stan developers to decide how easy it would be to maintain â€“ e.g. if Stan adds support for complex matrix exponential, this might be relatively simple to do for code that executes with an eigen backend but may be hard to do for other backends (e.g. would it need to be hand implemented for CUDA ? ).

My understanding is that as long as we have one working definition in /prim, the other folders are all optional (theyâ€™ll be used if they exist, but if they donâ€™t it falls back to the `prim` implementation)