[VA] SRC #015b - HP-15C & clones: COMPLEX Matrix Inverse up to 8x8

[Valentin Albillo]

 
Hi all,

This is my 1,000^th post in this New Forum (plus ~1,600 posts in the old one) since I joined back in 2015 and I want it to be a particularly noteworthy one, so here you are !

Welcome to this second part of my latest thread SRC #015 - HP-15C & clones: Big NxN Matrix Inverse & Determinant, which I created to commemorate the availability of the HP-15C Collector's Edition.

Note: This second part will also be included together with the first part in my future full article HP Article VA058 - Boldly Going - HP-15C CE Big Matrix Woes. Although I created both parts at the same time I eventually decided to split the whole into two separate threads lest it would be too long and unwieldy if left in one piece.

The first part dealt with computing the NxN real matrix inverse and determinant for large values of N (2 ≤ N ≤ 16,) in particular for N exceeding the firmware limitation of 8x8 matrices. Additional posts to that thread by Werner, J-F Garnier and myself further dealt with:

• Real matrix inversion and determinant up to 13x13 (my OP)
• Real linear systems up to 13x13,
• Complex matrix inversion up to 6x6
• Complex systems of equations up to 6x6,
• Computing the absolute value of the determinant of a complex matrix up to 6x6,

using various approaches. See that thread for the details, a must read. So far so good.

However, in the particular case of complex matrix inversion we can do better than 6x6, as we'll see now.

The Problem

The approach that Werner and J-F Garnier followed in Part 1 for inverting complex matrices is based in the technique featured in both the HP-15C CE Owner's Handbook (p.165) and the Advanced Functions Handbook (p.128,) consisting in converting the NxN complex matrix to a 2Nx2N real matrix, which is then dealt with by using the built-in functions for real matrices, as well as pre- and post-applying a number of ad-hoc transformations (MATRIX 2, MATRIX 3, Py,x, Cy,x.)

This is most fine and dandy but it has a number of serious problems, among them:

1. It's quite a cumbersome procedure, with up to 4 ad-hoc transformations to remember and apply and up to eight or nine manual steps to obtain the complex matrix inverse, let alone solving a system of complex equations. As Maximilian Hohmann said in this recent post (my bold):
    
   "[...] I have been into this calculator thing since 1976 or so. Yes, I have a "real" 15C and the "SE" and now the "CE", but, for example, the way complex matrices are dealt with is beyond my comprehension [...] I could work my [way] through the manual and follow the examples, but if I then need to perform such a calculation three months later I would have to start from scratch again [...] by today's standards, the way it is done is totally counter-intuitive [...]"
    
2. Converting the original NxN complex matrix to a 2Nx2N real matrix requires twice the memory, which seriously limits the size of the largest complex matrix that can be inverted within the available memory.
   
   For instance, an 8x8 complex matrix requires 128 registers ("regs" for short) just to be stored, but inverting it using HP's described conversion to a 2Nx2N real matrix would require twice that, i.e. 256 regs. This is extremely wasteful and matter of fact the HP-15C CE in 192-reg Mode can't invert complex matrices larger than 6x6 using HP's approach, as 7x7 and 8x8 require 196 and 256 regs, respectively, which simply aren't available. So, what to do ?
   

My Solution

Once more, an entirely new strategy is needed. The main problem is that conversion to a 2Nx2N real matrix uses 4x the memory that an NxN real matrix would need, instead of 2x, but implementing an LU-based procedure would result in a long program much slower than the 2Nx2N approach.

To sidestep this problem I use an approach which avoids wasting memory by not converting the NxN complex matrix to a 2Nx2N real one, while also using the built-in microcode real inversion instruction (1/x) for maximum speed, making sure it never has to invert a matrix larger than 8x8. Enter split complex matrices.

A split complex matrix is an NxN complex matrix M considered to be split into two real NxN matrices A, B, containing the real and imaginary parts, respectively, of the complex elements of M, like this:

       │ 5 + 3i   4 + 7i     8      1 + 2i   6 + 6i │
       │ 7 + 2i     5i     3 + 4i   8 +  i     1    │
   M = │ 8 + 4i   2 + 5i   6 + 7i   3 + 8i     5    │ = A + iB
       │ 6 +  i   4 + 2i   3 + 7i   4 + 2i   6 + 5i │
       │ 3 + 7i     i      8 + 7i   1 + 3i   2 + 4i │

