R help - Jul 2010 - Computing the power of a matrix efficiently and accurately

Good day,

I would like to know if there is an efficient and accurate method for
computing the power of a tridiagonal matrix in R? The "obvious" method
would
be eigen-decomposition, but I find that this is not very accurate, probably
because the dimensions of the matrix I am considering are large and the
eigenvectors are not well approximated. Here is an example of what I am
trying to do:

set.seed(1234)

m=500 # dimension of my matrix

X=diag(rnorm(500,20,2)) # set up tridiagonal matrix
for(i in 1:(m-1))
{
  X[i,i+1]=rnorm(1,5,1)
  X[i+1,i]=rnorm(1,5,1)
}

n=100 # the power of the matrix X i.e. X^n = X%*%X%*%X ... n times (n = 1000
in my application)

X_eigen=eigen(X)
P=X_eigen$vectors
P_inv=solve(P)
d=X_eigen$values

M=diag(m)

res_mult=c()
res_eigen=c()
for(i in 1:n)
{
# straightforward multiplication...
M = M%*%X
res_mult = rbind(res_mult, M[m,])

# eigen-decomposition...
res_eigen = rbind(res_eigen, (P%*%diag(d^i)%*%P_inv)[m,])
}

# Look how the diagonal element X[m,m] changes with successive
powers:
> res_mult[1,m]
[1] 23.45616
> res_eigen[1,m]
[1] 23.45616
> res_mult[n,m]
[1] 3.844556e+148
> res_eigen[n,m]
[1] 3.844556e+148

# Look how the off-diagonal element X[m,1] changes with successive
powers:
> res_mult[1,1]
[1] 0
> res_eigen[1,1]
[1] 1.942362e-16
> res_mult[n,1]
[1] 0
> res_eigen[n,1]
[1] 1.469650e+132


Notice how poorly the off-diagonal element is approximated by the
eigen-decomposition!! I am aware of the mtx.exp() function in the Malmig
package, although repeatedly calling this function is not efficient as I
need to extract the m-th row of the matrix raised to each power from 1:n.

Any suggestions? Could I possibly exploit the fact that X is tridiagonal in
my application?

Thanks!
Miguel

-- 
Miguel Lacerda
Department of Statistical Sciences
University of Cape Town

	[[alternative HTML version deleted]]