Re: MatrixPower[]

*To*: mathgroup at smc.vnet.net*Subject*: [mg28943] Re: [mg28919] MatrixPower[]*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Sat, 19 May 2001 22:28:01 -0400 (EDT)*References*: <200105180513.BAA06195@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Konstantin L Kouptsov wrote: > > 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. For nonnegative integer powers, yes. This limitation has been addressed in our development version. > 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. I am guessing this applies when/where the function can be identified with its Maclauren series? > 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. Here is another approach that may often be faster. Perturb the matrix in some sensible way (so that it takes on a full set of eigenvectors), take MatrixPower, and remove the perturbation with Limit. For example: In[8]:= mat = {{1,1,1},{0,2,0},{1,1,1}}; In[9]:= matp = mat; General::spell1: Possible spelling error: new symbol name "matp" is similar to existing symbol "mat". In[10]:= matp[[1,3]] += eps; In[11]:= mp = Limit[MatrixPower[matp,n], eps->0] n n n n 0 + 2 -1 + n -0 + 2 n Out[11]= {{-------, 2 n, --------}, {0, 2 , 0}, 2 2 n n n n -0 + 2 -1 + n 0 + 2 > {--------, 2 n, -------}} 2 2 The 0^n symbolic powers can of course be removed e.g. by Simplify[...,n>0]. There is another method that involves Z-transforms. Offhand I forget the details. I think these methods all share the characteristic that they can get into computationally stressful matrix algebra. Daniel Lichtblau Wolfram Research

**References**:**MatrixPower[]***From:*Konstantin L Kouptsov <klk206@nyu.edu>