Possible Method for MatrixExp



Well, probably a hard question, but I think it's better to cry out loud :).

I noticed that MatrixExp has a Method option when writing this answer, which is undocumented. Since this tutorial claims that MatrixExp uses variable-order Padé approximation, evaluating rational matrix functions using Paterson–Stockmeyer methods or Krylov subspace approximations, I guess there might exist non-Automatic Method(s), but I can't find it out. Sadly, techniques in this post seem not to help. Does anyone know? Or is Automatic really the only available option?


Posted 2015-05-18T09:11:33.550

Reputation: 44 878

I never found an answer to my own question of a similar nature. I wonder if the Method option is simply added for future extensibility.

– Mr.Wizard – 2015-05-18T09:18:39.957


This seems as good an opportunity as any to leave these two links.

– J. M.'s ennui – 2015-05-18T10:02:42.713



Experimentation based on the documentation you quoted led to two valid Method options:

MatrixExp[{{1.2, 5.6}, {3, 4}}, Method -> "Pade"]
{{346.557, 661.735}, {354.501, 677.425}}
MatrixExp[{{1.2, 5.6}, {3, 4}}, {1, 2}, Method -> "Krylov"]
{1670.03, 1709.35}

If "Krylov" is used for the single parameter syntax it complains:

MatrixExp[{{1.2, 5.6}, {3, 4}}, Method -> "Krylov"]

MatrixExp::novec: The method Krylov requires the specification of a vector. >>

"Pade" fails on exact or symbolic input:

MatrixExp[{{1, 2}, {3, 4}}, Method -> "Pade"]

MatrixExp::invmtd: Invalid method Pade for input with precision ∞. >>

Thanks to ilian we now have the last(?) piece of the puzzle:

MatrixExp[{{2, 0, 0}, {0, 1, -1}, {0, 1, 1}}, Method -> "BlockDecomposition"]
{{E^2, 0, 0}, {0, E Cos[1], -E Sin[1]}, {0, E Sin[1], E Cos[1]}}


Posted 2015-05-18T09:11:33.550

Reputation: 259 163


From the OP's quote, that would seem to be the only two settings. In particular, the "Krylov" setting is valuable for exponentiating sparse matrices, and is quite useful for the action of the exponential of a sparse matrix on a vector (that is, MatrixExp[A, v, Method -> "Krylov"] for computing $\exp(\mathbf A)\mathbf v$).

– J. M.'s ennui – 2015-05-18T09:34:28.740

@J. M. What method do you suppose is used for MatrixExp[{{1, 2}, {3, 4}}] as neither "Pade" nor "Krylov" is accepted? (I mean is there something that correlates to Automatic in this case.) – Mr.Wizard – 2015-05-18T09:36:11.850

Certainly, one does not use a Padé approximant for computing exact results. A look at the docs tells me that either Jordan decomposition or Putzer's algorithm is used.

– J. M.'s ennui – 2015-05-18T09:58:19.260

@J. M. Neither "Jordan" nor "Putzer" work either. I am just curious if there is a third valid setting or if whatever is used is only accessible via Automatic. (Oh, and thank you.) – Mr.Wizard – 2015-05-18T10:00:19.643

1I'm guessing MatrixExp[] is smart enough not to try to compute the eigenvalues of a symbolic matrix (imagine how many Root[] objects will that yield!), and will just use Putzer exclusively in that case. For matrices with exact entries, prolly a toss-up, and the precise details are known only to WRI devs. – J. M.'s ennui – 2015-05-18T10:06:03.850

1On that note, I'm surprised that there is no option to use the Schur decomposition for computing MatrixExp[] for inexact matrices, but a little bird told me MatrixFunction[] does use Schur. I guess one can use that if one suspects that Padé is returning something weird. – J. M.'s ennui – 2015-05-18T10:28:55.570

4The third valid setting can be accessed with Method -> "BlockDecomposition". – ilian – 2015-05-18T12:52:57.527

@ilian, interesting. That's the Schur-based method, I take it? – J. M.'s ennui – 2015-05-18T13:15:11.827

1More along the lines of splitting the matrix into blocks and applying the Putzer or Jordan decomposition to each. – ilian – 2015-05-18T13:39:10.563

@ilian Thanks! Added to my answer, with credit of course. – Mr.Wizard – 2015-05-18T13:41:23.187

@ilian, thanks for the detail! If it isn't too revealing, can you tell how MatrixExp[] picks between Jordan and Putzer when given a matrix with exact entries? – J. M.'s ennui – 2015-05-18T14:13:13.457

1Can't really go into more detail on the heuristic, but your guess above is pretty accurate in that all 'non-simple' cases tend to go to Putzer. – ilian – 2015-05-18T14:49:19.413

1Useful discussions and reponses. I have an inquiry. I am dealing with a very large sparse matrix and its matrix exponential on a vector. I am computing the solution using MatrixExp[mat, vector, Method -> "Krylov"]; // AbsoluteTiming, but it is still a bit time-consuming due to a very large size of my matrix. So, I want to know that is there a way to speed up this computation by imposing a tolerance or compiling so as to gain some speed? A result of lower accuracy is OK for me – Fazlollah – 2017-12-29T18:57:46.537

1@Fazlollah that is probably worth posting as a separate Question if you have not done so already. I at least don't have an answer for you. – Mr.Wizard – 2017-12-30T04:57:38.770