A bug and a related mini-challenge



#21

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.


#22

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

#23

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

#24

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.


#25

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.

#26

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.


#27

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.


#28

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...".


#29

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


#30

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".


#31

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.

#32

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


#33

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.


#34

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

#35

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


#36

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.

#37

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.


#38

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

#39

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


#40

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 292 09-13-2013, 04:27 PM
Last Post: BruceH
  Picture from the Mini-HHC (LA) Geir Isene 11 435 07-07-2013, 01:06 PM
Last Post: Jim Horn
  HP 41 Mcode related Questions Michael Fehlhammer 4 234 05-10-2013, 07:09 PM
Last Post: Michael Fehlhammer
  My birthday, so a little commemorative mini-challenge ! Valentin Albillo 43 1,324 03-07-2013, 03:44 AM
Last Post: Walter B
  WP 34S mini-challenge (B) Gerson W. Barbosa 17 577 12-27-2012, 04:39 PM
Last Post: Marcus von Cube, Germany
  Mini-challenge: HHC2012 RPL programming contest with larger input David Hayden 14 559 10-05-2012, 10:36 PM
Last Post: David Hayden
  Calculator O.T., instead HP related. Luiz C. Vieira (Brazil) 1 146 08-17-2012, 06:57 AM
Last Post: Frank Boehm (Germany)
  HP41 / SY41CL Mini-B USB Power Connector (Module) Matt Kernal 5 809 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 318 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 143 04-07-2012, 07:31 AM
Last Post: Luiz C. Vieira (Brazil)

Forum Jump: