I found a very odd thing. A matrix multiplied by its inverse matrix should be an identity matrix. But why the following thing happens? <a%*%solve(a) is not an identity matrix>> x%*%t(x)[,1] [,2] [,3] [,4] [,5] [1,] 108.16 58.24 32.24 66.56 225.68 [2,] 58.24 31.36 17.36 35.84 121.52 [3,] 32.24 17.36 9.61 19.84 67.27 [4,] 66.56 35.84 19.84 40.96 138.88 [5,] 225.68 121.52 67.27 138.88 470.89> a=x%*%t(x) > solve(a)[,1] [,2] [,3] [,4] [,5] [1,] -1.372649e+14 1.078492e+14 -2.553747e+14 1.126842e+14 4.120191e+13 [2,] -1.174558e+14 1.543860e+14 2.323143e+14 -1.375074e+14 2.381809e+13 [3,] -3.062129e+14 6.914056e+14 -1.656744e+15 7.007732e+14 -1.672741e+12 [4,] 1.761208e+14 -3.724017e+14 7.773835e+14 -2.156631e+14 -3.575351e+13 [5,] 8.789836e+13 -8.046912e+13 6.984285e+13 -5.502429e+13 -1.510937e+13> a%*%solve(a)[,1] [,2] [,3] [,4] [,5] [1,] 2.0410156 1.4433594 -4.169922 0.4003906 1.2170410 [2,] -2.5146484 1.9208984 -1.84912w=1 1.1879883 1.1203613 [3,] -1.5356445 -0.8549805 -1.008789 0.4116211 0.6218872 [4,] -0.5898438 -0.2539063 3.592773 2.0351563 0.7849121 [5,] 1.5000000 -1.1289063 -2.070313 -0.4003906 2.1386719>
On Sun, 30 Oct 2005, Robert wrote:> I found a very odd thing. > A matrix multiplied by its inverse matrix should be an > identity matrix.Well, an invertible matrix multiplied by its inverse... The matrix that you give (to two decimal places) is singular> solve(a)Error in solve.default(a) : system is computationally singular: reciprocal condition number = 3.25149e-20 Your matrix is slightly less singular, but the huge size of the elements in solve(a) shows that a%*%solve(a) will not be accurately computed. There are limits to what you can do with only 15 digits accuracy. -thomas> But why the following thing happens? > <a%*%solve(a) is not an identity matrix> >> x%*%t(x) > [,1] [,2] [,3] [,4] [,5] > [1,] 108.16 58.24 32.24 66.56 225.68 > [2,] 58.24 31.36 17.36 35.84 121.52 > [3,] 32.24 17.36 9.61 19.84 67.27 > [4,] 66.56 35.84 19.84 40.96 138.88 > [5,] 225.68 121.52 67.27 138.88 470.89 >> a=x%*%t(x) >> solve(a) > [,1] [,2] [,3] > [,4] [,5] > [1,] -1.372649e+14 1.078492e+14 -2.553747e+14 > 1.126842e+14 4.120191e+13 > [2,] -1.174558e+14 1.543860e+14 2.323143e+14 > -1.375074e+14 2.381809e+13 > [3,] -3.062129e+14 6.914056e+14 -1.656744e+15 > 7.007732e+14 -1.672741e+12 > [4,] 1.761208e+14 -3.724017e+14 7.773835e+14 > -2.156631e+14 -3.575351e+13 > [5,] 8.789836e+13 -8.046912e+13 6.984285e+13 > -5.502429e+13 -1.510937e+13 >> a%*%solve(a) > [,1] [,2] [,3] [,4] > [,5] > [1,] 2.0410156 1.4433594 -4.169922 0.4003906 > 1.2170410 > [2,] -2.5146484 1.9208984 -1.84912w=1 1.1879883 > 1.1203613 > [3,] -1.5356445 -0.8549805 -1.008789 0.4116211 > 0.6218872 > [4,] -0.5898438 -0.2539063 3.592773 2.0351563 > 0.7849121 > [5,] 1.5000000 -1.1289063 -2.070313 -0.4003906 > 2.1386719 >> > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html >Thomas Lumley Assoc. Professor, Biostatistics tlumley at u.washington.edu University of Washington, Seattle
On Sun, 30 Oct 2005, Robert wrote:> I found a very odd thing. > A matrix multiplied by its inverse matrix should be an > identity matrix. > But why the following thing happens? > <a%*%solve(a) is not an identity matrix> >> x%*%t(x) > [,1] [,2] [,3] [,4] [,5] > [1,] 108.16 58.24 32.24 66.56 225.68 > [2,] 58.24 31.36 17.36 35.84 121.52 > [3,] 32.24 17.36 9.61 19.84 67.27 > [4,] 66.56 35.84 19.84 40.96 138.88 > [5,] 225.68 121.52 67.27 138.88 470.89 >> a=x%*%t(x) >> solve(a) > [,1] [,2] [,3] > [,4] [,5] > [1,] -1.372649e+14 1.078492e+14 -2.553747e+14 > 1.126842e+14 4.120191e+13 > [2,] -1.174558e+14 1.543860e+14 2.323143e+14 > -1.375074e+14 2.381809e+13 > [3,] -3.062129e+14 6.914056e+14 -1.656744e+15 > 7.007732e+14 -1.672741e+12 > [4,] 1.761208e+14 -3.724017e+14 7.773835e+14 > -2.156631e+14 -3.575351e+13 > [5,] 8.789836e+13 -8.046912e+13 6.984285e+13 > -5.502429e+13 -1.510937e+13All those e+13, e+14, e+15 values didn't ring any alarms for you? What you are seeing is just the limits of accuracy of calculation on a computer. Yes, A times A^{-1} is the identity when calculations are infinitely accurate, but there are limits on accuracy when using a computer. Lots of other nice mathematical properties go out the window with floating point computation: it is not associative or distributive for starters. Best do a bit of reading on floating point computation: the Wiki article is a readily available starting point: http://en.wikipedia.org/wiki/Floating_point David Scott _________________________________________________________________ David Scott Department of Statistics, Tamaki Campus The University of Auckland, PB 92019 Auckland NEW ZEALAND Phone: +64 9 373 7599 ext 86830 Fax: +64 9 373 7000 Email: d.scott at auckland.ac.nz Graduate Officer, Department of Statistics