0
$\begingroup$

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
$\endgroup$
1
  • 3
    $\begingroup$ Why you try to re-implement MatrixExp? $\endgroup$
    – ybeltukov
    Dec 9, 2013 at 12:52

1 Answer 1

3
$\begingroup$

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.