Velocity Reviews > An error of matrix inversion using NumPy

An error of matrix inversion using NumPy

lancered
Guest
Posts: n/a

 04-04-2007
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.

=?ISO-8859-1?Q?BJ=F6rn_Lindqvist?=
Guest
Posts: n/a

 04-04-2007
On 4 Apr 2007 06:15:18 -0700, lancered <(E-Mail Removed)> wrote:
> 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().

Could it have something to do with floating point accuracy?

>>> r = matrix([[random.random() * 9999 for x in range(19)] for y in range(19)])
>>> allclose(linalg.inv(r) * r, identity(19))

True

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

--
mvh Björn

Dayong Wang
Guest
Posts: n/a

 04-04-2007
Here below is the matrix KK I used:

[[ 1939.33617572 -146.94170404 0. 0. 0.
0. 0. 0. 0.
-1172.61032101
0. 0. -193.69962687 -426.08452381 0.
0. 0. 0. 1792.39447168]
[ -146.94170404 5175.33392519 -442.430839 0. 0.
0. 0. 0. 0. 0.
-3409.58801135 0. 0. 0.
-767.46969697
-408.90367384 0. 0. 4585.96138215]
[ 0. -442.430839 4373.33685159 0. 0.
0. 0. 0. 0. 0.
0. -2354.70959362 0. 0. 0.
0. -855.36922061 -720.82719836 3930.90601259]
[ 0. 0. 0. 17064.73017917
-949.49987581 0. 0. 0. 0.
0. 0. 0. -16115.23030336 0.
0. 0. 0. 0.
16115.23030336]
[ 0. 0. 0. -949.49987581
11005.53604312 -1358.01000599 0. 0. 0.
0. 0. 0. 0.
-8698.02616132
0. 0. 0. 0.
8698.02616132]
[ 0. 0. 0. 0.
-1358.01000599
18322.57142994 -1495.29428718 0. 0. 0.
0. 0. 0. 0.
-15469.26713677
0. 0. 0. 15469.26713677]
[ 0. 0. 0. 0. 0.
-1495.29428718 12497.65812936 -1858.81899672 0. 0.
0. 0. 0. 0. 0.
-9143.54484546 0. 0. 9143.54484546]
[ 0. 0. 0. 0. 0.
0. -1858.81899672 20170.17739075 -1249.5298217 0.
0. 0. 0. 0. 0.
0. -17061.82857234 0. 17061.82857234]
[ 0. 0. 0. 0. 0.
0. 0. -1249.5298217 9476.04289846 0.
0. 0. 0. 0. 0.
0. 0. -8226.51307677 8226.51307677]
[ -1172.61032101 0. 0. 0. 0.
0. 0. 0. 0. 1500.8055591
-328.1952381 0. 0. 0. 0.
0. 0. 0. -1172.61032101]
[ 0. -3409.58801135 0. 0. 0.
0. 0. 0. 0. -328.1952381
4112.15248021 -374.36923077 0. 0. 0.
0. 0. 0. -3409.58801135]
[ 0. 0. -2354.70959362 0. 0.
0. 0. 0. 0. 0.
-374.36923077 2729.07882439 0. 0. 0.
0. 0. 0. -2354.70959362]
[ -193.69962687 0. 0. -16115.23030336 0.
0. 0. 0. 0. 0.
0. 0. 17726.91399397 -1417.98406375 0.
0. 0. 0. -16308.92993023]
[ -426.08452381 0. 0. 0.
-8698.02616132
0. 0. 0. 0. 0.
0. 0. -1417.98406375 12320.46305747
-1778.36830859 0. 0. 0.
-9124.11068513]
[ 0. -767.46969697 0. 0. 0.
-15469.26713677 0. 0. 0. 0.
0. 0. 0. -1778.36830859
19552.18019195 -1537.07504962 0. 0.
-16236.73683374]
[ 0. -408.90367384 0. 0. 0.
0. -9143.54484546 0. 0. 0.
0. 0. 0. 0.
-1537.07504962
12983.70625768 -1894.18268877 0. -9552.44851929]
[ 0. 0. -855.36922061 0. 0.
0. 0. -17061.82857234 0. 0.
0. 0. 0. 0. 0.
-1894.18268877 21039.17951514 -1227.79903343 -17917.19779295]
[ 0. 0. -720.82719836 0. 0.
0. 0. 0. -8226.51307677 0.
0. 0. 0. 0. 0.
0. -1227.79903343 10175.13930856 -8947.34027513]
[ 1792.39447168 4585.96138215 3930.90601259 16115.23030336
8698.02616132 15469.26713677 9143.54484546 17061.82857234
8226.51307677 -1172.61032101 -3409.58801135 -2354.70959362
-16308.92993023 -9124.11068513 -16236.73683374 -9552.44851929
-17917.19779295 -8947.34027513 85023.67196244]]

On 4/4/07, BJörn Lindqvist <(E-Mail Removed)> wrote:
> On 4 Apr 2007 06:15:18 -0700, lancered <(E-Mail Removed)> wrote:
> > 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().

>
> Could it have something to do with floating point accuracy?
>
> >>> r = matrix([[random.random() * 9999 for x in range(19)] for y in range(19)])
> >>> allclose(linalg.inv(r) * r, identity(19))

> True
>
> > 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.

>
>
> --
> mvh Björn
>

Lou Pecora
Guest
Posts: n/a

 04-04-2007
In article <(E-Mail Removed) om>,
"lancered" <(E-Mail Removed)> 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.
>
>

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.