Displaying 1 result from an estimated 1 matches for "x_eigen".
Did you mean:
m_eigen
2010 Jul 17
0
Computing the power of a matrix efficiently and accurately
...ple 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 dia...