free42s matrix inverse
#1

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.

#2

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

#3

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.

#4

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!

#5

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 :-)

#6

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



Possibly Related Threads…
Thread Author Replies Views Last Post
  AFTER HP-Prime update, Shift+Matrix CRASHES Joseph Ec 3 1,999 12-06-2013, 11:06 AM
Last Post: Joseph Ec
  HP Prime Matrix TERRIBLE bug and question uklo 19 4,959 11-25-2013, 12:10 PM
Last Post: Mic
  HP Prime: editing a matrix Alberto Candel 6 2,297 11-20-2013, 06:26 PM
Last Post: Helge Gabert
  Absolute Value and Matrix BruceTTT 5 2,055 11-11-2013, 11:52 PM
Last Post: Walter B
  WP-34S Matrix operations with routine-local registers? Tom Grydeland 1 1,155 09-04-2013, 10:46 AM
Last Post: Marcus von Cube, Germany
  Matrix Characteristic Polynomial - Reloaded. Ángel Martin 12 3,320 08-22-2013, 05:33 PM
Last Post: Thomas Klemm
  Matrix Richard Berler 3 1,390 08-18-2013, 06:24 PM
Last Post: Paul Dale
  Advantage/CCD Matrix Challenge Ángel Martin 1 1,141 08-09-2013, 06:22 PM
Last Post: Thomas Klemm
  [HP -Prime CAS] List, Matrix, Vector as one Array? CompSystems 0 972 07-26-2013, 05:22 PM
Last Post: CompSystems
  Inverse binomial Richard Berler 7 2,530 07-09-2013, 06:23 AM
Last Post: Marcus von Cube, Germany

Forum Jump: