Gaussian Elimination

Here we consider solving the following linear system with $M$ equations and $N$ unknowns:

$$\left\{\begin{array}{l} \sum_{j=1}^N a_{1j}x_j=a_{11}x_1+\cdots+a... ..._j+\cdots+a_{MN}x_N=b_M \cdots\cdots\cdots\cdots{\bf r}_M\\ \end{array}\right.$$ (84)

where the ith equation (or row) is labeled by ${\bf r}_i$ on the right hand side. This equation system can be expressed in matrix form as:

$${\bf A}{\bf x}=\left[\begin{array}{ccc} a_{11} & \cdots & a_{1N}\... ...ay}\right]= \left[\begin{array}{c} b_1\\ \vdots\\ b_M\end{array}\right]={\bf b}$$ (85)

where ${\bf A}$ is the coefficient matrix with $M$ rows and $N$ columns, and $a_{ij}$ is a component in the ith row and jth column.

When $M=N$ (the number of equations is the same as the number of unknowns), the methods of Gaussian elimination can be used to solve the equation ${\bf Ax}={\bf b}$. Specifically, we apply a sequence of elementary row operations listed below on both sides of the equations.

• Swap the order of two rows: ${\bf r}_i\leftrightarrow{\bf r}_j$;
• Multiply a row by a non-zero constant $c\ne 0$: $c{\bf r}_i\rightarrow {\bf r}_i$;
• Add a row multiplied by a non-zero constant $c\ne 0$ to another: ${\bf r}_i+c{\bf r}_j\rightarrow{\bf r}_i$ ($i\ne j$).

As the equations before and after each operation are equivalent, the equation system remains the same with the same solution.

These operations can be realized respectively by pre-multiplying the following matrices to both sides of the equations:

$${\bf T}_{ji}=\left[\begin{array}{ccccccc} 1 & & & & & & \\ & ... ...nd{array}\right] \begin{array}{c}\\ i \\ \\ j \\ \\ \end{array}$$ (86)

Now the equation system ${\bf Ax}={\bf b}$ can be solved by the following methods:

