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

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

```

• References:

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/