Неправильная обратная матрица приводит к R - PullRequest
Новая
       23

Неправильная обратная матрица приводит к R

2 голосов
/ 09 марта 2020

Я рассчитал обратную матрицу ( 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

Я понятия не имею, почему результаты отличается и был бы признателен за любую помощь.

1 Ответ

6 голосов
/ 09 марта 2020

Mathematica выполнит точные вычисления здесь, R выполнит вычисления с плавающей точкой. номер условия матрицы, которую вы пытаетесь инвертировать, очень велик 

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

дает оценку 1,04346e + 13. Как вы можете прочитать в Википедии ,

Как правило, если номер условия κ (A) = 10 ^ k, то вы можете потерять до k цифр Точность сверх того, что будет потеряно в численном методе из-за потери точности из арифметических c методов

Учитывая, что я действительно удивлен, что отличие от точного ответа таково: small как есть (например, проверка только первого значения): 

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

Возможно, вы уже знаете это, но вы также можете проверить каноническую SO R-ориентированную "почему эти числа не равны? FAQ

...