Add routine that computes the action of a matrix exponential on a vector.

Problem

In our application, we do not necessarily need the matrix exponential exp(A) but rather its action c = exp(A) * b on a vector. This is a quite common problem and this paper [1] presents an algorithm to compute the solution. Are there any efforts to implement this algorithm within Eigen? If not, would it be interesting to implement?

Solution approach

There are implementations that may serve as reference: [2] and [3] On this basis, we have implemented a solution that solves our concrete problem. However, it would be nice to offer a generalized solution, which takes into account different matrix classes etc. Since I am not familiar with the Eigen development process, could anyone outline the steps how to get code into Eigen?

IMHO this algorithm could be added to the unsupported/Eigen/MatrixFunctions module.

References

  1. Al-Mohy and Higham, SIAM J. Sci. Comp., 33(2), 2011, https://doi.org/10.1137/100788860
  2. https://github.com/higham/expmv
  3. https://github.com/stan-dev/math/blob/develop/stan/math/prim/fun/matrix_exp_action_handler.hpp
Edited by Michael Riesch