[R] identity matrix

[David Scott]

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