The $LU$ decomposition is closely related to gaussian elimination. It takes the original equation to be solved $A x = b$ and splits it up into two separate equations involving a unit lower triangular matrix $L$, and the row echelon matrix $U$:
\[ \begin{aligned} L y &= b \\ U x &= y \end{aligned} \]
where $A = LU$. The $L$ matrix is a unit lower triangular matrix and thus has ones on the diagonal, whereas $U$ is in row echelon form with pivot values in the leading coefficients of each row.
\[ A = \left( \begin{matrix} 1 & 0 & 0 & 0 \\ \ast & 1 & 0 & 0 \\ \ast & \ast & 1 & 0 \\ \ast & \ast & \ast & 1 \end{matrix}\right) \left(\begin{matrix} p_1 & \ast & \ast & \ast \\ 0 & p_2 & \ast & \ast \\ 0 & 0 & p_3 & \ast \\ 0 & 0 & 0 & p_4 \end{matrix}\right) \]
Thus, we have converted our original problem of solving $A x = b$ into two separate solves, first solving the equation $L y = b$, and then using the result $y$ to solve $U x = y$.
\[ \begin{aligned} \left(\begin{matrix} 1 & 0 & 0 & 0 \\ \ast & 1 & 0 & 0 \\ \ast & \ast & 1 & 0 \\ \ast & \ast & \ast & 1 \end{matrix}\right) \left(\begin{matrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{matrix}\right) &= \left(\begin{matrix} b_1 \\ b_2 \\ b_3 \\ b_4 \end{matrix}\right) \\ \left(\begin{matrix} p_1 & \ast & \ast & \ast \\ 0 & p_2 & \ast & \ast \\ 0 & 0 & p_3 & \ast \\ 0 & 0 & 0 & p_4 \end{matrix}\right) \left(\begin{matrix} x_1 \\ x_2 \\ x_3 \\ x_4 \end{matrix}\right) &= \left(\begin{matrix} y_1 \\ y_2 \\ y_3 \\ y_4 \end{matrix}\right) \end{aligned} \]
However, each of those solves is very cheap to compute, in this case for the 4x4 matrix shown above the solution of $L y = b$ only needs 6 multiplication and 6 additions, whereas $U x = y$ requires 4 divisions, 6 multiplications and 6 additions, leading to a total of 28 arithmetic operations, much fewer in comparison with the 62 operations required to solve the original equation $A x = b$. In general, $LU$ decomposition for an $n \times n$ matrix takes about $2 n^3 / 3$ flops, or floating point operations, to compute.
A relativly simple $LU$ algorithm can be described if we assume that no pivoting is required during a gaussian elimination. In this case, the gaussian elimination process is a sequence of $p$ linear operations $E_1, E_2, ..., E_p$, with each operation $E_i$ being a row replacement that adds a multiple of one row to another below it (i.e. $E_i$ is lower triangular). The final matrix after applying the sequence of row reductions is $U$ in row echelon form, that is:
\[ E_p \cdots E_2 E_1 A = U \]
Since we have $A = LU$, we can show that the sequence of operations $E_1, E_2, ..., E_p$ is also the sequence that reduces the matrix $L$ to an identity matrix:
\[ A = (E_p \cdots E_2 E_1)^{-1} U = LU, \]
therefore,
\[ L = (E_p \cdots E_2 E_1)^{-1}, \]
and,
\[ (E_p \cdots E_2 E_1) L = (E_p \cdots E_2 E_1) (E_p \cdots E_2 E_1)^{-1} = I \]
This implies how we can build up the matrix $L$. We choose values for $L$ such that the series of row operations $E_1, E_2, ..., E_p$ convert the matrix $L$ to the identity matrix. Since each $E_i$ is lower triangular, we know that both $(E_p \cdots E_2 E_1)$ and $(E_p \cdots E_2 E_1)^{-1}$ are also lower triangular.
For example, consider the following matrix
\[ A = \left(\begin{matrix} 3 & 2 & 1 & -3 \\ -6 & -2 & 1 & 5 \\ 3 & -4 & -7 & 2 \\ -9 & -6 & -1 & 15 \end{matrix}\right) \]
After three row reductions, $R_2 \mathrel{{+}{=}} 2 R_1$, $R_3 \mathrel{{+}{=}} -1 R_1$, and $R_3 \mathrel{{+}{=}} 3 R_1$, we have the following result:
\[ E_1 E_2 E_3 A = \left(\begin{matrix} 3 & 2 & 1 & -3 \\ 0 & 2 & * & * \\ 0 & -6 & * & * \\ 0 & 0 & * & * \end{matrix}\right) \]
To build the 1st column of $L$, we simply divide the 1st column of $A$ by the pivot value 3, giving
\[ L = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 1 & * & 1 & 0 \\ -3 & * & * & 1 \end{matrix}\right) \]
For the next column we do the same, using the new pivot value $A_{2,2} = 2$ in row 2 to reduce $A_{3,2}$ and $A_{4,2}$ to zero, and then dividing the column vector under the pivot $(-6, 0)^T$ by the pivot value 2 to obtain the next column of $L$.
Repeating this process for all the columns in $A$, we obtain the final factorisation. You can verify for yourself that repeating the same row operations we did to form $U$ to the matrix $L$ reduces it to the identity matrix.
\[ L = \left(\begin{matrix} 1 & 0 & 0 & 0 \\ -2 & 1 & 0 & 0 \\ 1 & -3 & 1 & 0 \\ -3 & 0 & 2 & 1 \end{matrix}\right) \]
\[ E_1 E_2 ... E_p A = U = \left(\begin{matrix} 3 & 2 & 1 & -3 \\ 0 & 2 & 3 & -1 \\ 0 & 0 & 1 & 2 \\ 0 & 0 & 0 & 2 \end{matrix}\right) \]
Of course, for any practial $LU$ factorisation we need to consider pivoting. Any matrix $A$ can be factorised into $PLU$, where $P$ is a permutation matrix, and $L$ and $U$ are defined as before. During the gaussian elimination steps we store an array of row indices $p_i$ indicating that row $i$ is interchanged with row $p_i$, and the resultant array of $p_i$ can be used to build the permutation matrix $P$ (It would be wasteful to store the entire martix $P$ so the array $p_i$ is stored instead).
Thus, the LU algorithm proceeds as follows:
In practice, most library implementation store $L$ and $U$ in the same matrix since they are lower and upper triangular respectivly.
It is often very benificial when solving linear systems to consider and take advantage of any special structure that the matrix $A$ might possesses. The $LDL$ decomposition is a varient on LU decomposition which is only applicable to a symmetric matrix $A$ (i.e. $A = A^T$). The advantage of using this decomposition is that it takes advantage of the redundent entries in the matrix to reduce the amount of computation to $n^3/3$, which is about a half that required for the $LU$ decomposition.
Take your gaussian elimination code that you wrote in the previous lesson and use it to
write an LU decomposition function that takes in a martix $A$, and returns $L$, $U$ and
the array $p_i$. You can check your answer using
scipy.linalg.lu_factor.
Hint: the permutation matrix is tricky to construct, so you might want to use the test
given in the documentation for lu_factor.
import numpy as np
import matplotlib.pylab as plt
import scipy
import scipy.linalg
import sys
def lu_decomposition(A):
m, n = A.shape
LU = np.copy(A)
pivots = np.empty(n, dtype=int)
# initialise the pivot row and column
h = 0
k = 0
while h < m and k < n:
# Find the k-th pivot:
pivots[k] = np.argmax(LU[h:, k]) + h
if LU[pivots[k], k] == 0:
# No pivot in this column, pass to next column
k = k+1
else:
# swap rows
LU[[h, pivots[k]], :] = LU[[pivots[k], h], :]
# Do for all rows below pivot:
for i in range(h+1, m):
f = LU[i, k] / LU[h, k]
# Store f as the new L column values
LU[i, k] = f
# Do for all remaining elements in current row:
for j in range(k + 1, n):
LU[i, j] = LU[i, j] - LU[h, j] * f
# Increase pivot row and column
h = h + 1
k = k + 1
return LU, pivots
def random_matrix(n):
R = np.random.rand(n, n)
A = np.zeros((n, n))
triu = np.triu_indices(n)
A[triu] = R[triu]
return A
def random_non_singular_matrix(n):
A = np.random.rand(n, n)
while np.linalg.cond(A) > 1/sys.float_info.epsilon:
A = np.random.rand(n, n)
return A
As = [
np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
random_non_singular_matrix(3),
random_non_singular_matrix(4),
random_non_singular_matrix(5),
random_non_singular_matrix(6),
]
def pivots_to_row_indices(pivots):
n = len(pivots)
indices = np.array(range(0, n))
for i, p in enumerate(pivots):
indices[i], indices[p] = indices[p], indices[i]
return indices
def calculate_L_mult_U(LU):
L = np.tril(LU)
np.fill_diagonal(L, 1)
U = np.triu(LU)
return L @ U
for A in As:
LU_scipy, pivots_scipy = scipy.linalg.lu_factor(A)
row_indices_scipy = pivots_to_row_indices(pivots_scipy)
LU_mine, pivots_mine = lu_decomposition(A)
row_indices_mine = pivots_to_row_indices(pivots_mine)
np.testing.assert_almost_equal(calculate_L_mult_U(LU_scipy), A[row_indices_scipy])
np.testing.assert_almost_equal(calculate_L_mult_U(LU_mine), A[row_indices_mine])