A function for calculating the eigenvalues of a matrix in R

Question

I want to write a function like eigen() to calculating eigenvalues and eigenvectors of an arbitary matrix. I wrote the following codes for calculation of eigenvalues and I need a function or method to solve the resulted linear equation.

eig <- function(x){
       if(nrow(x)!=ncol(x)) stop("dimension error")
          ff <- function(lambda){
                for(i in 1:nrow(x)) x[i,i] <- x[i,i] - lambda
                }
det(x)
}

I need to solve det(x)=0 that is a polynomial linear equation to find the values of lambda. Is there any way?

Answer 1

Here is one solution using uniroot.all:

library(rootSolve)
myeig <- function(mat){
  myeig1 <- function(lambda) {
    y = mat
    diag(y) = diag(mat) - lambda
    return(det(y))
  }

  myeig2 <- function(lambda){
    sapply(lambda, myeig1)
  }
  uniroot.all(myeig2, c(-10, 10))
}

R > x <- matrix(rnorm(9), 3)
R > eigen(x)$values
[1] -1.77461906 -1.21589769 -0.01010515
R > myeig(x)
[1] -1.77462211 -1.21589767 -0.01009019

Answer 2

Computing determinant is such a bad idea as it is not numerically stable. You can easily get Inf etc even for a moderately big matrix. I suggest reading the following answers (read them otherwise you have no idea what my code is doing):

• Are eigenvectors returned by R function eigen() wrong?
• eigenvectors when A-lx is singular with no solution

then use either of the following

NullSpace(A - diag(lambda, nrow(A)))
nullspace(A - diag(lambda, nrow(A)))

Answer 3

The solution from @liuminzhao won't work if there is two repeated eigenvalues. The function will fail to find the roots, because the characteristic polynomial of the matrix will not change sign (it is zero and does not cross the zero line), which is what rootSolve::uniroot.all() is doing when looking for roots. So you need another way to find a local minima (like optim()). Moreover, it will failed to determine the number of repeated eigenvalues.

A better way is to find the characteristic equation with, which is easily done with pracma::charpoly() and then using polyroot().

par <- pracma::charpoly(M) # find parameters of the CP of matrix M
par <- par[length(par):1]  # reverse order for polyroot()
roots <- Re(polyroot(par)) # keep real part of the polyroot()

The pracma::charpoly() is not too complicated in itself, see its source code, starting at line a1 <- a.