# speed up for computing MatrixExp

My code use two method to computing MatrixExp, I expected the second snippet code will faster than the first one, but in fact it a bit slower than that. How can I imporve my code?

nnx = 100;
pp = RandomReal[{-1, 1}, 2 {nnx, nnx}];
{n, q, num} = {2, 3, 16};
(*-------------------------*)
t1 = AbsoluteTime[];
{nn, nn} = Dimensions[pp];
a = Sum[((2^-num)^i MatrixPower[pp, i])/i!, {i, n*q}];
res1 = IdentityMatrix[nn] + Nest[(2 # + #.#) &, a, num];
AbsoluteTime[] - t1
(*-------------------------*)
t1 = AbsoluteTime[];
{nn, nn} = Dimensions[pp];
pp1 = 2^-num*pp;
QQ = MatrixPower[pp1, #] & /@ Range[q];
NN = MatrixPower[pp1, q*(# - 1)] & /@ Range[n];
a = Sum[NN[[nn]].QQ[[qq]]/(qq + (nn - 1)*q)!, {nn, n}, {qq, q}];
res2 = IdentityMatrix[nn] + Nest[(2 # + #.#) &, a, num];
AbsoluteTime[] - t1
(*-------------------------*)
res1 - res2 // Chop // Flatten // Union
res1 - MatrixExp[pp] // Chop // Flatten // Union

• Why you try to re-implement MatrixExp? Dec 9, 2013 at 12:52

A certain speedup is achieved when computing a with a NestList and Compiling .

MainEvaluate is called only for the IdentityMatrix but i think this is ok because it is not the time consuming part of the code. I also used the "InlineExternalDefinitions" -> True option but the gain was very small.

nnx = 200;
pp = RandomReal[{-1, 1}, 2 {nnx, nnx}];
{n, q, num} = {2, 3, 16};
{nn, nn} = Dimensions[pp];

snip = Compile[{},
Module[{a, e = 2.^(-num)},
a = Total@(NestList[#.pp &, pp, n*q - 1]*FoldList[#1*e/#2 &, e, Range[2, n*q]]);
IdentityMatrix[nn] + Nest[(2 # + #.#) &, a, num]
], CompilationOptions -> {"InlineExternalDefinitions" -> True}];


You can check timings:

t1 = AbsoluteTime[];
res3 = snip[];
AbsoluteTime[] - t1


You can check the correctnes with Chop[res1 - res3] === Table[0, {2 nnx}, {2 nnx}]

I experimented with compiling to C but there was no improvement.

• IdentityMatrix is not compilable. Maybe you gain another microsecond by using something which doesn't call MainEvaluate :-) Jul 14, 2014 at 1:40