Complex numbers and matrix inversion in Stan

Dear all.

I’m trying to figure our whether it is possible to use a model that includes complex matrices, and manipulations on such matrices in STAN. Specifically, I’m trying to fit a semi-markov model, which can be solved using Laplace transform which requires complex matrix manipulations.

Thanks!
Isaac

For the next Stan Math release, there is a lot of support for autodiffing complex numbers, but it is not exposed in the Stan language and there isn’t a concrete plan for doing so anytime soon. So, you would have to declare a function in the functions block of your Stan program and write your own C++ implementation, perhaps returning the real and imaginary parts separately.

Thanks for the prompt response. So can I just write my own C++ code in stan (using data types not included in stan)? What about using other libraries such as armadillo or eigen?

You can write your own C++ code that is called from Stan, in which case you can use any libraries you want but they will have to be templated to accept var types (which excludes Armadillo).

Ok. Is there any chance you can refer me to an example or two of something similar (ie standard cpp code within stan)? That would help me a lot.

Thanks!

Isaac

https://cran.r-project.org/web/packages/rstan/vignettes/external.html

That’s cool! Thanks. One last question if that’s ok - it seems like Stan is already using the Eigen library and types, is that correct? I have the function I need written using Eigen types already - can I just add this as a function? Or do I need to use more basic C++ types?

Isaac.

You can use Eigen, but you are going to have to do stuff like

Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>

instead of

Eigen::MatrixXd

Thanks Ben.

I’ve tried to implement it as you’ve suggested, but I get an error suggesting that the parser cannot understand this Eigen structure:

stanCode='
functions {
  Eigen::Matrix<std::complex, 4, 4>(Eigen::Matrix<std::complex, 4, 4>);
}
model {}

m1<-stan_model(model_code = stanCode)
'

I get the following error:


PARSER EXPECTED: "}"
Error in stanc(file = file, model_code = model_code, model_name = model_name,  : 
  failed to parse Stan model '65395727a452897c9f7949268dd2de08' due to the above error.
>

It is telling you the truth, in part because the complex number stuff is not released yet, in part because it is not planned to be released for the Stan language, and in part because that isn’t how you connect with external C++ code. The functions block of your Stan program should still look like it usually does, except you would need to pass the real and imaginary parts separately

functions {
  matrix[] foo(matrix real, matrix imag); // no definition
}
...

Then you need to write a C++ implementation that returns a std::array<Eigen::Matrix<T, -1, -1> > where the first element of the array is a matrix of real parts and the last element is a matrix of imaginary parts.

Oh. I see. So that means I can’t use std::complex or Eigen’s complex matrix inversion and multiplication schemes at all? Any chance I can use these types and functions within my function if I make it a separate hpp file? And then only return a matrix array as suggested?

You can do all that inside your C++ function. But the input needs to be real and the output needs to be real, which is why you need to separate the real and imaginary parts.

Ok. I see. That’s encouraging!

Thanks a lot!

Isaac.

Dear Ben.

I hope you’ll have the time for one more brief follow up question: I’ve implemented my function such that only what is inside the function will include Eigen types, whereas the function output and arguments are Stan matrices, vectors and real numbers. The only question now is whether there is anything special I need to do to convert a Stan matrix to an Eigen::matrix or an Eigen::array?

Thanks!

Isaac.

A matrix in the Stan language is an Eigen matrix, but there is nothing in the Stan language that is an Eigen array. Arrays in the Stan language are std::vectors.

I see. So if for example I have a function like the below it should work smoothly, right?

matrix foo(matrix a, real s){

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> b;

b=a*s;

return(b.inverse());

}

No. It would need to be something like this in Stan

functions {
  matrix[] foo(matrix[] a, real s);
}
...

and this in C++

template <typename T0__, typename T1__>
std::vector<Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, Eigen::Dynamic> >
foo(const std::vector<Eigen::Matrix<T0__, Eigen::Dynamic, Eigen::Dynamic> >& a,
      const T1__& s, std::ostream* pstream__) {
  Eigen::Matrix<std::complex<var>, -1, -1> b(a[0].rows(), a[0].cols();
  for (int i = 0; i < b.rows(); i++) for (int j = 0; j < b.cols(); j++) {
    b.coeffRef(i, j) = std::complex<var>(a[0].coeff(i, j), a[1].coeff(i,j)) * s;
  }
  auto b_inv = b.inverse().eval(); // Eigen::Matrix with the same type as b
  return {b_inv.real(), b_inv.imag()};
}

However, Stan’s complex var type is not being introduced until 2.23.

Ben,
Thanks for devoting the time to help me with this. It does seem more complicated than what I thought initially.
First - just to verify, that fact that “Stan’s complex var type is not being introduced until 2.23” means that I cannot use std::complex for now? (i.e., that the function you wrote above won’t work?).

I also had a few questions about the code, just to help me figure our how Stan interfaces with Cpp functions:
A) If I want to copy a Stan matrix to a regular Eigen::Matrix (i.e. var/double type), can I just use (or do I need to copy it element-by-element?):

Eigen::Matrix<var,2,2>=a;

B)
What is the meaning of defining the matrix argument with?:

