Numerical Analysis Home

## Systems of Linear Equations: Gaussian Elimination

A system of linear equations with $$n$$ variables $$x_1,\ldots,x_n$$ and $$n$$ equations can be written as follows: $\begin{eqnarray*} {\color{blue}\begin{array}{ccccccccc} a_{11}x_1&+&a_{12}x_2&+&\cdots &+&a_{1n}x_n&=&b_1\\ a_{21}x_1&+&a_{22}x_2&+&\cdots &+&a_{2n}x_n&=&b_2\\ \vdots&&\vdots&& &&\vdots&&\vdots\\ a_{n1}x_1&+&a_{n2}x_2&+&\cdots &+&a_{nn}x_n&=&b_n. \end{array}} \end{eqnarray*}$ The above linear system is equivalent to the matrix equation $$A\overrightarrow{x}=\overrightarrow{b}$$, where $A=\left[\begin{array}{cccc} a_{11}&a_{12}&\cdots &a_{1n}\\ a_{21}&a_{22}&\cdots &a_{2n}\\ \vdots&\vdots& &\vdots\\ a_{nn}&a_{n2}&\cdots &a_{nn} \end{array}\right], \overrightarrow{x}=\left[\begin{array}{c} x_1\\x_2\\ \vdots\\x_n\end{array} \right], \mbox{ and } \overrightarrow{b}=\left[\begin{array}{c} b_1\\b_2\\ \vdots\\b_n \end{array} \right].$ $$A$$ is called the coefficient matrix and $$[A\; \overrightarrow{b}]$$ is called the augmented matrix.

Example. $\begin{array}{rcrcrcr} 2x_1&+&x_2&&&=&7\\ -3x_1&&&+&x_3&=&-8\\ &&x_2&+&2x_3&=&-3 \end{array}.$ The above linear system is equivalent to $$A\overrightarrow{x}=\overrightarrow{b}$$, where $$A=\left[\begin{array}{rrr}2&1&0\\ -3&0&1\\ 0&1&2\end{array} \right]$$, $$\overrightarrow{x}=\left[\begin{array}{r}x_1\\x_2\\x_3 \end{array} \right]$$, $$\overrightarrow{b}=\left[\begin{array}{r}7\\-8\\-3 \end{array} \right]$$, and the augmented matrix is $$[A\; \overrightarrow{b}]= \left[\begin{array}{rrr|r}2&1&0&7\\ -3&0&1&-8\\ 0&1&2&-3\end{array} \right]$$.

Throughout this section assume
(1) $$A\overrightarrow{x}=\overrightarrow{b}$$ has a unique solution and
(2) $$a_{ii}\neq 0$$ for $$i=1,\ldots, n$$.

Techniques for solving $$A\overrightarrow{x}=\overrightarrow{b}$$:
There are multiple ways to solve $$A\overrightarrow{x}=\overrightarrow{b}$$. First we learn Gaussian elimination as illustrated in the following example:

Example. By Gaussian elimination find the unique solution of the following system of equations. $\begin{array}{rcrcrcr} 2x_1 &- &4x_2 &+ &2x_3 &=&2\\ x_1 &+ &x_2 &+ &2x_3 &=&0\\ -3x_1 &+ &8x_2 &- &3x_3 &=&-3 \end{array}$ Solution. \begin{align*} &\begin{array}{rrcrcrcr} E1: & 2x_1 &- &4x_2 &+ &2x_3 &=&2\\ E2: & x_1 &+ &x_2 &+ &2x_3 &=&0\\ E3: & -3x_1 &+ &8x_2 &- &3x_3 &=&-3\end{array}\\ &\\ \xrightarrow{\substack{-\frac{1}{2}E1+E2\\ \frac{3}{2}E1+E3}}\;\;\;\; & \begin{array}{rrcrcrcr} E1: & 2x_1 &- &4x_2 &+ &2x_3 &=&2\\ E2: & & &3x_2 &+ &x_3 &=&-1\\ E3: & & &2x_2 & & &=&0\end{array}\\ &\\ \xrightarrow{-\frac{2}{3}E2+E3}\;\;\;\; & \begin{array}{rrcrcrcr} E1: & 2x_1 &- &4x_2 &+ &2x_3 &=&2\\ E2: & & &3x_2 &+ &x_3 &=&-1\\ E3: & & & &- &\frac{2}{3}x_3 &=&\frac{2}{3}\end{array}\\ \end{align*} Note that if $$a_{ii}=0$$ for some $$i$$ and $$a_{ji}\neq 0$$ for some $$j > i$$, then we interchange $$Ei$$ and $$Ej$$. This is called partial pivoting. If $$a_{ji}=0$$ for all $$j\geq i$$, then there is no unique solution.
Now we do backward substitutions: \begin{align*} E3\implies &x_3=\frac{2/3}{-2/3}=-1\\ E2\implies &x_2=\frac{1}{3}(-1-x_3)=\frac{1}{3}(-1+1)=0\\ E1\implies &x_1=\frac{1}{2}(2+4x_2-2x_3)=\frac{1}{2}(2+4\cdot0-2(-1))=2 \end{align*} Thus the unique solution is $$(x_1,x_2,x_3)=\left(2,0,-1\right)$$.

