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

[Jitse Niesen <jitse@xxxxxxxxxxxxxxxxx>]

On Thu, 16 Feb 2012, Douglas Bates wrote:

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

Just out of interest, why?

The matrixFunction() routine in the unsupported modules uses a Schur 
decomposition. For how to use it to compute the exponential, see the 
example at 
http://eigen.tuxfamily.org/dox-devel/unsupported/group__MatrixFunctions__Module.html#matrixbase_matrixfunction

> 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.

Exactly. If A is not close to being non-diagonalizable, this approach 
works fine, and is in some cases even more accurate than the routine in 
the MatrixFunction module.

> 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.

What is x0? A random vector? It does surprise me that this test does not 
work. If A is non-diagonalizable, then the matrix of eigenvectors as 
computed by Eigen is singular so the system vecs * alpha = x0 does not 
have a solution for generic x0. As far as I can see, finite precision does 
not really change that (your test with sqrt(eps) takes that into account). 
So the only thing I can think of is that you're unlucky with the x0. 
Perhaps it would help to give a complete example with the matrix A?

There is another thing. How close A is to be non-diagonalizable depends 
not only on the condition number of the matrix of eigenvectors but also on 
the distance between the eigenvalues. As an extreme example, the 2-by-2 
identity matrix is perfectly conditioned but is very close to the 
non-diagonalizable matrix [ 1, 0.00001; 0, 1 ].

> Would it be better to check for a non-diagonalizable A using the T
> matrix from the Schur decomposition?

There is a general rule that whatever you do with the eigendecomposition 
can be done better with the Schur decomposition, so perhaps yes. But I 
don't see immediately how you can use the Schur decomposition here.

> 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)?

A simple way is to compute the condition number as the product of the 
norms of the matrix and its inverse. If you use the infinity-norm, that is 
perhaps not too costly, especially since your matrices are fairly small.

We do not have something like dtrcon as far as I know. It is actually a 
fairly important bit that is missing.


Jitse