where A (real parts) and B (imaginary parts) are:

       │ 5   4   8   1   6 │        │ 3   7   0   2   6 │
       │ 7   0   3   8   1 │        │ 2   5   4   1   0 |
   A = │ 8   2   6   3   5 │ ,  B = │ 4   5   7   8   0 │
       │ 6   4   3   4   6 │        │ 1   2   7   2   5 │
       │ 3   0   8   1   2 │        │ 7   1   7   3   4 │

and as it happens, there are efficient algorithms to invert such split complex matrices by interacting directly with their parts stored in real matrices, thus no complex arithmetic is ever needed as all computations are carried out in the real domain, with a noticeable speed advantage which is further increased because my routine spends most of its time executing microcode, not RPN user code.

Another important consideration is code size. Like the original one, the HP-15C CE has no way to save/load programs to/from mass storage (say by connecting to a laptop,) thus entering programs has to be done by manually typing them in, which becomes exceedingly bothersome, time-consuming and error-prone when the programs are long (say 100-200 steps or more,) and all you see are numeric keycodes.

It's one thing if you're using a single long program all the time but entirely another if you want to use your CE for various purposes requiring assorted programs, as eventually you'll be forced to erase one to make room for another and this will get utterly inconvenient, slow and annoying pretty soon.

To wit, a program might be a fine achievement and work well but if the user has to painstakingly key in a large number of program steps to use it, chances are it won't be used much if at all, not even to just check it out. Thus, I think it's important that my solution is very short, to maximize its potential use and minimize the burden of loading it into the calculator.


The implementation:  Complex Matrix Inversion up to 8x8

This 31-step (32-byte) routine will invert in place an NxN split complex matrix M for 1 ≤ N ≤ 8 (subject to available memory.) It takes no inputs but the caller (the user or another program) must have previously dimensioned and populated the two real matrices A, B with the corresponding parts of M's elements.

Once it returns, the original values in A, B will have been replaced with those of the complex inverse matrix, which the user or the caller program may then proceed to output or use as desired. In other words, this routine can be called from the keyboard or another program (in which case it could be directly embedded into it if called from a single location) but it doesn't perform any input or output operations, it just does the inversion in place.

Program listing

  LBL C           001- 42,21,13
  RCL MATRIX A    002- 45,16,11
  RCL MATRIX B    003- 45,16,12
  RESULT E        004- 42,26,15
   +              005-       40
  1/x             006-       15
  RCL MATRIX B    007- 45,16,12
  RCL MATRIX A    008- 45,16,11
  RESULT B        009- 42,26,12
   -              010-       30
  X<>Y            011-       34
  RESULT A        012- 42,26,11
   x              013-       20
  CHS             014-       16
  LASTX           015-    43 36
  RESULT E        016- 42,26,15
  1/x             017-       15
  X<>Y            018-       34
  RCL MATRIX B    019- 45,16,12
  MATRIX 6        020- 42,16, 6
  1/x             021-       15
  RCL MATRIX A    022- 45,16,11
  RESULT B        023- 42,26,12
   x              024-       20
  RESULT A        025- 42,26,11
   +              026-       40
  LASTX           027-    43 36
  RCL MATRIX E    028- 45,16,15
  RESULT B        029- 42,26,12
   -              030-       30
  RTN             031-    43 32

     Notes:

• This routine doesn't use or alter any numbered storage registers, including the three permanent index registers R0, R1 and RI.
   
• Also, it uses no flags, labels (save LBL C), branching, loops (simple or nested), logic tests of any kind or scalar operations and it runs sequentially from its first to its last step, executing each just once, which means it executes exactly 31 user-code instructions in all (not hundreds or thousands like other approaches,) so it runs very fast (0.93" to invert a 5x5 complex matrix, 2.4" for a 7x7 one.)
   
• The inversion is performed in place: once the process ends the elements of the complex inverse replace those of the original matrix (so reinverting the computed inverse would get back the original complex matrix) and for what is worth, on completion it leaves the A matrix (real parts) in Y and the B matrix (imaginary parts) in X.
   
