$LU$ factorization

Suppose now that $Ax=b$ admits a unique solution. How to find it? We are especially interested in constructing a general procedure, that will work no matter what the size of $A$ 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

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$.

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$.

Starting from the last equation we can now find $x_3=1$, replace in the second to obtain $-5x_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"

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.

2.$LU$-factorization

• We have introduced Gaussian elimination as a procedure to solve the linear system $Ax=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 $Ax=b$

with

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

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

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

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

Verify by carrying out the matrix multiplication

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

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$

• Recall the basic operation in row echelon reduction: constructing a linear combination of rows to form zeros beneath the main diagonal, e.g.

  ${\displaystyle 𝑨=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1m}\\ a_{21}& a_{22}& \dots & a_{2m}\\ a_{31}& a_{32}& \dots & a_{3m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \dots & a_{mm}\end{array}\right)\sim \left(\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1m}\\ 0& a_{22}-\frac{a_{21}}{a_{11}}{a}_{12}& \dots & a_{2m}-\frac{a_{21}}{a_{11}}{a}_{1m}\\ 0& a_{32}-\frac{a_{31}}{a_{11}}{a}_{12}& \dots & a_{3m}-\frac{a_{31}}{a_{11}}{a}_{1m}\\ ⋮& ⋮& \ddots & ⋮\\ 0& a_{m2}-\frac{a_{m1}}{a_{11}}{a}_{12}& \dots & a_{mm}-\frac{a_{m1}}{a_{11}}{a}_{1m}\end{array}\right)}$

• This can be stated as a matrix multiplication operation, with $l_{i1}=a_{i1}/a_{11}$

  ${\displaystyle \left(\begin{array}{ccccc}1& 0& 0& \dots & 0\\ -l_{21}& 1& 0& \dots & 0\\ -l_{31}& 0& 1& \dots & 0\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ -l_{m1}& 0& 0& \dots & 1\end{array}\right)\left(\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1m}\\ a_{21}& a_{22}& \dots & a_{2m}\\ a_{31}& a_{32}& \dots & a_{3m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \dots & a_{mm}\end{array}\right)=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \dots & a_{1m}\\ 0& a_{22}-l_{21}{a}_{12}& \dots & a_{2m}-l_{21}{a}_{1m}\\ 0& a_{32}-l_{31}{a}_{12}& \dots & a_{3m}-l_{31}{a}_{1m}\\ ⋮& ⋮& \ddots & ⋮\\ 0& a_{m2}-l_{m1}{a}_{12}& \dots & a_{mm}-l_{m1}{a}_{1m}\end{array}\right)}$

Definition. The matrix

with $l_{i,k}=a_{i,k}^{\left(k\right)}/a_{k,k}^{\left(k\right)}$, and ${𝑨}^{\left(k\right)}=\left(a_{i,j}^{\left(k\right)}\right)$ 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

  ${\displaystyle {𝑳}_{m-1}{𝑳}_{m-2}\dots {𝑳}_2{𝑳}_1𝑨=𝑼}$

• The inverse of a Gaussian multiplier is

  ${\displaystyle {𝑳}_k^{-1}=\left(\begin{array}{ccccc}1& \dots & 0& \dots & 1\\ 0& \ddots & 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 & ⋮& \ddots & ⋮\\ 0& \dots & l_{m,k}& \dots & 1\end{array}\right)=𝑰-({𝑳}_k-𝑰)}$

