Add a method to SelfAdjointEigenSolver for computing the matrix exponential

This commit is contained in:
Rasmus Munk Larsen
2025-10-05 15:06:04 +00:00
parent 32b0f386bc
commit 7eaf9ae68d
2 changed files with 27 additions and 0 deletions

View File

@@ -13,6 +13,7 @@
#include <limits>
#include <Eigen/Eigenvalues>
#include <Eigen/SparseCore>
#include <unsupported/Eigen/MatrixFunctions>
template <typename MatrixType>
void selfadjointeigensolver_essential_check(const MatrixType& m) {
@@ -135,11 +136,13 @@ void selfadjointeigensolver(const MatrixType& m) {
VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorExp());
eiSymmUninitialized.compute(symmA, false);
VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt());
VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorExp());
// test Tridiagonalization's methods
Tridiagonalization<MatrixType> tridiag(symmC);
@@ -167,6 +170,14 @@ void selfadjointeigensolver(const MatrixType& m) {
eiSymmTridiag.eigenvectors().real().transpose());
}
// Test matrix expponential from eigendecomposition.
// First scale to avoid overflow.
symmB = symmB / symmB.norm();
eiSymm.compute(symmB);
MatrixType expSymmB = eiSymm.operatorExp();
symmB = symmB.template selfadjointView<Lower>();
VERIFY_IS_APPROX(expSymmB, symmB.exp());
if (rows > 1 && rows < 20) {
// Test matrix with NaN
symmC(0, 0) = std::numeric_limits<typename MatrixType::RealScalar>::quiet_NaN();