• Gaussian elimination with back substitution

  The solution of ${\bf A}{\bf x}={\bf b}$ can be found in two phases:

  • Forward elimination

    For each of the $N-1$ columns ${\bf c}_j$ ( $j=1,\cdots,N-1$), eliminate all components $a_{ij}$ ( $i=j+1,\cdots,N$) below $a_{jj}$ on the diagonal:

    $${\bf T}_{ij}(-a_{ij}/a_{jj}):\;\;\;{\bf r}_i-\left(\frac{a_{ij}}{a_{jj}}\right){\bf r}_j \rightarrow {\bf r}_i$$ (87)

    At the end of this elimination process, the coefficient matrix becomes an upper triangular matrix with all components below the diagonal eliminated, $a_{ij}=0$ ($i>j$). Now the left-hand side of the ith equation has $N-i+1$ terms ( $i=1,\cdots,N$):

    $$\left[\begin{array}{ccccc} a_{11}& a_{12}& \cdots & a_{1,N-1} & a... ...] =\left[\begin{array}{l} b_1\\ b_2\\ \vdots \\ b_{N-1}\\ b_N\end{array}\right]$$ (88)

    This form of the coefficient matrix is called row echelon form (ref).

    As $O(N)$ multiplications are needed for each row operation ${\bf T}_{ij}$ ( $j=1,\cdots,N-1,\;\;i=j+1,\cdots,N$), the computational complexity of the forward elimination is $O(N^3)$.

  • Back substitution

    This is a bottom-up process. Starting with the Nth equation in the last row, we solve ith equation for $x_i$ for all $i=N,\cdots,1$:

    $$x_i=\frac{1}{a_{ii}}\left( b_i-\sum_{j=i+1}^N a_{ij} x_j\right)$$ (89)

  Example:

  $$\left\{ \begin{array}{rrrr} x_1+&2x_2+&3x_3&=14\cdots\cdots\cdots... ...}_2\\ 3x_1+&x_2+&2x_3&=11\cdots\cdots\cdots\cdots{\bf r}_3 \end{array} \right.$$ (90)

  • Forward elimination:

    We first consider column $j=1$. By operation:

    $${\bf T}_{21}(-a_{21}/a_{11})={\bf T}_{21}(-2): \:\;\;\; {\bf r}_2-2{\bf r}_1\rightarrow {\bf r}_2$$ (91)

    ${\bf r}_1$ is multiplied by $-a_{21}/a_{11}=-2/1=-2$:

    $$-2x_1-4x_2-6x_3=-28$$ (92)

    and added to ${\bf r}_2$:

    $$-x_2-5x_3=-17\;\;\;\;\;\; \;\;\;\;\;\;\;x_2+5x_3=17$$ (93)

    Also, by

    $${\bf T}_{31}(-a+{31}/a_{11})={\bf r}_{31}(-3):\;\;\; {\bf r}_3-3{\bf r}_1\rightarrow {\bf r}_3$$ (94)

    ${\bf r}_1$ is multiplied by $-a_{31}/a_{11}=-3/1=-3$:

    $$-3x_1-6x_2-9x_3=-42$$ (95)

    and added to ${\bf r}_3$:

    $$-5x_2-7x_3=-31\;\;\;\;\;\; \;\;\;\;\;\;5x_2+7x_3=31$$ (96)

    Now we have

    $$\left\{ \begin{array}{rrrr} x_1+&2x_2+&3x_3&=14\\ &x_2+&5x_3&=17\\ &5x_2+&7x_3&=31 \end{array} \right.$$ (97)

    We next consider column $j=2$. By

    $${\bf T}_{32}(-a_{32}/a_{22})={\bf T}_{32}(-5):\;\;\; {\bf r}_3-5{\bf r}_2\rightarrow {\bf r}_3$$ (98)

    ${\bf r}_2$ is multiplied by $-a_{32}/a_{22}=-5/1=-5$:

    $$-5x_2-25x_3=-85$$ (99)

    and added to ${\bf r}_3$:

    $$-18x_3=-54\;\;\;\;\;\; \;\;\;\;\;\;18x_3=54$$ (100)

    Now we have

    $$\left\{ \begin{array}{rrrr} x_1+&2x_2+&3x_3&=14\\ &x_2+&5x_3&=17\\ & &18x_3&=54 \end{array} \right.$$ (101)

  • Back substitution:

    Solving ${\bf r}_3$ we get $x_3=3$, substituting it into ${\bf r}_2$ we get $x_2=2$, substituting both $x_3=3$ and $x_2=2$ into ${\bf r}_1$ we get $x_1=1$.

• Gauss-Jordan elimination

  The method of Gauss-Jordan elimination can be used to find the inverse ${\bf A}^{-1}$ of the coefficient matrix ${\bf A}$. Once ${\bf A}^{-1}$ is available, we can simultaneously solve multiple linear systems if they are all based on the same coefficient matrix ${\bf A}$:

  $${\bf A}{\bf x}_i={\bf b}_i,\;\;\;\;\;\;\;\;\; {\bf x}_i={\bf A}^{-1}{\bf b}_i,\;\;\;\;\;\;\;\;\;(i=1,\cdots,m)$$ (102)

  To do so, we first construct the augmented matrix:

  $$\left[\begin{array}{ccccc}{\bf A}&{\bf b}_1&\cdots&{\bf b}_m&{\bf I}\end{array}\right]$$ (103)

  Pre-multiply ${\bf A}^{-1}$ (assumed available) to each of its components we would get
  

  $${\bf A}^{-1} \left[\begin{array}{ccccc}{\bf A}&{\bf b}_1&\cdots&{\bf b}_m&{\bf I}\end{array}\right] = \left[\begin{array}{ccccc}{\bf A}^{-1}{\bf A}&{\bf A}^{-1}{\bf b}_1&\cdots&{\bf A}^{-1}{\bf b}_m&{\bf A}^{-1}{\bf I}\end{array}\right]$$

  $$= \left[\begin{array}{ccccc}{\bf I}&{\bf x}_1&\cdots&{\bf x}_m&{\bf A}^{-1}\end{array}\right]$$ (104)

  
  However, as ${\bf A}^{-1}$ is unknown, the pre-multiplication by ${\bf A}^{-1}$ is actually carried be applying a sequence of elementary row operations ${\bf T}^{(i)}$ to the augmented matrix until eventually ${\bf A}$ on the left of the augmented matrix is converted to ${\bf I}$.

  $${\bf T}^{(k)}\cdots{\bf T}^{(1)}\left[\begin{array}{ccccc}{\bf A}... ...{array}{ccccc}{\bf I}&{\bf x}_1&\cdots&{\bf x}_m&{\bf A}^{-1}\end{array}\right]$$ (105)

  This conversion is done by the Gaussian elimination to eliminate all components below the diagonal, followed by a second phase to eliminate all components above the diagonal:
  • Convert ${\bf A}$ into the row echelon form by eliminating all $a_{ij}$ in equation ${\bf r}_i$ $(i=j+1,\cdots,N$) below the diagonal ( $j=1,\cdots,N-1$, left to right):

    $${\bf T}_{ij}(-a_{ij}/a_{jj}): \;\;\;{\bf r}_i-\frac{a_{ij}}{a_{jj}}{\bf r}_j \rightarrow {\bf r}_i$$ (106)

    By this phase of elimination all elements below the diagonal ($a_{ij}=0$ for all $i>j$) are eliminated, and ${\bf A}$ becomes an upper triangular matrix.

  • Convert the ref obtained into a diagonal matrix, called reduced row echelon form (rref), by eliminating all $a_{ij}$ in equation ${\bf r}_i$ ( $i=j-1,\cdots,1$) above the diagonal ( $j=N,\cdots,2$, right to left):

    $${\bf T}_{ij}(-a_{ij}/a_{jj}): \;\;\;{\bf r}_i-\frac{a_{ij}}{a_{jj}}{\bf r}_j \rightarrow {\bf r}_i$$ (107)

    By this phase of elimination all elements above the diagonal ($a_{ij}=0$ for all $i<j$) are eliminated, and ${\bf A}$ becomes a diagonal matrix, which can be easily turned into an identity matrix ${\bf I}$.

  Alternatively, we could combine the two phases to simply remove all components $a_{ij}$ in the jth column ( $i,j=1,\cdots,N,\;\;\;i\ne j$) by

  $${\bf T}_{ij}(-a_{ij}/a_{jj}): {\bf r}_i-\frac{a_{ij}}{a_{jj}}{\bf r}_j \rightarrow {\bf r}_i$$ (108)

  Note that as the process of the elementary row operations above is applied to all rows in the augmented matrix, it also generates at the same time the desired solutions ${\bf x}_i$ as well as the inverse matrix ${\bf A}^{-1}$.

  As each row operation ${\bf T}_{ij}$ ( $i,j=1,\cdots,N,\;\;i\ne j$) requires $O(N)$ multiplications, the total complexity for solving the linear system ${\bf Ax}={\bf b}$, the same for computing ${\bf A}^{-1}$, is $O(N^3)$.

  The Matlab code below implements the Gauss-Jordan method for finding the inverse of a given square matrix ${\bf A}$:

  function Ainv=GaussJordanInv(A) % Find inverse by Gaussian-Jordan elimination
      [m n]=size(A);
      n=2*m;
      A=[A eye(m)];           % construct augmented matrix
      Ainv=zeros(m);
      v=[1:m];                % row permutation record
      for k=1:m               % for each row of A  
          i=v(k);             % current row index
          p=abs(A(i,i));    
          ip=k;    
          for l=k+1:m         % find max in ith column below v(k)th row
              if abs(A(v(l), v(k)))>p          
                  p=abs(A(v(l), v(k)));            
                  ip=l;            
              end
          end
          temp=v(k);          % Switch the kth element and imaxth element in r-vector
          v(k)=v(ip);
          v(ip)=temp;
          for i=1:m           % Eliminate off-diagonal elements in v(k)th column 
              if i ~= k
                  c=A(v(i),k)/A(v(k),k);
                  for j=k:n
                      A(v(i),j)=A(v(i),j)-c*A(v(k),j);
                  end
              end
          end
      end
      for k=1:m
          i=v(k);
          [p ip]=max(abs(A(i,1:m)));
          for j=1:m
              Ainv(k,j)=A(i,m+j)/A(i,ip);
          end
      end
  end
  

  Example

  $${\bf A}{\bf x}=\left[ \begin{array}{ccc} 1 & 2 & 3\\ 2 & 3 & 1\\ ... ...end{array}\right] =\left[\begin{array}{c}14\\ 11\\ 11\end{array}\right]={\bf b}$$ (109)

  $$\left[\begin{array}{ccc}{\bf A}&{\bf b}&{\bf I}\end{array}\right]... ...0\\ 2 & 3 & 1 & 11 & 0 & 1 & 0\\ 3 & 1 & 2 & 11 & 0 & 0 & 1 \end{array}\right]$$ (110)

  • $j=1,\;\;i=2$: ${\bf T}_{21}(-2)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 2 & 3 & 14 & 1 & 0 ... ... 0 & -1 & -5 & -17 & -2 & 1 & 0\\ 3 & 1 & 2 & 11 & 0 & 0 & 1\end{array}\right]$$ (111)

  • $j=1,\;\;i=3$: ${\bf T}_{31}(-3)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 2 & 3 & 14 & 1 & 0 ... ... 0 & -1 & -5 & -17 & -2 & 1 & 0\\ 0 &-5 &-7 &-31 &-3 & 0 & 1\end{array}\right]$$ (112)

  • $j=2,\;\;i=3$: ${\bf T}_{32}(-5)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 2 & 3 & 14 & 1 & 0 ... ... 0 & -1 & -5 & -17 & -2 & 1 & 0\\ 0 & 0 & 18 &54 & 7 &-5 & 1\end{array}\right]$$ (113)

  • ${\bf T}_3(1/18)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 2 & 3 & 14 & 1 & 0 ... ...1 & -5 & -17 & -2 & 1 & 0\\ 0 & 0 & 1 & 3 & 7/18 &-5/18 &1/18\end{array}\right]$$ (114)

  • $j=3,\;\;i=2$: ${\bf T}_{23}(5)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 2 & 3 & 14 & 1 & 0 ... ...-2 & -1/18 &-7/18 & 5/18\\ 0 & 0 & 1 & 3 & 7/18 &-5/18 & 1/18\end{array}\right]$$ (115)

  • $j=3,\;\;i=1$: ${\bf T}_{13}(-3)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 2 & 0 & 5 & -3/18 &... ...-2 & -1/18 &-7/18 & 5/18\\ 0 & 0 & 1 & 3 & 7/18 &-5/18 & 1/18\end{array}\right]$$ (116)

  • $j=2,\;\;i=1$: ${\bf T}_{12}(2)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 0 & 0 & 1 & -5/18 &... ...-2 & -1/18 &-7/18 & 5/18\\ 0 & 0 & 1 & 3 & 7/18 &-5/18 & 1/18\end{array}\right]$$ (117)

  • ${\bf T}_2(-1)$

    $$\left[ \begin{array}{rrr\vert r\vert rrr} 1 & 0 & 0 & 1 & -5/18 &... ...& 2 & 1/18 & 7/18 &-5/18\\ 0 & 0 & 1 & 3 & 7/18 &-5/18 & 1/18\end{array}\right]$$ (118)

  Now we get both ${\bf x}$ and ${\bf A}^{-1}$ at the same time:

  $${\bf x}=\left[\begin{array}{c}1\\ 2\\ 3\end{array}\right],\;\;\;\... ...=\frac{1}{18}\left[\begin{array}{rrr}-5&1&7\\ 1&7&-5\\ 7&-5&1\end{array}\right]$$ (119)

Some issues related to Gaussian elimination methods:

• Pivoting

  In Gaussian elimination we always need to find the scaling factor by a division $a_{ij}/a_{jj}$, which cannot be carried out if $a_{jj}=0$. Also, even if $a_{jj}\ne 0$, its value may be too close to zero to be used in the division, as any error in $a_{ij}$ will be amplified when divided by $a_{jj}$ close to zero. To avoid such problems, we will always find the component with the largest absolute value in the current jth column. If $\vert a_{kj}\vert$ in the kth row is the maximum of all $\vert a_{ij}\vert\;i=1,\cdots,M$, then we use operation ${\bf T}_{kj}$ to swap row ${\bf r}_k$ with the current jth row ${\bf r}_j$. This new component $a_{jj}$ with the largest absolute value is called the pivot, and the process of finding such a component to use is called pivoting.

  It is possible that the pivot component in one of the columns is zero, indicating that the coefficient matrix is singular, its inverse does not exist.

• Computational complexity

  The complexity of the Gaussian elimination can be found based on the total number of operations:

  $${\bf r}_i-\left(\frac{a_{ij}}{a_{jj}}\right){\bf r}_j \rightarrow {\bf r}_i, \;\;\;\;\;\;(j=1,\cdots,N-1,\;i=j+1,\cdots,N)$$ (120)

  The division $a_{ij}/a_{jj}$ is carried for each of the components below the diagonal $a_{ij}$ for all $i>j$. Therefore the number of divisions is

  $$(N-1)+(N-2)+\cdots+1=N(N-1)/2=O(N^2)$$ (121)

  Having found the scaling factor $a_{ij}/a_{jj}$, we need to use it to multiply through out the entire row, therefore the number of multiplications will be $O(N^3)$. This is the overall complexity for the Gaussian elimination methods.