# Decomposition Methods


```
%matplotlib inline
```


```
import numpy as np
import numpy.linalg as la
```

In Gaussian elimination we performed row (or column) operations to reduce the linear system

\begin{equation}
  A {\bf x} = {\bf b}
\end{equation}

to a form that was easy to solve, essentially by transforming to an equivalent system where the matrix is triangular.

In *decomposition* methods we instead write the matrix as a product of other matrices in such a way that we have a set of easy-to-solve problems.

For example, suppose we can write

\begin{equation}
  A = L U
\end{equation}

where $L$ is a lower triangular matrix and $U$ upper triangular. In this case, the original linear system is equivalent to the two systems

\begin{equation}
  A {\bf x} = {\bf b} \qquad \Leftrightarrow \qquad \left\{ \begin{aligned} L {\bf y} & = {\bf b}, \\ U {\bf x} & = {\bf y}. \end{aligned} \right. 
\end{equation}

As both $L$ and $U$ are triangular, the two systems they define can easily be solved using forwards and backwards substitution respectively. Therefore, *if* we can decompose the matrix $A$ as the product of triangular matrices, the linear system is easily solved.

## Example

Consider the system

\begin{equation}
  \begin{pmatrix} 2 & 1 & -1 \\ 4 & 1 & 0 \\ -2 & -3 & 8 \end{pmatrix} {\bf x} = \begin{pmatrix} 0 \\ 6 \\ 12 \end{pmatrix}.
\end{equation}

We will see later that the matrix can be decomposed as

\begin{equation}
  \begin{pmatrix} 2 & 1 & -1 \\ 4 & 1 & 0 \\ -2 & -3 & 8 \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ -1 & 2 & 1 \end{pmatrix} \begin{pmatrix} 2 & 1 & -1 \\ 0 & -1 & 2 \\ 0 & 0 & 3 \end{pmatrix}.
\end{equation}

We can quickly check that.


```
A = np.array([[2.0, 1.0, -1.0] ,[4.0,  1.0, 0.0], [-2.0, -3.0, 8.0]])
L = np.array([[1.0, 0.0,  0.0] ,[2.0,  1.0, 0.0], [-1.0,  2.0, 1.0]])
U = np.array([[2.0, 1.0, -1.0] ,[0.0, -1.0, 2.0], [ 0.0,  0.0, 3.0]])
np.dot(L, U)
```




    array([[ 2.,  1., -1.],
           [ 4.,  1.,  0.],
           [-2., -3.,  8.]])



Then solve $L {\bf y} = {\bf b}$ by forward substitution:

\begin{equation}
  \begin{pmatrix}
    1 & 0 & 0 \\
    2 & 1 & 0 \\
   -1 & 2 & 1
  \end{pmatrix}
  \begin{pmatrix}
    y_1 \\ y_2 \\ y_3
  \end{pmatrix} =
  \begin{pmatrix}
    0 \\ 6 \\ -12
  \end{pmatrix}
  \Rightarrow \left\{
  \begin{aligned}
    y_1 & = 0, \\
    y_2 & = 6, \\
    y_3 & = -24.
  \end{aligned}
  \right.
\end{equation}

Finally solve $U {\bf x} = {\bf y}$ by backward substitution:

\begin{equation}
  \begin{pmatrix}
    2 & 1 & -1 \\
    0 & -1 & 2 \\
    0 & 0 & 3
  \end{pmatrix}
  \begin{pmatrix}
    x_1 \\ x_2 \\ x_3
  \end{pmatrix} =
  \begin{pmatrix}
    0 \\ 6 \\ -24
  \end{pmatrix}
  \Rightarrow \left\{
  \begin{aligned}
    x_1 & = 7, \\
    x_2 & = -22, \\
    x_3 & = -8.
  \end{aligned}
  \right.
\end{equation}

Again we can cross-check:


```
b = np.array([0.0, 6.0, -12.0])
la.solve(A, b)
```




    array([  7., -22.,  -8.])



## Uniqueness


  The factorisation is not unique: $L$ and $U$ together have $n^2+n$
  free coefficients; $A$ only has $n^2$.  Can freely choose $n$
  coefficients.   Two obvious choices are
  
1. Doolittle's factorisation: The diagonal entries of $L$ are 1.
2. Crout's factorisation: The diagonal entries of $U$ are 1.

  Then use the explicit matrix multiplication formula, and work from
  the top left corner. Each coefficient obeys
  
  \begin{equation}
    a_{i j} = \sum_{s=1}^n \ell_{i s} u_{s j} = \sum_{s=1}^{\min(i,j)}
     \ell_{i s} u_{s j}
  \end{equation}
  
  where the last equality holds because $L$ and $U$ are triangular:
  $\ell_{i s} = 0$ for $s > i$ and $u_{s j} = 0$ for $s > j$.
  

## Hand factorisation example


  The previous example had the matrix
  
  \begin{equation}
    A =
    \begin{pmatrix}
      2 & 1 & -1 \\
      4 & 1 & 0 \\
      -2 & -3 & 8
    \end{pmatrix}
  \end{equation}
  
  which we decompose using
  
  \begin{equation}
    a_{i j} = \sum_{s=1}^{\min(i, j)} \ell_{i s} u_{s j}.
  \end{equation}

  First coefficient, $i = j = 1$, obeys
  
  \begin{equation}
    a_{1 1}  = 2 = \ell_{1 1} u_{1 1}.
  \end{equation}
  
  Use our freedom to set some coefficients to choose $\ell_{1 1} = 1$, for example, which gives $u_{1 1} = a_{1 1} = 2$.
  
  
  Now consider the first row of $U$ ($i = 1$, $j$ free) and the first column of $L$ ($i$ free, $j = 1$). The first row of $U$ obeys
  
  \begin{equation}
    a_{1 j} = \ell_{1 1} u_{1 j} = u_{1 j}
  \end{equation}
  
  and the first column of $L$ obeys
  
  \begin{equation}
    a_{i 1} = \ell_{i 1} u_{1 1} = 2 \ell_{i 1}.
  \end{equation}
  
  Therefore we know that
  
  \begin{equation}
    L =
    \begin{pmatrix}
      1 & 0 & 0 \\
      \tfrac{4}{2} & ?? & 0 \\
      \tfrac{-2}{2} & ?? & ??
    \end{pmatrix}, \quad
    U =
    \begin{pmatrix}
      2 & 1 & -1 \\
      0 & ?? & ?? \\
      0 & 0 & ??
    \end{pmatrix}.
  \end{equation}

  Go to second row and column.
  
  \begin{equation}
    a_{2 2} = \ell_{2 1} u_{1 2} + \ell_{2 2} u_{2 2}.
  \end{equation}
  
  Already computed entries $\ell_{2 1} = 2$, $u_{1 2} = 1$. Again set $\ell_{2 2} = 1$, giving
  
  \begin{equation}
    a_{2 2} = 1 = 2 + u_{2 2} \Rightarrow u_{2 2} = -1.
  \end{equation}

  As in the previous step, consider the second row of $U$ and
  second column of $L$, finding
  
  \begin{align}
    u_{2 j}    & = a_{2 j} - \ell_{2 1} u_{1 j}, \\
    \ell_{i 2} & = \left(a_{i 2} - \ell_{i 1} u_{1 2} \right) / u_{2 2},
  \end{align}
  
  which implies that
  
  \begin{align}
    u_{2 3}    & = 0 - 2 \times (-1) = 2, \\
    \ell_{3 2} & = \left(-3 - (-1) \times 1\right) / (-1) = 2.
  \end{align}
  
  We now have the first two rows and columns:
  
  \begin{equation}
    L =
    \begin{pmatrix}
      1 & 0 & 0 \\
      2 & 1 & 0 \\
      -1 & 2 & ??
    \end{pmatrix}, \quad
    U =
    \begin{pmatrix}
      2 & 1 & -1 \\
      0 & -1 & 2 \\
      0 & 0 & ??
    \end{pmatrix}.
  \end{equation}
  
  Continuing as above, use last free choice to set $\ell_{3 3} =
  1$. Finally, as above, compute that $u_{3 3} = 3$. So we finally
  have
  
  \begin{equation}
    L =
    \begin{pmatrix}
      1 & 0 & 0 \\
      2 & 1 & 0 \\
      -1 & 2 & 1
    \end{pmatrix}, \quad
    U =
    \begin{pmatrix}
      2 & 1 & -1 \\
      0 & -1 & 2 \\
      0 & 0 & 3
    \end{pmatrix}.
  \end{equation}

## $LU$ Decomposition Algorithm

* Write down the explicit formula for the matrix coefficients.

  1. Work from the first row / column to the last.
  2. Look at the diagonal entry of $A$: use the freedom to choose the value of either $L$ or $U$'s diagonal entry.
  3. The appropriate row of $U$ and column of $L$ follows from the explicit formula as all other entries are known.
  
* If either $u_{k k}$ or $\ell_{k k}$ are zero the simple algorithm fails; however, an $LU$ decomposition may still exist.
* $U$ is the matrix that would be found using Gaussian elimination. $LU$ decomposition has the advantage that it holds for any ${\bf b}$ in $A {\bf x} = {\bf b}$, whereas with Gaussian elimination the whole algorithm needs repeating.
* As with Gaussian elimination, pivoting for accuracy is necessary for general matrices.

## Code example

Here is an example LU decomposition algorithm, taken from the solution to the worksheets. Note that, as it doesn't use pivoting, it may well fail.


```
def LU_decomposition(A):
    """Perform LU decomposition using the Doolittle factorisation."""
    
    L = np.zeros_like(A)
    U = np.zeros_like(A)
    N = np.size(A, 0)
    
    for k in range(N):
        L[k, k] = 1
        U[k, k] = (A[k, k] - np.dot(L[k, :k], U[:k, k])) / L[k, k]
        for j in range(k+1, N):
            U[k, j] = (A[k, j] - np.dot(L[k, :k], U[:k, j])) / L[k, k]
        for i in range(k+1, N):
            L[i, k] = (A[i, k] - np.dot(L[i, :k], U[:k, k])) / U[k, k]
    
    return L, U
```

Check this code against the example above.


```
L2, U2 = LU_decomposition(A)
print("Compare the code and hand answers for L:\n{},\n{}\nand for U:\n{},\n{}".format(L2, L, U2, U))
```

    Compare the code and hand answers for L:
    [[ 1.  0.  0.]
     [ 2.  1.  0.]
     [-1.  2.  1.]],
    [[ 1.  0.  0.]
     [ 2.  1.  0.]
     [-1.  2.  1.]]
    and for U:
    [[ 2.  1. -1.]
     [ 0. -1.  2.]
     [ 0.  0.  3.]],
    [[ 2.  1. -1.]
     [ 0. -1.  2.]
     [ 0.  0.  3.]]


## Conditions for factorisation

### Background


  We *assumed* an $LU$ decomposition is always possible. However,
  we then noted that the algorithm works only if both $u_{k k}$ and
  $\ell_{k k}$ are non-zero, which is not known in advance. 
  

  Two key concepts used in determining if a factorisation exists.
  
1.  The *principal minor* of order $k$ of a matrix $A$ is the (sub-) matrix
    \begin{equation*}
      \begin{pmatrix}
        a_{1 1} & \cdots & a_{1 k} \\
        \vdots & \ddots & \vdots \\
        a_{k 1} & \cdots & a_{k k}
      \end{pmatrix} .
    \end{equation*}
    
2. A matrix is *strictly diagonally dominant* if
    \begin{equation*}
      | a_{i i} | > \sum_{\substack{j = 1 \\ j \ne i}}^n | a_{i j} |,
        \quad (1 \le i \le n).
    \end{equation*}

### Theorems


  Summarize the two key theorems; see references given in the notes:

  **Theorem 1.1**: If all $n - 1$ leading principal minors of the $n
  \times n$ matrix $A$ are nonsingular, then an $LU$ decomposition of
  $A$ exists. 
            

  The relation between $LU$ decomposition and Gaussian elimination is
  exploited to prove this theorem. It also shows that, if the matrix
  is nonsingular, its rows can be permuted so that an $LU$
  decomposition of the permuted matrix can be found.
  
  **Theorem 1.2}**: Every strictly diagonally dominant matrix is
  nonsingular and has an $LU$ decomposition. 
  

  Diagonal dominance is a very simple condition to check. Most
  matrices involved in numerical methods for solving PDEs are
  diagonally dominant, leading to the practical utility of this
  theorem.


## Alternative decomposition methods

### Cholesky factorisation - example

We want to decompose the symmetric, positive definite matrix

  \begin{equation}
    A =
    \begin{pmatrix}
      2 & 1 \\
      1 & 3
    \end{pmatrix}
  \end{equation}
  
  as $A = L L^T$ where $L$ is lower triangular with positive diagonal
  elements. Again we write the coefficients out explicitly and
  look at the diagonal elements first:
  
  \begin{align}
                && A & = L L^T \\
    \Rightarrow &&
    \begin{pmatrix}
      2 & 1 \\
      1 & 3
    \end{pmatrix}
    & =
    \begin{pmatrix}
      \ell_{1 1}^2 & \ell_{1 1} \ell_{2 1} \\
      \ell_{1 1} \ell_{2 1} & \ell_{2 1}^2 +  \ell_{2 2}^2
    \end{pmatrix} \\ 
    \Rightarrow && &\left\{ 
      \begin{aligned}
        \ell_{1 1} & = \sqrt{2} \\
        \ell_{2 1} & = 1 / \sqrt{2} \\
        \ell_{2 2} & = \sqrt{3 - \ell_{2 1}^2} \\
                   & = \sqrt{5/2}
      \end{aligned}\right. .
  \end{align}


af02945e7a4d3e2442b0e52cf3de0dcb4dd2f7a3
ipynb
Jupyter Notebook
Lectures/03 - Decomposition Methods.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-02-12T19:51:52.000Z
2022-03-26T15:34:11.000Z
Lectures/03 - Decomposition Methods.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2017-05-24T19:49:52.000Z
2018-01-23T21:40:42.000Z
Lectures/03 - Decomposition Methods.ipynb
josh-gree/NumericalMethods
03cb91114b3f5eb1b56916920ad180d371fe5283
CC-BY-3.0
2015-01-05T13:30:47.000Z
2022-02-15T09:59:39.000Z
Qwen/Qwen-72B
1. YES
2. YES