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>