[eigen] Matrix exponential for evaluating a solution to a linear system of differential equations

[Douglas Bates <bates@xxxxxxxxxxxxx>]

As described in "Nineteen Dubious Ways of Computing the Exponential of
a Matrix" and in the "...: Twenty-five Years Later" follow-up paper
(SIAM Review, vol. 45, no. 1, 2003) there are many different
approaches to computing the matrix exponential.

I would appreciate recommendations on Eigen-based approaches to
determining the solution to linear systems of ordinary differential
equations using the matrix exponential.  In the cases I will consider,
the system matrix A is reasonably small, usually smaller than 5 by 5
and parameterized.  (They are from compartment models in
pharmacokinetics.)  For given values of the parameters we need to
evaluate e^{At} x_0 for several different values of t.  We also may
want to evaluate the sensitivities of the solution with respect to the
parameters.  The nature of the systems is such that there are no
complex eigenvalues for A.

I realize there is a matrix exponential function in the unsupported
MatrixFunctions but I would prefer to use a decomposition method if
possible.

So, assuming that A is diagonalizable, we can use the
Eigen::EigenSolver class and proceed happily.  However, we need to
detect if A is close to being non-diagonalizable.

My current approach is


        static double
precision(std::sqrt(std::numeric_limits<double>::epsilon()));
        Eigen::EigenSolver<Eigen::MatrixXd> es(A);
        if (es.eigenvalues().imag().squaredNorm() != 0)
	    throw std::invalid_argument("Complex eigenvalues");
        Eigen::VectorXd   vals(es.eigenvalues().real());
        if (es.eigenvectors().imag().squaredNorm() != 0)
	    throw std::invalid_argument("Complex eigenvectors");
	Eigen::MatrixXd   vecs(es.eigenvectors().real());
	Eigen::VectorXd  alpha(Eigen::FullPivLU<Eigen::MatrixXd>(vecs).solve(x0));
        bool nonsingular = (vecs * alpha).isApprox(x0, precision);

but I am confident that can be improved.  For one thing, even when A
is non-diagonalizable the matrix of eigenvectors is regarded as
nonsingular according to this test and that makes me uncomfortable.  I
could calculate the condition number from an Eigen::JacobiSVD but that
seems to be a very long way around.

Would it be better to check for a non-diagonalizable A using the T
matrix from the Schur decomposition?  Is there a simple way of getting
the condition number of a triangular matrix or an estimate of it in
Eigen (similar to dtrcon in Lapack)?