Gauss-Jordan method¶
We consider a \(n \times n\)-sized square matrix \(A_{ij}\). Assuming that \(A_{ij}\) is a regular matrix, our objective is to compute the inverse \(A_{ij}^{-1}\).
Starting from
\[A_{ik} X_{kj} = I_{ij},\]
the inverse matrix is obtained by
\[X_{ij}
=
A_{ik}^{-1}
I_{kj}
=
A_{ij}^{-1}.\]
Component-wise, starting from the initial state:
\[\begin{split}\newcommand{\x}[2]{a_{#1,#2}}
\newcommand{\y}[2]{b_{#1,#2}}
\begin{bmatrix}
\x{ 0}{ 0} & \x{ 0}{ 1} & \cdots & \x{ 0}{i-1} & \x{ 0}{i } & \x{ 0}{i+1} & \cdots & \x{ 0}{n-2} & \x{ 0}{n-1} \\
\x{ 1}{ 0} & \x{ 1}{ 1} & \cdots & \x{ 1}{i-1} & \x{ 1}{i } & \x{ 1}{i+1} & \cdots & \x{ 1}{n-2} & \x{ 1}{n-1} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
\x{i-1}{ 0} & \x{i-1}{ 1} & \cdots & \x{i-1}{i-1} & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} \\
\x{i }{ 0} & \x{i }{ 1} & \cdots & \x{i }{i-1} & \x{i }{i } & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} \\
\x{i+1}{ 0} & \x{i+1}{ 1} & \cdots & \x{i+1}{i-1} & \x{i+1}{i } & \x{i+1}{i+1} & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\x{n-2}{ 0} & \x{n-2}{ 1} & \cdots & \x{n-2}{i-1} & \x{n-2}{i } & \x{n-2}{i+1} & \cdots & \x{n-2}{n-2} & \x{n-2}{n-1} \\
\x{n-1}{ 0} & \x{n-1}{ 1} & \cdots & \x{n-1}{i-1} & \x{n-1}{i } & \x{n-1}{i+1} & \cdots & \x{n-1}{n-2} & \x{n-1}{n-1} \\
\end{bmatrix}
\begin{bmatrix}
1 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 \\
0 & 1 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & 0 & 0 & \cdots & 0 & 0 \\
0 & 0 & \cdots & 0 & 1 & 0 & \cdots & 0 & 0 \\
0 & 0 & \cdots & 0 & 0 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 1 & 0 \\
0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 1 \\
\end{bmatrix},\end{split}\]
and we aim at making the first and the second the identity and the inverse matrices, respectively.
Forward sweep¶
After forward-sweeping up to the \(i - 1\)-th row, we have
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} \\
0 & 0 & \cdots & 0 & \x{i }{i } & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} \\
0 & 0 & \cdots & 0 & \x{i+1}{i } & \x{i+1}{i+1} & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\end{bmatrix}
\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
\y{i-1}{ 0} & \y{i-1}{ 1} & \cdots & \y{i-1}{i-1} & 0 & 0 & \cdots & 0 & 0 \\
\y{i }{ 0} & \y{i }{ 1} & \cdots & \y{i }{i-1} & 1 & 0 & \cdots & 0 & 0 \\
\y{i+1}{ 0} & \y{i+1}{ 1} & \cdots & \y{i+1}{i-1} & 0 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
The \(i\)-th row is normalized (divided by \(\x{i}{i}\)) to obtain
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} \\
0 & 0 & \cdots & 0 & 1 & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} \\
0 & 0 & \cdots & 0 & \x{i+1}{i } & \x{i+1}{i+1} & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\end{bmatrix}
\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
\y{i-1}{ 0} & \y{i-1}{ 1} & \cdots & \y{i-1}{i-1} & 0 & 0 & \cdots & 0 & 0 \\
\y{i }{ 0} & \y{i }{ 1} & \cdots & \y{i }{i-1} & \y{i }{i } & 0 & \cdots & 0 & 0 \\
\y{i+1}{ 0} & \y{i+1}{ 1} & \cdots & \y{i+1}{i-1} & 0 & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
src/gauss_jordan.c¶
21const double f = 1. / arr[i * nitems + i];
22arr[i * nitems + i] = 1.;
23for (size_t j = i + 1; j < nitems; j++) {
24 arr[i * nitems + j] *= f;
25}
26for (size_t j = 0; j < i + 1; j++) {
27 inv[i * nitems + j] *= f;
28}
The new \(i\)-th row is utilized to eliminate the \(i\)-th column:
\[\begin{split}\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} \\
0 & 0 & \cdots & 0 & 1 & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} \\
0 & 0 & \cdots & 0 & 0 & \x{i+1}{i+1} & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\end{bmatrix}
\begin{bmatrix}
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
\y{i-1}{ 0} & \y{i-1}{ 1} & \cdots & \y{i-1}{i-1} & 0 & 0 & \cdots & 0 & 0 \\
\y{i }{ 0} & \y{i }{ 1} & \cdots & \y{i }{i-1} & \y{i }{i } & 0 & \cdots & 0 & 0 \\
\y{i+1}{ 0} & \y{i+1}{ 1} & \cdots & \y{i+1}{i-1} & \y{i+1}{i } & 1 & \cdots & 0 & 0 \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\end{bmatrix}.\end{split}\]
src/gauss_jordan.c¶
30for (size_t ii = i + 1; ii < nitems; ii++) {
31 const double f = arr[ii * nitems + i];
32 for (size_t j = i; j < nitems; j++) {
33 arr[ii * nitems + j] -= f * arr[i * nitems + j];
34 }
35 for (size_t j = 0; j < i + 1; j++) {
36 inv[ii * nitems + j] -= f * inv[i * nitems + j];
37 }
38}
Backward substitution¶
After the forward sweep, we have
\[\begin{split}\begin{bmatrix}
1 & \x{ 0}{ 1} & \cdots & \x{ 0}{i-1} & \x{ 0}{i } & \x{ 0}{i+1} & \cdots & \x{ 0}{n-2} & \x{ 0}{n-1} \\
0 & 1 & \cdots & \x{ 1}{i-1} & \x{ 1}{i } & \x{ 1}{i+1} & \cdots & \x{ 1}{n-2} & \x{ 1}{n-1} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
0 & 0 & \cdots & 1 & \x{i-1}{i } & \x{i-1}{i+1} & \cdots & \x{i-1}{n-2} & \x{i-1}{n-1} \\
0 & 0 & \cdots & 0 & 1 & \x{i }{i+1} & \cdots & \x{i }{n-2} & \x{i }{n-1} \\
0 & 0 & \cdots & 0 & 0 & 1 & \cdots & \x{i+1}{n-2} & \x{i+1}{n-1} \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 1 & \x{n-2}{n-1} \\
0 & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
\y{ 0}{ 0} & 0 & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 \\
\y{ 1}{ 0} & \y{ 1}{ 1} & \cdots & 0 & 0 & 0 & \cdots & 0 & 0 \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots & & \vdots & \vdots \\
\y{i-1}{ 0} & \y{i-1}{ 1} & \cdots & \y{i-1}{i-1} & 0 & 0 & \cdots & 0 & 0 \\
\y{i }{ 0} & \y{i }{ 1} & \cdots & \y{i }{i-1} & \y{i }{i } & 0 & \cdots & 0 & 0 \\
\y{i+1}{ 0} & \y{i+1}{ 1} & \cdots & \y{i+1}{i-1} & \y{i+1}{i } & \y{i+1}{i+1} & \cdots & 0 & 0 \\
\vdots & \vdots & & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\y{n-2}{ 0} & \y{n-2}{ 1} & \cdots & \y{n-2}{i-1} & \y{n-2}{i } & \y{n-2}{i+1} & \cdots & \y{n-2}{n-2} & 0 \\
\y{n-1}{ 0} & \y{n-1}{ 1} & \cdots & \y{n-1}{i-1} & \y{n-1}{i } & \y{n-1}{i+1} & \cdots & \y{n-1}{n-2} & \y{n-1}{n-1} \\
\end{bmatrix}.\end{split}\]
Starting from the bottom (\(n - 1\)-th) row, the upper-triangular elements are eliminated:
src/gauss_jordan.c¶
43for (size_t ii = 0; ii < i; ii++) {
44 const double f = arr[ii * nitems + i];
45 arr[ii * nitems + i] = 0.;
46 for (size_t j = 0; j < nitems; j++) {
47 inv[ii * nitems + j] -= f * inv[i * nitems + j];
48 }
49}