Matrix_exp for complex numbers

General
1

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.

2

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
3

Thank you for the update! I’m looking forward to the next release.

4

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?

5

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 ? ).

6

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)