Matrix Division on the 28C/S (and later)

[edryer]

There is a way to obtain the (approximate) result of a linear system of equations by using matrix transposition of the coefficients, multiplication and finally division... this is a quick handy way to perform this, and being new to Linear Algebra I assumed Matrix division was valid!

From the 28C/S Vectors and Matrices guide (see attachment) the formula (which was new to me), where A is the coefficient Matrix and B the results vector.

What exactly is the 28S doing in solving this system? (Which it does....)


Off topic:
Again, I have found yet ANOTHER way to enter matrices on the RPL machines....
3x3 for example:
1,2,3,4,5,6,7,8,9{3,3 ->ARR (via Softkey)

I have never found methods as fast as the RPL machines for entering and using Matrices, so many ways that all perform the same result! They put a LOT of thought into their work!

[John Keith]

AFAIK, the "/" key essentially inverts the matrix containing the equations then multiplies the inverted matrix by the array of values. The HP-48 manual does say that it uses 15 digit internal precision and is more accurate than inversion followed by multiplication.

[Werner]

No, it doesn't. It solves the system A*X=B, returning X, without having to calculate the inverse.
While mathematically solving A*X=B and calculating X=inv(A)*B are equivalent, numerically they are not.

Cheers, Werner

[edryer]

thank you both for your kind replies.

I have found reading thoroughly the 28C/S Manuals (all three!) and my 48G user manual helps me discover new things in Mathematics, the 28C/S Solutions Manuals are also invaluable!

I prefer the 28S over my old 48G and find it to be an excellent platform for such adventures being new to College mathematics in my forties, this was the dream machine of my youth so it seems fitting))

[J-F Garnier]

(02-25-2020 08:16 PM)John Keith Wrote:  AFAIK, the "/" key essentially inverts the matrix containing the equations then multiplies the inverted matrix by the array of values. The HP-48 manual does say that it uses 15 digit internal precision and is more accurate than inversion followed by multiplication.

(02-26-2020 08:44 AM)Werner Wrote:  No, it doesn't. It solves the system A*X=B, returning X, without having to calculate the inverse.
While mathematically solving A*X=B and calculating X=inv(A)*B are equivalent, numerically they are not.

Hi Werner, can you elaborate on this?
I always thought (and still think) that solving A*X=B is done by inverting A. It may exist other numerical methods, iterative or so, but I don't think HP is using such methods.
So I'm interested to get some justification of your point.

An indication: on Basic machines (Series 80, 75, 71), the MAT ... SYS statement is using matrix inversion, according to the manuals, and provides the determinant for the DETL function (last DET computed in INV or SYS statements).

J-F

[Werner]

Solving A*X=B for square A is done factorizing A into P*L*U, with P a permutation matrix and L and U lower and upper triangular matrices, respectively, of which one has only 1's on the diagonal (it depends on the algorithm used which one has the 1's). L and U are kept in A, overwriting it.
Then, A*X=B is solved by solving, in order:
(1) P*Z=B, effectively applying the same permutations to B as had been done to A
(2) L*Y=Z
(3) U*X=Y

Solving A*X=B as X=inv(A)*B requires a number of additional steps:
- invert A from A=P*L*U
1.inverting L: R=inv(L)
2.solving U*C=R to obtain C=inv(U)*inv(L)
3.applying P from the right: C=C*P
- multiply C and B
these steps require a lot more calculations, and are somewhat less accurate as a result.

BTW SYS on the 71B solves the system as I set out above, but it performs an additional 'iterative refinement' step to obtain a more accurate result:
Solving A*x=b does not give you x but x', an approximation of the true x.
Let x = x' + dx
Let's try and determine dx:
A*(x' + dx) = b
A*dx = b - A*x'
If we can calculate c = b-A*x' with greater accuracy, we can determine dx, solving A*dx = c.
In the 71B, b-A*x' is determined with 15-digit arithmetic (it's just a number of DOT products).
Alas, the 42S does not do this, but I have written a routine that will perform an iterative refinement step. It is only useful in the real 42S and emu42, not in Free42 or the DM42, as these do not carry out dot products in higher precision (sadly). I already posted it somewhere here before, I think.

Code:

                L       X       Y       Z       T
00 { 58-Byte Prgm }
01 LBL "/IR"            A       b
02 X<>Y
03 ENTER
04 TRANS                bt      b       A
05 RCL ST Z             A       bt      b       A
06 STO/ ST Z            A       bt      x       A
07 TRANS
08 EDIT
09 INSR                 0       bt      x       A
10 Rv                   bt      x       A       .
11 PUTM
12 CLX
13 RCLEL
14 EXITALL
15 TRANS                bA      x       A
16 X<>Y
17 LBL 02               x       bA      A
18 STO ST L    x        x       bA      A
19 EDIT
20 INSR
21 -1
22 EXITALL     x        -1x     bA      A
23 X<> ST L  -1x        x       bA      A
24 RCL ST Y  -1x        bA      x       bA      A
25 RCLx ST L            Ax-b    x       bA      A
26 RCL/ ST T           -dx      x       bA      A
27 -                    x        bA        A
28 RTN
29 GTO 02
30 END

Provide b ENTER A, then run "/IR".
It will solve the system and perform a single refinement step, leaving x on the
stack. If you're careful, you can examine x with EDIT and the arrow functions,
leaving the stack undisturbed. Exit x, and press R/S to carry out another
refinement step.
It will only work with vectors b (matrices with a single column).
Usually, a single refinement step suffices.

Cheers, Werner

[J-F Garnier]

Thanks, Werner, for your explanations. this clarifies your point and I agree with this description. Matrix inversion and system solving share the same LU decomposition, then differ in the result, either inverted matrix or solution of the system.

