2
$\begingroup$

For a square matrix $M$, consider $M^8$.

Does MatrixPower[M, 8] do seven matrix multiplications? Or is it smart enough to compute $((M^2)^2)^2$ that involves only three?

Or does Mathematica do something else?

$\endgroup$
4
  • $\begingroup$ The "Details" section of the documentation for MatrixPower states that "MatrixPower[m, n] effectively evaluates the product of a matrix with itself $n$ times." $\endgroup$
    – MarcoB
    Dec 7, 2016 at 6:51
  • 2
    $\begingroup$ @MarcoB I think a key word there is "effectively". It should give the same result, but I doubt that it does that since it wouldn't be the most efficient. I don't know what it does, but what makes sense to me is to switch methods based on whether the matrix is symbolic, exact numeric, inexact numeric or machine precision numeric, as well as based on the matrix size and the exponent. $\endgroup$
    – Szabolcs
    Dec 7, 2016 at 11:11
  • 2
    $\begingroup$ If you only want to know which method is faster, just benchmark it. $\endgroup$
    – Szabolcs
    Dec 7, 2016 at 11:13
  • 2
    $\begingroup$ This is just a guess, but for machine precision matrices a major concern could be numerical stability and hence a different method might be needed to ensure the answer is even correct. Now whether the 'nested power' method is as numerically stable as say diagonalization, I'm not sure... $\endgroup$
    – Greg Hurst
    Dec 9, 2016 at 18:18

1 Answer 1

6
$\begingroup$

Benchmark:

M = RandomReal[10, {10, 10}];
RepeatedTiming[M.M.M.M.M.M.M.M;, 10]
RepeatedTiming[MatrixPower[M, 8];, 10]
RepeatedTiming[temp = M.M;
 temp = temp.temp;
 temp.temp;, 10]

{7.7*10^-6, Null}

{3.1*10^-6, Null}

{3.60*10^-6, Null}

So it appears some optimization does take place, as the timing is much closer to the $((M^2)^2)^2$ method.

$\endgroup$
6
  • $\begingroup$ Repeat Feyre's benchmark with M = Array[x, {5, 5}]; and there's an even more dramatic difference. $\endgroup$
    – bill s
    Dec 7, 2016 at 13:53
  • $\begingroup$ Especially for high powers I could imagine it might be faster to diagonalize first if that's possible $\endgroup$ Dec 7, 2016 at 14:19
  • 3
    $\begingroup$ @JulesLamers The exact case e.g. integer or rational input would be ill suited for diagonalizing since that involves working in an algebraic extension field. For the inexact case this might make sense. For real input there will be fuzzy nonzero imaginary parts and offhand I don't know if they should be left alone or chopped off (probably the latter is my guess). $\endgroup$ Dec 7, 2016 at 15:57
  • 1
    $\begingroup$ @Jules, in addition to Daniel's comments, recall that not all matrices are diagonalizable; one can ostensibly use Schur before powering in the inexact case, tho. In the exact and symbolic cases, one could use Jordan, but then you hit the problem alluded in Daniel's first sentence. $\endgroup$ Dec 9, 2016 at 17:37
  • $\begingroup$ @J.M. Certainly, hence the "if that's possible" :) $\endgroup$ Dec 9, 2016 at 19:18

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.