inv ← {vals} ##.gauss_jordan mat            ⍝ Gauss-Jordan elimination.

Gauss-Jordan elimination is a classic algorithm, implemented the D-style.

NB: This  function  is  included  only  for interest as APL provides both matrix
inverse and matrix division as a primitive function: ⌹.

Further, the results of this D-function differ very slightly from those of APL's
primitive  function. This may be because APL's C-code maintains 80-bit floating-
point  accuracy during (much, if not all of) its calculation, whereas the inter-
mediate  arrays of this function use only 64 bits. Of course, the function could
be  adapted  to use, for example, rational →rats← arithmetic for perfect precis-
ion (see test script: ##.scripts.gauss_jordan).

Matrix  right  argument  [mat] represents the coefficients of a system of linear
equations and,  if present, left argument [vals] represents the equations' value
vector (right-hand side).

For example:

      x +   y + 2z = ¯3
    ¯1x + ¯2y + 3z = 14
     3x + ¯7y + 4z = ¯3

which may be written:

    A x = v
    ¯ ¯   ¯
where:
        ┌ 1  1 2┐   v = [¯3 14 ¯3]
    A = │¯1 ¯1 3│   ¯
    ¯   │ 3 ¯7 4│
        └       ┘

Then the solution vector x is:
                         ¯
    x = v ⌹ A               ⍝ using dyadic ⌹.
    ¯   ¯   ¯
or
    x = (⌹ A) × v           ⍝ using monadic ⌹ (where × is matrix product +.×).
    ¯      ¯  ¯ ¯                                    ¯
Gauss-Jordan uses three simple transformations, none of which changes the solut-
ion of the equations. For example, starting with:

[0] Initial equations:

      x +   y + 2z = ¯3
    ¯1x + ¯2y + 3z = 14
     3x + ¯7y + 4z = ¯3

[1] Exchange rows:

    ¯1x + ¯2y + 3z = 14     ⍝ swap first and second rows.
      x +   y + 2z = ¯3
     3x + ¯7y + 4z = ¯3

[2] Multiply each term in a row by the same factor:

    ¯1x + ¯2y + 3z = 14
     3x +  3y + 6z = ¯9     ⍝ multiply second row by 3.
     3x + ¯7y + 4z = ¯3

[3] Subtract one row from another:

    ¯1x +  ¯2y +  3z = ¯2
     3x +   3y +  6z = ¯9
     0x + ¯10y + ¯2z =  6   ⍝ subtract second row from third.

The algorithm applies these three steps repeatedly, until the coefficient matrix
is transformed into an identity matrix, giving:

    1x + 0y + 0z = ¯6
    0x + 1y + 0z = ¯1
    0x + 0y + 1z =  2

or:

     x    ·    · = ¯6
     ·    y    · = ¯1
     ·    ·    z =  2

Monadic case
------------
We  find  the explicit matrix inverse by starting with an identity matrix on the
right:

[0]                                         ⍝ initial matrix:
     1x +  1y + 2z = 1 0 0
    ¯1x + ¯2y + 3z = 0 1 0
     3x + ¯7y + 4z = 0 0 1
[1]                                         ⍝ exchange first and third rows:
     3x + ¯7y + 4z = 0 0 1
    ¯1x + ¯2y + 3z = 0 1 0
     1x +  1y + 2z = 1 0 0
[2]
     1x + ¯2.333y + 1.333z = 0 0 0.333      ⍝ multiply first row by ÷3:
    ¯1x +     ¯2y +     3z = 0 1 0
     1x +      1y +     2z = 1 0 0
[3]                                         ⍝ subtract first row from others:
     1x + ¯2.333y + 1.333z = 0 0  0.333
     0x + ¯4.333y + 4.333z = 0 1  0.333
     0x +  3.333y + 0.667z = 1 0 ¯0.333
...                                         ⍝ and so on, until:
     1x +  0y + 0z = 0.25 ¯0.346  0.135
     0x +  1y + 0z = 0.25 ¯0.038 ¯0.096
     0x +  0y + 1z = 0.25  0.192 ¯0.019

Choosing a pivot value
----------------------
In  a  finite-precision implementation, such as IEEE floating point, subtracting
relatively  large  but  commensurate  numbers  introduces  a  significant error.
Gauss-Jordan  reduces this as much as possible by choosing the "pivot" value for
transformation [2] to be the item with the largest absolute magnitude within the
remaining rows of the column under consideration.

       before                                                      after
    1  0  · · · ·                                               1  0  · · · ·
    0  1  · · · ·     In this example, working on the third     0  1  · · · ·
    0  0  3 · · ·     column,  the pivot (max abs) value is     0  0 ¯4 · · ·
    0  0 ¯2 · · ·     ¯4,  resulting in the exchange of the     0  0 ¯2 · · ·
    0  0 ¯4 · · ·     third and fifth rows.                     0  0  3 · · ·
    0  0  1 · · ·                                               0  0  1 · · ·

Notice that values in the first two rows of the third column are ignored.

Put simply,  the pivot value for column ⍺ is the item from rows ⍺ downwards with
the largest absolute value: ⌈/|⍺↓⍵[;⍺].

Illustration
------------
Here is an illustration of the steps in the Gauss-Jordan elimination of:

    ┌────────┐
    │4 8 4  0│
    │1 4 7  2│
    │1 5 4 ¯3│
    │1 3 0 ¯2│
    └────────┘

Append an identity matrix:

    ┌────────┬───────┐
    │4 8 4  0│1 0 0 0│ append identity matrix.
    │1 4 7  2│0 1 0 0│
    │1 5 4 ¯3│0 0 1 0│
    │1 3 0 ¯2│0 0 0 1│
    └────────┴───────┘

Eliminate off-diagonal values from first column:
    ┌┬───────┬───────┐
    │4 8 4  0│1 0 0 0│ first col; pivot value is 4.        (4)· · ·
    │1 4 7  2│0 1 0 0│                                      1 · · ·
    │1 5 4 ¯3│0 0 1 0│                                      1 · · ·
    │1 3 0 ¯2│0 0 0 1│                                      1 · · ·
    └┴───────┴───────┘
    ┌────────┬──────────┐
    ├1─2─1──0┼0.25─0─0─0┤ divide first row by pivot value.
    │1 4 7  2│0    1 0 0│
    │1 5 4 ¯3│0    0 1 0│
    │1 3 0 ¯2│0    0 0 1│
    └────────┴──────────┘
    ┌─────────┬───────────┐
    │1 2  1  0│ 0.25 0 0 0│ subtract first row from remaining rows.
    │0 2  6  2│¯0.25 1 0 0│ to leave 0s in all off-diagonal rows
    │0 3  3 ¯3│¯0.25 0 1 0│ of first column.
    │0 1 ¯1 ¯2│¯0.25 0 0 1│
    └─────────┴───────────┘

Eliminate off-diagonal values from second column:
    ┌──┬──────┬───────────┐
    │1 2  1  0│ 0.25 0 0 0│                                 · · · ·
    │0 2  6  2│¯0.25 1 0 0│                                 · 2 · ·
    │0 3  3 ¯3│¯0.25 0 1 0│ second col; pivot value is 3.   ·(3)· ·
    │0 1 ¯1 ¯2│¯0.25 0 0 1│                                 · 1 · ·
    └──┴──────┴───────────┘
    ┌─────────┬───────────┐
    │1 2  1  0│ 0.25 0 0 0│                                 · · · ·
    ├0─3──3─¯3┼¯0.25─0─1─0┤ swap second row with            ·(3)· ·
    ├0─2──6──2┼¯0.25─1─0─0┤ pivot value row.                · 2 · ·
    │0 1 ¯1 ¯2│¯0.25 0 0 1│                                 · 1 · ·
    └─────────┴───────────┘
    ┌─────────┬───────────────────┐
    │1 2  1  0│ 0.25    0 0      0│
    ├0─1──1─¯1┼¯0.08333─0─0.3333─0┤ divide second row by pivot value.
    │0 2  6  2│¯0.25    1 0      0│
    │0 1 ¯1 ¯2│¯0.25    0 0      1│
    └─────────┴───────────────────┘
    ┌─────────┬────────────────────┐
    │1 0 ¯1  2│ 0.4167  0 ¯0.6667 0│ subtract multiples of second row from
    │0 1  1 ¯1│¯0.08333 0  0.3333 0│ remaining rows to leave 0s in all
    │0 0  4  4│¯0.08333 1 ¯0.6667 0│ off-diagonal rows of second column.
    │0 0 ¯2 ¯1│¯0.1667  0 ¯0.3333 1│
    └─────────┴────────────────────┘

Eliminate off-diagonal values from third column:
    ┌─────┬───┬────────────────────┐
    │1 0 ¯1  2│ 0.4167  0 ¯0.6667 0│                                · · · ·
    │0 1  1 ¯1│¯0.08333 0  0.3333 0│                                · · · ·
    │0 0  4  4│¯0.08333 1 ¯0.6667 0│ third col; pivot value is 4.   · ·(4)·
    │0 0 ¯2 ¯1│¯0.1667  0 ¯0.3333 1│                                · ·¯2 ·
    └─────┴───┴────────────────────┘
    ┌─────────┬───────────────────────┐
    │1 0 ¯1  2│ 0.4167  0    ¯0.6667 0│
    │0 1  1 ¯1│¯0.08333 0     0.3333 0│
    ├0─0──1──1┼¯0.02083─0.25─¯0.1667─0┤ divide third row by pivot value.
    │0 0 ¯2 ¯1│¯0.1667  0    ¯0.3333 1│
    └─────────┴───────────────────────┘
    ┌────────┬────────────────────────┐
    │1 0 0  3│ 0.3958   0.25 ¯0.8333 0│ subtract multiples of third row from
    │0 1 0 ¯2│¯0.0625  ¯0.25  0.5    0│ remaining rows to leave 0s in all
    │0 0 1  1│¯0.02083  0.25 ¯0.1667 0│ off-diagonal rows of third column.
    │0 0 0  1│¯0.2083   0.5  ¯0.6667 1│
    └────────┴────────────────────────┘

Eliminate off-diagonal values from fourth column:
    ┌───────┬┬────────────────────────┐
    │1 0 0  3│ 0.3958   0.25 ¯0.8333 0│                                 · · · ·
    │0 1 0 ¯2│¯0.0625  ¯0.25  0.5    0│                                 · · · ·
    │0 0 1  1│¯0.02083  0.25 ¯0.1667 0│                                 · · · ·
    │0 0 0  1│¯0.2083   0.5  ¯0.6667 1│ fourth col; pivot value is 1.   · · ·(1)
    └───────┴┴────────────────────────┘
    ┌───────┬────────────────────────┐
    │1 0 0 0│ 1.021  ¯1.25  1.167  ¯3│ subtract multiples of fourth row
    │0 1 0 0│¯0.4792  0.75 ¯0.8333  2│ from remaining rows to leave 0s in all
    │0 0 1 0│ 0.1875 ¯0.25  0.5    ¯1│ off-diagonal rows of fourth column.
    │0 0 0 1│¯0.2083  0.5  ¯0.6667  1│
    └───────┴────────────────────────┘

Drop leading identity matrix for inverse of original matrix:
            ┌────────────────────────┐
            │ 1.021  ¯1.25  1.167  ¯3│
            │¯0.4792  0.75 ¯0.8333  2│
            │ 0.1875 ¯0.25  0.5    ¯1│
            │¯0.2083  0.5  ¯0.6667  1│
            └────────────────────────┘

(muse:

    As supplied,  the function is offensive to the functional programming purist
    because it uses destructive assignment:

        mat←⍵                               ⍝ name matrix for updating.
        mat[⍺ p;]←mat[p ⍺;]                 ⍝ exchange ⍺th and pth rows.
        mat[⍺;]÷←mat[⍺;⍺]                   ⍝ reduce col diagonal to 1.

    For fun, we could recast it in a (non-destructive) purer form:

    gauss_jordan←{⎕IO←0                       ⍝ Gauss-Jordan elimination.

        elim←{                                ⍝ elimination of row/col ⍺.
            mask←⍺=⍳⊃⍴⍵                       ⍝ selection mask for ⍺th row.
            pivt←{⍵⍳⌈/⍵}|⍵[;⍺]×∨\mask         ⍝ position of pivot row.
            mat1←⍺ pivt swap ⍵                ⍝ exchanged ⍺th and pth rows.
            mat2←mat1÷[0]⍵[pivt;⍺]*mask       ⍝ ⍺th row divided by pivot.
            mat2-(mat2[;⍺]×~mask)∘.×mat2[⍺;]  ⍝ remaining rows reduced by pvt.
        }

        swap←{                                ⍝ matrix ⍵ with rows ⍺ exchanged.
            i←⍳⊃⍴⍵                            ⍝ index vector for rows.
            y n←1 0=⊂i∊⍺                      ⍝ mask and not-mask for swap rows.
            ⍵[(n\n/i)+y\⌽y/i;]                ⍝ merge with reversed rows.
        }

        0::⎕SIGNAL ⎕EN                        ⍝ pass back error to caller.

        ⍺←=/↑⍳⍴⍵                              ⍝ id matrix for monadic case.

        (⍴⍺)⍴(0 1×⍴⍵)↓↑elim/(⌽⍳⌊/⍴⍵),⊂⍵,⍺     ⍝ elimination/ ··· 2 1 0 (⍵,⍺)

        See also →at←
    }

    Further,  we could remove square-bracket-indexing and represent our matrices
    as vectors-of-row-vectors:

    gauss_jordan←{⎕IO←0                     ⍝ Gauss-Jordan elimination.

        elim←{                              ⍝ elimination of row/col ⍺.
            mask←⍺=⍳⍴⍵                      ⍝ selection mask for ⍺th row/col.
            pivt←{⍵⍳⌈/⍵}|(⍺⊃¨⍵)×∨\mask      ⍝ position of pivot row.
            mat1←⍺ pivt swap ⍵              ⍝ exchanged ⍺th and pivot'th rows.
            mat2←mat1÷(pivt ⍺⊃⍵)*mask       ⍝ ⍺th row divided by pivot.
            mat2-((⍺⊃¨mat2)×~mask)×⊂⍺⊃mat2  ⍝ remaining rows reduced by pivot.
        }

        swap←{                              ⍝ swap rows ⍺ in matrix ⍵.
            i←⍳⍴⍵                           ⍝ index vector for rows.
            y n←1 0=⊂i∊⍺                    ⍝ mask and not-mask for swap rows.
            ((n\n/i)+y\⌽y/i)⊃¨⊂⍵            ⍝ merge with reversed rows.
        }

        0::⎕SIGNAL ⎕EN                      ⍝ pass back error to caller.

        ⍺←=/↑⍳⍴⍵                            ⍝ id matrix for monadic case.

        (⍴⍺)⍴(0 1×⍴⍵)↓↑↑elim/(⌽⍳⌊/⍴⍵),⊂↓⍵,⍺ ⍝ elimination/ ··· 2 1 0 (⍵,⍺)
    }

    The  clinically  obsessive would now be in a position to refine the function
    still further by replacing local names with additional inner D-functions and
    operators.

    Pursuing  this  approach  to  the edge of reason, we wind up with a function
    that  is a single expression of functions and operators, with no assignments
    or  guards.  Here  is  the  code for the monadic case (the dyadic equivalent
    replaces  (=/↑⍳⍴⍵)  with  ⍺ in the first line). The expression is strung out
    over a number of lines with some white space:

    {
    ·   (=/↑⍳⍴⍵){
    ·   ·   (⍴⍺)⍴(0 1×⍴⍵)↓↑↑{
    ·   ·   ·   ⍺(
    ·   ·   ·   ·   (⍺=⍳⍴⍵){
    ·   ·   ·   ·   ·   ⍺(
    ·   ·   ·   ·   ·   ·   ⍺⍺{
    ·   ·   ·   ·   ·   ·   ·   ⍵-(
    ·   ·   ·   ·   ·   ·   ·   ·   (⍺⊃¨⍵)×~⍺⍺
    ·   ·   ·   ·   ·   ·   ·   )×⊂⍺⊃⍵
    ·   ·   ·   ·   ·   ·   }
    ·   ·   ·   ·   ·   )⍺(
    ·   ·   ·   ·   ·   ·   (
    ·   ·   ·   ·   ·   ·   ·   {⍵⍳⌈/⍵}|(⍺⊃¨⍵)×∨\⍺⍺
    ·   ·   ·   ·   ·   ·   ){
    ·   ·   ·   ·   ·   ·   ·   (
    ·   ·   ·   ·   ·   ·   ·   ·   ⍺ ⍺⍺{
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ⍺(
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   (⍳⍴⍵){
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ↑⍵{
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   (
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   (⍵\⍵/⍵⍵)+⍺\⌽⍺/⍵⍵
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   )⊃¨⊂⍺⍺
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   }⍺⍺/1 0=⊂⍺⍺∊⍺
    ·   ·   ·   ·   ·   ·   ·   ·   ·   ·   }
    ·   ·   ·   ·   ·   ·   ·   ·   ·   )⍵
    ·   ·   ·   ·   ·   ·   ·   ·   }⍵
    ·   ·   ·   ·   ·   ·   ·   )÷(⍺⍺ ⍺⊃⍵)*⍵⍵
    ·   ·   ·   ·   ·   ·   }⍺⍺
    ·   ·   ·   ·   ·   )⍵
    ·   ·   ·   ·   }
    ·   ·   ·   )⍵
    ·   ·   }/(⌽⍳⌊/⍴⍵),⊂↓⍵,⍺
    ·   }⍵
    }

    We  can check that the coding is correct by pasting the lines into character
    variable gj_src, say, and then executing:

        ⎕io ⎕pp ← 0 4                                   ⍝ condition environment.

        (⍎gj_src~⎕tc,'·') 3 3⍴ 1 2 3, ¯3 1 5, 2 4 ¯1    ⍝ apply function.
     0.4286 ¯0.2857 ¯0.1429
    ¯0.1429  0.1429  0.2857
     0.2857  0      ¯0.1429
)

Examples:

    ⎕←mat←4 4⍴ 4 8 4 0, 1 4 7 2, 1 5 4 ¯3, 1 3 0 ¯2
4 8 4  0
1 4 7  2
1 5 4 ¯3
1 3 0 ¯2

    ⎕pp←4                           ⍝ display numbers to 4 sig-figs.

    gauss_jordan mat                ⍝ matrix inverse.
 1.021  ¯1.25  1.167  ¯3
¯0.4792  0.75 ¯0.8333  2
 0.1875 ¯0.25  0.5    ¯1
¯0.2083  0.5  ¯0.6667  1

    1 2 3 4 gauss_jordan mat        ⍝ matrix divide.
¯9.979 6.521 ¯2.813 2.792

    (⌹mat)≡gauss_jordan mat         ⍝ check against primitive matrix inverse.
1

    (1 2 3 4⌹mat)≡1 2 3 4 gauss_jordan mat      ⍝ ditto matrix divide.
1

⍝ Using an order-⍵ Hilbert matrix, we can see the slight variation between
⍝ primitive ⌹ and the D-function, as errors accumulate.

    hil←{÷1+∘.+⍨(⍳⍵)-⎕IO}           ⍝ order-⍵ Hilbert matrix.

⍝ The results are identical up to hil 5:

    ↑{(⌹⍵)(gauss_jordan ⍵)}∘hil¨ 0 to 5
┌───────────────────────────────────┬───────────────────────────────────┐
│0                                  │0                                  │
├───────────────────────────────────┼───────────────────────────────────┤
│1                                  │1                                  │
├───────────────────────────────────┼───────────────────────────────────┤
│ 4 ¯6                              │ 4 ¯6                              │
│¯6 12                              │¯6 12                              │
├───────────────────────────────────┼───────────────────────────────────┤
│  9  ¯36   30                      │  9  ¯36   30                      │
│¯36  192 ¯180                      │¯36  192 ¯180                      │
│ 30 ¯180  180                      │ 30 ¯180  180                      │
├───────────────────────────────────┼───────────────────────────────────┤
│  16  ¯120   240  ¯140             │  16  ¯120   240  ¯140             │
│¯120  1200 ¯2700  1680             │¯120  1200 ¯2700  1680             │
│ 240 ¯2700  6480 ¯4200             │ 240 ¯2700  6480 ¯4200             │
│¯140  1680 ¯4200  2800             │¯140  1680 ¯4200  2800             │
├───────────────────────────────────┼───────────────────────────────────┤
│   25   ¯300    1050   ¯1400    630│   25   ¯300    1050   ¯1400    630│
│ ¯300   4800  ¯18900   26880 ¯12600│ ¯300   4800  ¯18900   26880 ¯12600│
│ 1050 ¯18900   79380 ¯117600  56700│ 1050 ¯18900   79380 ¯117600  56700│
│¯1400  26880 ¯117600  179200 ¯88200│¯1400  26880 ¯117600  179200 ¯88200│
│  630 ¯12600   56700  ¯88200  44100│  630 ¯12600   56700  ¯88200  44100│
└───────────────────────────────────┴───────────────────────────────────┘

⍝ However, at hil 6, owing to error accumulation,
⍝ small differences begin to appear:

    ⍪ {(⌹⍵)(gauss_jordan ⍵)}hil 6
┌────────────────────────────────────────────────────────────────────────────────────────────┐
│   36            ¯630.0000001     3360          ¯7560.000001     7560.000001    ¯2772.000001│
│ ¯630.0000001   14700           ¯88200.00001   211680         ¯220500           83160.00002 │
│ 3360          ¯88200.00001     564480.0001  ¯1411200         1512000         ¯582120.0001  │
│¯7560.000001   211680         ¯1411200        3628800.001    ¯3969000.001     1552320       │
│ 7560.000001  ¯220500          1512000       ¯3969000.001     4410000.001    ¯1746360       │
│¯2772.000001    83160.00002    ¯582120.0001   1552320        ¯1746360          698544.0001  │
├────────────────────────────────────────────────────────────────────────────────────────────┤
│   36           ¯630           3360          ¯7560     7560.000001    ¯2772                 │
│ ¯630          14700         ¯88200.00001   211680  ¯220500           83160.00001           │
│ 3360         ¯88200.00001   564480       ¯1411200  1512000         ¯582120                 │
│¯7560         211680       ¯1411200        3628800 ¯3969000         1552320                 │
│ 7560.000001 ¯220500        1512000       ¯3969000  4410000        ¯1746360                 │
│¯2772          83160.00001  ¯582120        1552320 ¯1746360          698544.0001            │
└────────────────────────────────────────────────────────────────────────────────────────────┘

See also: rats det Cholesky at

Back to: contents

Back to: Workspaces