R 中错误的矩阵求逆结果

Wrong matrix inverse results in R

我计算了矩阵的逆(I-Q)(I是单位矩阵)在 R 和 Mathematica 中,但与理论结果相比,R 给了我错误的结果。我附上了R和Mathematica中的代码,你可以看到结果是不同的。

R 中的代码:

> Q <- matrix(c(25/26, 1/26, 0,    0,    0,    0,    0,    0, 
+               24/26, 1/26, 1/26, 0,    0,    0,    0,    0, 
+               25/26, 0,    0,    1/26, 0,    0,    0,    0, 
+               24/26, 1/26, 0,    0,    1/26, 0,    0,    0, 
+               24/26, 0,    0,    1/26, 0,    1/26, 0,    0, 
+               24/26, 1/26, 0,    0,    0,    0,    1/26, 0,
+               24/26, 1/26, 0,    0,    0,    0,    0,    1/26, 
+               24/26, 1/26, 0,    0,    0,    0,    0,    0), 8, 8, byrow = T)
> N <- solve(diag(8) - Q)
> N
             [,1]       [,2]      [,3]     [,4]     [,5]    [,6]     [,7]     [,8]
[1,] 200487454281 8019974160 308460545 11881443 456978.6 17576.1 676.0038 26.00015
[2,] 200487454255 8019974160 308460545 11881443 456978.6 17576.1 676.0038 26.00015
[3,] 200487453631 8019974134 308460545 11881443 456978.6 17576.1 676.0038 26.00015
[4,] 200487437380 8019973484 308460519 11881443 456978.6 17576.1 676.0038 26.00015
[5,] 200487014904 8019956584 308459869 11881417 456978.6 17576.1 676.0038 26.00015
[6,] 200476047392 8019517857 308442995 11880767 456952.6 17576.1 676.0038 26.00015
[7,] 200190875205 8008110293 308004242 11863867 456302.6 17550.1 676.0038 26.00015
[8,] 192776398346 7711513615 296596678 11424465 439402.5 16900.1 650.0037 26.00014
> N %*% rep(1, 8)
             [,1]
[1,] 208828245685
[2,] 208828245659
[3,] 208828245009
[4,] 208828228083
[5,] 208827788031
[6,] 208816364242
[7,] 208519328162
[8,] 200796390082

Mathematica 中的代码:

Q := {{25/26, 1/26, 0, 0, 0, 0, 0, 0},
  {24/26, 1/26, 1/26, 0, 0, 0, 0, 0},
  {25/26, 0, 0, 1/26, 0, 0, 0, 0},
  {24/26, 1/26, 0, 0, 1/26, 0, 0, 0},
  {24/26, 0, 0, 1/26, 0, 1/26, 0, 0},
  {24/26, 1/26, 0, 0, 0, 0, 1/26, 0},
  {24/26, 1/26, 0, 0, 0, 0, 0, 1/26},
  {24/26, 1/26, 0, 0, 0, 0, 0, 0}}


Inverse[DiagonalMatrix[{1, 1, 1, 1, 1, 1, 1, 1}] - Q]
{{200486320346,8019928800,308458800,11881376,456976,17576,676,26},
{200486320320,8019928800,308458800,11881376,456976,17576,676,26},
{200486319696,8019928774,308458800,11881376,456976,17576,676,26},
{200486303446,8019928124,308458774,11881376,456976,17576,676,26},
{200485880972,8019911224,308458124,11881350,456976,17576,676,26},
{200474913522,8019472500,308441250,11880700,456950,17576,676,26},
{200189742948,8008065000,308002500,11863800,456300,17550,676,26},
{192775308024,7711470000,296595000,11424400,439400,16900,650,26}}


Inverse[DiagonalMatrix[{1, 1, 1, 1, 1, 1, 1, 1}] - Q].{1, 1, 1, 1, 1,1, 1, 1}
(208827064576
208827064550
208827063900
208827046974
208826606924
208815183200
208518148800
200795254400
)

Mathematica 结果与理论结果相符,因为 (I-Q) 的倒数的行总和应该是

208827064576
208827064550
208827063900
208827046974
208826606924
208815183200
208518148800
200795254400

我不知道为什么结果不同,希望得到任何帮助。

Mathematica 将在此处进行精确计算,R 将进行浮点计算。您要反转的矩阵的条件数 非常大

Matrix::condest(diag(8)-Q)

给出的估计值为 1.04346e+13。您可以在 Wikipedia

上阅读

As a rule of thumb, if the condition number κ ( A ) = 10^k , then you may lose up to k digits of accuracy on top of what would be lost to the numerical method due to loss of precision from arithmetic methods

鉴于此,我真的很惊讶与确切答案的差异 small 因为它是(例如,只检查第一个值):

x1 <- 208828245685
x2 <- 208827064576
all.equal(x1,x2)
## [1] "Mean relative difference: 5.655887e-06"

您可能已经知道这一点,但您也可以查看规范的 SO R-oriented "why are these numbers not equal?" FAQ