An error of matrix inversion using NumPy

Discussion in 'Python' started by lancered, Apr 4, 2007.

  1. lancered

    lancered Guest

    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!
    lancered, Apr 4, 2007
    #1
    1. Advertising

  2. On 4 Apr 2007 06:15:18 -0700, lancered <> 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.


    Please post it.

    --
    mvh Björn
    =?ISO-8859-1?Q?BJ=F6rn_Lindqvist?=, Apr 4, 2007
    #2
    1. Advertising

  3. lancered

    Dayong Wang Guest

    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 <> wrote:
    > On 4 Apr 2007 06:15:18 -0700, lancered <> 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.

    >
    > Please post it.
    >
    > --
    > mvh Björn
    >
    Dayong Wang, Apr 4, 2007
    #3
  4. lancered

    Lou Pecora Guest

    In article <>,
    "lancered" <> 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.
    Lou Pecora, Apr 4, 2007
    #4
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Robert M. Gary

    Matrix inversion algorithm examples

    Robert M. Gary, Feb 10, 2006, in forum: Java
    Replies:
    11
    Views:
    34,582
    ersid
    Oct 29, 2008
  2. Jan-Hendrik Huehne

    matrix inversion

    Jan-Hendrik Huehne, Feb 7, 2004, in forum: C++
    Replies:
    2
    Views:
    634
    Victor Bazarov
    Feb 7, 2004
  3. Rafal 'Raf256' Maj

    [ot?] matrix inversion

    Rafal 'Raf256' Maj, Jan 14, 2004, in forum: C Programming
    Replies:
    87
    Views:
    1,856
    Dan Pop
    Jan 19, 2004
  4. lancered
    Replies:
    6
    Views:
    753
    Travis Oliphant
    Apr 4, 2007
  5. Sprechen sie von C++

    Matrix inversion

    Sprechen sie von C++, Jun 19, 2010, in forum: C++
    Replies:
    14
    Views:
    835
    Francesco S. Carta
    Jun 20, 2010
Loading...

Share This Page