The Museum of HP Calculators

HP Forum Archive 14

[ Return to Index | Top of Index ]

free42s matrix inverse
Message #1 Posted by hugh steers on 11 Nov 2004, 8:44 p.m.

hi there,

ive been going through the 42s owners manual trying stuff on free42. its great! now i'm even more impressed that every little detail appears to be there.

ive gotten onto trying out the matrices and for fun, ive feed it some nasty ones. there appears to be a problem with inverse that’s not down to floating point error nor the method used (ie lu-decomposition & back substitution. i'm reading the code).

i'm trying this

3 1 2 1 1e-040 1e-040 2 1e-040 -1e-040

ok, so its nearly singular (this example is treated in the 15c advanced function handbook p105). i tried it with my own matrix code which gives the right answer using 1eee754 doubles as,

-2e-040 3 -1 3 -4e+040 2e+040 -1 2e+040 -1e+040

but free42 is wildly out. i tied to debug it. i don’t think its your use of “tiny” for the pivot because that happens to be calculated as 1e-20 (the same as i use). so i suspect the backsubstitution. my version doesn’t have an inner loop. is there a bug here?

and many thanks for making all the fix releases. much appreciated.

      
Re: free42s matrix inverse
Message #2 Posted by whuy on 12 Nov 2004, 2:06 a.m.,
in response to message #1 by hugh steers

Well, that's my point. If the matrix is singular to calculator precision, the result you get is effectively 'anything and nothing'. A slight change - algorithmic, order of operations - may return completely different answers, that's what a badly conditioned problem is all about. The result is meaningless, and what is worse - you are not told that it is, there is no 'singularity indicator'.

Werner

            
Re: free42s matrix inverse
Message #3 Posted by hugh steers on 12 Nov 2004, 7:05 a.m.,
in response to message #2 by whuy

yes indeed, i completely agree. although ive found that this example is within the limits of 1eee754, albeit out of range for the 15c.

this came up whilst testing the matrix implementation. i'm thinking that free42s should be able to cope with this case. what im doing is comparing the results of my own matrix library with free42's answers. both use the same underlying number representation.

      
Re: free42s matrix inverse
Message #4 Posted by Thomas Okken on 12 Nov 2004, 6:41 p.m.,
in response to message #1 by hugh steers

I'll look into it this weekend. BTW, you mention your code does give the correct answer. Does your code agree with mine as far as the LU decomposition is concerned?
I'll just step through the code with the NR original for reference, but maybe we can narrow this down a bit beforehand.
Thanks!

P.S. I'm posting another bug fix release (1.0.8) as I write this, but it does NOT fix this particular bug yet, just a couple of minor ones.

Again, thanks for the bug report!

            
Re: free42s matrix inverse
Message #5 Posted by hugh steers on 12 Nov 2004, 7:28 p.m.,
in response to message #4 by Thomas Okken

hi there,

ive done some more investigations using different inversion techniques as an experiment. firstly, it looks like i was too hasty in my accusation of a bug in free42 and Whuy was right in this respect. what started it, is that inverting the matrix,

3                    1                    2
1                    1e-040               1e-040
2                    1e-040               -1e-040

gave the exact correct answer in my code when using lu decomposition then backsubstitution. however, examples using numbers less far out than 10^40 wind up with worse answers. this leads me to think that the result was a fluke. here’s what i get with the above under different inversion algorithms (all use 1eee754 doubles):

original matrix:
3                    1                    2
1                    1e-040               1e-040
2                    1e-040               -1e-040
LU Inverse Matrix:
-2e-040              3                    -1
3                    -4e+040              2e+040
-1                   2e+040               -1e+040
SVD Inverse Matrix:
2.77555756156289e-017 0.2                  0.4
0.2                  -0.12                -0.24
0.4                  -0.24                -0.48
Gauss Inverse Matrix:
0                    0                    0.5
0                    0                    -0.5
0.5                  0.5                  -0.75

as you can see, the lu inverse is spot on but the other methods, like svd (which is sometimes better) are way off. i now think this is a fluke and your code is different because its different and i just got lucky.

taking 10^-40 down to 10^-10, say, i get a different picture:

original matrix:
3                    1                    2
1                    1e-010               1e-010
2                    1e-010               -1e-010
LU Inverse Matrix:
-2.0000000012e-010   3.0000000018         -1.0000000006
3.0000000018         -40000000027         20000000009
-1.0000000006        20000000009          -10000000003
SVD Inverse Matrix:
-1.99999246971831e-010 2.99999065204657     -0.999995325723285
2.99998375179818     -39999820941.8902    19999910466.4451
-0.999991875599092   19999910466.4451     -9999955231.72258
Gauss Inverse Matrix:
-2.00000016548074e-010 3                    -0.9999999997
3                    -39999996690.3854    19999998340.6927
-0.9999999997        19999998340.6927     -9999999168.84636

now the lu isn’t spot on anymore and the others are still worse. trying again, in range at 10^-5

original matrix:
3                    1                    2
1                    1e-005               1e-005
2                    1e-005               -1e-005
LU Inverse Matrix:
-2.00012000720043e-005 3.00018001080065     -1.00006000360022
3.00018001080065     -400027.001620097    200009.000540032
-1.00006000360022    200009.000540032     -100003.000180011
SVD Inverse Matrix:
-2.00012000714772e-005 3.00018001076159     -1.00006000358069
3.00018001074169     -400027.001613791    200009.000536879
-1.00006000357074    200009.00053688      -100003.000178434
Gauss Inverse Matrix:
-2.00012000719087e-005 3.00018001078885     -1.00006000359432
3.00018001078885     -400027.001617928    200009.000538948
-1.00006000359432    200009.000538948     -100003.000179468

we’re seeing about 9 correct digits for lud & svd. so the picture is more like what you'd expect, except that 10^-40 just happend to work for lud :-)

                  
Re: free42s matrix inverse
Message #6 Posted by Thomas Okken on 13 Nov 2004, 8:52 p.m.,
in response to message #5 by hugh steers

I'm glad to hear it! Having copied the code straight out of Numerical Recipes, I was kinda hoping it would be OK. :-D
Still intrigued, though. My first gut feeling was that with numbers 40 powers of 10 apart, you wouldn't have a chance in hell of getting any result -- I'm thinking with IEEE-754 having just about 16 decimal digits, you'd be looking at an end result with -24 significant digits. When I tried e-10 instead of e-40, my results were consistent with that thought.
Oh, well, I'm still intrigued but I guess this issue does not count as an actual bug, then. So I can enjoy a quiet Sunday!

Cheers,

- Thomas


[ Return to Index | Top of Index ]

Go back to the main exhibit hall