LU decomposition, without pivoting

The LU decomposition is a factorization of a matrix $A$ into

\begin{equation*} A = LU \end{equation*}

where $L$ is lower triangular and $U$ is upper-triangular. Usually we also require the diagonal elements of $L$ to be 1. For simplicity of exposition we will assume $A$ is square, $m \times m$.

The main use of the LU decomposition is in the solution of $Ax=b$ problems, where $A$ is a known square matrix and $b$ is a known right-hand-side vector and $x$ is an unknown vector. We'll get into the details of that in the next lecture. For now let's concentrate on computing the LU decomposition $A=LU$ given a matrix $A$.

LU algorithm: the big picture

The LU decomposition is computed as as a series of left multiplications of $A$ by lower-triangular matrices $L_1, L_2, \ldots, L_{m-1}$ to form an upper-triangular matrix $U$.

\begin{equation*} L_{m-1} L_{m-2} \ldots L_3 L_2 L_1 A = U \end{equation*}

Each $L_j$ multiplication introduces zeros on the subdiagonal elements of the $j$th column of the product. Let $L = (L_{m-1} L_{m-2} \ldots L_3 L_2 L_1)^{-1}$. It turns out that $L$ is lower-triangular as well. Now we multiply from left by $L$ to get the LU factorization.

\begin{equation*} A = LU \end{equation*}

LU algorithm, a 4 x 4 example

The easiest way to see how this works is by example. Let's start with the following matrix $A$.

\begin{equation*} A = \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ -6 & -8 & 1 & 0 \\ 8 & 9 & 6 & -4 \\ -2 & 2 & -17 & 7 \end{array}\right] \end{equation*}

Step 1. We now multiply $A$ from left by a lower-triangular matrix $L_1$ that is designed to put zeros in the subdiagonal elements of the 1st column.

\begin{align*} \left[\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ 3 & 1 & 0 & 0 \\ -4 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 \end{array}\right] \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ -6 & -8 & 1 & 0 \\ 8 & 9 & 6 & -4 \\ -2 & 2 & -17 & 7 \end{array}\right] &= \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ 0 & 1 & -2 & 3 \\ 0 & -3 & 10 & -8 \\ 0 & 5 & -18 & 8 \end{array}\right]\\ L_1~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~A~~~~~~~~~~~~~~~~~~~~~&=~~~~~~~~~~~~~L_1 A \end{align*}

How did we know to put -2, 4, 3 on the subdiagonal of the 1st column of $L_1$? Well, [-2 1 0 0] in the 2nd row of $L_1$

Some text

\begin{align*} \left[\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 3 & 1 & 0 \\ 0 & -5 & 0 & 1 \end{array}\right] \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ 0 & 1 & -2 & 3 \\ 0 & -3 & 10 & -8 \\ 0 & 5 & -18 & 8 \end{array}\right] &= \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ 0 & 1 & -2 & 3 \\ 0 & 0 & 4 & 1 \\ 0 & 0 & -8 & -7 \end{array}\right] \\ L_2~~~~~~~~~~~~~~~~~~~~~~~~~L_1 A~~~~~~~~~~~~~~~~~~~~~&=~~~~~~~~~~~L_2 L_1 A \end{align*}

Some text

\begin{align*} \left[\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 2 & 1 \end{array}\right] \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ 0 & 1 & -2 & 3 \\ 0 & 0 & 4 & 1 \\ 0 & 0 & -8 & -7 \end{array}\right] &= \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ 0 & 1 & -2 & 3 \\ 0 & 0 & 4 & 1 \\ 0 & 0 & 0 & -5 \end{array}\right] \\ L_3~~~~~~~~~~~~~~~~~~~~~~~L_2\, L_1 \,A~~~~~~~~~~~~~~~~~~~~~&=~~~~~~~~~L_3\, L_2\, L_1\,A = U \end{align*}

Thus we have $L_3 L_2 L_1 A = U$. Let $L = (L_3 L_2 L_1)^{-1}$. Then multiply both sides of $L_3 L_2 L_1 A = U$ by $L = (L_3 L_2 L_1)^{-1}$ to get $A = LU$. It is not hard to show that $L$ is lower-triangular and has subdiagonal elements given by negating the corresponding elements of $L_1, L_2,$ and $L_3$. I.e.

\begin{align*} L = (L_3 L_2 L_1)^{-1} = \left[\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 4 & -5 & 1 & 0 \\ -3 & 3 & 2 & 1 \end{array}\right] \end{align*}

Putting all the pieces together, we arrive at an LU decomposition

\begin{align*} \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ -6 & -8 & 1 & 0 \\ 8 & 9 & 6 & -4 \\ -2 & 2 & -17 & -7 \end{array}\right] &= \left[\begin{array}{rrrr} 1 & 0 & 0 & 0 \\ -3 & 1 & 0 & 0 \\ 4 & -3 & 1 & 0 \\ -1 & 5 & -2 & 1 \end{array}\right] \left[\begin{array}{rrrr} 2 & 3 & -1 & 1 \\ 0 & 1 & -2 & 3 \\ 0 & 0 & 4 & 1 \\ 0 & 0 & 0 & -5 \end{array}\right] \\ A ~~~~~~~~~~~~~~~~~~~~~~ &= ~~~~~~~~~~~~~~~~~L~~~~~~~~~~~~~~~~~~~~~~~~~~U \end{align*}

