1

How does Mathematica compute MatrixPower[m,n]?

The documentation states it is equivalent to Dot @@ ConstantArray[m, {n}], but this construction is limited to $n\leq 4096$ and simple Timing exercises show it is slower.

This answer to a similar question mentions JordanDecomposition, but I don't know if that is what Mathematica actually uses.

In Linear Algebra classes, it is often taught that diagonalization allows for fast exponentiation. But is that actually how Mathematica does it?

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
GregH
  • 1,909
  • 12
  • 25
  • 1
  • (I'm still curious to know the answer, though) – Michael Seifert Feb 24 '22 at 18:41
  • 1
    The Internals document's second paragraph states "... knowing about the internals of the Wolfram System may be of intellectual interest ...", to which I reply "Yes. Yes it is." I'm a pure mathematician at heart, and am less interested in implementation than others. But I don't like telling students "This is a great way to accomplish XYZ" when, in reality, no one actually does. So I wonder if MMA uses diagonalization; it looks like it does use the generalized form via JordanDecomposition, but I'm not positive. – GregH Feb 25 '22 at 13:31
  • For comparison, matrix_power in numpy.linang appears to simply perform repeated matrix multiplication. It gains some efficiency by (effectively) repeatedly squaring $M$ to get matrices of the form $M^{2n}$ and then using the binary decomposition of the desired exponent to express it as a product of matrices of that form. – Michael Seifert Feb 25 '22 at 14:40
  • 1
    ... So if Mathematica uses a similar algorithm then it wouldn't be a good example. However, this Wolfram page notes that "MatrixExp uses Putzer's method or Jordan decomposition" for exact numerical results, so maybe you can mention that to your students. – Michael Seifert Feb 25 '22 at 14:59
  • 2
    To modify @Michael's statement a little: Jordan or Putzer is used in the exact case, and the Schur decomposition is used in the inexact case for noninteger powers. – J. M.'s missing motivation May 25 '22 at 16:57
  • MatrixPower and MatrixExp do not work the same way in general. – Daniel Lichtblau May 25 '22 at 18:17
  • 1
    @DanielLichtblau: Makes sense. My comment was mainly directed towards the OP's desire to be able to tell his students about a "real-world" implementation of Jordan decomposition. – Michael Seifert May 25 '22 at 18:20

0 Answers0