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