# Abrupt increase of evaluation time of MatrixExp

I test MatrixExp with different matrix sizes. When its dimension is smaller than 1435*1435, the time consumed is not much. But when it goes to 1435, the time increases a lot.

AbsoluteTiming[
MatrixExp[RandomReal[{0, 1}, {1433, 1433}]]
] // First


0.76363

AbsoluteTiming[
MatrixExp[RandomReal[{0, 1}, {1434, 1434}]]
] // First


0.744058

AbsoluteTiming[
MatrixExp[RandomReal[{0, 1}, {1435, 1435}]]
] // First


7.308

Is there anything special about 1435? I found it on MMA 14.0 version. Maybe it is reproducible on your machine but with some other value instead of 1435.

• in V 14.1 on windows 10 with 128 GB RAM, I get 0.36 for the 1433 size and a wapping 157 seconds for the 1435 size. This could be related to which version on Intel MKL is used if any by Mathematica. I had similar sudden change of timing when doing Matrix rank. see cpu-timing-for-matrix-rank-calculation At matrix size over 2500, even by just one, a dramatic speed increase was seen. In your case it is sudden speed decrease. Commented Aug 9 at 3:19
• Finding what implementation Mathematica uses for MatrixExp will help. screen shot i.sstatic.net/L9935bdr.png Commented Aug 9 at 3:25

I think the problem is that around a dimension of 1434 or 1435, the values of MatrixExp overflow and become arbitrary precision. Perhaps the variation in the length of the computation from one to a couple hundred seconds depends on when in the process the overflow occurs. Here is the value where the overflow occurs:

$MaxMachineNumber Log[%] (* 1.79769*10^308 709.783 *)  Here are a couple of moderately slow 1433 x 1433 matrix-exp examples: SeedRandom[43]; mat = RandomReal[{0, 1}, {1433, 1433}]; AbsoluteTiming[MatrixExp[mat]] // {First[#], Max[#]} & SeedRandom[46]; mat = RandomReal[{0, 1}, {1433, 1433}]; AbsoluteTiming[MatrixExp[mat]] // {First[#], Max[#]} & (* {1.0524, 1.9525398297099231548474866201172180477115.954589770191005*^308} {5.24424, 2.050737336427998218875455784857492478815.954589770191005*^308} *)  Here is a relatively fast 2000 x 2000 example: SeedRandom[43]; (* any seed will do, I think *) mat = RandomReal[{0, 1/10}, {2000, 2000}]; AbsoluteTiming[MatrixExp[mat]] // {First[#], Max[#]} & (* {0.797571, 1.5334042747017165*^40} *)  The slower ones are arbitrary-precision and the faster are machine-precision. Sometimes the arbitrary-precision examples are only a little slower. The slowest example I have observed was a 1435 x 1435 matrix at 360 sec. This 5000 x 5000 example takes around 11 sec. even though it stays at machine precision. It's a fairly large matrix, though. SeedRandom[1]; mat = RandomReal[{0, 1/10}, {5000, 5000}]; AbsoluteTiming[MatrixExp[mat]] // {First[#], Max[#]} & (* {10.7786, 7.99613*10^104} *)  Note the max entry of the matrix exponential of an $$n \times n$$ matrix with uniformly distributed entries between $$0$$ and $$b$$ seems to be around $$\exp(b n/2)/n$$. Compare with the previous $$n=5000,\ b=1/10$$ example: Exp[1/10 * 5000 / 2.] / 5000 (* 7.49291*10^104 *)  I don't know if this asymptotic approximation is the case — surely it would be known if it is. It has just seemed to be the case in the examples I have run. It also gives an idea why $$n=1434,\ b=1$$ is around the the machine-to-arbitrary precision jump in timing: Exp[1434/2.*1]/1434$MaxMachineNumber
(*
1.708418923521*10^308
1.79769*10^308
*)

• Maybe I should have pointed out that machine-precision arithmetic is carried out by the hardware and arbitrary-precision is carried by libraries and tracks precision. Commented Aug 9 at 13:50
• Thank you. This is very informative. How did you approximate the max entry value by exp(bn/2)/n?
– Ming
Commented Aug 10 at 3:11
• @Ming That's a good question. As I meant to imply, I tested several examples and guessed. One expects a relationship between the matrix exp and $\exp|\lambda_{\text{max}}|$, where $|\lambda_{\text{max}}|$ is the magnitude of the largest eigenvalue (or largest real part, but we're talking heuristics here, not mathematical proof). This is connected to matrix norms, which are connected to row or column sums, which in turn should average $bn/2$. So I compared the max entry with $\exp(bn/2)$. By eye it seemed off by $1/n$. So I recalculated, and it seemed close enough to use for explanation. Commented Aug 10 at 11:48

It becomes much faster when you use the form with vector. From help

r = RandomReal[{0, 1}, {1435, 1435}];
v = RandomReal[{0, 1}, 1435];
AbsoluteTiming[ MatrixExp[r, v]] // First

(* 0.0041711 *)


Compared to 147 second without the v second argument.

This does not explain the sudden change in timing, but something to consider using.

Help says

ie. using v or not.

(Longer comment) Usually, functions use some heuristics to determine which method to use internally. This heuristics are rarely explained, and there can be somehow unexpected jumps in performance when the method switches. You can see that something different is indeed happening if you turn on the messages for unpacking:

On["Packing"];

MatrixExp[RandomReal[{0, 1}, {1433, 1433}]];

MatrixExp[RandomReal[{0, 1}, {1434, 1434}]];
(* DeveloperFromPackedArray::punpack1: Unpacking array with dimensions {1434,1434}. *)

Off["Packing"];


Furthermore, you can compare the outputs:

MatrixExp[RandomReal[{0, 1}, {1433, 1433}]][[1, 1]] // FullForm
(* 1.3919352140584585*^308 *)

MatrixExp[RandomReal[{0, 1}, {1434, 1434}]][[1, 1]] // FullForm
(* 2.3758085965346514001391833509916824644215.954589770191005*^308 *)


Note that I don't think there is anything special with this particular dimensions, $$1434\times 1434$$, but some other properties of the matrix are checked in the heuristics. You can see this if you run the code multiple times:

Table[MatrixExp[RandomReal[{0, 1}, {1434, 1434}]] // AbsoluteTiming //
First, {i, 10}]

(* {5.5406, 4.83421, 0.803066, 3.96151, 5.5701, 0.590244,
0.580606, 6.36402, 5.58893, 5.86998} *)