Redoing the 4 x 4 example in Julia

The Julia syntax for entering a matrix and assigning it to a variable is pretty similar to Matlab.

In [1]:
A = [2 3 -1 1 ; -6 -8 1 0 ; 8 9 6 -4 ; -2 2 -17 7]
Out[1]:
4×4 Array{Int64,2}:
  2   3   -1   1
 -6  -8    1   0
  8   9    6  -4
 -2   2  -17   7

Here's a cool Juliaism: Variable names can use Unicode characters, e.g. Greek letters, subscripts, superscripts, etc. (This is true for function names and operators, as well). You can get $L_1$ as a variable name by typing L\_1<TAB> (where <TAB> means hitting the tab button).

In [2]:
L₁ = [1 0 0 0 ; 3 1 0 0; -4 0 1 0 ; 1 0 0 1]
Out[2]:
4×4 Array{Int64,2}:
  1  0  0  0
  3  1  0  0
 -4  0  1  0
  1  0  0  1

Now compute the product of $L_1$ and $A$ and assign it to $X$.

In [3]:
X = L₁ * A
Out[3]:
4×4 Array{Int64,2}:
 2   3   -1   1
 0   1   -2   3
 0  -3   10  -8
 0   5  -18   8

$X$ now has value $X= L_1 A$.

In [4]:
L₂ = [1 0 0 0 ; 0 1 0 0 ; 0 3 1 0 ; 0 -5 0 1]
Out[4]:
4×4 Array{Int64,2}:
 1   0  0  0
 0   1  0  0
 0   3  1  0
 0  -5  0  1
In [5]:
X = L₂ * X
Out[5]:
4×4 Array{Int64,2}:
 2  3  -1   1
 0  1  -2   3
 0  0   4   1
 0  0  -8  -7

$X$ now has value $X= L_2 L_1 A$.

In [6]:
L₃ = [1 0 0 0 ; 0 1 0 0 ; 0 0 1 0 ; 0 0 2 1]
Out[6]:
4×4 Array{Int64,2}:
 1  0  0  0
 0  1  0  0
 0  0  1  0
 0  0  2  1
In [7]:
X = L₃ * X
Out[7]:
4×4 Array{Int64,2}:
 2  3  -1   1
 0  1  -2   3
 0  0   4   1
 0  0   0  -5

$X$ now has value $X = U = L_3 L_2 L_1 A$. We've reduced $X$ to upper-tridiagonal form!

In [8]:
U = X
Out[8]:
4×4 Array{Int64,2}:
 2  3  -1   1
 0  1  -2   3
 0  0   4   1
 0  0   0  -5
In [9]:
L = [1 0 0 0 ; -3 1 0 0 ; 4 -3 1 0 ; -1 5 -2 1]
Out[9]:
4×4 Array{Int64,2}:
  1   0   0  0
 -3   1   0  0
  4  -3   1  0
 -1   5  -2  1

Now let's verify that the $L$ matrix obtained by negating the entries of the $L_j$'s and putting 1's on the diagonal is really the inverse of $L_3 L_2 L_1$.

In [10]:
L₃L₂L₁ = L₃ * L₂ * L₁
Out[10]:
4×4 Array{Int64,2}:
  1  0  0  0
  3  1  0  0
  5  3  1  0
 -4  1  2  1
In [11]:
L * L₃L₂L₁
Out[11]:
4×4 Array{Int64,2}:
 1  0  0  0
 0  1  0  0
 0  0  1  0
 0  0  0  1

Yes! Finally, let's verify that $LU$ equals $A$.

In [12]:
L*U
Out[12]:
4×4 Array{Int64,2}:
  2   3   -1   1
 -6  -8    1   0
  8   9    6  -4
 -2   2  -17   7
In [13]:
A
Out[13]:
4×4 Array{Int64,2}:
  2   3   -1   1
 -6  -8    1   0
  8   9    6  -4
 -2   2  -17   7

Yes!

The 4 x 4 case, symbolically

Step j=1. Let's review the 4 x 4 case from a more general perspective. The first multiplication is of this form

\begin{align*} \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ -\ell_{21} & 1 & 0 & 0 \\ -\ell_{31} & 0 & 1 & 0 \\ -\ell_{41} & 0 & 0 & 1 \\ \end{array}\right] \left[\begin{array}{cccc} x_{11} & x_{12} & x_{13} & x_{14} \\ x_{21} & x_{22} & x_{23} & x_{24} \\ x_{31} & x_{32} & x_{33} & x_{34} \\ x_{41} & x_{42} & x_{43} & x_{44} \end{array}\right] &= \left[\begin{array}{cccc} x_{11} & x_{12} & x_{13} & x_{14} \\ 0 & * & * & * \\ 0 & * & * & * \\ 0 & * & * & * \end{array}\right] \\ L_1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ A ~~~~~~~~~~~~~~~~~~ &= ~~~~~~~~~~~~~~~ L_1 A \end{align*}

