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
MatrixExp
? $\endgroup$