I was aware of the refinements of the HP71 MAT ... SYS statement.
These refinements are not done in the RPL "matrix division" operation, but the 28S and later 48 series provide a RSD command to compute the residual.

J-F

[John Keith]

I'm surprised that the RPL machines do not do the iterative refinement that the HP-71 Math ROM does. I recall reading that most of the math routines in the RPL calculators were carried over from the HP-71.

[rprosperi]

(02-27-2020 06:35 PM)John Keith Wrote:  I'm surprised that the RPL machines do not do the iterative refinement that the HP-71 Math ROM does. I recall reading that most of the math routines in the RPL calculators were carried over from the HP-71.

I'd guess it's because on the RPL machines, matrices are their own new data type so the old algorithm/code probably would not apply, without as much re-work as doing it from scratch in the new format. But only guessing...

[J-F Garnier]

To go back to the OT:

(02-25-2020 03:26 PM)edryer Wrote:  From the 28C/S Vectors and Matrices guide (see attachment) the formula (which was new to me), where A is the coefficient Matrix and B the results vector.

The HP-71B Math ROM provides an efficient transpose-and-multiply statement for this purpose, so the equivalent sequence on the HP-71B is:
MAT B=TRN(A)*B
MAT A=TRN(A)*A
MAT X=SYS(A,B)

J-F

[Paul Dale]

This is the least squares solution that was ingrained into my memory during my statistics classes: \( (X^TX)^{-1}X^TY \)

[Thomas Okken]

(02-27-2020 07:02 AM)Werner Wrote:  Alas, the 42S does not do this, but I have written a routine that will perform an iterative refinement step. It is only useful in the real 42S and emu42, not in Free42 or the DM42, as these do not carry out dot products in higher precision (sadly). I already posted it somewhere here before, I think.

Werner, I have a two-week vacation coming up, and I might look into higher-precision DOT. I have the links from your old post on the topic, but could you also provide a pointer on how to create test cases for this?

[Werner]

Hi Thomas.
Do you want me to give you examples of DOT products, SYS solves etc that would show the effects of the increased accuracy?
Werner

[Thomas Okken]

(03-02-2020 12:26 PM)Werner Wrote:  Hi Thomas.
Do you want me to give you examples of DOT products, SYS solves etc that would show the effects of the increased accuracy?
Werner

Hi Werner,

Yes, the DOT products in particular. The references you cited look detailed enough that I should be able to write the code from them, but since I'm not too familiar with the numerical issues involved, I could use some help with test cases. (Of course I'd also be happy to provide a test build for you to play with once I get to that point, if you like.)

I haven't really thought about the SIMQ refinement technique, since it's not built-in 42S functionality, but it does sound like something I should look into as well, and maybe put up on my 42S programs page.

Thanks!

Thomas

[Werner]

No, what I meant was that since LU-decomp uses DOT products as well, some accuracy gains can be expected there as well, even without an iterative refinement step. Test cases for DOT products shouldn’t be too hard, for Sys solve (and MATINV) a bit harder.
Cheers, Werner

[Albert Chan]

A rule of thumb for LU decomposition method, "Numerical Methods in Engineering" by Salvadori and Baron, 1952, page 29

Quote:The number of figures lost in the computations varies from system to system, but statistically is of the order 0.3n.
Thus a system of 10 equations should be solved carrying 3 more significant figures than those required in the roots

The book is old, and above might not hold now. Example, this is particularly funny

Quote:A system of 10 equations can easily be solved in approximately 10 hours.

[Thomas Okken]

(03-02-2020 05:48 PM)Werner Wrote:  No, what I meant was that since LU-decomp uses DOT products as well, some accuracy gains can be expected there as well, even without an iterative refinement step. Test cases for DOT products shouldn’t be too hard, for Sys solve (and MATINV) a bit harder.

The LU decomposition in Free42 is not implemented using the DOT function, though, so having it benefit from a more accurate dot product algorithm would be a separate exercise. The same thing is true of all the other functions that perform dot products under the hood, like matrix-matrix multiplication.

[Werner]

okay, one thing at a time then.
I can give a hand converting the existing routines to DOT-based ones.

I think using the FMA to determine the correction term will allow results as if 68 digits were used throughout, no?
So

Code:

@ dot test 1: input a, result = DOT([a 1 -1],[1 1 1])
>LBL "DT1" 
 3
 1
 NEWMAT
 EDIT
 Rv
 \->
 1
 \->
 -1
 EXITALL
 ENTER
 SIGN
 ABS
 DOT
 END

should return 'a' then, for 1E-67<=a<1E-33
Currently, of course, Free42 returns 1E-33 for 5E-34<a<1E-33 and 0 for smaller a.
(incidentally, I notice that Free42 performs round-to-even even in dot products - something the 42S does not do, due to the 15 intermediate digits ;-)

Is this the kind of thing you would like to see? I'll think of some more involved examples.

Cheers, Werner

[Thomas Okken]

Yes, that's perfect, thanks! A simple test that I can throw at the new DOT, with expected results, that's what I'm looking for.

I'll write the new version of DOT next week, I think. I'll have to add FMA to the "phloat" functions, using bid128_fma() in the decimal version, and (I think) std::fma() in the binary version, and then implement the functions from your citations. It looks pretty straightforward.

[Werner]

Please allow your internal dot function to be able to access elements that are NOT adjacent (as you'll always have with user DOT), but that are a constant distance apart. That way, all other functions (DET, SOLVE(/), INVMAT, matrix multiply) can be rewritten making use of that.
And then we have the COMPLEX cases to consider as well..
Werner