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

Question

I have no trouble implementing a code to find the biggest eigenvalue, and corresponding eigenvector of a matrix using the power method.

What I have more trouble with, is thinking of a code that can output all eigenvalues and eigenvectors of a given matrix at once. I am able to do it manually on a small matrix, but can't seem to properly generalize it. I suspect it can be done in a beautiful way with some recursion but I'd need some help on that.

EDIT: Also I don't have trouble finding all the eigenvalues either, it's the eigenvectors that cause me trouble

Here would be the code that does it manually:

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

The function to get one 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)

}

Now the example by hand:

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

I am using the answer to the following question to do this:

Answer to: Power method for finding all eigenvectors

How to make a clean function that does everything at one out of that?

I have to admit that it was just for the sake of trying to implement my own thing for practice. I know of eigen() indeed. Maybe the power method isn't even meant to be used to get everything, and I should just use it to get just the biggest one. (I know of QR and stuff like that) — Joel H

Zac Zac · Accepted Answer · 2018-04-11T15:39:51

A recursive function like this seems to work:

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

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

1 Answers