Linear Systems of Equations and Gauss Elimination

Linear Systems of Equations (LSE’s)

We are going to handle equations of type Ax = b in this chapter. Let’s see how linear systems of equations can be represented as a matrix and a vector. Let our LSE in $n$ unknowns ($x_i$) and $m$ equations (the rows) be as in the following:

$$\begin{aligned} 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 \\ a_{m1}x_1 + a_{m2}x_2 + \cdots + a_{mn}x_n &= b_m, \end{aligned}$$

This suspiciously looks similar to a matrix vector multiplication! If it is not clear to you how, notice that each row $i$ of coefficients $a_{i1}\ldots a_{in}$ gets multiplied with our unknown variables $x_i$. Seeing the matrix-vector picture might help:

$$A\textbf{x} = \textbf{b} : \quad \underbrace{ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} \end{bmatrix} }_{A,\, m \times n} \; \underbrace{ \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix} }_{\textbf{x} \in \mathbb{R}^n} = \underbrace{ \begin{pmatrix} b_1 \\ b_2 \\ \vdots \\ b_m \end{pmatrix} }_{\textbf{b} \in \mathbb{R}^m}.$$

This way of writing linear systems of equations lets us argue about the properties of the equations using the properties of the matrix $A$. For example, we know that $A\textbf{x}$ calculates a linear combination of $A$‘s columns where the coefficients are the elements of the vector $\textbf{x}$. This means we are trying to write $\textbf{b}$ as a linear combination of $A$‘s columns. This is possible if and only if $\textbf{b}\in C(A)$, ie if $\textbf{b}$ is in the column space of the matrix $A$.

Gauss Elimination

We can represent linear systems of equations in form of $A\textbf{x} = \textbf{b}$ using matrix vector multiplication. But how about solving the LSE? Our first algorithm is Gauss Elimination, which actually is the classic algorithm for solving linear systems of equations. It consists of two steps: elimination and back substitution. For this week we assume that our coefficient matrix is $m \times m$, which means that we have a linear system of equations in $m$ equations and $m$ variables. This is the nice case, we will deal with different number of equations and variables later.

If you find a step difficult to understand, go to the subsection below 2.3 Example: Solving $A\textbf{x}=\textbf{b}$ using Gaussian Elimination to see the theory in practice.

Elimination

We want to generate an upper triangular matrix $U$ using our initial matrix $A$ and the vector $\textbf{b}$. We want to transform the equation $A\textbf{x} = \textbf{b}$ into another equation $U\textbf{x} = \textbf{c}$ where $U$ is an upper triangular matrix and $\textbf{x}$, the solutions are the same for both systems. In other words we want to create equivalent systems.

It is going to become obvious why we want an upper triangular matrix when we handle back substitution.

There are two main row operations to perform the elimination:

• row subtraction
• row exchange

These operations do not change the solution of the system when applied on both $A$ and $\textbf{b}$. This is why we usually concatenate $A$ and $\textbf{b}$ when we perform these elimination steps as follows:

$$\begin{bmatrix} & \vert & \\ A& \vert & \textbf{b} \\ & | & \end{bmatrix}$$

Both operations can be represented as linear transformations!

In the following there are the matrices that represent these linear transformations in the general case:

$$E_{ij} = \begin{bmatrix} \ddots & & & & \\ & 1 & & & \\ & & \ddots & & \\ & -c & & 1 & \\ & & & & \ddots \\ & \substack{\uparrow \\ j} & & \substack{\uparrow \\ i} & \end{bmatrix} \begin{matrix} \vphantom{\ddots} \\ \leftarrow j \\ \vphantom{\ddots} \\ \leftarrow i \\ \vphantom{\ddots} \\ \vphantom{\substack{\uparrow \\ j}} \end{matrix}$$

(a) Elimination: all diagonal entries are $1$, row $i$, col $j$ is $-c$ and everything else is 0

