CODEWITHSUNDEEP
reigenvectoreigenvaluelinear-equation
Mahmoud
Mahmoud

Reputation: 864

A function for calculating the eigenvalues of a matrix in R

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?

Upvotes: 3

Views: 6281

Answers (3)

POC
POC

Reputation: 278

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.

Upvotes: 0

RandomUser
RandomUser

Reputation: 11

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

then use either of the following

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

Upvotes: 1

liuminzhao
liuminzhao

Reputation: 2455

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

Upvotes: 3

Related Questions