0

我需要solve一个列表中的一千多个矩阵。但是，我得到了错误Lapack routine dgesv: system is exactly singular。我的问题是我的输入数据是非奇异矩阵，但是根据我需要对矩阵进行的计算，其中一些矩阵变得奇异。然而，我的数据子集的可重复示例是不可能的，因为它会很长（我已经尝试过）。这是我的问题的一个基本示例（经过一些计算，A 将是一个矩阵，而我需要做的下一个计算是 R）：

A=matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,0,1), nrow=4)
R = solve(diag(4)-A)

你有想法如何“解决”这个问题，也许是其他功能？或者如何非常非常轻微地变换奇异矩阵，以免产生完全不同的结果？谢谢

编辑：根据@Roman Susi，我包括我必须做的函数（所有计算）：

function(Z, p) {
  imp <- as.vector(cbind(imp=rowSums(Z)))
  exp <- as.vector(t(cbind(exp=colSums(Z))))
  x = p + imp
  ac = p + imp - exp  
  einsdurchx = 1/as.vector(x)    
  einsdurchx[is.infinite(einsdurchx)] <- 0
  A = Z %*% diag(einsdurchx)
  R = solve(diag(length(p))-A) %*% diag(p)    
  C = ac * einsdurchx
  R_bar = diag(as.vector(C)) %*% R
  rR_bar = round(R_bar)
  return(rR_bar)
}

问题出在函数计算的第 8 行solve(diag(length(p))-A)。Z在这里，我可以为and提供新的示例数据p，但是在此示例中它工作正常，因为我无法重新创建一个会带来错误的示例：

p = c(200, 1000, 100, 10)
Z = matrix(
  c(0,0,100,200,0,0,0,0,50,350,0,50,50,200,200,0),
  nrow = 4, 
  ncol = 4,
  byrow = T)

所以，根据@Roman Susi 的问题是：有没有办法改变之前的计算，以便det(diag(length(p))-A)永远不会得到 0 来solve得到方程？我希望你能明白我想要什么:) 想法，谢谢。Edit2：也许问得更容易：如何避免这个函数中的奇异性（至少在第 8 行之前）？

4

2 回答 2

3

MASS 包中的广义逆ginv可以处理奇异矩阵，但必须确定它是否对您的问题有意义。

 library(MASS)
 ginv(diag(4) - A)

给予：

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0
[4,]    0    0    0    0

ibdreg包里也有这个Ginv功能。

于 2015-11-06T19:37:57.097 回答
0

R 的 QR 分解函数可能有你的答案。它们提供了一种稳健地求解线性方程的方法。QR 分解不提供逆分解，而是矩阵分解，通常可以在使用逆分解的地方使用。

对于矩形矩阵，QR 分解可用于找到最小二乘拟合。对于正方形，（近）奇异矩阵，qr()检测这个（近）奇异性，qr.coef()然后可以用于获得没有任何错误但可能有一些 NA（可以转换为零）的解。

于 2016-12-10T08:25:48.073 回答