Re: [R] identity matrix

From: David Scott <d.scott_at_auckland.ac.nz>
Date: Mon 31 Oct 2005 - 08:29:52 EST

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@auckland.ac.nz


Graduate Officer, Department of Statistics



R-help@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 Received on Mon Oct 31 08:36:44 2005

This archive was generated by hypermail 2.1.8 : Mon 31 Oct 2005 - 15:15:28 EST