Householder

# Householder Transformation for the HP-41

Overview

1°)  Householder Algorithm
2°)  A Modified Algorithm

-The 1st program uses the Householder method to triangularize one - or several - linear system(s)
-Then "TRS" is called to compute the solutions.

-In the 2nd program, the method is modified to diagonalize the principal matrix - if it is possible - so that the diagonal only contains "1"s
-No subroutine is needed in this case.

-In both cases, the elements are to be stored in column order, the right-hand side first.

B = A.X

1°)  Householder Algorithm

-The Householder method uses orthogonal symmetries to gradually triangularize a matrix.

-If ( a1 , ...... , an ) is the first column of the matrix A,  the other columns  ( x1 , ...... , xn )  are replaced by  ( y1 , ...... , yn )  defined as

( y1 , ...... , yn ) = ( x1 , ...... , xn ) - C [ a1 + sgn(a1) ( a12 + ....... + an2 )1/2 , a2 , ...... , an ]

where   C = [ SUMi=1,...,n ai xi + x1 sgn(a1) ( SUMi=1,...,n ai2 )1/2 ] / [ SUMi=1,...,n ai2 + a1 sgn(a1) ( SUMi=1,...,n ai2 )1/2 ]

-The first column becomes  ( A1 , 0 , .......... , 0 )   with   A1 = -sgn(a1) ( SUMi=1,...,n ai2 )1/2
-Then, the same process is repeated with the 2nd column, without taking into account the 1st column and the 1st raw, and so on...

-Actually, the program listed below starts with the last column and the last element, and then it continues with the next-to-last column and so on...

Data Registers:               R00 = Det A                               ( Registers R01 thru Rn.m are to be initialized before executing "HLS" )

•  R01 = b1           •  Rn+1 = a1,1      .....................      •  Rnm-n+1 = a1,n
•  R02 = b2           •  Rn+2 = a2,1      .....................      •  Rnm-n = a2,n

................            ...............................................................................

•  Rnn = bn           •  R2n = an,1         .....................      •  Rnm = an,n

>>>  At the end  x1 , .......... , xn   are in registers R01 to Rnn  and  the square matrix has become a lower triangular matrix.

•  Of course, this program may be used to solve several linear systems with the same "principal" matrix,
especially if you want to invert the square-matrix itself.

Flag:   F00

CF 00 = we solve the linear system(s) except if det A = 0
SF 00 = we stop when the array has only zeros above the "diagonal"

F00 is automatically set if m <= n  ( overdetermined systems )

Subroutine:    "TRS"  ( cf "Triangular Systems for the HP-41" )

-Lines 51 & 143 are three-byte GTOs
-Line 72 is a TEXT0 or any other NOP instruction ( STO X .... )

 01  LBL "HLS"    02  X<=Y?   03  SF 00   04  X<>Y   05  .1   06  %   07  +   08  STO N   09  ST* Y   10  FRC   11  -   12  STO O   13  1   14  STO 00   15  LASTX   16  %   17  RCL Z   18  INT   19  +   20  LBL 01   21  STO M   22  INT   23  RCL O   24  FRC   25  +   26   E-6   27  RCL O   28  INT   29  RCL N   30  INT   31  -   32  STO P 33  CLX   34  LBL 02   35  RCL IND Z   36  X^2   37  +   38  DSE Z   39  GTO 02   40  SQRT   41  RCL IND M   42  SIGN   43  *   44  STO Q   45  ABS   46  X<=Y?   47  CLX   48  X=0?   49  STO 00   50  X=0?   51  GTO 07   52  LBL 03   53  RCL P   54  X=0?   55  GTO 06   56  RCL O   57  INT   58  ST- P   59  RCL M   60  INT   61  ST+ P   62  RCL O   63  FRC   64  ST+ Y 65  CLX   66  LBL 04   67  RCL IND Y   68  RCL IND P   69  *   70  +   71  DSE P   72  ""   73  DSE Y   74  GTO 04   75  RCL Q   76  /   77  RCL N   78  RCL P   79  +   80  RCL M   81  +   82  INT   83  RCL O   84  INT   85  -   86  STO P   87  X<>Y   88  RCL IND Y   89  +   90  RCL IND M   91  STO Z   92  RCL Q   93  ST+ IND M   94  +   95  /   96  RCL M 97  INT   98  RCL O   99  FRC 100  + 101  X<>Y 102  SIGN 103  LBL 05 104  CLX 105  RCL IND Y 106  LASTX 107  * 108  ST- IND P 109  DSE P 110  CLX 111  DSE Y 112  GTO 05 113  R^ 114  STO IND M 115  GTO 03 116  LBL 06 117  RCL O 118  FRC 119  RCL M 120  INT 121  PI 122  INT 123  10^X 124  ST* Z 125  / 126  + 127  ISG X 128  CLRGX 129  RCL Q 130  ST* 00 131  CHS 132  STO IND M 133  LBL 07 134  RCL N 135  ST- O 136  RCL O 137  RCL M 138  1 139  - 140  X<=Y? 141  CLX 142  DSE X 143  GTO 01 144  RCL M 145  RCL N 146  INT 147  / 148  INT 149  RCL N 150  INT 151  RCL 00 152  CLA 153  FC? 00 154  X=0? 155  RTN 156  X<> Z 157  XEQ "TRS" 158  RCL 00 159  END

