Abrupt increase of evaluation time of MatrixExp

Question

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.

Answer 1

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.95253982970992315484748662011721804771`15.954589770191005*^308}
{5.24424, 2.0507373364279982188754557848574924788`15.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
*)

Answer 2

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

So your example becomes

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.

Answer 3

(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}]];
(* Developer`FromPackedArray::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.37580859653465140013918335099168246442`15.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} *)