3

我想写一个函数eigen()来计算任意矩阵的特征值和特征向量。我写了以下代码来计算特征值，我需要一个函数或方法来求解得到的线性方程。

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

我需要解决det(x)=0一个多项式线性方程来找到 的值lambda。有什么办法吗？

4

3 回答 3

3

这是使用的一种解决方案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
于 2013-04-21T14:37:39.547 回答
1

计算行列式是个坏主意，因为它在数值上不稳定。Inf即使对于中等大小的矩阵，您也可以轻松获得等。我建议阅读以下答案（阅读它们，否则你不知道我的代码在做什么）：

然后使用以下任一

NullSpace(A - diag(lambda, nrow(A)))
nullspace(A - diag(lambda, nrow(A)))
于 2018-10-16T23:11:59.757 回答
0

如果有两个重复的特征值，@liuminzhao 的解决方案将不起作用。该函数将无法找到根，因为矩阵的特征多项式不会改变符号（它为零并且不穿过零线），这就是rootSolve::uniroot.all()寻找根时所做的事情。因此，您需要另一种方法来找到局部最小值（如optim()）。此外，它将无法确定重复特征值的数量。

更好的方法是用 找到特征方程，这很容易完成，pracma::charpoly()然后使用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()

本身并pracma::charpoly()不太复杂，请参阅它的源代码，从第 1 行开始a1 <- a。

于 2022-01-10T19:10:11.693 回答