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

[ Thread Index | Date Index | More Archives ]

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

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

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


Mail converted by MHonArc 2.6.19+