• From $\left({𝑳}_{m-1}{𝑳}_{m-2}\dots {𝑳}_2{𝑳}_1\right)𝑨=𝑼$ obtain

  ${\displaystyle 𝑨={\left({𝑳}_{m-1}{𝑳}_{m-2}\dots {𝑳}_2{𝑳}_1\right)}^{-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

  ${\displaystyle 𝑳=\left(\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\\ ⋮& ⋮& ⋮& \ddots & ⋮& ⋮\\ 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}\right)}$

We will show that this indeed possible if $Ax=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

and obtain

or

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 $Ax=b$

1. Find the factorization $LU=A$

2. Insert the factorization into $Ax=b$ to obtain $(LU)x=L(Ux)=Ly=b$, where the notation $y=Ux$ has been introduced. The system

   ${\displaystyle Ly=b}$

   is easy to solve by forward substitution to find $y$ for given $b$

3. Finally find $x$ by backward substitution solution of

   ${\displaystyle Ux=y}$

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

3.Inverse matrix

By analogy to the simple scalar equation $ax=b$ with solution $x=a^{-1}b$ when $a\ne 0$, we are interested in writing the solution to a linear system $Ax=b$ as $x=A^{-1}b$ for $A\in {ℝ}^{m\times m}$, $x\in {ℝ}^m$. Recall that solving $Ax=b=Ib$ 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

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 $AB=I$. Introducing the column notation for $B,I$ leads to

and identification of each column in the equality states

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

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$

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

to obtain

4.Determinants

• $𝑨\in {ℝ}^{m\times m}$ a square matrix, $\det \left(𝑨\right)\in ℝ$ is the oriented volume enclosed by the column vectors of $𝑨$ (a parallelipiped)

• Geometric interpretation of determinants

• Determinant calculation rules

• Algebraic definition of a determinant

Definition. The determinant of a square matrix $𝑨=\left(\begin{array}{lll}{𝒂}_1& \dots & {𝒂}_m\end{array}\right)\in {ℝ}^{m\times m}$ is a real number

giving the (oriented) volume of the parallelepiped spanned by matrix column vectors.

• $m=2$

  ${\displaystyle 𝑨=\left(\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right),\det \left(𝑨\right)=\left|\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right|}$

• $m=3$

  ${\displaystyle 𝑨=\left(\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right),\det \left(𝑨\right)=\left|\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right|}$

• Computation of a determinant with $m=2$

  ${\displaystyle \left|\begin{array}{cc}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right|=a_{11}{a}_{22}-a_{12}{a}_{21}}$

• Computation of a determinant with $m=3$

  ${\displaystyle \begin{array}{rcl}\left|\begin{array}{ccc}{a}_{11}& a_{12}& a_{13}\\ a_{21}& a_{22}& a_{23}\\ a_{31}& a_{32}& a_{33}\end{array}\right|& =& a_{11}{a}_{22}{a}_{33}+a_{21}{a}_{32}{a}_{13}+a_{31}{a}_{12}{a}_{23}\\ & & -a_{13}{a}_{22}{a}_{31}-a_{23}{a}_{32}{a}_{11}-a_{33}{a}_{12}{a}_{21}\end{array}}$

• Where do these determinant computation rules come from? Two viewpoints

  • Geometric viewpoint: determinants express parallelepiped volumes

  • Algebraic viewpoint: determinants are computed from all possible products that can be formed from choosing a factor from each row and each column

• $m=2$

  Figure 1.

• In two dimensions a “parallelepiped” becomes a parallelogram with area given as

  ${\displaystyle \left(\mathrm{Area}\right)=(\mathrm{Length}\mathrm{of}\mathrm{Base})\times (\mathrm{Length}\mathrm{of}\mathrm{Height})}$

• Take ${𝒂}_1$ as the base, with length $b=||{𝒂}_1||$. Vector ${𝒂}_1$ is at angle ${\phi }_1$ to $x_1$-axis, ${𝒂}_2$ is at angle ${\phi }_2$ to $x_2$-axis, and the angle between ${𝒂}_1$, ${𝒂}_2$ is $\theta ={\phi }_2-{\phi }_1$. The height has length

  ${\displaystyle h=||{𝒂}_2||\sin \theta =||{𝒂}_2||\sin .(.{\phi }_2-{\phi }_1)=||{𝒂}_2||(\sin {\phi }_2\cos {\phi }_1-\sin {\phi }_1\cos {\phi }_2)}$

• Use $\cos {\phi }_1=a_{11}/||{𝒂}_1||$, $\sin {\phi }_1=a_{12}/||{𝒂}_1||$, $\cos {\phi }_2=a_{21}/||{𝒂}_2||$, $\sin {\phi }_2=a_{22}/||{𝒂}_2||$

  ${\displaystyle \left(\mathrm{Area}\right)=||{𝒂}_1||||{𝒂}_2||(\sin {\phi }_2\cos {\phi }_1-\sin {\phi }_1\cos {\phi }_2)=a_{11}{a}_{22}-a_{12}{a}_{21}}$

• The geometric interpretation of a determinant as an oriented volume is useful in establishing rules for calculation with determinants:

  • Determinant of matrix with repeated columns is zero (since two edges of the parallelepiped are identical). Example for $m=3$

    ${\displaystyle \Delta =\left|\begin{array}{ccc}a& a& u\\ b& b& v\\ c& c& w\end{array}\right|=abw+bcu+cav-ubc-vca-wab=0}$

    This is more easily seen using the column notation

    ${\displaystyle \Delta =\det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_1& {𝒂}_3& \dots \end{array}\right)=0}$

  • Determinant of matrix with linearly dependent columns is zero (since one edge lies in the 'hyperplane' formed by all the others)

