3

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?

Mahmoud
  • 844
  • 1
  • 10
  • 17

3 Answers3

3

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
liuminzhao
  • 2,385
  • 17
  • 28
  • you have the right tools, but can you write one `myeig` function that only takes a matrix as input and returns the eigenvalues? – flodel Apr 21 '13 at 14:47
  • Thank you for your advice. I wrapped up and updated it in my answer. – liuminzhao Apr 21 '13 at 18:05
  • Thank you liuminzhao. For determining eigenvectors I should solve another system of linear equations for every eigenvalue. This is the codes I wrote using your help: `for(i in 1:nrow(mat)){ solve(mat-values[i]*diag(nrow(mat)),rep(0,nrow(mat)))}` But this returs just zero vector that isn't a eigrnvector. How can I find this answer? – Mahmoud Apr 21 '13 at 19:28
  • Can `uniroot.all()` find the complex roots in `myeig2` in above example? If no, which function can do such? – Mahmoud Apr 23 '13 at 14:35
  • @Mahmood I don't think it can do that. [polyroot](http://stat.ethz.ch/R-manual/R-patched/library/base/html/polyroot.html) can find complex roots, but I think you need to get the polynomial coefficients first. – liuminzhao Apr 23 '13 at 14:44
  • This solution does not work for all matrices, e.g., `R <- matrix(c(1, .5, .5, .5, 1, .5 , .5, .5, 1), 3 , 3)` – POC Jan 09 '22 at 21:31
1

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)))
RandomUser
  • 11
  • 1
0

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.

POC
  • 268
  • 1
  • 2
  • 7