用于查找所有特征值和特征向量的幂法代码(在 R 中)

Code for power method to find all eigenvalues and eigenvectors ( in R)

我可以毫不费力地实现一个代码来找到最大的特征值,以及使用幂法的矩阵的相应特征向量。

我比较麻烦的是想一个可以一次输出给定矩阵的所有特征值和特征向量的代码。我可以在一个小矩阵上手动完成,但似乎无法正确概括它。 我怀疑它可以通过一些递归以一种漂亮的方式完成,但我需要一些帮助。

编辑:我也不难找到所有的特征值,是特征向量给我带来了麻烦

下面是手动执行此操作的代码:

#Example matrix
test.set=matrix(0,4,4)
test.set[1,]=c(-2,2,1,5)
test.set[2,]=c(2,5,8,8)
test.set[3,]=c(4,2,6,3)
test.set[4,]=c(5,-2,4,9)

获得一个的函数Eigenvalue/Eigenvector

#Power method simple : return biggest egeinvalue and corresponding eigenvector
power_method_simple=function(A,n_rep) {


  #Initialize with a random column of the matrix
  b_0=A[,sample(1:ncol(A),size=1)]

  for (k in 1:n_rep) {

  b_0=A%*%b_0
  b_0_norm=sqrt(t(b_0)%*%b_0)
  b_0=b_0/b_0_norm[1,1]
  print(b_0)
  }

eigenvalue=(t(b_0)%*%A%*%b_0)[1,1]
res_list=list(b_0,eigenvalue)
names(res_list)=c("vector","eigenvalue")
return(res_list)

}

现在手写例子:

#################
#By hand


#First
res1=power_method_simple(test.set,n_rep=1000)

first_ev=res1$vector
first_value=res1$eigenvalue


#Second
residual_matrix1=test.set-first_value*first_ev%*%t(first_ev)

res2=power_method_simple(residual_matrix1,n_rep=1000)

second_value=res2$eigenvalue
second_ev=(second_value-first_value)*res2$vector + first_value*((t(first_ev)%*%res2$vector)[1,1])*first_ev
second_ev=second_ev/sqrt((t(second_ev)%*%second_ev)[1,1])

#Third
residual_matrix2=residual_matrix1-second_value*res2$vector%*%t(res2$vector)


res3=power_method_simple(residual_matrix2,n_rep=1000)
third_value=res3$eigenvalue

u3=(third_value-second_value)*res3$vector + second_value*((t(res2$vector)%*%res3$vector)[1,1])*res2$vector
u3=u3/sqrt((t(u3)%*%u3)[1,1])

third_ev=(third_value-first_value)*u3 + first_value*((t(first_ev)%*%u3)[1,1])*first_ev
third_ev=third_ev/sqrt((t(third_ev)%*%third_ev)[1,1])


#I works for first three
print(eigen(test.set)$vector)
print(cbind(first_ev,second_ev,third_ev))

我正在使用以下问题的答案来执行此操作:

Answer to: Power method for finding all eigenvectors

如何制作一个干净的函数来一次完成所有事情?

像这样的递归函数似乎有效:

find_vals=function(matrix, values=list(), vectors=list(), evs=list(), n=nrow(matrix)) {
  if (n<1) return(list(values, evs))
  res=power_method_simple(matrix,n_rep=1000)
  curr_val = res$eigenvalue
  res_v = res$vector
  i = nrow(matrix) - n + 1
  values[i] = curr_val
  vectors[[i]] = res_v
  if (i == 1) {
    evs[[i]] = res_v
  } else {
    curr_v = vectors[[i]]
    for (k in (i-1):1) {
      curr_v = (values[[i]] - values[[k]])*curr_v + values[[k]]*((t(vectors[[k]])%*%curr_v)[1,1])*vectors[[k]]
      curr_v=curr_v/sqrt((t(curr_v)%*%curr_v)[1,1])
    }
    evs[[i]] = curr_v
  }
  matrix=matrix-curr_val*res_v%*%t(res_v)
  return (find_vals(matrix, values, vectors, evs, n-1))
}