Gaussian elimination steps:

1. For each $$i=1,\ldots,n-1$$, $$Ej=Ej-\frac{a_{ji}}{a_{ii}}Ei$$, $$j=i+1,\ldots,n$$.
2. $$x_n=\frac{b_n}{a_{nn}}$$ and $$x_i=\frac{1}{a_{ii}}\Bigg(b_i-\displaystyle\sum_{j=i+1}^n a_{ij}x_j\Bigg)$$, $$i=n-1,n-2,\ldots,1$$. (backward substitution)

Operation Counts:
In step 1, for each $$i=1,\ldots,n-1$$, we do $$(n-i)$$ divisions $$\frac{a_{ji}}{a_{ii}}$$, $$j=i+1,\ldots,n$$. Then for each $$j=i+1,\ldots,n$$, to get $$\frac{a_{ji}}{a_{ii}}Ei$$, we do $$n-i+1$$ multiplication when we multiply $$n-i+1$$ numbers $$a_{i,i+1},\ldots,a_{i,n},b_i$$ in $$Ei$$ by $$\frac{a_{ji}}{a_{ii}}$$. So for each $$i=1,\ldots,n-1$$, the total number of multiplications/divisions is $(n-i)+(n-i)(n-i+1)=(n-i)(n-i+2).$ For a fixed $$i=1,\ldots,n-1$$, we do $$(n-i+1)$$ subtraction for $$Ej=Ej-\frac{a_{ji}}{a_{ii}}Ei$$ for each $$j=i+1,\ldots,n$$, i.e., total $$(n-i)(n-i+1)$$ subtractions.
So the total number of multiplications/divisions is $\sum_{i=1}^{n-1}(n-i)(n-i+2)=\frac{2n^3+3n^2-5n}{6}.$ The total number of subtractions is $\sum_{i=1}^{n-1}(n-i)(n-i+1)=\frac{n^3-n}{3}.$ Note that in a computer the time for a multiplication/division is greater than that of an addition/subtraction. So the total number of operations in step 1 is $$\displaystyle \frac{2n^3+3n^2-5n}{6}+\frac{n^3-n}{3}=\frac{4n^3+3n^2-7n}{6}$$ which is $$O(n^3)$$. Similarly we can show the total number of operations in step 2 is $$\displaystyle \frac{n^2+n}{2}+\frac{n^2-n}{2}=n^2$$. Thus total number of operations in Gaussian elimination is $$\displaystyle \frac{4n^3+3n^2-7n}{6}+n^2=\frac{4n^3+9n^2-7n}{6}$$ which is $$O(n^3)$$.

