Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

MatrixPower[]

• To: mathgroup at smc.vnet.net
• Subject: [mg28919] MatrixPower[]
• From: Konstantin L Kouptsov <klk206 at nyu.edu>
• Date: Fri, 18 May 2001 01:13:26 -0400 (EDT)
• Organization: New York University
• Sender: owner-wri-mathgroup at wolfram.com

Hello, Mathgroup,

I needed to calculate the n-th power of a given matrix M.
Usually it can be done with the MatrixPower[M,n] function.
But for the symbolic power it becomes a bit tricky.

Take M={{1,1,1},{0,1,0},{1,1,1}},
then
MatrixPower[M,n] gives
  {{2^(-1+n),-1+2^n,2^(-1+n)},
   {0, 1,0},
   {2^(-1+ n),-1+2^n,2^(-1+n)}}

Now, if you change the matrix to
M={{1,1,1},{0,2,0},{1,1,1}}, then the MatrixPower function
gives an error

  MatrixPower::"zvec": "Cannot compute MatrixPower because the
  argument has a zero eigenvector."

which seems irrelevant, since the power of a square matrix always
exist.

The reason for this error is in the way Mathematica calculates the
power. Here is the explicit form of the used algorithm

  {eval, evec} = Eigensystem[M0];
  s = Transpose[evec];
  A = DiagonalMatrix[eval];
  M[n_] := s.MatrixPower[A, n].Inverse[s]

which of course does not work if s has no inverse.

To circumvent the problem one can use the Hamilton-Cailey theorem
(which can be found in the Korn,Korn mathematical handbook).
The theorem gives a recipe to calculate any function of a matrix
using Vandermonde matrix built on its eigenvalues and a finite
number of powers of the matrix itself.

Below is the implementation of this theorem. There is one thing to watch
out: it the matrix has equal eigenvalues, the Vandermonde matrix is
zero, so you get indeterminate 0/0. The way the theorem can be applied 
is that you can resolve indeterminancy
by cancelling out zeroes in numerator and denominator, as in

	(2-2)(7-5)    1
        __________  = _

	(2-2)(9-3)    3

which must be done in mathematically correct way (remember calculus?)

I used the function ReplaceEqual[] to discern the similar eigevalues
in the following manner:

	ReplaceEqual[{1,1,2,3,4,4,5}] gives

		{{1,1+i1,2,3,4,4+i2,5},{i1,i2}}

i.e. shift similar eigenvalues to make them different, and at the end
calculate the result by setting i1,i2 to zero. This is done by the following
code

ReplaceEqual[x_List] := Module[{i = 1, a, k, ls1 = {}, ls2, s, x1},
    x1 = Sort[x]; a = First[x1]; ls2 = List[a];
    For[k = 2, k <= Length[x1], k++,
      If[x1[[k]] == a,
        s = ToExpression["i" <> ToString[i++]];
        AppendTo[ls2, x1[[k]] + s]; AppendTo[ls1, s],
        AppendTo[ls2, x1[[k]]]; a = x1[[k]]
        ]
      ];
    {ls2, ls1}
    ]


Now the Vandermonde matrix determinant can be calculated in a rather 
efficient way


VandermondeMatrix[eval_List] := Module[{x, y},
    x = Table[1, {Length[eval]}];
    y = (eval^# &) /@ Range[1, Length[eval] - 1];
    PrependTo[y, x]
    ]

but we also need the determinant of the Vandermonde matrix with one of 
the columns replaced by {F[eval[[1]]],F[eval[[2]]],...}

VandermondeMatrixReplaced[eval_List,f_, n_] := Module[{w, fs},
    w = VandermondeMatrix[eval]; (* or pull this call out for efficiency *)
    fs = Map[f, eval];
    ReplacePart[w, fs, n]
    ]

Calculation of this deteminant was a problem which caused my recent messages to
MathGroup. Rasmus Debitsch have suggested using the row echelon form
of the matrix. For smaller matrices simple Det[] works.

Finally we arrive at the HamiltonCailey[] function which computes any
function of a square matrix:

HamiltonCailey[A_, F_] := Module[{eig, n, w, wr, x},
    {eig, vars} = ReplaceEqual[Eigenvalues[A]];
    n = Dim[A]; (* dimension of the n x n matrix A *)
    w = VandermondeMatrix[eig];
    wr = (1/Det[w]) Sum[
          Det[VandermondeMatrixReplaced[eig, F, n - k + 1]]* 
            MatrixPower[A, n - k], {k, 1, n}]
    ]

Do not forget to apply the Limit[] to the result. There is still a lot
to be said about the above method, but you get the idea.

For the matrix in above example
M={{1,1,1},{0,2,0},{1,1,1}}
you write

Limit[HamiltonCailey[M,#^n&]//Simplify,i1->0]

to get the answer

{{2^(-1+n),n 2^(-1+n),2^(-1+n)},
 {0,2^n,0},
 {2^(-1+n),n 2^(-1+n),2^(-1+n)}}

Konstantin.

PS. I wish to thank Rasmus Debitsch and all who responded to my post.

• Follow-Ups:
  • Re: MatrixPower[]
    • From: Daniel Lichtblau <danl@wolfram.com>