$$P_{jk} = \begin{bmatrix} \ddots & & & & \\ & 0 & & 1 & \\ & & \ddots & & \\ & 1 & & 0 & \\ & & & & \ddots \\ & \substack{\uparrow \\ j} & & \substack{\uparrow \\ k} & \end{bmatrix} \begin{matrix} \vphantom{\ddots} \\ \leftarrow j \\ \vphantom{\ddots} \\ \leftarrow k \\ \vphantom{\ddots} \\ \vphantom{\substack{\uparrow \\ j}} \end{matrix}$$

(b) Permutation: all diagonal entries except $j$th and $k$th are equal to 1. $j$th and $k$th diagonal entries are $0$ and row $j$ col $k$ and row $k$ col $j$ are 1.

The permutation matrix $P_{jk}$ swaps rows $j$ and $k$ and $E_{ij}$ represents the elimination matrix. When multiplied with the original matrix $A$ as $E_{ij}A$, it subtracts $c$ times row $j$ from row $i$. To better convince yourself to this you can think of the matrix multiplication $E_{ij}A$ as linear combinations of the rows of $A$ using the row vectors of $E_{ij}$ as the coefficients. The $i$-th row of $E_{ij}A$ is $-c$ times the $j$th row of $A$ plus $1$ times the $i$th row of $A$.

When do we use the row operations?

We try to generate an upper triangular matrix with nonzero elements on the diagonal. So if we have a nonzero element under the diagonal, we subtract one row from the other to get rid of it. We exchange rows when we get a zero element on the diagonal, $i.e.$ when we have a pivot that is equal to 0.

$$\boxed{ \text{The diagonal entry of the current column is called the \textbf{pivot}.} }$$

After completing the elimination we are left with the equation $U\textbf{x} = \textbf{c}$ that has the same solutions as our initial equation $A\textbf{x} = \textbf{b}$. We are going to prove this but first let’s continue solving the LSE with back substitution:

Back Substitution

After we get to an upper triangular matrix, we have our equation suitable for back substitution. This step is not different than what you would do in high school to solve a linear system of equations. Let’s see it in a $3 \times 3$ example. Here we name the matrix $U$ and the right hand side vector $\textbf{c}$ to stress that they are transformed and not the same as the initial $A$ and $\textbf{b}$ at this point. As mentioned in previous paragraph, the solution is still the same.

$$\underbrace{\begin{bmatrix} 3 & 2 & 1 \\ \color{green}0 & 4 & 3 \\ \color{green}0 & \color{green}0 & 5 \end{bmatrix}}_{U\in \R^{m\times m}} \underbrace{\begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}}_{\textbf{x}\in \R^m} = \underbrace{\begin{bmatrix} 11 \\ 13 \\ 15 \end{bmatrix}}_{\textbf{c}\in \R^m}$$

And the back substitution in the direction of the arrow on the right:

$$\left. \begin{array}{c|r|r|r} \text{equation} & \text{before substitution} & \text{after substitution} & \text{solution} \\ \hline 1 & 3x_1 + 2x_2 + x_3 = 11 & 3x_1 + 5 = 11 & x_1 = 2 \\ 2 & 4x_2 + 3x_3 = 13 & 4x_2 + 9 = 13 & x_2 = 1 \\ 3 & 5x_3 = 15 & & x_3 = 3 \end{array} \quad \right\uparrow$$

Using the variables, for which we already have calculated the value, we can calculate the values of other variables. A way to formalize this is:

$$x_i = \frac{c_i - \sum_{j=i+1}^{m} u_{ij} x_j }{u_{ii}}$$

Since we know all $x_j$ here, we can calculate $x_i$. Note that this only works if the diagonal entry $u_{ii}$ is not zero.

Example: Solving $A\textbf{x} = \textbf{b}$ using Gauss Elimination

We are given the system of linear equations represented as:

$$\begin{bmatrix} 0 & 2 & 1 & \vert & 1 \\ 2 & -3 & 4 & \vert & -2 \\ 1 & 1 & -2 & \vert & 0 \end{bmatrix}$$

We wish to solve the system $A\textbf{x} = \textbf{b}$. First, we will perform Gaussian elimination by applying elimination matrices and a permutation matrix.

