48/50 Linear System Solver

[Tim Wessman]

Have someone with a complaint about the linear system solver on the 48 series (50g specifically).

x - 5y + 4z = -3
2x - 7y + 3z = -2
-2x + y + 7z =-1

Now he knows these are dependent, but the concern is that instead of giving "no solutions" it returns

x = -0.04869...
y = 0.25937...
z = -0.221295...

In looking at the old 48 manual, it talks about the linear solver in 18-11 where it talks about 3 different cases.

[attachment=294]

I am curious as to where exactly these numbers are coming from. Anyone know which case, or what they are supposed to represent? Or is this an artifact of some lower level matrix routines?

The output of the LINSOLVE command is interesting, where it returns a '1=0' as the Z solution. :-)

[Manolo Sobrino]

Yes, this happens too with emu48 for the 48g. I guess it's the LU decomposition. You can try to decompose the coefficient matrix with the 50g and check by yourself that there is a nasty 1E-14 for the (3,3) element in the matrix L. It should be zero as it comes from a singular matrix (flag 54 doesn't seem to help BTW). The linear solver should check first that the matrix is invertible IMO, at least for small dimension matrices and/or integer coefficients, these potential side effects can cause quite a bit of trouble.

[Dieter]

(02-13-2014 01:30 AM)Tim Wessman Wrote:  Have someone with a complaint about the linear system solver on the 48 series (50g specifically).

The system is exactly determined (three variables in three equations), but its determinant is zero. Which is a quite challenging case for iterative numerical algorithms like a LU-decomposition: due to accumulated roundoff errors the result may be close to, but not exactly zero. Which makes the calculator think a valid solution can be calculated, and so it continues.. with a meaningless result.

I do not have a 48-series calculator here. What does it return if you try to evaluate the determinant? Something very close to (but not equal to) zero?

BTW: that's why I like simple direct methods like Cramer's rule for small linear systems with 2, 3 or maybe even 4 variables. ;-)

Dieter

[Marcus von Cube]

(02-16-2014 06:00 PM)Dieter Wrote:  The system is exactly determined (three variables in three equations), but its determinant is zero.

If I try to solve the system manually I get 0 = 5 which is obviously false irrespective of the values of x, y or z. So the system has no solution. A determinant of zero can say "no solution" or "a linear set of solutions" (depending on the right hand side) but it certainly does not say "exactly one solution".

Any calculator should error out on this equation. If I understand the excerpt from the 48G manual correctly, the result of LINSOLV occasionally is a least squares approximation. This is OK if the coefficients come from real world data and a solution is expected to exist but does not because of limited accuracy of the input.

[Manolo Sobrino]

To put it bluntly, the problem is that the HP 48/50 can find easily that the determinant is zero, yet the linear solver just calls the LU decomposition routine without checking the coefficient matrix first.

[Thomas Klemm]

(02-16-2014 06:00 PM)Dieter Wrote:  What does it return if you try to evaluate the determinant? Something very close to (but not equal to) zero?

Exactly 0. The RANK of the matrix gives 2.
When flag -22 is cleared (Infinite -> error) I get an error (Infinite Result) when I use / or 1/x to solve the equation. However I get Tim's result when the flag -22 is set or when I use LSQ or the linear solver. But when you calculate the residual with RSD you will notice that the computed solution isn't close to an actual solution. Furthermore the condition number is huge: 9.8e15. The approximate number of accurate digits is thus 0. When the flag -22 is cleared I get the same error: Infinite Result.

We get a slightly different result when the last equation is replaced by 4*Ⅰ-3*Ⅱ:
-2x + y + 7z =-6

Now we get a real solution which can be checked with the residual.

Cheers
Thomas

[Thomas Klemm]

(02-16-2014 06:49 PM)Marcus von Cube Wrote:  

(02-16-2014 06:00 PM)Dieter Wrote:  The system is exactly determined (three variables in three equations), but its determinant is zero.

A determinant of zero can say "no solution" or "a linear set of solutions" (depending on the right hand side) but it certainly does not say "exactly one solution".

Exactly! Who said so?

[gestrickland]

I don't consider that the calculator is giving a "meaningless result", although it is certainly the case that the offered result does not satisfy the equation with zero or negligible residuals.

To start with a simple illustration, consider the two-dimensional case with the equations x - y = 1 and x - y = -1 . If you put this system into the linear solver you get the "solution" of x = 0, y = 0 . Obviously this does not solve the equations in the usual sense, but it, and any other point midway between the two lines defined by the equations, will give smaller residuals than other points. Among all these midway points, the calculator chooses the one with the smallest distance from the origin.

