identity matrix

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

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

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
All 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