• Separating sums in a column (similar for rows)

  ${\displaystyle \det \left(\begin{array}{cccc}{𝒂}_1+{𝒃}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)=\det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)+\det \left(\begin{array}{cccc}{𝒃}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)}$

  with ${𝒂}_i,{𝒃}_1\in {ℝ}^m$

• Scalar product in a column (similar for rows)

  ${\displaystyle \det \left(\begin{array}{cccc}\alpha {𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)=\alpha \det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)}$

  with $\alpha \in ℝ$

• Linear combinations of columns (similar for rows)

  ${\displaystyle \det \left(\begin{array}{cccc}{𝒂}_1& {𝒂}_2& \dots & {𝒂}_m\end{array}\right)=\det \left(\begin{array}{cccc}{𝒂}_1& \alpha {𝒂}_1+{𝒂}_2& \dots & {𝒂}_m\end{array}\right)}$

  with $\alpha \in ℝ$.

• A determinant of size $m$ can be expressed as a sum of determinants of size $m-1$ by expansion along a row or column

  ${\displaystyle \begin{array}{rcl}\left|\begin{array}{ccccc}{a}_{11}& a_{12}& a_{13}& \dots & a_{1m}\\ a_{21}& a_{22}& a_{23}& \dots & a_{2m}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& a_{m3}& \dots & a_{mm}\end{array}\right|& =& a_{11}\left|\begin{array}{cccc}{a}_{22}& a_{23}& \dots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m2}& a_{m3}& \dots & a_{mm}\end{array}\right|-\\ & & a_{12}\left|\begin{array}{cccc}{a}_{21}& a_{23}& \dots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m3}& \dots & a_{mm}\end{array}\right|+\\ & & a_{13}\left|\begin{array}{ccccc}{a}_{21}& a_{22}& a_{24}& \dots & a_{2m}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& a_{m4}& \dots & a_{mm}\end{array}\right|-\\ & & \dots \\ & & +{(-1)}^{m+1}{a}_{1m}\left|\begin{array}{cccc}{a}_{21}& a_{23}& \dots & a_{2,m-1}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m3}& \dots & a_{m,m-1}\end{array}\right|\end{array}}$

• The formal definition of a determinant

  ${\displaystyle \det A=\sum _{\sigma \in \Sigma }\nu \left(\sigma \right)a_{1i_1}{a}_{2i_2}\dots a_{m{i}_m}}$

  requires $mm!$ operations, a number that rapidly increases with $m$

• A more economical determinant is to use row and column combinations to create zeros and then reduce the size of the determinant, an algorithm reminiscent of Gauss elimination for systems

  Example:

  ${\displaystyle \left|\begin{array}{ccc}1& 2& 3\\ -1& 0& 1\\ -2& -1& 4\end{array}\right|=\left|\begin{array}{ccc}1& 2& 3\\ 0& 2& 4\\ 0& 3& 10\end{array}\right|=\left|\begin{array}{cc}2& 4\\ 3& 10\end{array}\right|=20-12=8}$

  The first equality comes from linear combinations of rows, i.e. row 1 is added to row 2, and row 1 multiplied by 2 is added to row 3. These linear combinations maintain the value of the determinant. The second equality comes from expansion along the first column

4.1.Cross product

• Consider $u,v\in {ℝ}^3$. We've introduced the idea of a scalar product

  ${\displaystyle u\cdot v=u^{T}v=u_1{v}_1+u_2{v}_2+u_3{v}_3}$

  in which from two vectors one obtains a scalar

• We've also introduced the idea of an exterior product

  ${\displaystyle u{v}^T=\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right)\left(\begin{array}{ccc}{v}_1& v_2& v_3\end{array}\right)=\left(\begin{array}{ccc}{u}_1{v}_1& u_1{v}_2& u_1{v}_3\\ u_2{v}_1& u_2{v}_2& u_2{v}_3\\ u_3{v}_1& u_3{v}_2& u_3{v}_3\end{array}\right)}$

  in which a matrix is obtained from two vectors

• Another product of two vectors is also useful, the cross product, most conveniently expressed in determinant-like form

  ${\displaystyle u\times v=\left|\begin{array}{ccc}{𝒆}_1& {𝒆}_2& {𝒆}_3\\ u_1& u_2& u_3\\ v_1& v_2& v_3\end{array}\right|=(u_2{v}_3-v_2{u}_3){𝒆}_1+(u_3{v}_1-v_3{u}_1){𝒆}_2+(u_1{v}_2-v_1{u}_2){𝒆}_3}$