• Step 1: Apply a Permutation Matrix to Swap Rows

  Since the first pivot is zero, we swap row 1 with row 2 using a permutation matrix $P$:

  $$P = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix}$$

  After applying the permutation, the augmented matrix becomes:

  $$\begin{bmatrix} 2 & -3 & 4 & \vert & -2 \\ 0 & 2 & 1 & \vert & 1 \\ 1 & 1 & -2 & \vert & 0 \end{bmatrix}$$

• Step 2: First Elimination Matrix $E_{31}$

  We eliminate the entry below the first pivot (2 in the $(1,1)$ position) by subtracting half of row 1 from row 3. The elimination matrix $E_{31}$ is:

  $$E_{31} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -\frac{1}{2} & 0 & 1 \end{bmatrix}$$

  After applying $E_{31}$, the matrix becomes:

  $$\begin{bmatrix} 2 & -3 & 4 & \vert & -2 \\ 0 & 2 & 1 & \vert & 1 \\ 0 & \frac{5}{2} & -4 & \vert & 1 \end{bmatrix}$$

• Step 3: Second Elimination Matrix $E_{32}$

  We now eliminate the entry below the second pivot (2 in the $(2,2)$ position) by subtracting $\frac{5}{4}$ of row 2 from row 3. The elimination matrix $E_{32}$ is:

  $$E_{32} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & -\frac{5}{4} & 1 \end{bmatrix}$$

  After applying $E_{32}$, the augmented matrix becomes:

  $$\underbrace{\begin{bmatrix} 2 & -3 & 4 & \vert & -2 \\ 0 & 2 & 1 & \vert & 1 \\ 0 & 0 & -\frac{21}{4} & \vert & -\frac{1}{4} \end{bmatrix}}_{U \quad \quad \quad\quad | \quad \textbf{c}}$$

Final Upper Triangular Form

We have reduced the system to an upper triangular form. The final system is:

$$U\textbf{x} = \textbf{c}$$

where:

$$U = \begin{bmatrix} 2 & -3 & 4 \\ 0 & 2 & 1 \\ 0 & 0 & -\frac{21}{4} \end{bmatrix} \quad \text{and} \quad \textbf{c} = \begin{bmatrix} -2 \\ 1 \\ -\frac{1}{4} \end{bmatrix}$$

We can now solve for $\textbf{x}$ by back-substitution.

From the third equation:

$$-\tfrac{21}{4}x_3 = -\tfrac{1}{4} \quad\Rightarrow\quad x_3 = \frac{-\tfrac{1}{4}}{-\tfrac{21}{4}} = \tfrac{1}{21}.$$

From the second equation using $x_3 = \frac1{21}$:

$$2x_2 + x_3 = 1 \quad\Rightarrow\quad 2x_2 = 1 - x_3 = 1 - \tfrac{1}{21} = \tfrac{20}{21} \quad\Rightarrow\quad x_2 = \tfrac{10}{21}.$$

From the first equation using $x_2 = \frac{10}{21}$ and $x_3 = \frac1{21}$:

$$2x_1 - 3x_2 + 4x_3 = -2 \quad\Rightarrow\quad 2x_1 - 3\left(\tfrac{10}{21}\right) + 4\left(\tfrac{1}{21}\right) = -2.$$ $$2x_1 - \tfrac{30}{21} + \tfrac{4}{21} = -2 \quad\Rightarrow\quad 2x_1 - \tfrac{26}{21} = -2 \quad\Rightarrow\quad 2x_1 = -2 + \tfrac{26}{21} \quad\Rightarrow\quad x_1 = -\tfrac{8}{21}.$$

And we are done! The vector $\textbf{x}\in \R^3$ that satisies $A\textbf{x} = \textbf{b}$ is:

$$\textbf{x} = \begin{pmatrix} -\tfrac{8}{21} \\ \tfrac{10}{21} \\ \tfrac{1}{21} \end{pmatrix}.$$