$$LU$$-factorization:
If $$M=[A\; \overrightarrow{b}]$$, the first step transforms it into $$M'=[U\; \overrightarrow{b'}]$$ where $$U$$ is an upper triangular matrix with nonzero diagonals. This is obtained by premultiplying $$M$$ by elementary matrices. For example, premultiplying $$M$$ by the elementary matrix $$E_1=\left[\begin{array}{rrr}1&0&0\\ -\frac{1}{2}&1&0\\ 0&0&1\end{array} \right]$$ gives $E_1M=\left[\begin{array}{rrr}1&0&0\\ -\frac{1}{2}&1&0\\ 0&0&1\end{array} \right] \left[\begin{array}{rrr|r}2&-4&2&2\\ 1&1&2&0\\ -3&8&-3&-3\end{array} \right]= \left[\begin{array}{rrr|r}2&-4&2&2\\ 0&3&1&-1\\ -3&8&-3&-3\end{array} \right].$ Similarly using the elementary matrices $$E_2=\left[\begin{array}{rrr}1&0&0\\ 0&1&0\\ \frac{3}{2}&0&1\end{array} \right]$$ and $$E_3=\left[\begin{array}{rrr}1&0&0\\ 0&1&0\\ 0&-\frac{2}{3}&1\end{array} \right]$$, we get $E_3E_2E_1M=E_3E_2E_1[A|\overrightarrow{b}]= \left[\begin{array}{rrr|r}2 & -4 & 2 & 2\\ 0 & 3 & 1 & -1\\ 0 & 0& -2/3 & 2/3\end{array} \right]=[U|\overrightarrow{b'}],$ which implies $$E_3E_2E_1A=U$$ an upper-triangular matrix with nonzero diagonals. Now $E_3E_2E_1A=U\implies A=(E_3E_2E_1)^{-1}U=LU= \left[\begin{array}{rrr} 1 & 0 & 0\\ 1/2 & 1 & 0\\ -3/2 & 2/3 & 1\end{array} \right] \left[\begin{array}{rrr} 2 & -4 & 2\\ 0 & 3 & 1\\ 0 & 0& -2/3 \end{array} \right],$ where $$L=[l_{ij}]$$ is a lower-triangular matrix with $$l_{ji}=a_{ji}^{(i)}/a_{ii}^{(i)}$$. So we have the $$LU$$-factorization of $$A$$: $$A=LU$$. Then $$A\overrightarrow{x}=\overrightarrow{b} \implies LU\overrightarrow{x}=\overrightarrow{b}$$. Now we do steps to find $$\overrightarrow{x}$$:

1. From $$L\overrightarrow{y}=\overrightarrow{b}$$ solve for $$\overrightarrow{y}$$ by forward substitution

2. From $$U\overrightarrow{x}=\overrightarrow{y}$$ solve for $$\overrightarrow{x}$$ by backward substitution

Example. Solve $$A\overrightarrow{x}=[2\;\; 0\; -3]^T$$, where $A=\left[\begin{array}{rrr} 1 & 0 & 0\\ 1/2 & 1 & 0\\ -3/2 & 2/3 & 1\end{array} \right] \left[\begin{array}{rrr} 2 & -4 & 2\\ 0 & 3 & 1\\ 0 & 0& -2/3 \end{array} \right].$ Solution. With $$A=LU$$, we solve for $$\overrightarrow{y}$$ from $$L\overrightarrow{y}=\overrightarrow{b}$$ by forward substitution: \begin{align*} L\overrightarrow{y}=\overrightarrow{b} &\implies \left[\begin{array}{rrr} 1 & 0 & 0\\ 1/2 & 1 & 0\\ -3/2 & 2/3 & 1\end{array} \right] \left[\begin{array}{r}y_1\\y_2\\y_3 \end{array} \right]=\left[\begin{array}{r}2\\0\\-3 \end{array} \right]\\ &\implies \left\lbrace\begin{array}{rcrcrcr} y_1 & & & & &=&2\\ y_1/2 &+ &y_2 & & &=&0\\ -3y_1/2 &+ &2y_2/3 &+ &y_3 &=&-3 \end{array} \right.\\ &\implies \left\lbrace\begin{array}{rl} y_1 &=2\\ y_2 &=0-y_1/2=0-2/2=-1\\ y_3 &=-3+3y_1/2-2y_2/3=-3+3\cdot 2/2-2(-1)/3=2/3 \end{array} \right. \end{align*} Now we solve for $$\overrightarrow{x}$$ from $$U\overrightarrow{x}=\overrightarrow{y}$$ by forward substitution: \begin{align*} U\overrightarrow{x}=\overrightarrow{y} &\implies \left[\begin{array}{rrr} 2 & -4 & 2\\ 0 & 3 & 1\\ 0 & 0& -2/3 \end{array} \right] \left[\begin{array}{r}x_1\\x_2\\x_3 \end{array} \right]=\left[\begin{array}{r}2\\-1\\2/3 \end{array} \right]\\ &\implies \left\lbrace\begin{array}{rcrcrcr} 2x_1 &- &4x_2 &+ &2x_3 &=&2\\ & &3x_2 &+ &x_3 &=&-1\\ & & &- &2x_3/3 &=&2/3 \end{array} \right.\\ &\implies \left\lbrace\begin{array}{rl} x_3 &=\frac{2/3}{-2/3}=-1\\ x_2 &=\frac{1}{3}(-1-x_3)=\frac{1}{3}(-1+1)=0\\ x_1 &=\frac{1}{2}(2+4x_2-2x_3)=\frac{1}{2}(2+4\cdot0-2(-1))=2 \end{array} \right. \end{align*} Thus the unique solution is $$(x_1,x_2,x_3)=\left(2,0,-1\right)$$.

Advantage of $$LU$$-factorization:
If we need to solve $$A\overrightarrow{x}=\overrightarrow{b_1}$$, $$A\overrightarrow{x}=\overrightarrow{b_2}$$,...,$$A\overrightarrow{x}=\overrightarrow{b_k}$$, then Gaussian elimination solves each problem separately. But once we know the $$LU$$-factorization of $$A=LU$$, the solution $$\overrightarrow{x_i}$$ of $$A\overrightarrow{x}=\overrightarrow{b_i}$$ is obtained just from forward and backward substitution in $$L\overrightarrow{y_i}=\overrightarrow{b_i} \text{ and } U\overrightarrow{x_i}=\overrightarrow{y_i}$$ respectively.

Last edited