▼
Posts: 305
Threads: 17
Joined: Jun 2007
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 minichallenge 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.
▼
Posts: 320
Threads: 59
Joined: Dec 2006
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
Posts: 120
Threads: 9
Joined: Aug 2005
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
Posts: 320
Threads: 59
Joined: Dec 2006
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 MoorePenrose inverse. This type of inverse exists even for singular matrices, and also gives the leastsquares solution. Now, how to get the HP50 to find this inverse??
Edited: 28 Mar 2008, 7:26 p.m.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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.
Posts: 305
Threads: 17
Joined: Jun 2007
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 MoorePenrose inverse. This type of inverse exists even for singular matrices, and also gives the leastsquares 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.
▼
Posts: 320
Threads: 59
Joined: Dec 2006
Thanks Roger. With your suggestions, I wrote a program for the HP50g which calculates the pseudoinverse 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 PseudoInverse 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.
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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...".
▼
Posts: 320
Threads: 59
Joined: Dec 2006
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
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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".
▼
Posts: 556
Threads: 9
Joined: Jul 2007
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.
Posts: 163
Threads: 7
Joined: Jul 2007
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
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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.
▼
Posts: 163
Threads: 7
Joined: Jul 2007
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 premultiplies 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
Posts: 163
Threads: 7
Joined: Jul 2007
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.01890419166484E14 1.9873061521461
And unfortunately the negligibility test is
Smax*1e14 >= 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.50000000000000E14 1.9873061521461
which just passes the negligibility test, and we get rank three..
Perhaps using 1e14 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 1e14 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
▼
Posts: 305
Threads: 17
Joined: Jun 2007
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.
Posts: 305
Threads: 17
Joined: Jun 2007
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.01890419166484E14 1.9873061521461
And unfortunately the negligibility test is
Smax*1e14 >= 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.74289198711E15 ]
It's strange that a problem with only small integers as input data would run into such difficulties.
▼
Posts: 163
Threads: 7
Joined: Jul 2007
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
Posts: 163
Threads: 7
Joined: Jul 2007
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 1e14 as tolerance is the issue here. Since we're only working with *estimates* for the singular values, that tolerance is a bit too strict. 1e13 or 5e12 would've been more appropriate. And alas, we have no way of providing that ourselves.
Cheers, Werner
▼
Posts: 305
Threads: 17
Joined: Jun 2007
Werner, would you mind sending me a private email with an email address I can use to communicate with you?
