Search All of the Math Forum:
Views expressed in these public forums are not endorsed by
NCTM or The Math Forum.



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


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, TwentyFive 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 secondorder 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 nondiagonal element non null (1) > or > 2] two diagonal elements nonnull (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 selfcorrect.
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 matrixmatrix multiplications). Again the matrices are 3x3 so this isn't all that bad. In method (3) we have the simplest operations (just matrixmatrix 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



