Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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
share|improve this question
3  
Why you try to re-implement MatrixExp? –  ybeltukov Dec 9 '13 at 12:52

1 Answer 1

up vote 3 down vote accepted

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.

share|improve this answer
1  
IdentityMatrix is not compilable. Maybe you gain another microsecond by using something which doesn't call MainEvaluate :-) –  halirutan Jul 14 at 1:40

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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