Finding gradient of ComplexEigenSolver calculation




I have written some C++ code to evaluate the likelihood of some time-series data in the frequency domain.

One of the steps in the likelihood calculation uses functions in Eigen’s ComplexEigenSolver.

I would like to find the gradient (and the Hessian) of my likelihood function with respect to model parameters (theta). Is this possible using the Stan Math library? If not, do you have any advice on what would be the best way of going about this?

Here is a snippet of my code,

double like(const vector<double> & y_f, const vector<double> & theta)
  MatrixXcd J;
  // evaluate Jacobian
  jac(theta, J);

  // evaluate eigen-decomposition of Jacobian
  ComplexEigenSolver<MatrixXcd> ces(J);
  VectorXcd EigVals = ces.eigenvalues();
  MatrixXcd R = ces.eigenvectors();
  MatrixXcd L = R.inverse();
  // more calculations using R, L, EigVals to calculate ll ...
  return ll;


The Stan Math Library does not know how to do autodiff for complex numbers. You have to do the derivatives analytically using only double values.