We found a solution! However this is not always the case. Depending on the particular LSE there might be infinitely many solutions or no solutions at all and our algorithm might fail in some cases. To be more precise, we get stuck if we can’t make the pivot nonzero using row exchanges. You can see the pdf on my website What does failure actually mean in Gaussian Elimination? to gain a little more insight on what might happen when the algorithm fails.

Row Operations: the general picture

There are important lemmas that formalize our intuition about the algorithm and actually prove that the solution to the transformed LSE $U\textbf{x} = \textbf{c}$ is the same as the solution for $A\textbf{x}=\textbf{b}$. The most important one is:

Lemma 3.2 justifies that Gauss elimination with back substitution finds the solution of the actual LSE and does not change it with row subtractions and permutations. Note that both row subtractions and permutations are repreented with invertible matrices. To invert them simply add the row back or swap the two rows back!

You can use Lemma 3.3 to prove Lemma 3.4

The row space is the column space of the transpose by its definition. Hence we need to show $\textbf{C}(A^\top) = \textbf{C}(A^\top M^\top)$. Note that $M^\top$ is also invertible (by Lemma 2.59 $\sim$ “transpose of invertible matrix is invertible”). Now we argue about the equality of two sets by showing if a vector is an element of one, then it is also an element of the other:

For $\textbf{C}(A^\top) \subseteq \textbf{C}(A^\top M^\top)$:

$$\begin{aligned} \textbf{v} \in \textbf{C}(A^\top) &\Rightarrow \textbf{v} = A^\top \textbf{x} \text{ for some $\textbf{x}\in \R^m$} \\ &\Rightarrow \textbf{v} = A^\top M^\top \textbf{x}' \text{ for some $\textbf{x}'\in \R^m$ } \quad (\textbf{x}' = (M^\top)^{-1} \textbf{x}) \\ & \Rightarrow \textbf{v} \in \textbf{C}(A^\top M^\top) \end{aligned}$$

For $\textbf{C}(A^\top M^\top) \subseteq \textbf{C}(A^\top)$ :

$$\begin{aligned} \textbf{v} \in \textbf{C}(A^\top M^\top) &\Rightarrow \textbf{v} = A^\top M^\top \textbf{x} \text{ for some $\textbf{x}\in \R^m$ } \\ &\Rightarrow \textbf{v} = A^\top \textbf{x}' \text{ for some $\textbf{x}'\in \R^m$} \quad (\textbf{x}' = M^\top \textbf{x}) \\ & \Rightarrow \textbf{v} \in \textbf{C}(A^\top) \end{aligned}$$

$\square$

The key element that enables this proof is that $M^\top$ is invertible. This is why we could find $\textbf{x}'$ for arbitrary $\textbf{x}$ in the first part of the proof such that $\textbf{x} = M^\top \textbf{x}'$.

This proof is written in a more compact way in the lecture notes and I recommend to read the proofs of all lemmas in this section from there.

One must be careful as the same property does NOT hold for the column space. In general $C(A)$ is not equal to $C(MA)$. Here is a counterexample:

$$A = \begin{bmatrix} 0 & 0 \\ 1 & 0 \end{bmatrix}, \quad M = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}$$

$M$ is invertible, in fact it is self-inverse ie. $M^2 = I$. We have

$$MA = \begin{bmatrix} 1 & 0 \\ 0 & 0 \end{bmatrix}$$

Now look at $\textbf{C}(A)$ and $\textbf{C}(MA)$:

$$\begin{aligned} \textbf{C}(A) &= \{\lambda \cdot \begin{pmatrix} 0 \\ 1 \end{pmatrix} \mid \lambda \in \R\} \\ \textbf{C}(MA) &= \{\lambda \cdot \begin{pmatrix} 1 \\ 0 \end{pmatrix} \mid \lambda \in \R\} \end{aligned}$$

which are different from each other. $\begin{pmatrix} 0 & 1 \end{pmatrix}^\top \not \in \textbf{C}(MA)$ and $\begin{pmatrix} 1 & 0 \end{pmatrix}^\top \not \in \textbf{C}(A)$.

The column spaces of $A$ and $MA$ are not the same, however their ranks are the same. In particular we can make the following statement:

Finally, we can say when the elimination fails and when it succeeds. This is the case when the columns of $A$ are linearly independent. This makes sense since $m$ linearly independent vectors can span the whole $\R^m$ so there exists a unique $\textbf{x} \in \R^m$ for every $\textbf{b} \in \R^m$ such that $A\textbf{x} = \textbf{b}$, given that the columns of $A$ are linearly independent. The following theorem formalizes this. You can find a beautiful proof in the lecture notes.

Runtime

It is always good to know the runtime of an algorithm you use. Here it is for Gauss Elimination. You can find the full runtime analysis in the lecture notes.

For $m$ equations and $m$ unknowns:

• Gauss Elimination: $O(m^3)$
• Back Substitution: $O(m^2)$
• Gauss Elimination with Back Substitution: $O(m^3)$
• Calculating the inverse of an invertible $m\times m$ matrix as in Algorithm 3 from the lecture notes: $O(m^3)$

Algorithm: Calculating the Inverse

We can use Gauss elimination with back substitution to calculate the inverse of a matrix. Let’s put the problem of calculating the inverse in terms of solving LSE’s. We want to find the inverse of a matrix $A \in \R^{m\times m}$. Call $A^{-1} = B$ for sake of simplicity in notation. We look for $B \in \R^{m\times m}$ with columns $\textbf{b}_1,\ldots,\textbf{b}_m \in \R^m$ such that:

$$AB = I$$

Taking a closer look into the matrix multiplication, we can formulate this expression as follows:

$$\begin{bmatrix} \vert & \vert & &\vert \\ A\textbf{b}_1 & A\textbf{b}_2 & \cdots& A\textbf{b}_m \\ \vert & \vert & &\vert \end{bmatrix} = \begin{bmatrix} \vert & \vert & &\vert \\ \textbf{e}_1 & \textbf{e}_2& \cdots &\textbf{e}_m \\ \vert & \vert & &\vert \\ \end{bmatrix}$$

The equality of two matrices means their columns must match. So we want the columns of $B$ to satisfy:

$$\begin{aligned} A\textbf{b}_1 =& \textbf{e}_1 \\ A\textbf{b}_2 =& \textbf{e}_2 \\ \vdots \\ A\textbf{b}_m =& \textbf{e}_m \end{aligned}$$

There we have $m$ linear systems of equations, all of which we can try to solve using Gauss elimination with back substitution. If Gauss elimination fails we know by Theorem 3.7 and Theorem 3.8 that the matrix $A$ has linearly dependent columns and hence is not invertible. Therefore we can assume that Gauss elimination succeeds for the following computations.

Applying Gauss elimination transforms A into an upper triangular form and all $\textbf{e}_i$‘s are also modified during the process. For upper triangular $U$ we have:

$$\begin{aligned} U\textbf{b}_1 =& \textbf{c}_1 \\ U\textbf{b}_2 =& \textbf{c}_2 \\ \vdots \\ U\textbf{b}_m =& \textbf{c}_m \end{aligned}$$

Now we can do back substitution to find the individual columns $\textbf{b}_i$. We will have calculated the inverse once we have all the columns.

Runtime Optimization: $m$ right hand sides

If you try to do the Gauss elimination (that has runtime $O(m^3)$) step from scratch for all $m$ equations we get a runtime of $O(m^4)$. We do not want this. You probably agree that calculating the $U$ from scratch is not necessary since $U$ is determined only by the matrix $A$. What we can do to improve this is to apply the algorithm for $m$ right hand side vectors.

We usually write the matrix in augmented form as $[A|\textbf{y}]$ when we try to solve one linear system of equation $A\textbf{x} = \textbf{y}$. We do the same here, but for all right hand sides. To continue with the example from above we have:

$$\begin{bmatrix} \,| & | & & | & \big| & | & | & & |\, \\ \mathbf{a}_1 & \mathbf{a}_2 & \cdots & \mathbf{a}_m & \bigg| & \mathbf{e}_1 & \mathbf{e}_2 & \cdots & \mathbf{e}_m \\ \,| & | & & | & \big| & | & | & & |\, \end{bmatrix}.$$

