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.

• Prev by Date: Re: Creating graph with only a few data points
• Next by Date: Request for contour plotting algorithm
• Previous by thread: Re: Creating graph with only a few data points
• Next by thread: Re: MatrixPower[]