EigenSolver question

[Alexios Galanos]

I am trying to compute the absolute eigenvalues of a square non-symmetric matrix using the Eigen::EigenSolver function:

template<class Type>
vector<Type> esolver(matrix<Type>A){
Eigen::EigenSolver<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>> eigs(A);
Eigen::MatrixXcd lambda(eigs.eigenvalues());
vector<Type> tmp = lambda.cwiseAbs();
return tmp;
};

Matrix A is a is formed from a combination of PARAMETERS and DATA, and relates to a constraint I am trying to impose for system stability.

I seem to have read a similar question on the GitHub issues page posted a long time ago, but haven’t seen any solutions to this (https://github.com/kaskr/adcomp/issues/144).

The problem I think appears to be casting from CppAD::AD<double> to double. I’ve managed to make this work using RTMB, but would prefer to not use it as I am already importing TMB in my package. I also took a look at the tmbutils/matexp.hpp code which also makes use of EigenSolver, but couldn’t figure out exactly how to adjust this for my purpose.

Would appreciate any suggestions.

Thanks,

Alexios

[Kasper Kristensen]

It's possible to make your example compile (you just need to use the right types...) *but I do not recommend it* :

Eigen::EigenSolver<Eigen::Matrix<Type, Eigen::Dynamic, Eigen::Dynamic>> eigs(A);
Eigen::Matrix<std::complex<Type>, Eigen::Dynamic, 1> lambda(eigs.eigenvalues());

vector<Type> tmp = lambda.cwiseAbs();

While this appears to work, here's why it actually doesn't:

1. Internally, the general eigen solver uses if/else statements. The AD tape is then valid for the A matrix used when the AD tape was generated. However, when you evaluate the same tape for a different matrix, you'll get wrong eigen values (unless the new matrix is sufficiently close to the initial matrix in some sense...).
2. According to the C++ standard, std::complex<Type> is 'undefined behaviour'. So, while std::complex<Type> might work as expected for most compilers, I wouldn't feel confident that it works on all of the CRAN flavors.

By contrast, in RTMB, the general eigen decompostion is an atomic function. It's therefore safe to use. In theory it would be straight forward to implement the same atomic in normal TMB if it wasn't for (2).

[Alexios Galanos]

Thanks Kaspar,

Appreciate the explanation. I'll use RTMB.

Cheers,

Alexios