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
我计算了矩阵的逆(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