相关疑难解决方法(0)

我试图在R中计算P任意N x J矩阵的投影矩阵S:

P = S (S'S) ^ -1 S'

我一直试图用以下功能执行此操作:

P <- function(S){
  output <- S %*% solve(t(S) %*% S) %*% t(S)
  return(output)
}

但是当我使用它时,我得到的错误看起来像这样:

# Error in solve.default(t(S) %*% S, t(S), tol = 1e-07) : 
#  system is computationally singular: reciprocal condition number = 2.26005e-28

我认为这是数值下溢和/或不稳定的结果在许多地方一样讨论R-帮助和这里,但我没有足够的经验使用SVD或QR分解来解决这个问题,要不然就把这个现有的代码到行动.我也尝试了建议的代码,即将solve写成一个系统:

output <- S %*% solve (t(S) %*% S, t(S), tol=1e-7)

但它仍然无效.任何建议,将不胜感激.

我很确定我的矩阵应该是可逆的并且没有任何共线性,只是因为我尝试用正交虚拟变量矩阵进行测试,但它仍然不起作用.

另外,我想将它应用于相当大的矩阵,所以我正在寻找一个简洁的通用解决方案.

我需要以下对角线:

diag(X %*% solve(A) %*% t(X))

其中A是满秩矩阵,X是矩形矩阵.这两个A和X稀疏.

我知道发现矩阵的逆是坏的,除非你真的需要它.但是,我无法看到如何重写公式,以便用两个参数solve(A)替换solve,这样线性系统就可以在没有显式反转的情况下得到解决.那可能吗？