|Re: [eigen] Matrix exponential for evaluating a solution to a linear system of differential equations|
[ Thread Index |
| More lists.tuxfamily.org/eigen 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
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
if (es.eigenvalues().imag().squaredNorm() != 0)
throw std::invalid_argument("Complex eigenvalues");
if (es.eigenvectors().imag().squaredNorm() != 0)
throw std::invalid_argument("Complex eigenvectors");
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.