当一些重复时计算R中的特征值的函数

Question

按照这个post，更具体地说是这个answer，我正在尝试在R 中从头开始构建一个eigen() 函数来计算特征值和特征向量。但是，我偶然发现了重复特征值的问题，而上述解决方案对此不起作用。

有没有一种简单而聪明的方法来解决这个问题？我知道它来自solveroot() 函数，因为在重复特征值的情况下，它不会改变符号。我想不出什么是特征值以及重复多少次的解决方法。

这是我目前拥有的。

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

此代码适用于没有重复特征值的矩阵。

谢谢，

Answer 1

行列式在两个重复的特征值处接触但不交叉 0。 （类似于 x^2 永远不会为负，但两个根都为 0。）AFAIK，uniroot 不适合找到这些根。请改用optim：

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