12

While trying to compute the Matrix Exponential of an $n \times n$ array I decided to take advantage of a Python function called scipy.linalg.expm().

According to the documentation, this function adopts the Padé's Approximant to perform a calculation that, otherwise, is said to require a lot of resources and is still a matter of discussion in the mathematics community (many of you might be aware of this).

The documentation itself references a work by Al-Mohy and Higham (doi: 10.1137/04061101X), in which I was not able to find any indication regarding the Time Complexity (either T or Theta would work) of this procedure.

Has anyone come across the same problem? What is your suggestion for determining said Time Complexity?

Rodrigo de Azevedo
  • 1
FaCoffee
  • 241
  • 5
    What's the time complexity of the scalar exponential? – user541686 May 17 '16 at 19:59
  • 1
    @Mehrdad -- $O(M(n) \log n)$, where $M$ is the cost of multiplication and $n$ is the number of bits of accuracy. See http://mathoverflow.net/questions/19946/ – Steve Huntsman May 18 '16 at 11:39
  • @SteveHuntsman: Is this for integer exponentials or does it work for anything? The reason I asked was that I could swear I'd come across a statement that said it is unknown how many bits of precision you'd need to compute an exponential to $n$ bits of accuracy, but I can't find it right now... – user541686 May 18 '16 at 17:28
  • @Mehrdad - It's for anything. A PDF of the relevant paper is linked to in my answer to the question linked above. – Steve Huntsman May 18 '16 at 18:08
  • 1
    @SteveHuntsman: Ah, what I was having trouble reconciling it with was the Table-Maker's Dilemma... could you explain why that isn't a problem here? As you can see, the article says, "Nobody knows how much it would cost to compute $y^w$ correctly rounded for every two floating-point arguments at which it does not over/underflow."... – user541686 May 18 '16 at 19:11
  • @Mehrdad- You can have lots of bits that are wrong while still having high accuracy (and precision). Consider the approximation $0.999999 \approx 1.000000$ in decimal by analogy. – Steve Huntsman May 18 '16 at 20:27
  • @Rodrigo editing a dozen old questions at a time drives new question off the front page. Please limit your retagging efforts to three or four questions a day. – Gerry Myerson Apr 05 '20 at 12:32
  • 1
    @GerryMyerson My apologies. I had a little time to kill and got carried away. – Rodrigo de Azevedo Apr 05 '20 at 14:35

2 Answers2

13

The exponential of a matrix is not, in general, computable explicitly. Strictly speaking, the Complexity of Matrix Exponential does not exist. Instead, you have an approximate calculation and you are asking about its complexity. Now, a few comments :

  • there is not just one approximate exponential, but actually a sequence of approximations. Therefore, the complexity must involve the tolerance parameter.
  • One cannot overlook the stability of the approximation. At least two ingredients are at stake. Either the diagonalisation of the matrix can be ill-conditionned (the matrix is far from normal). Or (even worse) its spectrum is moderately large. Say that $30$ and $-30$ are eigenvalues, then the magnitude of $\exp A$ is $e^{30}$, a large number. But $e^{-30}$ is also an (extremely small) eigenvalue, and the action of the approximate $\exp A$ on the corresponding eigenvector will be totally eroneous. More generally, a moderate condition number of $A$ causes huge errors in the modes associated with the smallest eigenvalues. These phenomena are well known by numerical analysts when dealing with ODEs. And they must choose a version of the Euler scheme adapted to the matrix itself.
Denis Serre
  • 51,599
12

Moler's paper "Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later" contains the following extracts:

In estimating the time required by matrix computations it is traditional to estimate the number of multiplications and then employ some factor to account for the other operations. We suggest making this slightly more precise by defining a basic floating point operation, or “flop”, to be the time required for a particular computer system to execute the FORTRAN statement

$A(I; J) = A(I; J) + T^*A(I;K)$.

This involves one floating point multiplication, one floating point addition, a few subscript and index calculations, and a few storage references. We can then say, for example, that Gaussian elimination requires $n^3/3$ flops to solve an $n$-by-$n$ linear system $Ax = b$....

About $qn^3$ flops are required to evaluate [a Padé approximant]...

The scaling and squaring method of section 3 and some of the matrix decomposition methods of section 6 require on the order of $10$ to $20n^3$ flops and they obtain higher accuracies than those obtained with $200 n^3$ or more flops for the o.d.e. solvers...

Most of the computational cost [in block diagonal methods] lies in obtaining the basic Schur decomposition. Although this cost varies somewhat from matrix to matrix because of the iterative nature of the QR algorithm, a good average figure is $15 n^3$ flops, including the further reduction to block diagonal form.

and there are a few more nuggets buried in there.

Based purely on this (and without reading the paper you cite but noting that its title is "A new scaling and squaring algorithm for the matrix exponential", for which see also here and http://eprints.ma.man.ac.uk/1442/), I would guess that Al-Mohy and Higham's scaling and squaring method is in the ballpark of $Cn^3$ for $C \approx 15$.

Steve Huntsman
  • 15,258
  • The complexity of scaling and squaring depends on how many times one has to square, so ultimately it depends on the matrix at hand. – Federico Poloni May 17 '16 at 21:37
  • 4
    That would be funny to put in programming documentation: "This function is based on the paper 'Nineteen dubious ways to compute the exponential of a matrix' by Moler. Way 13 was used." – Christopher King May 17 '16 at 22:10
  • Well, it'd be fun for sure. Another funny thing is to have the Python documentation for each function reporting the relevant estimated time complexity, although sometimes it can be determined by looking at the source code. – FaCoffee May 18 '16 at 08:10