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

[ Thread Index |
Date 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
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