In this equation we are using $x_{ij}$ to represent potentially nonzero numbers, and $*$ to represent an entry in the right-hand-side matrix that changed from its correspinding $x_{ij}$ value on the left-hand-side. The $[1 \quad 0 \quad 0 \quad 0]$ first row of $L_1$ means the top row doesn't change, so we have the same top row $[x_{11} \quad x_{12} \quad x_{13} \quad x_{14}]$ on both sides.

We want to multiply from the left by an $L_1$ matrix to put zeros in the marked spots in the right-hand-side matrix. What value do we assign to $\ell_{21}$? Well, the $0$ in the $i,j = 2,1$ spot on the right arises from the product of the 2nd row of $L_1$ and the 1st column of $x$'s on the left. Writing that out algebraically, we get

\begin{equation*} -\ell_{21} x_{11} + 1 \cdot x_{21} = 0 \end{equation*}

Solving for $\ell_{21}$ gives

\begin{equation*} \ell_{21} = \frac{x_{12}}{x_{11}} \end{equation*}

The other $\ell_{i1}$'s work out similarly:

\begin{equation*} \ell_{i1} = \frac{x_{i1}}{x_{11}} \end{equation*}

Step j=2

The second lower-triangular matrix $L_2$ introduces zeros in the subdiagonal elements of the 2nd column. Here the $x$'s in the left-hand matrix r represent whatever values appear on the right-hand-side of the above $L_1 A$ multiplication.

\begin{align*} \left[\begin{array}{cccc} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & -\ell_{32} & 0 & 1 & 0 \\ 0 & -\ell_{42} & 0 & 0 & 1 \\ \end{array}\right] \left[\begin{array}{cccc} x_{11} & x_{12} & x_{13} & x_{14} \\ 0 & x_{22} & x_{23} & x_{24} \\ 0 & x_{32} & x_{33} & x_{34} \\ 0 & x_{42} & x_{43} & x_{44} \end{array}\right] &= \left[\begin{array}{cccc} x_{11} & x_{12} & x_{13} & x_{14} \\ 0 & x_{22} & x_{23} & x_{24} \\ 0 & 0 & * & * \\ 0 & 0 & * & * \end{array}\right] \\ L_2 ~~~~~~~~~~~~~~~~~~~~~~~~~~~ L_1 A ~~~~~~~~~~~~~~~~~~~~~~~ &= ~~~~~~~~~~ L_2 L_1 A \end{align*}

Similarly to the $j=1$ step, solving this system for $\ell_{i2}$ gives

\begin{equation*} \ell_{i2} = \frac{x_{i2}}{x_{22}} \end{equation*}

The general case

In general, we proceed stepping through columns $j=1$ to $j=m-1$, and choosing multipliers $\ell_{ij}$ for rows $i$ from $j+1$ to $m$ according to

\begin{equation*} \ell_{ij} = \frac{x_{ij}}{x_{jj}} \end{equation*}

based on the current entries $x_{ij}$ of the matrix being reduced to upper-triangular form. Multiplying $L_j$ times $X = L_{j-1} L_{j-2} \ldots L_2 L_1 A$ to form a new $X' = L_j L_{j-1} L_{j-2} \ldots L_2 L_1 A$ is equivalent to updating the $x_{ij}$ according to the formula

\begin{equation*} x'_{ik} = x_{ik} - \ell_{ij} x_{jk} \end{equation*}

for $i$ from $j+1$ to $m$ and $k$ from $j$ to $m$. After the stepping through this procedure from $j=1$ to $j=m-1$, the multipliers $\ell_{ij}$ form the subdiagonal elements of $L$ and the upper-tridiagonal matrix $U=X$, forming the LU decomposition $A=LU$.

Pseudo code for LU decomposition without pivoting.

Pseudo code is an algorithm written out as a series of instructions in no particular programming language, using mathematical notation. Your first job for HW3 is to write an LU decomposition in Julia by translating the following pseudo code into Julia. You will have to add the necessary Julia syntax for declaring a Julia function, specifying its input arguments, checking that inputs are of proper form, doing the calculation, and then returning the proper outputs. You will have to translate some English-language expressions and some mathematical notation into proper Julia syntax.

In [14]:
# Input: m x m matrix A
# Outputs: lower-triangular matrix L and upper-triangular matrix U such that A = LU

# let X = A
# let L be an m x m matrix of zeros
#
# for j = 1 to m         (loop over cols of X)
#   ℓⱼⱼ = 1
#   for i = j+1 to m     (loop over rows of jth col of L)
#     ℓᵢⱼ = xᵢⱼ ÷ xⱼⱼ
#     for k = j to m     (update ith row of X)
#        xᵢₖ = xᵢₖ - ℓᵢⱼ xⱼₖ
#     end
#   end
# end
# let U = X

Hints for pseudo-code translation

  • In Julia you get the size of a matrix via (m,n) = size(A).
  • Mathematical notation like $\ell_{ij}$ for an element of matrix $L$ must be translated to L[i,j] in Julia.
  • English "1 to m" is 1:m in Julia
In [ ]: