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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
