The Museum of HP Calculators

HP Forum Archive 18

[ Return to Index | Top of Index ]

A bug and a related mini-challenge
Message #1 Posted by Rodger Rosenbaum on 26 Mar 2008, 9:19 p.m.

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.

      
Re: A bug and a related mini-challenge
Message #2 Posted by Chuck on 27 Mar 2008, 1:04 a.m.,
in response to message #1 by Rodger Rosenbaum

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
Message #3 Posted by Chris Dean on 28 Mar 2008, 12:51 p.m.,
in response to message #1 by Rodger Rosenbaum

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

      
Re: A bug and a related mini-challenge
Message #4 Posted by Chuck on 28 Mar 2008, 7:06 p.m.,
in response to message #1 by Rodger Rosenbaum

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
Message #5 Posted by Rodger Rosenbaum on 28 Mar 2008, 9:49 p.m.,
in response to message #4 by Chuck

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
Message #6 Posted by Rodger Rosenbaum on 28 Mar 2008, 9:55 p.m.,
in response to message #4 by Chuck

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.

                  
Re: A bug and a related mini-challenge (solution)
Message #7 Posted by Chuck on 29 Mar 2008, 3:18 p.m.,
in response to message #6 by Rodger Rosenbaum

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.

                        
Re: A bug and a related mini-challenge (solution)
Message #8 Posted by Rodger Rosenbaum on 29 Mar 2008, 5:53 p.m.,
in response to message #7 by Chuck

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)
Message #9 Posted by Chuck on 30 Mar 2008, 12:41 a.m.,
in response to message #8 by Rodger Rosenbaum

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)
Message #10 Posted by Rodger Rosenbaum on 30 Mar 2008, 4:46 a.m.,
in response to message #9 by Chuck

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)
Message #11 Posted by reth on 30 Mar 2008, 5:04 a.m.,
in response to message #10 by Rodger Rosenbaum

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.

                              
Re: A bug and a related mini-challenge (solution)
Message #12 Posted by Werner on 31 Mar 2008, 8:06 a.m.,
in response to message #8 by Rodger Rosenbaum

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

                                    
Re: A bug and a related mini-challenge (solution)
Message #13 Posted by Rodger Rosenbaum on 31 Mar 2008, 8:48 a.m.,
in response to message #12 by Werner

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)
Message #14 Posted by Werner on 2 Apr 2008, 4:35 a.m.,
in response to message #13 by Rodger Rosenbaum

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

                                    
Re: A bug and a related mini-challenge (solution)
Message #15 Posted by Werner on 2 Apr 2008, 2:36 p.m.,
in response to message #12 by Werner

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

                                          
Re: A bug and a related mini-challenge (solution)
Message #16 Posted by Rodger Rosenbaum on 2 Apr 2008, 5:57 p.m.,
in response to message #15 by Werner

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)
Message #17 Posted by Rodger Rosenbaum on 2 Apr 2008, 9:37 p.m.,
in response to message #15 by Werner

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.

                                                
Re: A bug and a related mini-challenge (solution)
Message #18 Posted by Werner on 3 Apr 2008, 2:50 a.m.,
in response to message #17 by Rodger Rosenbaum

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

                                          
Re: A bug and a related mini-challenge (solution)
Message #19 Posted by Werner on 3 Apr 2008, 6:55 a.m.,
in response to message #15 by Werner

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)
Message #20 Posted by Rodger Rosenbaum on 3 Apr 2008, 8:50 a.m.,
in response to message #19 by Werner

Werner, would you mind sending me a private email with an email address I can use to communicate with you?


[ Return to Index | Top of Index ]

Go back to the main exhibit hall