查找所有特征值和特征向量的幂方法代码（在 R 中）

Question

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

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

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

以下是手动执行的代码：

#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)

获取一个特征值/特征向量的函数

#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

如何制作一个干净的函数，可以同时完成所有操作？

Answer 1

这样的递归函数似乎可以工作：

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))
}