10
$\begingroup$

I recently asked a question along the same lines for skew-Hermitian matrices. Inspired by the success of that question, and after banging my head against a wall for a couple of hours, I'm looking at the matrix exponential of real asymmetric matrices. The route to finding the eigenvalues and eigenvectors seems rather convoluted, and I'm afraid I've gotten lost.

Background: Some time ago I asked this question on the theoretical physics SE. The result allows me to phrase master equations as real asymmetric matrices. In the time-independent case, the master equation is solved by exponentiating this matrix. In the time-dependent case it will require integration. I'm only concerned with time-independence at the moment.

Upon looking at the various subroutines I think I should be calling (?gehrd, ?orghr, ?hseqr...) it is unclear if it would be simpler to cast the matrix from real*8 to complex*16 and proceed with the complex double versions of these routines, or stick with real*8 and take the hit of doubling the number of my arrays and making a complex matrix of them later.

So, which routines should I be calling (and in what order), and should I use the real double versions or the complex double versions? Below is an attempt at doing this with real double versions. I've become stuck finding the eigenvalues and eigenvectors of L*t.

function time_indep_master(s,L,t)
  ! s is the length of a side of L, which is square.
  ! L is a real*8, asymmetric square matrix.
  ! t is a real*8 value corresponding to time.
  ! This function (will) compute expm(L*t).

  integer, intent(in)    :: s
  real*8,  intent(in)    :: L(s,s), t
  real*8                 :: tau(s-1), work(s), wr(s), wi(s), vl
  real*8, dimension(s,s) :: time_indep_master, A, H, vr
  integer                :: info, m, ifaill(2*s), ifailr(2*s)
  logical                :: sel(s)

  A = L*t
  sel = .true.

  call dgehrd(s,1,s,A,s,tau,work,s,info)
  H = A
  call dorghr(s,1,s,A,s,tau,work,s,info)
  call dhseqr('e','v',s,1,s,H,s,wr,wi,A,s,work,s,info)
  call dhsein('r','q','n',sel,H,s,wr,wi,vl,1,vr,s,2*s,m,work,ifaill,ifailr,info)

  ! Confused now...

end function
Improve this question
$\endgroup$

3 Answers 3

Reset to default
8
$\begingroup$

I would first think really hard about whether or not the matrix is really completely arbitrary: Is there any transformation that would make it Hermitian? Does the physics guarantee that the matrix should be diagonalizable (with a reasonably conditioned eigenvector matrix)?

If it turns out that there really isn't any symmetry to exploit, then you should start by reading Nineteen Dubious Ways to Compute the Matrix Exponential, which is the standard reference (and is written by the author of MATLAB and the coauthor of G&vL).

Improve this answer
$\endgroup$
7
  • 1
    $\begingroup$ The best I can do is turn it into a block diagonal matrix with asymmetric blocks. That in itself is very interesting though. Most of those blocks are $2\times2$, and I can just solve those analytically. There is a remaining $4\times4$ block with no symmetries to exploit though. $\endgroup$
    – qubyte
    Feb 14, 2012 at 7:45
  • 1
    $\begingroup$ I like this answer; the unsymmetric case possesses enough pitfalls that it is worth considering if there might be a formulation of your problem that leads to symmetric matrices instead of unsymmetric ones. $\endgroup$
    – J. M.
    Feb 14, 2012 at 11:11
  • $\begingroup$ @MarkS.Everitt: You seem to be almost there...how big are the matrices? ~36 x 36 again? $\endgroup$
    – Jack Poulson
    Feb 14, 2012 at 13:34
  • $\begingroup$ In this case $16\times16$, but there is the possibility of it going up to $36\times36$. $\endgroup$
    – qubyte
    Feb 14, 2012 at 13:45
  • 2
    $\begingroup$ @MarkS.Everitt: So your problem is effectively now just how to exponentiate 4x4 matrices. This is sufficiently small for asymptotic analysis to be irrelevant, so the answer will depend entirely on the values. I can't really say anymore unless you translate your linked physics post into linear algebra (what is a superoperator?!?). $\endgroup$
    – Jack Poulson
    Feb 14, 2012 at 17:47
7
$\begingroup$

To build on what Jack has said, the standard approach that seems to be used in software (like EXPOKIT, mentioned in your earlier question) is scaling-and-squaring followed by Padé approximation (Methods 2 and 3) or Krylov subspace methods (Method 20). In particular, if you're looking at exponential integrators, you'll want to consider the Krylov subspace methods, and look at papers on exponential integrators (some references are mentioned along with Method 20 in the Moler & van Loan paper).

If you're hell bent on using eigenvectors, consider using triangular systems of eigenvectors (Method 15); since your matrix may be nondiagonalizable, this approach might not be best, but it's better than trying to calculate the eigenvectors and eigenvalues directly (i.e., Method 14).

Reduction to Hessenberg form is not a good idea (Method 13).

It's not apparent to me whether you would be better served with real or complex arithmetic, since Fortran complex arithmetic is fast, but might overflow/underflow (see "How much better are Fortran compilers really?").

You can safely ignore Methods 5-7 (ODE solver-based methods are inefficient), Methods 8-13 (expensive), Method 14 (calculating eigenvectors of large matrices is hard without special structure and prone to numerical error in ill-conditioned cases), and Method 16 (calculating the Jordan decomposition of a matrix is numerically unstable). Methods 17-19 are trickier to implement; in particular, Methods 17 and 18 would require more reading. Method 1 is a fall-back option for scaling-and-squaring if Padé approximations don't work well.

Edit #1: Based on the comments in response to Jack's answer, block diagonalizing seems like an option, in which case, something like Method 18 (block-triangular diagonalization) is a very good method to use. I hesitated to recommend it at the beginning because your question didn't specify this structure, but if you have a transformation that block diagonalizes your matrix, it takes most of the complexity out of the approach. You'll just want to make sure that you use G. W. Stewart's trick of decomposing each block diagonal matrix $B_{j}$ into

$$B_j = \gamma_j I + E_j,$$

where $\gamma_{j}$ is the average of the eigenvalues of the $j$th block diagonal matrix. This decomposition will make $E_{j}$ nearly nilpotent, which will enhance the accuracy of your matrix exponential calculations. This trick is discussed on page 26 of the version of the Moler & van Loan paper that Jack linked to.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ Upvoted, but the folks that implemented LAPACK are not naive about complex arithmetic, especially given all of the time Kahan put into analyzing it. However, he should still compute in real arithmetic if it saves flops; the cost of the conversion later is just $O(n^2)$ versus all of the $O(n^3)$ dense linear algebra. $\endgroup$
    – Jack Poulson
    Feb 14, 2012 at 4:54
  • $\begingroup$ No doubt they know what they're doing; I'm not worried about LAPACK's implementation. I'm more surprised about the Fortran compiler behavior. $\endgroup$
    – Geoff Oxberry
    Feb 14, 2012 at 5:05
  • 2
    $\begingroup$ Yeah, the compiler might be more of a problem than well-written LAPACK. It can be disconcerting to find that your program was failing only because the implementations for absolute value and division used by the compiler were botched... $\endgroup$
    – J. M.
    Feb 14, 2012 at 11:14
-1
$\begingroup$

I have a simple Fortran subroutine that computes the exponent of an arbitrary matrix. I checked it up against Matlab command and it is fine. It is based on scaling and squaring. I wrote it a few years back.

I would have like to fins another subroutine, like the ones I download from gams.nist.gov. But no luck yet.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.