( 244 bytes / SIZE  1+n.m )

 STACK INPUTS OUTPUTS Y n / X m Det A

where  n = number of raws
and   m = number of columns

Example1:    The instructions are similar to those of "LS3"  ( cf "Linear & Non-Linear Systems for the HP-41" )

-3 = 2 x + 3 y - 4 z              R01   R04   R07   R10
21 = 4 x - 5 y + 7 z      =     R02   R05   R08   R11
38 = 4 x + 2 y + 6 z            R03   R06   R09   R12

CF 00
3   ENTER^
4   XEQ "HLS"   >>>>   Det =  -188                       ---Execution time = 26s---

R01 = x = 2              R04   R07   R10          -3.6778         0            0
R02 = y = 3    and    R05   R08   R11    =      1.8181    5.0864        0                   rounded to 4D
R03 = z = 4              R06   R09   R12          -4.3782    3.4826   -10.0499

R01 = -7.3556   R02 = 18.8953   R03 = -38.5079

Example2:    To invert the same matrix, store

R01   R04   R07   R10   R13   R16         1    0    0    2    3   -4
R02   R05   R08   R11   R14   R17   =    0    1    0    4   -5   7
R03   R06   R09   R12   R15   R18         0    0    1    4    2    6

CF 00
3   ENTER^
6   XEQ "HLS"   >>>>   Det =  -188                      ---Execution time = 47s---

R01   R04   R07         11/47    13/94    -1/188
R02   R05   R08   =    -1/47    -7/47      15/94          with small roundoff-errors as usual
R03   R06   R09         -7/47    -2/47      11/94

Note:

-With n = 7 & m = 8, Execution time = 2mn39s

2°)  A Modified Algorithm

-A similar method may be used to diagonalize a matrix:

-After using the formulae mentionned above, the last raw of the matrix is divided by Am So, the last column is now ( 0 , ........... , 0 , 1 )
-And we continue in the same way with the other complete columns.

Data Registers:               R00 = Det A                             ( Registers R01 thru Rn.m are to be initialized before executing "HLS" )

•  R01 = b1           •  Rn+1 = a1,1      .....................      •  Rnm-n+1 = a1,n
•  R02 = b2           •  Rn+2 = a2,1      .....................      •  Rnm-n = a2,n

................            ...............................................................................

•  Rnn = bn           •  R2n = an,1         .....................      •  Rnm = an,n

>>>  At the end  x1 , .......... , xn   are in registers R01 to Rnn  and  the square matrix - if it's not singular - has become a Unit matrix.

•  This program may also be used to solve several linear systems with the same "principal" matrix,
especially if you want to invert the square-matrix itself.

Flags: /
Subroutines: /

