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$
  • 3
    $\begingroup$ Why you try to re-implement MatrixExp? $\endgroup$ – ybeltukov Dec 9 '13 at 12:52
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
    $\begingroup$ IdentityMatrix is not compilable. Maybe you gain another microsecond by using something which doesn't call MainEvaluate :-) $\endgroup$ – halirutan Jul 14 '14 at 1:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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