Math Forum Discussions

Chip Eastham

Posts: 2,341
Registered: 12/13/04

Re: Computation of the matrix exponential
Posted: May 25, 2010 7:33 PM

Plain Text

On May 25, 4:49 am, luca <luca.frammol...@gmail.com> wrote:
> On 25 Mag, 07:59, Chip Eastham <hardm...@gmail.com> wrote:
>
>
>
> > On May 24, 8:37 pm, luca <luca.frammol...@gmail.com> wrote:
>
> > > Hi,
>
> > > i have the following problem: given a 3x3 real matrix, compute exp(A).
>
> > > I need a really fast way to do this. I have searched a bit with
> > > google, but it seems to me that
> > > computing the matrix exponential is not so simple, at least if your
> > > matrix does not have a special
> > > structure (for example A=diagonal matrix).
>
> > > I have found a simple method that use the diagonalization of A. If A
> > > has 3 distinct eigenvalues, than compute
> > > A=PDP^-1, where P is the matrix of the eigenvectors, D is a diagonal
> > > matrix (whose diagonal elements are
> > > the eigenvalues of A). Than, exp(A) = P exp(D) P^-1. Since P^-1 is
> > > fast enough and exp(D) is simple
> > > to compute, this should be a fast method.
>
> > > But, the problem is: i am not sure that the matrix A will always have
> > > 3 distinct eigenvalues...what happens
> > > if this does not happen? Can i use that formula even if 2 (or all
> > > three) eigenvalues are equal?
>
> > > Are there any other ways to compute exp(A) in a fast way?
>
> > > Thank you,
> > > Luca
>
> > You'll probably find the "updated" discussion of matrix
> > exponentiation by Cleve Moler and Charles Van Loan to be
> > helpful:
>
> > [Nineteen Dubious Ways to Compute the Exponential of a
> > Matrix, Twenty-Five Years Later]http://www.cs.cornell.edu/cv/ResearchPDF/19ways+.pdf
>
> > Your statement that you need to compute the matrix
> > exponential "really fast" suggests that you intend
> > to do this repeatedly. A bit more information about
> > the reason for doing this might lead to a more focused
> > evaluation of algorithms.
>
> > regards, chip- Nascondi testo citato
>
> > - Mostra testo citato -
>
> Hi Chip,
>
> thank you for your replay. Yes, i need to compute the matrix
> exponential during an iterative algorithm.
>
> The algorithm is a second-order minimization one.
>
> Suppose that you have a vector x (of size 8 x 1), computed during the
> iterations of the algorithm.
> When || x || <= eps, with eps a small real number, there is
> convergence.
>
> The method is used to find a homography (3x3 matrix). In practice, one
> starts with an estimate of the
> homography H', compute the vector x (i will not give other details on
> this since it would
> require too much space), update the homography using H' = H' *
> exp(A(x)).
>
> What is A? A is a 3x3 matrix of the form:
>
> A(x) = sum_i=1..8 (x_i*A_i)
>
> where the A_i are 3x3 matrices. They are chosen to be a basis for the
> Lie algebra sl(3). sl(3) is the algebra associated to the group SL(3).
> The projective transformation H is in SL(3).
> The papers says that :"the exponential map is a homeomorphism between
> a neighborhood of I in SL(3) and a neighborhood of the null matrix 0
> of sl(3)".
>
> The A_i are constant matrices. They have trace null. The matrix A_i
> will have:
>
> 1] only one non-diagonal element non null (1)
> or
> 2] two diagonal elements non-null (1 and -1)
>
> Since the methods is such that the vector x will tends to have smaller
> and smaller components, if we start with a x sufficiently near to the
> null vector,
> we can use the above parameterization of H.
>
> This is all.
> If you are interested i can give you the paper of reference. I am
> working in the field of Computer Vision.
>
> Luca

Thanks, Luca, for taking the time to explain the
setup to your question in some detail. For the
purpose of discussing computation of exp(A) it
is important to know that the majority of these
cases will involve a small ||x|| and hence a
correspondingly small A. [The fact that trace(A)
is zero means the only way an eigenvalue could
have multiplicity 3 is to be zero, which would
simplify the Taylor series to a low degree
polynomial as previously sketched.]

Since the algorithm amounts to a sequence of
updates to H', one can optimistically hope that
errors in the initial choice of H' and in the
computation of exp(A(x)) will self-correct.

Let's consider three possible matrix exponential
methods:

(1) Find the Jordan canonical form of A and the
matrix exponential of the corresponding blocks

(2) Pade (rational) approximation of the matrix
exponential (possibly with scaling & squaring
when A is not yet "small")

(3) Taylor (power) series (when A is "small")

Method (1) is probably the most expensive because
of the need to compute the Jordan canonical form
for each different A, but since A is 3x3 this is
probably not too expensive, at least for some
fixed number of initial steps when A is not yet
small. Method (2) corresponds to using rational
approximations, rather than polynomial ones, to
the exponential function, so has a certain expense
for computing at least one matrix inverse (as well
as some matrix-matrix multiplications). Again the
matrices are 3x3 so this isn't all that bad. In
method (3) we have the simplest operations (just
matrix-matrix multiplies and adds) but the Taylor
series may require several more terms than a
Pade (rational) approximation would for the same
accuracy.

I'll send you some links to details about these
alternatives.

regards, chip