-Lines 49 & 143 are three-byte GTOs

 01  LBL "HLS"    02  X<>Y   03  .1   04  %   05  +   06  STO N   07  ST* Y   08  FRC   09  -   10  STO O   11  1   12  STO 00   13  LASTX   14  %   15  RCL Z   16  INT   17  +   18  LBL 01   19  STO M   20  INT   21  RCL O   22  FRC   23  +   24   E-6   25  RCL O   26  INT   27  RCL N   28  INT   29  -   30  STO P 31  CLX   32  LBL 02   33  RCL IND Z   34  X^2   35  +   36  DSE Z   37  GTO 02   38  SQRT   39  RCL IND M   40  SIGN   41  *   42  STO Q   43  ABS   44  X<=Y?   45  CLX   46  X=0?   47  STO 00   48  X=0?   49  GTO 08   50  LBL 03   51  RCL P   52  X=0?   53  GTO 06   54  RCL O   55  INT   56  ST- P   57  RCL M   58  INT   59  ST+ P   60  RCL O 61  FRC   62  ST+ Y   63  CLX   64  LBL 04   65  RCL IND Y   66  RCL IND P   67  *   68  +   69  DSE P   70  ""   71  DSE Y   72  GTO 04   73  RCL Q   74  /   75  RCL N   76  INT   77  ST+ P   78  CLX   79  RCL P   80  RCL M   81  +   82  INT   83  RCL O   84  INT   85  -   86  X<>Y   87  RCL IND Y   88  +   89  RCL IND M   90  STO Z 91  RCL Q   92  ST+ IND M   93  +   94  /   95  RCL O   96  X<>Y   97  SIGN   98  LBL 05   99  CLX 100  RCL IND Y 101  LASTX 102  * 103  ST- IND P 104  DSE P 105  CLX 106  DSE Y 107  GTO 05 108  R^ 109  STO IND M 110  GTO 03 111  LBL 06 112  RCL O 113  INT 114  LASTX 115  FRC 116  PI 117  INT 118  10^X 119  ST/ Z 120  * 121  + 122  ISG X 123  CLRGX 124  RCL M 125  RCL Q 126  ST* 00 127  CHS 128  STO IND Y 129  LBL 07 130  ST/ IND Y 131  DSE Y 132  GTO 07 133  LBL 08 134  RCL N 135  ST- O 136  RCL O 137  RCL M 138  1 139  - 140  X<=Y? 141  CLX 142  DSE X 143  GTO 01 144  RCL 00 145  CLA 146  END

( 223 bytes / SIZE 1+n.m )

 STACK INPUTS OUTPUTS T / n.nnn Z / / Y n / X m Det A

where  n = number of raws
and   m = number of columns

Example1:

-3 = 2 x + 3 y - 4 z              R01   R04   R07   R10
21 = 4 x - 5 y + 7 z      =     R02   R05   R08   R11
38 = 4 x + 2 y + 6 z            R03   R06   R09   R12

3   ENTER^
4   XEQ "HLS"   >>>>   Det =  -188                       ---Execution time = 25s---

R01 = x = 2              R04   R07   R10           1    0    0
R02 = y = 3    and    R05   R08   R11    =     0    1    0
R03 = z = 4              R06   R09   R12           0    0    1

Example2:    To invert the same matrix, store

R01   R04   R07   R10   R13   R16         1    0    0    2    3   -4
R02   R05   R08   R11   R14   R17   =    0    1    0    4   -5   7
R03   R06   R09   R12   R15   R18         0    0    1    4    2    6

3   ENTER^
6   XEQ "HLS"   >>>>   Det =  -188                       ---Execution time = 42s---

R01   R04   R07         11/47    13/94    -1/188
R02   R05   R08   =    -1/47    -7/47      15/94          with small roundoff-errors...
R03   R06   R09         -7/47    -2/47      11/94

Notes:

-With n = 7 & m = 8, Execution time = 2mn50s
-The Householder method is very sound and no pivoting is needed.
-On the other hand, it's slower and roundoff-errors may be slightly larger than when using Gaussian elimination.
-"LS3" inverses the matrix above in 25 seconds instead of 42 !