In the three dimensional case differing pathologies are possible. In Tim's case it seems that the three planes defined by the three equations intersect in three parallel distinct lines. Since these never meet in a single point, there is no exact solution. In order to have something a bit easier to think about, let us consider an analogous case where the three equations are y + z = 1, y = 0, z = 0 . Here we have three planes, and their three lines of intersection are parallel. If we put this system into the linear solver, we get the "solution" x = 0, y = 0.3333..., z = 0.3333... . Again, this is not a solution in the usual sense, but it, and any other point with the same y and z values will give smaller residuals than all other points. Again, the calculator offers the one closest to the origin.

I agree that it would be nice if the calculator gave warning for this kind of a situation. I don't know whether there might be cases where this type of "solution" might be of practical value, but it is not meaningless.

[Werner]

There are a few issues here..
for me, my flag settings are as follows:
-20 Clear (underflow -> 0)
-21 Set (Overflow -> error)
-22 Clear (Infinite -> error)
Then there's another flag that influences the outcome: -54 Use tiny element.

With flag -54 Clear (Tiny element ->0), the 'stack solve' (/) returns Infinite Result, and the interactive Linear Solver will return the minimum norm Least Squares solution, because the system is underdetermined.
In this mode, the determinant is returned as 0 exactly, because there's a check to see if the matrix contains integers only, and with flag -54 clear, the calculator knows the determinant must be integer and will round the result to the nearest integer.
While solving systems, elements less than 1e-14 (relative size) will be set to zero.

With flag -54 Set (use tiny element), both methods will return

[ -5.41666..e14 -2.0833..e14 -1.25..e14]

The determinant now is calculated as -1.2e-13, due to roundoff. So the system is considered NOT singular, and a result can be produced.
In any case, it pays to compute the condition number (COND) that returns 9.8e15.
As a rule of thumb, the exponent signifies the number of wrong digits in your result, indicating here that the matrix is singular to working precision.

Werner


In any case

[Manolo Sobrino]

So that was the minimum norm solution! I assumed it was just garbage coming from the LU. I should have calculated it through (don't really do dynamical systems... didn't see that coming.)

I'm used to 20 clear, 21 clear, 22 check. With flag 54 checked I get for the determinant +1.2E-13, and the same solution as Werner. With flag 54 cleared, the determinant is zero and I get Tim's solution.

Looks like the linear solver is a powerful beast indeed, although giving that kind of answer without a warning is kind of mean. If I asked for a Coke, I'm not expecting a Pepsi.