Now you can apply row subtractions and row swaps to $A$ and at the same time to all $m$ right hand side vectors. In the end you have:

$$\begin{bmatrix} \,| & | & & | & \big| & | & | & & |\, \\ \mathbf{u}_1 & \mathbf{u}_2 & \cdots & \mathbf{u}_m & \bigg| & \mathbf{c}_1 & \mathbf{c}_2 & \cdots & \mathbf{c}_m \\ \,| & | & & | & \big| & | & | & & |\, \end{bmatrix}.$$

where $u_i$‘s are the columns of the upper triangular matrix $U$ and $c_j$‘s are the right hand side vectors after elimination. This operation takes $O(m^3)$ time since we apply the elimination to all right hand sides at the same time. To convince yourself to this fact you can read section 3.2.6 from the lecture notes.

Example

Let $A = \begin{bmatrix} 1&2&3\\ 0&1&4\\ 5&6&0 \end{bmatrix}$ then we have:

$$\begin{aligned} \left[ \begin{array}{ccc|ccc} 1&2&3 & 1&0&0\\ 0&1&4 & 0&1&0\\ 5&6&0 & 0&0&1 \end{array} \right] &\overset{R_3\leftarrow R_3-5R_1}{\longrightarrow} \left[ \begin{array}{ccc|ccc} 1&2&3 & 1&0&0\\ 0&1&4 & 0&1&0\\ 0&-4&-15 & -5&0&1 \end{array} \right]\\ &\overset{R_3\leftarrow R_3+4R_2}{\longrightarrow} \left[ \begin{array}{ccc|ccc} 1&2&3 & 1&0&0\\ 0&1&4 & 0&1&0\\ 0&0&1 & -5&4&1 \end{array} \right] \end{aligned}$$

Now we have our upper triangular matrix $U$ on the left side and the vectors $\textbf{c}_j$ from before on the right side. We can do back substitution 3 times to get the inverse:

$$\begin{aligned} \begin{bmatrix} 1&2&3\\ 0&1&4\\ 0&0&1 \end{bmatrix} \textbf{b}_1 &= \begin{pmatrix} 1 \\ 0 \\ -5 \end{pmatrix} \\ \begin{bmatrix} 1&2&3\\ 0&1&4\\ 0&0&1 \end{bmatrix} \textbf{b}_2 &= \begin{pmatrix} 0 \\ 1 \\ 4 \end{pmatrix} \\ \begin{bmatrix} 1&2&3\\ 0&1&4\\ 0&0&1 \end{bmatrix} \textbf{b}_3 &= \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix} \end{aligned}$$

The back substitutions yield:

$$B = \begin{bmatrix} \vert & \vert & \vert \\ \textbf{b}_1 & \textbf{b}_2 &\textbf{b}_3 \\ \vert & \vert &\vert \\ \end{bmatrix} = A^{-1}= \begin{bmatrix} -24 & 18 & 5\\ 20 & -15 & -4\\ -5 & 4 & 1 \end{bmatrix}.$$

Hints

1. Solved In Class
2. A is invertible if and only if $A\textbf{x} \ne \textbf{0}$ for all $\textbf{x} \ne \textbf{0}$. Use it and show that one of the entries of the vector $A\textbf{x}$ can not be equal to 0.
3. This one is very similar to the in class exercise. Make a linear system of equations out of the given data points. Remember, the variables you want to find are $a,b,c$.
4. Try to represent $AA^{-1} = I$ as three separate systems of equations. Be careful with division. You can divide by a variable $\mu$ with the constraint that $\mu \neq 0$ and afterwards you must specify what happens when that variable is zero. In other words make a small case distinction.
5. You can first guess the inverse and then verify that $AA^{-1} = I$ by direct computation or calculate the inverse using Gauss elimination with back substitution. b) Assume that the pattern from $A$ holds in general. Try to generalize it by writing the columns of the inverse using unit vectors. Then show that the matrix you described is actually the inverse.
6. Try to write the equations using matrices as in $W = VM$.