48/50 Linear System Solver

02132014, 01:30 AM
(This post was last modified: 02132014 01:34 AM by Tim Wessman.)
Post: #1




48/50 Linear System Solver
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 1811 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. :) TW Although I work for HP, the views and opinions I post here are my own. 

02132014, 04:01 PM
Post: #2




RE: 48/50 Linear System Solver
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 1E14 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.


02162014, 06:00 PM
(This post was last modified: 02162014 06:01 PM by Dieter.)
Post: #3




RE: 48/50 Linear System Solver
(02132014 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 LUdecomposition: 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 48series 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 

02162014, 06:49 PM
Post: #4




RE: 48/50 Linear System Solver
(02162014 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. Marcus von Cube Wehrheim, Germany http://www.mvcsys.de http://wp34s.sf.net http://mvcsys.de/doc/basiccompare.html 

02162014, 06:54 PM
Post: #5




RE: 48/50 Linear System Solver
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.


02162014, 11:40 PM
Post: #6




RE: 48/50 Linear System Solver
(02162014 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 

02162014, 11:54 PM
Post: #7




RE: 48/50 Linear System Solver
(02162014 06:49 PM)Marcus von Cube Wrote:Exactly! Who said so?(02162014 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". 

02172014, 07:19 AM
Post: #8




RE: 48/50 Linear System Solver
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 twodimensional 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. 

02172014, 08:18 AM
Post: #9




RE: 48/50 Linear System Solver
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 1e14 (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.2e13, 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 

02172014, 02:32 PM
Post: #10




RE: 48/50 Linear System Solver
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.2E13, 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.) 

01142015, 05:35 PM
(This post was last modified: 01142015 06:53 PM by Han.)
Post: #11




RE: 48/50 Linear System Solver
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 pseudoinverse. Code: << 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 Axb, 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 Axb. However, I seem to get different answers when I manually compute Axb as compared to using RSD. Is this the case for anyone else? What am I overlooking? Graph 3D  QPI  SolveSys 

01142015, 07:59 PM
(This post was last modified: 01152015 11:17 AM by Gilles.)
Post: #12




RE: 48/50 Linear System Solver
By the way, in exact mode with the 50G :
Code:
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' ] } 

01152015, 12:17 PM
Post: #13




RE: 48/50 Linear System Solver
(01142015 05:35 PM)Han Wrote: Let A = [[12.25, 1.25, 4.5],[1.25,12.25,4.5],[4.5,4.5,3]] 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.485e12, 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 

01152015, 05:35 PM
(This post was last modified: 01152015 06:32 PM by Han.)
Post: #14




RE: 48/50 Linear System Solver
(01152015 12:17 PM)Werner Wrote: As I said before: 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 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 pseudoinverse (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. Graph 3D  QPI  SolveSys 

01162015, 09:52 AM
Post: #15




RE: 48/50 Linear System Solver  
01162015, 01:44 PM
Post: #16




RE: 48/50 Linear System Solver
(01162015 09:52 AM)Werner Wrote:(01152015 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]]. 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 Graph 3D  QPI  SolveSys 

01162015, 02:15 PM
Post: #17




RE: 48/50 Linear System Solver
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 

01162015, 04:18 PM
Post: #18




RE: 48/50 Linear System Solver
(01162015 02:15 PM)Werner Wrote: Ah OK. Aha! Thank you so much for your patience, Werner! Graph 3D  QPI  SolveSys 

« Next Oldest  Next Newest »

User(s) browsing this thread: 1 Guest(s)