0

I used %*% to multiple a matrix and its inverse. I don't get the identity matrix. What am i missing?

D 

    [,1] [,2] [,3]
[1,] 1 2 3

[2,] 4 2 1

[3,] 2 2 0

solve(D)

       [,1]       [,2]       [,3]
[1,] -0.1428571 0.4285714 -0.2857143

[2,] 0.1428571 -0.4285714 0.7857143

[3,] 0.2857143 0.1428571 -0.4285714

D %*% solve(D)

          [,1]          [,2]          [,3]
[1,] 1.000000e+00 0.000000e+00 -2.220446e-16

[2,] -5.551115e-17 1.000000e+00 0.000000e+00

[3,] -1.110223e-16 -1.110223e-16 1.000000e+00
Tchotchke
  • 3,061
  • 3
  • 22
  • 37
user2956302
  • 1
  • 3
    That is the identity matrix, if you squint. `all.equal( diag(3), D%*%solve(D) ) # TRUE` – Frank Oct 16 '15 at 19:11
  • Not as obvious, but this is an example of the [R FAQ 7.31](https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f). – r2evans Oct 17 '15 at 05:42

2 Answers2

5

You're not getting back exactly zero for the off-diagonals because of floating-point precision errors.

You can see that this really is the identity matrix if you round:

round(D %*% solve(D))
C_Z_
  • 7,427
  • 5
  • 44
  • 81
0

As Frank mentioned, it's precision errors. Usually e-16 ish numbers and smaller are a good indicator this happens. Also consider

> 10/3-3-1/3
[1] 1.665335e-16

Clearly, we'd consider this to be 0.

In addition to r2evans's answer, the answer to this question has a far more detail about any language. Why are these numbers not equal?

Community
  • 1
  • 1
mateen
  • 83
  • 5