A function to compute the eigenvalues in R when some are repeated

Question

Following this post and more specifically this answer, I am trying to build from scratch an eigen() function in R to compute eigenvalues and eigenvectors. However, I have stumbled accross the problems of repeated eigenvalues for which the aforementionned solution does not work.

Is there a simple and clever way to get around this issus? I understand that it comes from the solveroot() function as in the case of repeated eigenvalues, it does not change sign. I can not think of a workaround though on what is the eigenvalues and how many times it is repeated.

Here is what I have for the moment.

# a matrix with two pairs of repeated eigenvalues
mat <- matrix(c( 1, .5,  0,  0,
                .5,  1,  0,  0,
                 0,  0,  1, .5,
                 0,  0, .5,  1), 4, 4)
svd(mat)

# compute the determinant
eig <- function(lambda, mat) {
  y <-  mat - diag(lambda, nrow(mat))
  return(det(y))
}

# A function to find potential 0 with uniroot
solveroot <- function(f, lower, upper, n = 100, tol = .Machine$double.eps^.5, ...){
  xseq <- seq(lower, upper, length = n + 1)
  mod <- sapply(X = xseq, FUN = f, ...)
  Equi <- xseq[which(mod == 0)]    # if any already = 0
  ss   <- mod[1:n] * mod[2:(n+1)]  # check sign
  ii   <- which(ss < 0) 
  for (i in ii){
    Equi <- c(Equi, uniroot(f = f,
                            lower = xseq[i],
                            upper = xseq[i + 1],
                            tol = tol,
                            ...)$root)
  }
  return(Equi)
}
# find the eigenvalues
solveroot(f = eig, lower = 0.0001, upper = sum(diag(mat)), mat = mat) 
# no eigenvalues are returned.

This code works for matrix which do not have repeated eigenvalues.

Thank you,

Answer 1

The determinant touches, but does not cross, 0 at the two repeated eigenvalues. (Similar to how x^2 is never negative, but has both roots at 0.) AFAIK, uniroot is not suitable for finding those roots. Use optim instead:

mat <- matrix(c( 1, .5,  0,  0,
                 .5,  1,  0,  0,
                 0,  0,  1, .5,
                 0,  0, .5,  1), 4, 4)
svd(mat)$d
#> [1] 1.5 1.5 0.5 0.5

# compute the determinant
eig <- function(lambda, mat) {
  y <-  mat - diag(lambda, nrow(mat))
  return(det(y))
}

# A function to find potential 0 with uniroot
solveroot <- function(f, lower, upper, n = 100L, tol = .Machine$double.eps^.5, ...){
  xseq <- seq(lower, upper, length = n + 1L)
  mod <- sapply(X = xseq, FUN = f, ...)
  Equi <- xseq[which(mod == 0)]    # if any already = 0
  ss   <- mod[1:n] * mod[2:(n+1L)]  # check sign
  ii   <- which(ss < 0) 
  for (i in ii){
    Equi <- c(Equi, uniroot(f = f,
                            lower = xseq[i],
                            upper = xseq[i + 1],
                            tol = tol,
                            ...)$root)
  }
  
  mod <- abs(mod)
  ii <- which(mod[2:n] <= mod[1:(n-1)] & mod[2:n] <= mod[3:(n+1)] & ss[1:(n-1)] > 0 & ss[2:n] > 0) + 1L
  for (i in ii){
    sol <- optim(par = xseq[i],
                 fn = f,
                 lower = xseq[i-1],
                 upper = xseq[i+1],
                 method = "L-BFGS-B",
                 ...)
    if (sol$value <= tol) {Equi <- c(Equi, sol$par)}
  }
  return(Equi)
}
# find the eigenvalues
solveroot(f = eig, lower = 0.0001, upper = sum(diag(mat)), mat = mat)
#> [1] 0.500001 1.499999