• The approach used requires 3NxN regs instead of the 4NxN regs required by the conversion to a real 2Nx2N matrix. The resulting NxN regs saved can be used to allow processing larger matrices or for other purposes. After the inversion, yet another batch of NxN regs can be freed, as discussed next.
   
• Once the inversion is over, you can redimension auxiliary matrix E to 0x0, which frees NxN regs for other uses such as additional data or code (up to 7NxN program steps.) For instance, to solve an NxN system of N complex equations you could dimension the constant and solutions matrices to be Nx2 each, which would still leave enough memory for the program code to matrix-multiply the inverse matrix and the constant matrix to produce the complex solution matrix.
   
  Hint: To redimension a matrix (say E) to 0x0 you don't need to execute 0, ENTER, DIM E, just 0, DIM E or CLX, DIM E will do no matter what's in the Y register (even if it holds a huge or negative number or even a matrix descriptor).

     Requirements:

• The maximum size NxN complex matrix M you can invert depends on the memory available in your physical or virtual device, as per this table:
  
         M    A,B,E   Regs  +Prog    Comments
      -------------------------------------------------------------
        1x1    1x1      3      8        -
        2x2    2x2     12     17        -
        3x3    3x3     27     32        -
        4x4    4x4     48     53    Max. size w/ 15C/64 but see (*)
        5x5    5x5     75     80    ditto CE/96
        6x6    6x6    108    113        -
        7x7    7x7    147    152    ditto CE/192 but see (**)
        8x8    8x8    192    197    ditto DM15/M1B
  
  so e.g. if you want to invert an 8x8 complex matrix you need 197 regs available, which includes matrices A, B (and E,) plus 5 regs to hold the routine itself. Auxiliary matrix E also needs to be NxN because the HP-15C can't multiply two NxN matrices in place, a third NxN matrix (E) is needed to store the result.
   
  Implementation note: A 7x7 real matrix occupies 49 regs so you can have only three such matrices simultaneously at any given time because using a fourth matrix would require 196 regs, exceeding the 192 regs available in the CE/192.
  
  To overcome this limitation, I've used a trick to make do with only three matrices. If there was room for a fourth matrix (which is the case when using a CE/192 and dealing with 6x6 or smaller complex matrices,) a version of my routine limited to those sizes would be ~30% faster.
   
• The routine works as-is for complex matrices up to and including 8x8. Larger complex matrices wouldn't work because of both the 8x8 firmware limitation for (real) inversion and memory constraints. As of Sept. 2023 no physical or virtual devices exist which provide more than 229 regs because the HP-15C has a RAM limit of 256 regs and 27 of those are used internally by the firmware.
  
  This means that matrices 9x9 or larger can't be tackled using this routine on any currently existing device. However, it might be possible that some future patch would override those restrictions, in which case this routine would invert complex matrices larger than 8x8 without any modifications.
   
• The matrix A+B must be invertible, i.e. det(A+B) # 0. For most real-life uses this will be the case but if this condition isn't met there are slight variations to this routine that would work Ok in those cases.
  
  

     (*) Observation re the original HP-15C and 4x4 complex matrices:

• The original HP-15C could invert a 4x4 complex matrix but it took all 64 regs available and was a complicated, completely manual process as there wasn't any memory left for program code. On the other hand, running my routine is a fast, automated process that leaves out all the drudgery (transformations, etc.) and still leaves 11 regs (up to 77 extra program steps) free for additional code or data.
  
  Furthermore, now it seems possible to solve a system of 4 complex equations with equal ease, like this (see Worked Example below for basic details):
   
  • Initialize and store the data in the A, B matrices.
     
  • Call my complex inversion routine. When it ends, the complex inverse is stored in A, B and there's still 11 regs free. Optionally, output the inverse by the usual procedure (RCL A, B in USER mode) either now or after the system is solved, as the solving procedure doesn't affect the inverse.
     
  • Get rid of auxiliary matrix E (redimension it to 0x0.) This frees another 16 regs, so there's now 27 regs available for what follows.
     
  • Dimension both the constant matrix (say C) and the solution matrix (say D) to be 4x2 and populate the constant matrix. This will leave 11 regs (up to 77 program steps) still available for the matrix-multiplication code, which is left as a fairly easy exercise for the reader (just a loop which multiplies each row of the inverse by the constant matrix using complex arithmetic.)
     
  • Run said matrix-multiplication code, which should use the inverse's data in A, B to multiply the inverse and the constant matrix C, storing the result in the complex solution matrix D.
     
  • Output or otherwise make use of the complex solution matrix.
  
  All these steps can be coded as a single program which inverts the complex matrix, gets rid of E to make room, dimensions the constant and solution matrices and stops for the user to populate the constant matrix; once done, the user presses R/S to continue execution and the solution matrix is computed for the user to output or other code to use it.
  
  The inverse matrix is left undisturbed so other systems having the same original matrix and different constant matrices can be solved outright using the same inverse matrix. Moreover, the original matrix can be recovered if desired (negligible rounding errors aside,) by ensuring that the inverse remains intact and there's at least 16 regs available, then reinverting the inverse matrix with a GSB C.
  
  This would completely automate the task of inverting a 4x4 complex matrix or solving a system of 4 complex equations on an original HP-15C, with no complicated transformations or manual steps whatsoever.
  
  

     (*) Observation re the HP-15C CE and 8x8 complex matrices:

• Though there's not enough room in the HP-15C CE in 192-regs Mode to invert an 8x8 complex matrix M by running the routine featured here (197 regs would be needed; the routine itself wouldn't fit) there's just enough room for the split matrices A, B and the auxiliary matrix E (192 regs in all,) so the 8x8 complex matrix can be inverted in a pinch if the user executes the 29 program instructions manually in sequential order. The procedure would be like this:
   
  • Initialize and store the data in the A, B matrices.
     
  • Carefully execute manually the 29 instructions from  002 RCL MATRIX A  to  030 -  in sequential order. Assuming you're reasonably proficient with the HP-15C, this should take ~ 3 min. or less.
     
  • Output or otherwise make use of the inverted complex matrix data in A, B. Now you can redimension auxiliary matrix E to 0x0 to free 64 regs (up to 448 program steps) for further processing using the inverse matrix just computed (e.g. solving a system of 8 complex equations.)
    
  
  This procedure is perfectly workable and reasonably fast if you want to invert an 8x8 complex matrix using a CE/192 but you should be very careful as most errors keying in an instruction could mean having to restart from the beginning, including re-storing the elements anew in A, B.

Worked example