(Please, don't "the user should know better" me... I'm probably using a calculator to solve a 3x3 system because I'm in a hurry, if the system is singular just give it to me straight.)

[Han]

This was a question that I stumbled upon while working on the solver an equation library program. Solve for x in Ax=b if:

Let A = [[12.25, 1.25, 4.5],[1.25,12.25,4.5],[4.5,4.5,3]]
Let b = [[-166.2],[157.2],[-3]]

Note that A is singular. On the HP48, with flag -22 SET (infinite -> MAXREAL), I get:

b A / -> x=[[-25.0333...],[4.3666...],[30]]
b A LSQ -> x=[[-14.972727...],[14.42727...],[-.1818...]]

The b A LSQ method gives a solution with smaller norm than the b A / method. I am able to compute the same solution using SVD decomposition on A and applying the pseudo-inverse.

Code:

<<
A SVD
OBJ-> DROP DROP @ get nonzero eigenvalues
{} + + INV @ turn into a list and compute reciprocals
0 + EVAL ->V3 { 3 3 } DIAG-> @reconstruct inverse 3x3 diagonal matrix
SWAP TRN SWAP * SWAP TRN * @ compute the pseudo inverse
b *
>>

However, I am not sure how the HP48 got the answer for b A / ... Comparing the norms of the residuals, b A / gives a smaller norm for Ax-b, but clearly the norm of the solution is larger than that from b A LSQ. So what algorithm does b A / apply? (Would really like to not have to break out Jazz and decompile some ROM routines just to figure out the algorithm...)

EDIT: Additionally, RSD according the the advanced user's guide is supposed to compute Ax-b. However, I seem to get different answers when I manually compute Ax-b as compared to using RSD. Is this the case for anyone else? What am I overlooking?

[Gilles]

By the way, in exact mode with the 50G :

Code:


[
 'x - 5*y + 4*z = -3'
 '2*x - 7*y + 3*z = -2'
 '-2*x + y + 7*z =-1'
]

['x' 'y' 'z']

SOLVE

returns

{}

I sometimes used SOLVE instead of LINSOLVE. As far I know It's not documented, but works in most cases even with things like :

[
'x^2- 5*y - 4*z = 3'
'2*x - 7*y + 3*z = -2'
'-2*x + y + 7*z =-1'
]

['x' 'y' 'z']

SOLVE

{ [ 'x=-((-37+2*√911)/26)' 'y=-((-513+20*√911)/676)' 'z=-((-105+12*√911)/676)' ] [ 'x=(37+2*√911)/26' 'y=(513+20*√911)/676' 'z=(105+12*√911)/676' ] }

[Werner]

(01-14-2015 05:35 PM)Han Wrote:  Let A = [[12.25, 1.25, 4.5],[1.25,12.25,4.5],[4.5,4.5,3]]
Let b = [[-166.2],[157.2],[-3]]

Note that A is singular. On the HP48, with flag -22 SET (infinite -> MAXREAL), I get:

b A / -> x=[[-25.0333...],[4.3666...],[30]]
b A LSQ -> x=[[-14.972727...],[14.42727...],[-.1818...]]

As I said before:
with flag -54 SET, both b A / and the linear solver application return x=[[-25.0333...],[4.3666...],[30]]. In this case, the 48 does not know A is singular - the computed determinant is 1.485e-12, and a result can be calculated using the normal LU decomposition and subsequent solve. Of course, COND returns 3.e15 meaning the result is meaningless.
with flag -54 CLEAR, b A / returns Infinite Result, and the solver app returns the Least Squares solution.

Cheers, Werner

[Han]

(01-15-2015 12:17 PM)Werner Wrote:  As I said before:
with flag -54 SET, both b A / and the linear solver application return x=[[-25.0333...],[4.3666...],[30]]. In this case, the 48 does not know A is singular - the computed determinant is 1.485e-12, and a result can be calculated using the normal LU decomposition and subsequent solve. Of course, COND returns 3.e15 meaning the result is meaningless.
with flag -54 CLEAR, b A / returns Infinite Result, and the solver app returns the Least Squares solution.

Cheers, Werner

Thank you, Werner. I don't know how I glossed over your original post and didn't see your comments -- probably had tunnel vision and too focused on something else. Your explanation is very helpful.

Edit: So I did the LU decomposition of the matrix, and got a different solution: [[84.9666...],[114.3666...],[-300]]. In original example posted by Tim:

Quote:With flag -54 Set (use tiny element), both methods will return

[ -5.41666..e14 -2.0833..e14 -1.25..e14]

Using LU, I got the same result. However, for my particular example, LU doesn't give the same answer as the / operator (under the same flag settings) which suggests that the HP48 doesn't seem to always use LU.

Just to be clear, I'm not so much after meaningful solutions as I am about understanding how the HP48 solves linear systems using the / operator. So far I've tried QR, LU, and pseudo-inverse (least squares) and none match. QR gives [ 139.27, -181.42, 5.43 ], LU gives something with large norm (useless answer), and least squares I already listed.

[Werner]

(01-15-2015 05:35 PM)Han Wrote:  Edit: So I did the LU decomposition of the matrix, and got a different solution: [[84.9666...],[114.3666...],[-300]].

How did you arrive at this?
And the '/' operator does LU-decomposition.

Werner

[Han]

(01-16-2015 09:52 AM)Werner Wrote:  

(01-15-2015 05:35 PM)Han Wrote:  Edit: So I did the LU decomposition of the matrix, and got a different solution: [[84.9666...],[114.3666...],[-300]].

How did you arrive at this?
And the '/' operator does LU-decomposition.

Werner

My stack was:

2: [[-166.2],[157.2],[-3]]
1: [[12.25,1.25,4.5],[1.25,12.25,4.5],[4.5,4.5,3]]

I typed in the following in the command line; remove the comments after the @'s

Code:

LU @ b L U P
 4 ROLL * @ L U P*b
ROT INV SWAP * @ solve Ly=P * b by y=L^-1 * P *b
SWAP INV SWAP * @ solve Ux = y by x = U^-1 y

[Werner]

Ah OK.
Two remarks here:
- the condition number of the matrix is 3e15, meaning that very small perturbations in the input lead to huge changes in the calculated solution - exactly what we see here.
- And then, to solve Ly = Pb or Uy=x, simply use / instead of INV SWAP *. Faster and more accurate, but of course not in this case.
(now you get [-1621.7 -1592.3 4820 ], equally meaningless, but that's what the condition number tells you)

Calculating Ax=b with / vs with LU and then manually doing the triangular solves makes a difference: the former uses 15 digits throughout, while the second has the intermediate amounts that make up L and U rounded to 12 digits. And that small difference is enough to completely change the calculated solution.

Werner

[Han]

(01-16-2015 02:15 PM)Werner Wrote:  Ah OK.
Two remarks here:
- the condition number of the matrix is 3e15, meaning that very small perturbations in the input lead to huge changes in the calculated solution - exactly what we see here.
- And then, to solve Ly = Pb or Uy=x, simply use / instead of INV SWAP *. Faster and more accurate, but of course not in this case.
(now you get [-1621.7 -1592.3 4820 ], equally meaningless, but that's what the condition number tells you)

Calculating Ax=b with / vs with LU and then manually doing the triangular solves makes a difference: the former uses 15 digits throughout, while the second has the intermediate amounts that make up L and U rounded to 12 digits. And that small difference is enough to completely change the calculated solution.

Werner

Aha! Thank you so much for your patience, Werner!