[42S] iterative refinement

[Werner]

0.Abstract

The 71B performed an iterative refinement step as part of its MAT SYS routine,
thus delivering more accurate results than the 42S, while deploying the same
processor and precision. The routine presented here will perform the same thing
on the 42S (and EMU42). While it will run on Free42, it will have no (or little)
effect, as Free42 does not use higher precision in matrix products. The 42S
uses 15 digits intermediate precision to evaluate dot products, (as does the
71B) and that we can use to improve the result obtained by b A / (the way to
solve a system A*x=b on the stack).


1.Iterative refinement

Solving A*x = b does not give you x but x', an approximation of the true x.
Let x = x' + dx
We'll try and determine dx:

A*(x' + dx) = b
A*dx = b - A*x'
dx = inv(A)*(b-A*x')

The multiplication by the inverse is just a shorthand way of saying 'solve'.
We need to be able to calculate b - A*x' with greater accuracy - like the RSD
function in the 48. However, the 42S does not have an RSD function and will
round the result of A*x' to 12 digits before subtracting it from b - which
defeats the purpose as the two will be very close and cancellation will occur.
If however we place b as an extra column in A and add a -1 to x':

Code:

  b1 a11 a12 a13 a14   -1
  b2 a21 a22 a23 a24 * x1
  b3 a31 a32 a33 a34   x2
  b4 a41 a42 a43 a44   x3
                       x4

so

Code:

 x := x' - inv(A)*( [ b A ] * [ -1  ] )
                              [  x' ]

We can then calculate A*x'-b in a single operation, and make use of the 42S's
intermediate 15-digits. Adding a column can be done adding a row to the transpose with INSR.


2.The Program

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

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

3.Example

A small matrix with a high condition number is:

[[ 12969 8648 ]
[ 2161 1441 ]]

Determinant=1, cond. number is about 2.5e8

exact inverse is

[[ 1441 -8648 ]
[ -2161 12969 ]]

calculated inverse on a 48G using 15 digits is
[[ 1440.99995021 -8647.99970121 ]
[ -2160.99992534 12968.9995519 ]]

on a 42S with 12 digits:
[[ 1440.98630787 -8647.91782819 ]
[ -2160.97946655 12968.8767708 ]]

We lost 8 digits, and have 12-8=4 digits left

after 1 refinement (calculated by column):
[[ 1441.0000431 -8647.99880302 ]
[ -2161.00006463 12968.9982049 ]]

after 2:
[[ 1440.99994236 -8648.00008641 ]
[ -2160.99991356 12969.0001296 ]]

The first column did not improve further, but the second did, and the result is
now about as accurate as the one from the 48G. Since we evaluate Ax-b with
15-digit accuracy, and the condition number tells us we will lose 8 digits, we
expect no further improvement. Usually a single refinement suffices.

Werner