A bug and a related mini-challenge « Next Oldest | Next Newest »

 ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 03-26-2008, 09:19 PM 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 ] [ 1 0 -1 0 1 0 ] [ 0 1 1 1 0 0 ] [ 0 0 0 1 1 1 ]] ``` and B: ```[[ 4 ] [ 8 ] [ 12 ] [ 16 ]] ``` LSQ gives a solution: ```[[6.99541284404 ] [ 6.15137614679 ] [ 2.84403669725 ] [ 3.00458715596 ] [ 3.84862385321 ] [ 9.14678999083 ]] ``` 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. ▼ Chuck Senior Member Posts: 320 Threads: 59 Joined: Dec 2006 03-27-2008, 01:04 AM 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 Chris Dean Member Posts: 120 Threads: 9 Joined: Aug 2005 03-28-2008, 12:51 PM 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 ``` For the column maxtrix x=(A, B, C, D, E, F)T the four equations can be derived A + B - F = 4 A - C + E = 8 B + C + D = 12 D + E + F = 16 Euclidean length = Sqrt(A^2 + B^2 + C^2 + D^2 + E^2 + F^2) I wonder if problems are caused when determining the least Euclidean value because adding together the first and last equation and then separately adding the second and third equation both yield the single equation A + B + D + E = 20 In this simple type of problem dismissing the matrix approach maybe seems okay but it is no contender to matrix methods when trying to obtain speedy results! ``` Regards Chris Dean Chuck Senior Member Posts: 320 Threads: 59 Joined: Dec 2006 03-28-2008, 07:06 PM 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. ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 03-28-2008, 09:49 PM 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. Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 03-28-2008, 09:55 PM Quote: 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?? 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. ▼ Chuck Senior Member Posts: 320 Threads: 59 Joined: Dec 2006 03-29-2008, 03:18 PM 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 << -> A << A SVD 10 RND 'V' STO 1 V SIZE EVAL FOR I V I GET DUP IF 0 =/= the "=/=" is not equal THEN INV END V SWAP I SWAP PUT 'V' STO NEXT V A SIZE DIAG-> SWAP * * TRAN 10 RND >> >> ``` 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. ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 03-29-2008, 05:53 PM 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...". ▼ Chuck Senior Member Posts: 320 Threads: 59 Joined: Dec 2006 03-30-2008, 12:41 AM 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 ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 03-30-2008, 04:46 AM 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". ▼ Reth Senior Member Posts: 556 Threads: 9 Joined: Jul 2007 03-30-2008, 05:04 AM I wonder how many people would be interested in the topic (worldwide) and can HP48 express the ratio with decent precision :) cheers, Reth Edited: 30 Mar 2008, 6:48 a.m. Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 03-31-2008, 08:06 AM I disagree! I'm very interested as well, as always when you post something, Rodger. I have no idea how the 48/49 have implemented LSQ, but I'll find out. Cheers, Werner ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 03-31-2008, 08:48 AM 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. ▼ Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 04-02-2008, 04:35 AM Hi Rodger, no I responded to your message #8 to Chuck. I looked at the code, it performs a QR decomposition to turn A into an upper trapezoidal R, and then determines the rank before postmultiplying by another orthogonal matrix U to obtain a triangular matrix, then solves the triangular system and pre-multiplies the result by U^H. It's in the rank determination that things go awry, your original matrix is diagnosed as rank four, while the same matrix with the first two columns reversed is considered (correctly) as rank three. It'll take me a while to pinpoint the error - it's quite a complicated part of the code with incremental estimates of the min and max singular values.. and not in Fortran, but in SysRPL. Cheers, Werner Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 04-02-2008, 02:36 PM 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. With A and B as in your original post, the successive estimates for max and min singular values are: Smin Smax 1x1 1.41421356237309 1.41421356237309 2x2 1.41421356237309 1.41421356237309 3x3 0.84807051216015 1.96630585393211 4x4 3.01890419166484E-14 1.9873061521461 And unfortunately the negligibility test is Smax*1e-14 >= Smin which is false in all four cases - so the rank is determined to be four.. With the first two columns swapped, we get: min max 1x1 1.41421356237309 1.41421356237309 2x2 1.41421356237309 1.41421356237309 3x3 0.84807051216015 1.96630585393211 4x4 0.50000000000000E-14 1.9873061521461 which just passes the negligibility test, and we get rank three.. Perhaps using 1e-14 as tolerance was cutting it a bit close, as internally only 15 digits are used - while double precision Fortran is equivalent to almost 16 decimal digits (many of the routines read as if they were based on linpack or lapack, and I *think* they use a 1e-14 tolerance in there as well). And a single digit can be lost pretty easily. ``` I'll go on investigating the differences in the two cases anyway. Cheers, Werner ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 04-02-2008, 05:57 PM 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. Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 04-02-2008, 09:37 PM Quote: ```With A and B as in your original post, the successive estimates for max and min singular values are: Smin Smax 1x1 1.41421356237309 1.41421356237309 2x2 1.41421356237309 1.41421356237309 3x3 0.84807051216015 1.96630585393211 4x4 3.01890419166484E-14 1.9873061521461 And unfortunately the negligibility test is Smax*1e-14 >= Smin which is false in all four cases - so the rank is determined to be four.. ``` 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. ▼ Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 04-03-2008, 02:50 AM Hi Rodger, they represent the submatrix the estimates are based on. The rank determination is based on successive estimates for the singular values, starting with the 1x1 top left submatrix of the upper trapezoidal R (from A*P=Q*R decomposition), then the 2x2 submatrix etc. Apparently, swapping the first two columns in A changes P, Q and R, and not just P. Cheers, Werner Werner Member Posts: 163 Threads: 7 Joined: Jul 2007 04-03-2008, 06:55 AM 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 ▼ Rodger Rosenbaum Senior Member Posts: 305 Threads: 17 Joined: Jun 2007 04-03-2008, 08:50 AM Werner, would you mind sending me a private email with an email address I can use to communicate with you?

 Possibly Related Threads... Thread Author Replies Views Last Post HPCC Mini Conference 2013 hugh steers 6 441 09-13-2013, 04:27 PM Last Post: BruceH Picture from the Mini-HHC (LA) Geir Isene 11 664 07-07-2013, 01:06 PM Last Post: Jim Horn HP 41 Mcode related Questions Michael Fehlhammer 4 348 05-10-2013, 07:09 PM Last Post: Michael Fehlhammer My birthday, so a little commemorative mini-challenge ! Valentin Albillo 43 2,035 03-07-2013, 03:44 AM Last Post: Walter B WP 34S mini-challenge (B) Gerson W. Barbosa 17 838 12-27-2012, 04:39 PM Last Post: Marcus von Cube, Germany Mini-challenge: HHC2012 RPL programming contest with larger input David Hayden 14 824 10-05-2012, 10:36 PM Last Post: David Hayden Calculator O.T., instead HP related. Luiz C. Vieira (Brazil) 1 209 08-17-2012, 06:57 AM Last Post: Frank Boehm (Germany) HP41 / SY41CL Mini-B USB Power Connector (Module) Matt Kernal 5 1,100 07-08-2012, 06:23 PM Last Post: Diego Diaz OT and TAS related, but ... What Unisonic Models are shown in this auction picture? Gene Wright 6 467 07-04-2012, 12:53 PM Last Post: Guenter Schink HP65 related: I know it belongs in the adds, but I need advice. Luiz C. Vieira (Brazil) 2 218 04-07-2012, 07:31 AM Last Post: Luiz C. Vieira (Brazil)

Forum Jump: