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←{⎕ML ⎕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←{ ⍝ 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 (⍵,⍺)
}
Further, we could remove square-bracket-indexing and represent our matrices
as vectors-of-row-vectors:
gauss_jordan←{⎕ML ⎕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 ⎕ml) ⎕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:
0 1 0 disp ↑{(⌹⍵)(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:
0 1 0 disp ,[⍬] {(⌹⍵)(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
Back to: contents
Back to: Dyalog APL
Trouble seeing APL font?