The following example can be run on the HP-15C CE (Collector's Edition) in its 96- or 192-regs Mode, as well as in any physical/virtual clone with at least 80 regs allocatable.

Note: The usefulness of this relatively small 5x5 example is twofold: first, to get to know how to use the routine and get comfortable using it and second, to ascertain that you've loaded it correctly into program memory by running it and checking the results it produces.

Invert the following 5x5 complex matrix M: (matrix taken from this J-F's post)

       │ 5 + 3i   4 + 7i     8      1 + 2i   6 + 6i │
       │ 7 + 2i     5i     3 + 4i   8 +  i     1    │
   M = │ 8 + 4i   2 + 5i   6 + 7i   3 + 8i     5    │
       │ 6 +  i   4 + 2i   3 + 7i   4 + 2i   6 + 5i │
       │ 3 + 7i     i      8 + 7i   1 + 3i   2 + 4i │

- The A (real parts) and B (imaginary parts) matrices are:

       │ 5   4   8   1   6 │       │ 3   7   0   2   6 │
       │ 7   0   3   8   1 │       │ 2   5   4   1   0 |
   A = │ 8   2   6   3   5 │ , B = │ 4   5   7   8   0 │
       │ 6   4   3   4   6 │       │ 1   2   7   2   5 │
       │ 3   0   8   1   2 │       │ 7   1   7   3   4 │

- Ensure disabled complex stack and set 4 decimal places:

   CF 8, FIX 4 

- Allocate all memory to matrices, initialize them and dimension the real/imag parts A, B to 5x5:

   1, DIM (i), MATRIX 0
   5, ENTER, DIM A, DIM B

- Store the real parts of M's elements into matrix A:

   USER, MATRIX 1

   5 STO A, 4 STO A, 8 STO A, 1 STO A, 6 STO A,
   7 STO A, 0 STO A, 3 STO A, 8 STO A, 1 STO A,
   8 STO A, 2 STO A, 6 STO A, 3 STO A, 5 STO A,
   6 STO A, 4 STO A, 3 STO A, 4 STO A, 6 STO A,
   3 STO A, 0 STO A, 8 STO A, 1 STO A, 2 STO A

- Store their imaginary parts into matrix B:

   3 STO B, 7 STO B, 0 STO B, 2 STO B, 6 STO B,
   2 STO B, 5 STO B, 4 STO B, 1 STO B, 0 STO B,
   4 STO B, 5 STO B, 7 STO B, 8 STO B, 0 STO B,
   1 STO B, 2 STO B, 7 STO B, 2 STO B, 5 STO B,
   7 STO B, 1 STO B, 7 STO B, 3 STO B, 4 STO B

- Compute in-place the complex inverse matrix M':

   GSB C -> B     5 5   {in 0.93"}

- Recall the inverse matrix elements from A and B:

   Real parts:

      RCL A: -0.1197, RCL A: -0.0070, RCL A:  0.0635, RCL A:  0.0322, RCL A:  0.0285,
      RCL A:  0.0247, RCL A: -0.0206, RCL A:  0.0232, RCL A: -0.0363, RCL A:  0.0638,
      RCL A:  0.0277, RCL A: -0.0416, RCL A:  0.0183, RCL A: -0.0348, RCL A:  0.0439,
      RCL A:  0.0397, RCL A:  0.0886, RCL A: -0.0975, RCL A:  0.0319, RCL A: -0.0036,
      RCL A:  0.0587, RCL A:  0.0388, RCL A: -0.0482, RCL A:  0.0509, RCL A: -0.0549

   Imaginary parts:

      RCL B: -0.0115, RCL B: -0.1044, RCL B:  0.0663, RCL B:  0.0940, RCL B: -0.1395,
      RCL B: -0.1038, RCL B: -0.0652, RCL B: -0.0403, RCL B:  0.0796, RCL B:  0.0945,
      RCL B:  0.0430, RCL B:  0.0005, RCL B: -0.0459, RCL B: -0.0329, RCL B:  0.0005,
      RCL B: -0.0081, RCL B:  0.1066, RCL B: -0.0765, RCL B: -0.0271, RCL B:  0.0520,
      RCL B: -0.0181, RCL B:  0.0712, RCL B:  0.0676, RCL B: -0.1183, RCL B:  0.0114

and so the 5x5 complex inverse matrix M' is:

       │ -0.1197-0.0115i -0.0070-0.1044i  0.0635+0.0663i  0.0322+0.0940i  0.0285-0.1395i │
       │  0.0247-0.1038i -0.0206-0.0652i  0.0232-0.0403i -0.0363+0.0796i  0.0638+0.0945i │
  M' = │  0.0277+0.0430i -0.0416+0.0005i  0.0183-0.0459i -0.0348-0.0329i  0.0439+0.0005i │
       │  0.0397-0.0081i  0.0886+0.1066i -0.0975-0.0765i  0.0319-0.0271i -0.0036+0.0520i │
       │  0.0587-0.0181i  0.0388+0.0712i -0.0482+0.0676i  0.0509-0.1183i -0.0549+0.0114i │

Notes:

• If in doubt, a simple way to check the inverse's correction is to reinvert it once you've written down its elements, which should still be undisturbed in memory. Simply run the routine again:
  
     C (in USER mode) or GSB C (out of USER mode) -> B     N N
  
  and you'll get the original matrix back (ignoring negligible rounding errors,) which can be verified by using RCL in USER mode, as we did in the Worked Example above.
   
• Once you're done with using the routine you might want to free memory by resizing all matrices A, B and the auxiliary matrix E to 0x0 as well as resetting the row/col indexes and reallocating the default numbered storage registers R0-R.9, like this:
  
     MATRIX 0, MATRIX 1, 19, DIM (i)
  

That's all, hope you enjoyed it and found it useful for your own purposes. Comments welcome.

V.