<typename boost::math::tools::promote_args<T0__, T1__>::type>

Instead of just:

<double>

or

<var>

Again - thanks a lot for the help! I hope at the end of this I will be able to contribute a function for fitting semi-markov models in Stan…
Isaac.

Ok, so I’ve tried to just write my function (which is a bit more complex that what I wrote before, however it should only return a simple matrix, and it’s arguments are simple matrices and a real complex number). Here is what I wrote in Stan:

 matrix foo(real sReal, real sImag, matrix P, matrix lambda, real Ua, real Ub);

and here is what I wrote in a separate cpp file:

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> foo(const double sReal,const double sImag,
                                                           const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> P,
                                                           const Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> lambda,
                                                           const double Ua,
                                                           const double Ub,
                                                           std::ostream* pstream__) {
  std::complex<double> s;
  s.real(sReal);
  s.imag(sImag);
  
  Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp1 =P.array()*(lambda.array()/(lambda.array()+s));
  Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> Id=tmp1;
  Id.setIdentity();
  Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp2=(Id-tmp1).inverse(); 
  Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> tmp3=tmp2.diagonal().asDiagonal().inverse();
  std::complex<double> LTuni = (std::exp(-s*Ua)-std::exp(-s*Ub))/(s*(Ub-Ua)); 
  return((tmp1*tmp2*tmp3*LTuni).real());
}

I then wrote a stupid model with a transformed parameters block that just gets the output of foo for a given set of arguments, and it seems to work - i.e., the results match the real results of this function. Does that mean that I am good to go with implementing the other parts of the function (which do not include complex matrix inversion or complex numbers) in Stan and the model itself? Or am I missing something? I am not sure how exactly does this work with the fact that std::complex should not yet work in Stan…

One issue that I’ve noticed is that I can call this CPP function from the model/transformed parameter/generated quantities block, but if I’m trying to call it from another Stan function - it doesn’t work. For example, I created another Stan function that simply calls foo:

matrix test(real sReal, real sImag, matrix P, matrix lambda, real Ua, real Ub){
    matrix[cols(P),cols(P)] a;
    a=foo(sReal,sImag,P,lambda,Ua,Ub);
    return(a);

  }

I get the following error message:

file125c5eca5307.o:file125c5eca5307.cpp:(.text$_ZN60model125c42415a5b_c54d1236bbbfb99ce518497ac68d1123_namespace4testIddddddEEN5Eigen6MatrixIN5boost4math5tools12promote_argsIT_T0_T1_T2_NS6_IT3_T4_ffffE4typeEfE4typeELin1ELin1ELi0ELin1ELin1EEERKS7_RKS8_RKNS2_IS9_Lin1ELin1ELi0ELin1ELin1EEERKNS2_ISA_Lin1ELin1ELi0ELin1ELin1EEERKSB_RKSC_PSo[_ZN60model125c42415a5b_c54d1236bbbfb99ce518497ac68d1123_namespace4testIddddddEEN5Eigen6MatrixIN5boost4math5tools12promote_argsIT_T0_T1_T2_NS6_IT3_T4_ffffE4typeEfE4typeELin1ELin1ELi0ELin1ELin1EEERKS7_RKS8_RKNS2_IS9_Lin1ELin1ELi0ELin1ELin1EEERKNS2_ISA_Lin1ELin1ELi0ELin1ELin1EEERKSB_RKSC_PSo]+0x2d6): undefined reference to `Eigen::Matrix<boost::math::tools::promote_args<double, double, double, doub

What am I missing here?
BW,
Isaac.

Dear Ben and other users.
I am assuming you are quite busy, but I wanted to try and ask these last questions again. Specifically - why do I get an error message when calling my external C++ function (with Eigen::Matrix input and output) from a Stan function, but not from the transformed parameters/model/generated quantities block.
And second - assuming that we will solve the first issue, which will allow me to use complex numbers (double type) in the likelihood function (that will return a real scalar, however), or will it be a problem because the var type still doesn’t allow for complex numbers.

BW,
Isaac