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│
└────────────────────────┘
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