Wrong matrix inverse results in R

Question

I calculated the inverse of a matrix (I-Q) (I is the identity matrix) in both R and Mathematica, but R gives me wrong results compared with the theoretical results. I have attached the code in R and Mathematica, and you can see the results are different.

Code in 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

Code in 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 results match the theoretical result as the row sums of the inverse of (I-Q) should be

208827064576
208827064550
208827063900
208827046974
208826606924
208815183200
208518148800
200795254400

I have no idea why the results differ and would appreciate any help.

Answer 1

Mathematica will do exact computations here, R will do floating-point computation. The condition number of the matrix you're trying to invert is very large

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

gives an estimate of 1.04346e+13. As you can read on 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

Given that, I'm actually surprised that the difference from the exact answer is as small as it is (e.g., checking just the first value):

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

You may already know this, but you can also check the canonical SO R-oriented "why are these numbers not equal?" FAQ