A bug and a related mini-challenge - Printable Version +- HP Forums (https://archived.hpcalc.org/museumforum) +-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html) +--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html) +--- Thread: A bug and a related mini-challenge (/thread-134903.html) |
A bug and a related mini-challenge - Rodger Rosenbaum - 03-26-2008 I recently came across an unexpected bug in the HP48G (and descendants). The LSQ command is supposed to return the solution with minimum Euclidean length when given an underdetermined system (with an infinite number of solutions). Given system A*x=B, with A:
[[ 1 1 0 0 0 -1 ] and B:
[[ 4 ] LSQ gives a solution:
[[6.99541284404 ] with a Euclidean length of 14.2255745922. But this is not the minimum length solution. For example, [ 6 2 1 9 3 4 ]T is a solution and has Euclidean length 12.124355653, but there IS a solution with minimum possible length. The mini-challenge is to find it, using your HP48G, HP49G+ or HP50G. Don't use MATLAB or something similar until you give it a shot on your HP calculator.
I was very surprised by this bug. I've used LSQ a lot, and never before seen such behavior, and I have no idea what particular circumstance leads to the error. Maybe Werner can find out.
Re: A bug and a related mini-challenge - Chuck - 03-27-2008 I gave it the old fashioned (by hand) Least Squares method, namely x = (Trans[A] A)^-1 Trans[A] B. The problem is... the determinant of (Trans[A] A) = 0, which means it has no inverse, thus the least squares method fails (I forget if they call this a critical matrix???) or gives erroneous results. Other "minimal" methods are needed. I'll play with it some more.
CHUCK
Re: A bug and a related mini-challenge - Chris Dean - 03-28-2008 Programmatically using 6 loops (sledgehammer approach) and four simultaneous equations I found that the minimum Euclidean length (for integer values) was 10.95445. I did this labouriously on my HP17BII+. I will not spoil the challenge for anyone else by supplying the 6 values. The equations I used were
Regards
Chris Dean
Re: A bug and a related mini-challenge - Chuck - 03-28-2008 After hacking at it with the HP, I threw Mathematica at it, which produced a startling nice solution of {3,4,1,7,6,3} with a length of 10.95 (So I assume this is the same solution Chris achieved). So, two questions now... 1) How was LeastSquares implemented inside Mathematica, and 2) How to get the HP to find this nice answer. CHUCK [Addendum]
I may have discovered what Mathematica is doing. I believe it is finding a pseudoinverse for the matrix Transpose[A].A, otherwise known as the Moore-Penrose inverse. This type of inverse exists even for singular matrices, and also gives the least-squares solution. Now, how to get the HP50 to find this inverse?? Edited: 28 Mar 2008, 7:26 p.m.
Re: A bug and a related mini-challenge - Rodger Rosenbaum - 03-28-2008 You and Chris both got the correct minimum length. If you simply swap the first and second columns in the original A matrix, the LSQ command will give the correct solution. LSQ just doesn't like that original A matrix. A little program to find the pseudoinverse is: PSEU << DUP SIZE OBJ-> DROP DUP2 IF > THEN SWAP DROP IDN SWAP TRN LSQ TRN ELSE DROP IDN SWAP LSQ END >>
If you apply it to the original A matrix you will get a bad result (because it's using LSQ), but if you swap the first and second columns of that matrix, then PSEU will give you a correct pseudoinverse, which can be used to premultiply the B matrix to give the correct solution.
Re: A bug and a related mini-challenge - Rodger Rosenbaum - 03-28-2008 Quote:You don't need to find the pseudoinverse of Transpose[A].A All you need to do is find the pseudoinverse of A and use that to premultiply B.
Forming Transpose[A].A is unnecessary and squares the condition number.
Re: A bug and a related mini-challenge (solution) - Chuck - 03-29-2008 Thanks Roger. With your suggestions, I wrote a program for the HP50g which calculates the pseudo-inverse using SVD (singular value decomposition) instead of LSQ, which means it works without switching any columns or rows. I haven't programmed the 50g for quite some time, so the program is probably very inefficient. Here it is:
PSEUI This will return the Pseudo-Inverse matrix. Then, premultiplying by B will give the desired result, [3 4 1 7 6 3]. I'm sure the program can be cleaned up and improved a bit, but I'm satisfied I even was able to do it. Thanks for the challenge; I learned quite a bit more about linear algebra and HP50g commands. CHUCK
Edited: 29 Mar 2008, 3:22 p.m.
Re: A bug and a related mini-challenge (solution) - Rodger Rosenbaum - 03-29-2008 Very good! I think you may be the only person on this forum that I've gotten interested enough to do something with the SVD on their own. In your program, if you use TRN instead of TRAN, compatibility with the HP48G will be maintained (as long as complex numbers aren't involved). Also, you said "...premultiplying by B will give the desired result...". I think that there is an implied "pseudoinverse(A)" in there; in other words you're saying "...premultiplying pseudoinverse(A) by B will give the desired result...", but that is backwards. You need to premultiply B by pseudoinverse(A) to get the desired result.
So, you should have said "...postmultiplying by B will give the desired result...".
Re: A bug and a related mini-challenge (solution) - Chuck - 03-30-2008 Thanks Rodger. Your syntax makes sense ( I was never good at English, hence my going into mathematics). Multiply "PseudoInverse(A) by B" makes the most sense, now that I look at it. I still need to experiment with my program a bit more to see what downfalls it may have. All in all, this was a great challenge to expand my linear algebra knowledge a bit.
CHUCK
Re: A bug and a related mini-challenge (solution) - Rodger Rosenbaum - 03-30-2008 You have to be really careful here. Saying: 'Multiply "PseudoInverse(A) by B"...' is ambiguous. Consider the following two statements: (1) Premultiply Pseudoinverse(A) by B and: (2) Postmultiply Pseudoinverse(A) by B The first one means B * Pseudoinverse(A) while the second means Pseudoinverse(A) * B and since matrix multiplication is not commutative, they can (and usually will) give different results even if defined. In the case of the problem we've been considering, the first multiplication can't be carried out because the matrices in that order aren't conformable for multiplication.
So, to avoid ambiguity, if you're describing a matix multiplication in english rather than with mathematical symbology, be sure to use "premultiply" or "postmultiply".
Re: A bug and a related mini-challenge (solution) - Reth - 03-30-2008 I wonder how many people would be interested in the topic (worldwide) and can HP48 express the ratio with decent precision :) Edited: 30 Mar 2008, 6:48 a.m.
Re: A bug and a related mini-challenge (solution) - Werner - 03-31-2008 I disagree!
Cheers, Werner
Re: A bug and a related mini-challenge (solution) - Rodger Rosenbaum - 03-31-2008 Hi, Werner. I hoped I would attract your attention! You noticed the last thing I said in message #1, I guess.
I remember that you found a tiny bug in the 48G's eigen function for a particular matrix. This anomaly is similarly peculiar.
Re: A bug and a related mini-challenge (solution) - Werner - 04-02-2008 Hi Rodger,
I looked at the code, it performs a QR decomposition to turn A into an upper trapezoidal R, and then determines the rank before
Cheers, Werner
Re: A bug and a related mini-challenge (solution) - Werner - 04-02-2008 Update #1: Mhh well. It doesn't seem to be a bug after all.. just the limits of accuracy and tolerance that have been reached.I'll go on investigating the differences in the two cases anyway.
Cheers, Re: A bug and a related mini-challenge (solution) - Rodger Rosenbaum - 04-02-2008 I have a book, "Linear Algebra Investigations with the HP48G/GX", by Donald LaTorre, and in the acknowledgement section he says: "I must express my deep appreciation to...Paul McClellan of the software design team at Hewlett Packard for his splendid development of calculator versions of the LAPACK Fortran code for the matrix operations that are built in to the HP48G/GX."
Hence the similarity to LAPACK routines.
Re: A bug and a related mini-challenge (solution) - Rodger Rosenbaum - 04-02-2008 Quote: What do the 4 lines labeled 1x1, 2x2, etc., represent? If I calculate the singular values of the original A matrix with the SVD command, I get: [ 2 2 2 5.74289198711E-15 ]
It's strange that a problem with only small integers as input data would run into such difficulties.
Re: A bug and a related mini-challenge (solution) - Werner - 04-03-2008 Hi Rodger,
Cheers, Werner
Re: A bug and a related mini-challenge (solution) - Werner - 04-03-2008 OK, I looked it up. The routine is a blueprint of LAPACK's DGELSX (or newer DGESLY) routine (http://www.netlib.org/lapack/double/dgelsx.f), but there, the tolerance is an *input* parameter.. To me, using 1e-14 as tolerance is the issue here. Since we're only working with *estimates* for the singular values, that tolerance is a bit too strict. 1e-13 or 5e-12 would've been more appropriate. And alas, we have no way of providing that ourselves.
Cheers, Werner
Re: A bug and a related mini-challenge (solution) - Rodger Rosenbaum - 04-03-2008 Werner, would you mind sending me a private email with an email address I can use to communicate with you?
|