Solve a Lyapunov / Sylvester equation: include custom c++ function using Eigen library -- possible?

Modeling
#1

I’d like to solve a Lyapanov equation where AX + XA’= Q within my stan model.
I have a c++ function using Eigen that can do this, and I have tried to include it in the stan model, but I just get the seemingly standard " Compilation ERROR, function(s)/method(s) not created!" message. My c++ is not great, I know the code works within R / Rcppeigen, but I’m also not at all sure whether Stan will allow this sort of more elaborate / library dependent c++ , and if so, which if any includes are needed.

library(rstan)
mc='
functions{
  matrix syl(matrix A, matrix Q);
}
data{
int d;
matrix[d,d] S;
matrix[d,d] A;
}
transformed parameters{
matrix[d,d] C = syl(A,S);
}
'

s=stan_model(model_code = mc, model_name = "external", allow_undefined = TRUE,
  includes = paste0('\n#include "', 
    file.path(getwd(), 'syl.hpp'), '"\n'))

The c++ is:

Eigen::MatrixXd syl(const Eigen::MatrixXd & A, const Eigen::MatrixXd & Q, std::ostream* pstream){

const MatrixXd B = A.transpose();

Eigen::ComplexSchur<Eigen::MatrixXd> SchurA(A);
Eigen::MatrixXcd R = SchurA.matrixT();
Eigen::MatrixXcd U = SchurA.matrixU();

Eigen::ComplexSchur<Eigen::MatrixXd> SchurB(B);
Eigen::MatrixXcd S = SchurB.matrixT();
Eigen::MatrixXcd V = SchurB.matrixU();

Eigen::MatrixXcd F = (U.adjoint() * Q) * V;

Eigen::MatrixXcd Y =
Eigen::internal::matrix_function_solve_triangular_sylvester(R, S, F);

Eigen::MatrixXcd X = (U * Y) * V.adjoint();

Eigen::MatrixXcd Q_calc = A * X + X * B;

return Q_calc.real();
}
Schur decomposition ready?
#2

Hi, I’m sure you have read this page

https://mc-stan.org/rstan/articles/external.html

Remember to add ostream for the function arguments.

This could also help

https://pystan.readthedocs.io/en/latest/external_cpp.html

3 Likes
#3

Thanks, there seems to be some discrepancy as to whether the argument is pstream or pstream__ though? Edit: it doesn’t matter.

#4

Ok, the non created function seems to be due to Eigen::internal::matrix_function_solve_triangular_sylvester
Does the fact that it’s an unsupported function in Eigen mean that it can’t be found? Including:
#include <unsupported/Eigen/MatrixFunctions>
didn’t help.

#5

I don’t know if you need to add some flags so compiler can see that function?

Maybe some -isystem/-I?

#6

edited:
I couldn’t work it out so just copied / adapted that function. I’m now stuck on the autograd, not sure how to pass the vars through… code below works but all the commented out lines are needed and things break.
Any hand holding would be appreciated!

template <typename T1>
Eigen::Matrix<T1, Eigen::Dynamic, Eigen::Dynamic> sylcpp(const Eigen::Matrix<T1, Eigen::Dynamic, Eigen::Dynamic> & A,
  const Eigen::Matrix<T1, Eigen::Dynamic, Eigen::Dynamic>  & C){
  
  using matrix_t = Eigen::Matrix<T1, Eigen::Dynamic, Eigen::Dynamic>;
  matrix_t B = A.transpose();

  // Eigen::ComplexSchur<matrix_t> SchurA(A);
  // Eigen::MatrixXcd R = SchurA.matrixT();
  // Eigen::MatrixXcd U = SchurA.matrixU();
  // Eigen::ComplexSchur<Eigen::MatrixXd> SchurB(B);
  // Eigen::MatrixXcd S = SchurB.matrixT();
  // Eigen::MatrixXcd V = SchurB.matrixU();
  // Eigen::MatrixXcd F = (U.adjoint() * C) * V;
  // Eigen::MatrixXcd Y = matrix_function_solve_triangular_sylvester(R, S, F);
  // Eigen::MatrixXcd X = (U * Y) * V.adjoint();
  // return X.real();
  return(A*C+C*A);
}


template <typename T0__, typename T1__>
Eigen::Matrix<typename boost::math::tools::promote_args<T0__, T1__>::type, Eigen::Dynamic, Eigen::Dynamic>
syl(const Eigen::Matrix<T0__, Eigen::Dynamic, Eigen::Dynamic>& A,
         const Eigen::Matrix<T1__, Eigen::Dynamic, Eigen::Dynamic>& C, std::ostream* pstream__){

  Eigen::Matrix<T0__, Eigen::Dynamic,Eigen::Dynamic> X = sylcpp(A,C);
  return(X);
}
#7

I think the blogpost has a good example.

I can not test these right now, but I can try to help later this week.

Maybe our C++ folks can help cc @wds15 @syclik ?

1 Like
#8

Thanks. The blogpost code didn’t work for me, maybe I misunderstood something, but the code above is working fine up to simple functions like in the blog, but I can’t manage to pass the var through all the eigen functions. Possibly a complex number issue in there as well and I will have to find a way to do it with reals.

#9

Have you seen

https://cran.r-project.org/web/packages/StanHeaders/vignettes/stanmath.html

1 Like