An error of matrix inversion using NumPy

Lou Pecora pecora at anvil.nrl.navy.mil
Wed Apr 4 09:52:51 EDT 2007


In article <1175692518.282887.3040 at y66g2000hsf.googlegroups.com>,
 "lancered" <wangday at gmail.com> wrote:

> Hi dear all,
>         I am using Python2.4.2+NumPy1.0.1 to deal with a parameter
> estimation problem with the least square methods. During the
> calculations, I use  NumPy package to deal with matrix operations,
> mostly matrix inversion and trasposition. The dimentions of the
> matrices used are about 29x19,19x19 and 29x29.
> 
>         During the calculation, I noticed an apparent error of
> inverion of a 19x19 matrix. Denote this matrix as KK, U=KK^ -1, I
> found the product of U and KK is not equivalent to unit matrix! This
> apparently violate the definition of inversion. The inversion is
> through the function linalg.inv().
> 
>        I have checked that det(KK)=-1.2E+40. At first, I thought the
> error may be caused by such a large determinant, so I scaled it as
> LL=KK/100, then invert LL. Since det(LL)=11.5 and all its elements are
> within -180.0 to 850.0, this seems easier. But the result is still not
> correct, the product of LL^-1 thus obtained and LL still not unit
> matrix ... At the same time, the inversion results of some 29x19
> matrices are correct.
> 
>        So,  can you tell me what goes wrong?  Is this a bug in
> Numpy.linalg? How to deal with this situation?  If you need, I can
> post the matrix I used below, but it is so long,so not at the moment.
> 
>        Thanks in advance!
> 

Changing the overall scaling won't help unless the numbers themselves 
are near the double floating point boundary (perhaps).  The real number 
you want to calculate is the condition number of the matrix. Roughly the 
ratio of the largest to the smallest eigenvalue.  If that's large, then 
you're attempt at inversion is dealing with differences between very 
large and very small numbers and/or very small differences between two 
large numbers which probably goes beyond the number of digits (bits) the 
machine can provide to represent floating point numbers.  No rescaling 
of the original matrix will change that.  Do a google on "condition 
number".

-- Lou Pecora  (my views are my own) REMOVE THIS to email me.



More information about the Python-list mailing list