MATH661

Lecture 10: Square Linear Systems

1.Gaussian elimination and row echelon reduction

Suppose now that $𝑨 𝒙 = 𝒃$ admits a unique solution with $𝑨 \in ℝ^{m \times m}$ , $𝒙, 𝒃 \in ℝ^{m}$ , known as a square system of linear equations where the number of unknowns equals that of the equations. How to find it? We are especially interested in constructing a general procedure, that will work no matter what the size of $𝑨$ might be. This means we seek an algorithm that precisely specifies the steps that lead to the solution, and that we can program a computing device to carry out automatically. One such algorithm is Gaussian elimination.

Consider the system

{\begin{array}{ccc} x_{1} + 2 x_{2} - x_{3} & = & 2 \\ 2 x_{1} - x_{2} + x_{3} & = & 2 \\ 3 x_{1} - x_{2} - x_{3} & = & 1 \end{array} .

The idea is to combine equations such that we have one fewer unknown in each equation. Ask: with what number should the first equation be multiplied in order to eliminate $x_{1}$ from sum of equation 1 and equation 2? This number is called a Gaussian multiplier, and is in this case $- 2$ . Repeat the question for eliminating $x_{1}$ from third equation, with multiplier $- 3$ .

{\begin{array}{ccc} x_{1} + 2 x_{2} - x_{3} & = & 2 \\ 2 x_{1} - x_{2} + x_{3} & = & 2 \\ 3 x_{1} - x_{2} - x_{3} & = & 1 \end{array} . \Rightarrow {\begin{array}{ccc} x_{1} + 2 x_{2} - x_{3} & = & 2 \\ - 5 x_{2} + 3 x_{3} & = & - 2 \\ - 7 x_{2} + 2 x_{3} & = & - 5 \end{array} .

Now, ask: with what number should the second equation be multiplied to eliminate $x_{2}$ from sum of second and third equations. The multiplier is in this case $- 7 / 5$ .

{\begin{array}{ccc} x_{1} + 2 x_{2} - x_{3} & = & 2 \\ - 5 x_{2} + 3 x_{3} & = & - 2 \\ - 7 x_{2} + 2 x_{3} & = & - 5 \end{array} . \Rightarrow {\begin{array}{ccc} x_{1} + 2 x_{2} - x_{3} & = & 2 \\ - 5 x_{2} + 3 x_{3} & = & - 2 \\ - \frac{11}{5} x_{3} & = & - \frac{11}{5} \end{array} .

Starting from the last equation we can now find $x_{3} = 1$ , replace in the second to obtain $- 5 x_{2} = - 5$ , hence $x_{2} = 1$ , and finally replace in the first equation to obtain $x_{1} = 1$ .

The above operations only involve coefficients. A more compact notation is therefore to work with what is known as the "bordered matrix"

(\begin{array}{cccc} 1 & 2 & - 1 & 2 \\ 2 & - 1 & 1 & 2 \\ 3 & - 1 & - 1 & 1 \end{array}) \sim (\begin{array}{cccc} 1 & 2 & - 1 & 2 \\ 0 & - 5 & 3 & - 2 \\ 0 & - 7 & 2 & - 5 \end{array}) \sim (\begin{array}{cccc} 1 & 2 & - 1 & 2 \\ 0 & - 5 & 3 & - 2 \\ 0 & 0 & - \frac{11}{5} & - \frac{11}{5} \end{array})

Once the above triangular form has been obtain, the solution is found by back substitution, in which we seek to form the identity matrix in the first 3 columns, and the solution is obtained in the last column.

(\begin{array}{cccc} 1 & 2 & - 1 & 2 \\ 0 & - 5 & 3 & - 2 \\ 0 & 0 & - \frac{11}{5} & - \frac{11}{5} \end{array}) \sim (\begin{array}{cccc} 1 & 2 & - 1 & 2 \\ 0 & - 5 & 3 & - 2 \\ 0 & 0 & 1 & 1 \end{array}) \sim (\begin{array}{cccc} 1 & 0 & 0 & 1 \\ 0 & 1 & 0 & 1 \\ 0 & 0 & 1 & 1 \end{array})

2. $L U$ -factorization

We have introduced Gaussian elimination as a procedure to solve the linear system $A x = b$ ("find coordinates of vector $b$ in terms of column vectors of matrix $A$ "), $x, b \in ℝ^{m}, A \in ℝ^{m \times m}$
We now reinterpret Gaussian elimination as a sequence of matrix multiplications applied to $A$ to obtain a simpler, upper triangular form.

2.1.Example for $m = 3$

Consider the system $A x = b$

{\begin{array}{ccc} x_{1} + 2 x_{2} - x_{3} & = & 2 \\ 2 x_{1} - x_{2} + x_{3} & = & 2 \\ 3 x_{1} - x_{2} - x_{3} & = & 1 \end{array} .

with

A = (\begin{array}{ccc} 1 & 2 & - 1 \\ 2 & - 1 & 1 \\ 3 & - 1 & - 1 \end{array}), b = (\begin{array}{c} 2 \\ 2 \\ 1 \end{array})

We ask if there is a matrix $L_{1}$ that could be multiplied with $A$ to produce a result $L_{1} A$ with zeros under the main diagonal in the first column. First, gain insight by considering multiplication by the identity matrix, which leaves $A$ unchanged

(\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}) (\begin{array}{ccc} 1 & 2 & - 1 \\ 2 & - 1 & 1 \\ 3 & - 1 & - 1 \end{array}) = (\begin{array}{ccc} 1 & 2 & - 1 \\ 2 & - 1 & 1 \\ 3 & - 1 & - 1 \end{array})

In the first stage of Gaussian multiplication, the first line remains unchanged, so we deduce that $L_{1}$ should have the same first line as the identity matrix

L_{1} = (\begin{array}{ccc} 1 & 0 & 0 \\ ? & ? & ? \\ ? & ? & ? \end{array})

(\begin{array}{ccc} 1 & 0 & 0 \\ ? & ? & ? \\ ? & ? & ? \end{array}) (\begin{array}{ccc} 1 & 2 & - 1 \\ 2 & - 1 & 1 \\ 3 & - 1 & - 1 \end{array}) = (\begin{array}{ccc} 1 & 2 & - 1 \\ 0 & - 5 & 3 \\ 0 & - 7 & 2 \end{array})

Next, recall the way Gaussian multipliers were determined: find number to multiply first line so that added to second, third lines a zero is obtained. This leads to the form

L_{1} = (\begin{array}{ccc} 1 & 0 & 0 \\ ? & 1 & 0 \\ ? & 0 & 1 \end{array})

Finally, identify the missing entries with the Gaussian multipliers to determine

L_{1} = (\begin{array}{ccc} 1 & 0 & 0 \\ - 2 & 1 & 0 \\ - 3 & 0 & 1 \end{array})

Verify by carrying out the matrix multiplication

L_{1} A = (\begin{array}{ccc} 1 & 0 & 0 \\ - 2 & 1 & 0 \\ - 3 & 0 & 1 \end{array}) (\begin{array}{ccc} 1 & 2 & - 1 \\ 2 & - 1 & 1 \\ 3 & - 1 & - 1 \end{array}) = (\begin{array}{ccc} 1 & 2 & - 1 \\ 0 & - 5 & 3 \\ 0 & - 7 & 2 \end{array})

Repeat the above reasoning to come up with a second matrix $L_{2}$ that forms a zero under the main diagonal in the second column

L_{2} = (\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & - 7 / 5 & 1 \end{array})

L_{2} L_{1} A = (\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & - 7 / 5 & 1 \end{array}) (\begin{array}{ccc} 1 & 2 & - 1 \\ 0 & - 5 & 3 \\ 0 & - 7 & 2 \end{array}) = (\begin{array}{ccc} 1 & 2 & - 1 \\ 0 & - 5 & 3 \\ 0 & 0 & - 11 / 5 \end{array}) = U

We have obtained a matrix with zero entries under the main diagonal (an upper triangular matrix) by a sequence of matrix multiplications.

2.2.General $m$ case

From the above, we assume that we can form a sequence of multiplier matrices such that the result is an upper triangular matrix $U$

L_{m - 1} . . . L_{2} L_{1} A = U

Recall the basic operation in row echelon reduction: constructing a linear combination of rows to form zeros beneath the main diagonal, e.g.
$𝑨 = (\begin{array}{cccc} a_{11} & a_{12} & \dots & a_{1 m} \\ a_{21} & a_{22} & \dots & a_{2 m} \\ a_{31} & a_{32} & \dots & a_{3 m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{m 1} & a_{m 2} & \dots & a_{m m} \end{array}) \sim (\begin{array}{cccc} a_{11} & a_{12} & \dots & a_{1 m} \\ 0 & a_{22} - \frac{a_{21}}{a_{11}} a_{12} & \dots & a_{2 m} - \frac{a_{21}}{a_{11}} a_{1 m} \\ 0 & a_{32} - \frac{a_{31}}{a_{11}} a_{12} & \dots & a_{3 m} - \frac{a_{31}}{a_{11}} a_{1 m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & a_{m 2} - \frac{a_{m 1}}{a_{11}} a_{12} & \dots & a_{m m} - \frac{a_{m 1}}{a_{11}} a_{1 m} \end{array})$
This can be stated as a matrix multiplication operation, with $l_{i 1} = a_{i 1} / a_{11}$
$(\begin{array}{ccccc} 1 & 0 & 0 & \dots & 0 \\ - l_{21} & 1 & 0 & \dots & 0 \\ - l_{31} & 0 & 1 & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ - l_{m 1} & 0 & 0 & \dots & 1 \end{array}) (\begin{array}{cccc} a_{11} & a_{12} & \dots & a_{1 m} \\ a_{21} & a_{22} & \dots & a_{2 m} \\ a_{31} & a_{32} & \dots & a_{3 m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ a_{m 1} & a_{m 2} & \dots & a_{m m} \end{array}) = (\begin{array}{cccc} a_{11} & a_{12} & \dots & a_{1 m} \\ 0 & a_{22} - l_{21} a_{12} & \dots & a_{2 m} - l_{21} a_{1 m} \\ 0 & a_{32} - l_{31} a_{12} & \dots & a_{3 m} - l_{31} a_{1 m} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & a_{m 2} - l_{m 1} a_{12} & \dots & a_{m m} - l_{m 1} a_{1 m} \end{array})$

Definition. The matrix

$𝑳_{k} = (\begin{array}{ccccc} 1 & \dots & 0 & \dots & 1 \\ 0 & ⋱ & 0 & \dots & 0 \\ 0 & \dots & 1 & \dots & 0 \\ 0 & \dots & - l_{k + 1, k} & \dots & 0 \\ 0 & \dots & - l_{k + 2, k} & \dots & 0 \\ ⋮ & \dots & ⋮ & ⋱ & ⋮ \\ 0 & \dots & - l_{m, k} & \dots & 1 \end{array})$

with $l_{i, k} = a_{i, k}^{(k)} / a_{k, k}^{(k)}$ , and $𝑨^{(k)} = (a_{i, j}^{(k)})$ the matrix obtained after step $k$ of row echelon reduction (or, equivalently, Gaussian elimination) is called a Gaussian multiplier matrix.

For $𝑨 \in ℝ^{m \times m}$ nonsingular, the successive steps in row echelon reduction (or Gaussian elimination) correspond to successive multiplications on the left by Gaussian multiplier matrices
$𝑳_{m - 1} 𝑳_{m - 2} \dots 𝑳_{2} 𝑳_{1} 𝑨 = 𝑼$
The inverse of a Gaussian multiplier is
$𝑳_{k}^{- 1} = (\begin{array}{ccccc} 1 & \dots & 0 & \dots & 1 \\ 0 & ⋱ & 0 & \dots & 0 \\ 0 & \dots & 1 & \dots & 0 \\ 0 & \dots & l_{k + 1, k} & \dots & 0 \\ 0 & \dots & l_{k + 2, k} & \dots & 0 \\ ⋮ & \dots & ⋮ & ⋱ & ⋮ \\ 0 & \dots & l_{m, k} & \dots & 1 \end{array}) = 𝑰 - (𝑳_{k} - 𝑰)$

From $(𝑳_{m - 1} 𝑳_{m - 2} \dots 𝑳_{2} 𝑳_{1}) 𝑨 = 𝑼$ obtain
$𝑨 = {(𝑳_{m - 1} 𝑳_{m - 2} \dots 𝑳_{2} 𝑳_{1})}^{- 1} 𝑼 = 𝑳_{1}^{- 1} 𝑳_{2}^{- 1} \cdot \dots \cdot 𝑳_{m - 1}^{- 1} 𝑼 = 𝑳 𝑼$
Due to the simple form of $𝑳_{k}^{- 1}$ the matrix $𝑳$ is easily obtained as
$𝑳 = (\begin{array}{cccccc} 1 & 0 & 0 & \dots & 0 & 0 \\ l_{2, 1} & 1 & 0 & \dots & 0 & 0 \\ l_{3, 1} & l_{3, 2} & 1 & \dots & 0 & 0 \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ l_{m - 1, 1} & l_{m - 1, 2} & l_{m - 1, 3} & \dots & 1 & 0 \\ l_{m, 1} & l_{m, 2} & l_{m, 3} & \dots & l_{m, m - 1} & 1 \end{array})$

We will show that this indeed possible if $A x = b$ admits a unique solution. Furthermore, the product of lower triangular matrices is lower triangular, and the inverse of a lower triangular matrix is lower triangular (same applies for upper triangular matrices). Introduce the notation

L^{- 1} = L_{m - 1} . . . L_{2} L_{1}

and obtain

L^{- 1} A = U

A = L U

The above result permits a basic insight into Gaussian elimination: the procedure depends on "factoring" the matrix $A$ into two "simpler" matrices $L, U$ . The idea of representing a matrix as a product of simpler matrices is fundamental to linear algebra, and we will come across it repeatedly.

For now, the factorization allows us to devise the following general approach to solving $A x = b$

Find the factorization $L U = A$
Insert the factorization into $A x = b$ to obtain $(L U) x = L (U x) = L y = b$ , where the notation $y = U x$ has been introduced. The system
$L y = b$
is easy to solve by forward substitution to find $y$ for given $b$
Finally find $x$ by backward substitution solution of
$U x = y$

Algorithm Gauss elimination without pivoting

for $s = 1$ to $m - 1$

for $i = s + 1$ to $m$

$t = - a_{i s} / a_{s s}$

for $j = s + 1$ to $m$

$a_{i j} = a_{i j} + t \cdot a_{s j}$

$b_{i} = b_{i} + t \cdot b_{s}$

for $s = m$ downto 1

$x_{s} = b_{s} / a_{s s}$

for $i = 1$ to $s - 1$

$b_{i} = b_{i} - a_{i s} \cdot x_{s}$

return $x$

Algorithm Gauss elimination with partial pivoting

$p = 1 : m$ (initialize row permutation vector)

for $s = 1$ to $m - 1$

piv = $abs (a_{p (s), s})$

for $i = s + 1$ to $m$

mag = $abs (a_{p (i), s})$

if $mag > piv$ then

$piv = mag; k = p (s); p (s) = p (i); p (i) = k$

if $piv < ϵ$ then break(“Singular matrix”)

$t = - a_{p (i) s} / a_{p (s) s}$

for $j = s + 1$ to $m$

$a_{p (i) j} = a_{p (i) j} + t \cdot a_{p (s) j}$

$b_{p (i)} = b_{p (i)} + t \cdot b_{p (s)}$

for $s = m$ downto 1

$x_{s} = b_{p (s)} / a_{p (s) s}$

for $i = 1$ to $s - 1$

$b_{p (i)} = b_{p (i)} - a_{p (i) s} \cdot x_{s}$

return $x$

Given $𝑨 \in ℝ^{m \times n}$

Singular value decomposition	Gram-Schmidt	Lower-upper
Transformation of coordinates	$𝑨 𝒙 = 𝒃$
$𝑼 𝚺 𝑽^{T} = 𝑨$	$𝑸 𝑹 = 𝑨$	$𝑳 𝑼 = 𝑨$
$(𝑼 𝚺 𝑽^{T}) 𝒙 = 𝒃 \Rightarrow 𝑼 𝒚 = 𝒃 \Rightarrow 𝒚 = 𝑼^{T} 𝒃$	$(𝑸 𝑹) 𝒙 = 𝒃 \Rightarrow 𝑸 𝒚 = 𝒃$ , $𝒚 = 𝑸^{T} 𝒃$	$(𝑳 𝑼) 𝒙 = 𝒃 \Rightarrow 𝑳 𝒚 = 𝒃 (forward sub to find)$ y
$𝚺 𝒛 = 𝒚 \Rightarrow 𝒛 = 𝚺^{+} 𝒚$	$𝑹 𝒙 = 𝒚 (back sub to find . 𝒙$ )	$𝑼 𝒙 = 𝒚$ (back sub to find $𝒙$ )
$𝑽^{T} 𝒙 = 𝒛 \Rightarrow 𝒙 = 𝑽 𝒛$

3.Matrix inverse

By analogy to the simple scalar equation $a x = b$ with solution $x = a^{- 1} b$ when $a \neq 0$ , we are interested in writing the solution to a linear system $A x = b$ as $x = A^{- 1} b$ for $A \in ℝ^{m \times m}$ , $x \in ℝ^{m}$ . Recall that solving $A x = b = I b$ corresponds to expressing the vector $b$ as a linear combination of the columns of $A$ . This can only be done if the columns of $A$ form a basis for $ℝ^{m}$ , in which case we say that $A$ is non-singular.

Definition 1. For matrix $A \in ℝ^{m \times m}$ non-singular the inverse matrix is denoted by $A^{- 1}$ and satisfies the properties

$A A^{- 1} = A^{- 1} A = I$

3.1.Gauss-Jordan algorithm

Computation of the inverse $A^{- 1}$ can be carried out by repeated use of Gauss elimination. Denote the inverse by $B = A^{- 1}$ for a moment and consider the inverse matrix property $A B = I$ . Introducing the column notation for $B, I$ leads to

A (\begin{array}{ccc} B_{1} & . . . & B_{m} \end{array}) = (\begin{array}{ccc} e_{1} & . . . & e_{m} \end{array})

and identification of each column in the equality states

A B_{k} = e_{k}, k = 1, 2, . ., m

with $e_{k}$ the column unit vector with zero components everywhere except for a $1$ in row $k$ . To find the inverse we need to simultaneously solve the $m$ linear systems given above.

Gauss-Jordan algorithm example. Consider

A = (\begin{array}{ccc} 1 & 2 & 3 \\ - 1 & 3 & 1 \\ 2 & - 1 & - 2 \end{array})

To find the inverse we solve the systems $A B_{1} = e_{1}, A B_{2} = e_{2}, A B_{3} = e_{3}$ . This can be done simultaneously by working with the matrix $A$ bordered by $I$

(A | I) = (\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ - 1 & 1 & 1 & 0 & 1 & 0 \\ 2 & 4 & - 2 & 0 & 0 & 1 \end{array})

Carry out now operations involving linear row combinations and permutations to bring the left side to $I$

(\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ - 1 & 1 & 1 & 0 & 1 & 0 \\ 2 & 4 & - 2 & 0 & 0 & 1 \end{array}) \sim (\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 2 & 1 & 1 & 1 & 0 \\ 0 & 2 & - 2 & - 2 & 0 & 1 \end{array}) \sim (\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 2 & 1 & 1 & 1 & 0 \\ 0 & 0 & - 3 & - 3 & - 1 & 1 \end{array}) \sim

\sim (\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 2 & 1 & 1 & 1 & 0 \\ 0 & 0 & 1 & 1 & \frac{1}{3} & - \frac{1}{3} \end{array}) \sim (\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 2 & 0 & 0 & \frac{2}{3} & \frac{1}{3} \\ 0 & 0 & 1 & 1 & \frac{1}{3} & - \frac{1}{3} \end{array}) \sim (\begin{array}{cccccc} 1 & 1 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & \frac{1}{3} & \frac{1}{6} \\ 0 & 0 & 1 & 1 & \frac{1}{3} & - \frac{1}{3} \end{array}) \sim

(\begin{array}{cccccc} 1 & 0 & 0 & 1 & - \frac{1}{3} & - \frac{1}{6} \\ 0 & 1 & 0 & 0 & \frac{1}{3} & \frac{1}{6} \\ 0 & 0 & 1 & 1 & \frac{1}{3} & - \frac{1}{3} \end{array})

to obtain

A^{- 1} = (\begin{array}{ccc} 1 & - \frac{1}{3} & - \frac{1}{6} \\ 0 & \frac{1}{3} & \frac{1}{6} \\ 1 & \frac{1}{3} & - \frac{1}{3} \end{array})