<jupyter_start><jupyter_text>Gaussian elimination and LU decompositionSay we want to compute the solution of$$Ax = b$$for the vector $x$. We learn how to do this by transforming it to the problem of solving$$U x = y$$where $U$ is an upper-triangular matrix obtained by performing Gaussian elimiantion on $A$ and $y$ is obtained by performing the same operations on $b$. We can then use back substitution to solve $Ux=y$ more easily than solving $Ax=b$ directly.This approach is directly related to the LU decomposition of a matrix, where we wish to factor a matrix $A$ into a product of a lower triangular matrix $L$ and an upper triangular matrix $U$ to give $A = LU$. To understand how to compute the LU decomposition of a matrix, let us start by reminding ourselves of how to do Gaussian elimination. Gaussian elimination by handTo start, consider the following 3x3 matrix$$ A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 10 \end{bmatrix}$$1. Use Gaussian elimination to transform this by hand to an upper triangular matrix $U$ (in row echelon form). Record each elementary row operation you perform along the way.2. Apply the same sequence of row operations to the vector$$b = \begin{bmatrix} 1 \\ 2 \\ 3 \end{bmatrix}$$to obtain the transformed vector $y$.3. Use back substitution to solve $U x = y$.  SolutionUsing the standard Gaussian elimination algorithm, we would perform the following steps:1. Subtract 4x(row 1) from row 2. This leaves us with $$\begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & -6 \\ 7 & 8 & 10  \end{bmatrix}$$2. Subtract 7x(row 1) from row 3. This leaves us with $$\begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & -6 \\ 0 & -6 & -11  \end{bmatrix}$$3. Subtract 2x(row 2) from row 3. This leaves us with $$\begin{bmatrix} 1 & 2 & 3 \\ 0 & -3 & -6 \\ 0 & 0 & 1 \end{bmatrix}$$We now have an upper-triangular matrix $U$. Applying the same sequence of operations to $b$:1. Subtract 4x(row 1) from row 2. This leaves us with $$\begin{bmatrix} 1 \\ -2 \\ 3\end{bmatrix}$$2. Subtract 7x(row 1) from row 3. This leaves us with $$\begin{bmatrix} 1 \\ -2 \\ -4 \end{bmatrix}$$3. Subtract 2x(row 2) from row 3. This leaves us with $$\begin{bmatrix} 1 \\ -2 \\ 0 \end{bmatrix}$$Finally, we use backsubstitution to solve $Ux = y$ for x. Starting with the last entry$$ x_n = 0 / 1 = 0$$$$ x_{n-1} = \frac{-2 - (-6)(0)}{-3} = \frac23$$$$ x_{n-2} = \frac{1 - (2)(\frac23) - (3)(0)}{1} = -\frac13$$so we have the solution$$x = \begin{bmatrix} -\frac13 \\ \frac23 \\ 0\end{bmatrix}$$  Gaussian elimination in PythonWe will now transform the previous algorithm into Python code. First of all we define the matrix $A$ and the vector $b$.<jupyter_code>import numpy as np

A = np.array([[1,2,3],[4,5,6],[7,8,10]])
b = np.array([1,2,3])
n = 3<jupyter_output><empty_output><jupyter_text>Now perform Gaussian elimination and store the result in a matrix $U$ and a vector $y$. We keep track of the multiplication factors for each step in a matrix $L$.<jupyter_code>U = np.array(A)
y = np.array(b)
L = np.identity(n)

# loop over columns and transform to get zeros below the pivot
for k in range(0,n):
    # loop over all rows below the pivot
    for i in range(k + 1, n):
        # Store the multiplication factors in the matrix L
        L[i,k] = U[i,k] / U[k,k]
        
        # Subtract a multiple of row k from row i
        # for j in range(k, n):
        #    U[i,j] = U[i,j] - L[i,k] * U[k,j]
        U[i,:] = U[i,:] - L[i,k] * U[k,:]
        y[i] = y[i] - L[i,k] * y[k]
U
y<jupyter_output><empty_output><jupyter_text>If we consider how many operations this took, there are: ($n$ iterations of the outer loop) x ($n-(k+1)$ iterations of the inner loop) x (n multiplications for the subtraction). This means we require $\mathcal{O}(n^3)$ operations for the Gaussian elimination step. Let us now solve for $x$ using back substitution on $U x = y$.<jupyter_code>x = np.zeros(U.shape[1])

# Start with the last entry in x
x[n-1] = y[n-1] / U[n-1,n-1]

# Iterate upwards from the second last entry to the first entry
for i in range(n-2,-1,-1):
    # Subtract all of the previously computed values from y, then divide by U[i,i]
    #sum = 0
    #for j in range(i+1,n):
    #    sum += U[i,j] * x[j]
    x[i] = (y[i] - U[i,i+1:n]@x[i+1:n])/U[i,i]
x<jupyter_output><empty_output><jupyter_text>We can check that our original matrix is given by $A=LU$:<jupyter_code>L@U<jupyter_output><empty_output><jupyter_text>Gaussian elimination by matrix multiplicationWe could consider each of the steps in Gaussian elimination in terms of multiplication on the left by a sequence of *elementary elimination matrices*. These come in three forms:1. Multiplying row $i$ by a scalar $c$: $\mathbf{r}_i \to c\, \mathbf{r}_i$. This is equivalent to pre-multiplying by a matrix with $1$'s along the diagonal and c in the $i$-th diagonal,$$E_1(i, c) = \begin{bmatrix}  1 &        &   &   &   &        &   \\    & \ddots &   &   &   &        &   \\    &        & 1 &   &   &        &   \\    &        &   & c &   &        &   \\    &        &   &   & 1 &        &   \\    &        &   &   &   & \ddots &   \\    &        &   &   &   &        & 1\end{bmatrix}$$Note that the inverse is given by $E_1(c)^{-1} = E_1(c^{-1})$.2. Add a multiple $c$ of row $j$ to row $i$: $\mathbf{r}_i \to \mathbf{r}_i  + c\, \mathbf{r}_j$. This is equivalent to premultiplying by a matrix with $1$'s along the diagonal and $c$ in $(i, j)$-th entry:$$E_2(i,j,c) = \begin{bmatrix}  1 &        &   &        &   &        &   \\    & \ddots &   &        &   &        &   \\    &        & 1 &        &   &        &   \\    &        &   & \ddots &   &        &   \\    &        & c &        & 1 &        &   \\    &        &   &        &   & \ddots &   \\    &        &   &        &   &        & 1\end{bmatrix}$$In this case the inverse is given by $E_2(c)^{-1} = E_2(-c)$.3. Interchanging rows $i$ and $j$: $\mathbf{r}_i \leftrightarrow \mathbf{r}_j$. This is equivalent to pre-multiplying by a matrix which is the identity with rows $i$ and $j$ swapped: $$E_3(i,j) = \begin{bmatrix}  1 &        &   &        &   &        &   \\    & \ddots &   &        &   &        &   \\    &        & 0 &        & 1 &        &   \\    &        &   & \ddots &   &        &   \\    &        & 1 &        & 0 &        &   \\    &        &   &        &   & \ddots &   \\    &        &   &        &   &        & 1\end{bmatrix}$$In this case the $E_3$ is a permutation matrix and it is its own inverse $E_3^{-1} = E_3$. Let's work out the sequence of elimination matrices we need to perform the Gaussian elimination from the previous example. First, we define Python functions produce each of the three types of elimination matrix:<jupyter_code>def E1(i, c):
    e1 = np.identity(n)
    e1[i, i] = c
    return e1

def E2(i, j, c):
    e2 = np.identity(n)
    e2[i, j] = c
    return e2

def E3(i, j):
    e3 = np.identity(n)
    e3[i, i] = 0
    e3[j, j] = 0
    e3[i, j] = 1
    e3[j, i] = 1
    return e3<jupyter_output><empty_output><jupyter_text>Now, we can see that the Gaussian elimination steps correspond to$$ U = E_2(2,1,-2) E_2(2,0,-7) E_2(1,0,-4) A$$<jupyter_code>E2(1,0,-4)@A
E2(2,0,-7)@E2(1,0,-4)@A
E2(2,1,-2)@E2(2,0,-7)@E2(1,0,-4)@A<jupyter_output><empty_output><jupyter_text>We therefore have$$\begin{aligned}A &= [E_2(2,1,-2) E_2(2,0,-7) E_2(1,0,-4)]^{-1} U \\  &= E_2(1,0,-4)^{-1} E_2(2,0,-7)^{-1} E_2(2,1,-2)^{-1} U \\  &= E_2(1,0,4) E_2(2,0,7) E_2(2,1,2) U \\  &= L U\end{aligned}$$so we have $L$ in terms of elementry elimination matrices.<jupyter_code>E2(1,0,4)@E2(2,0,7)@E2(2,1,2)
L<jupyter_output><empty_output><jupyter_text>LU decomposition and rank-1 matricesIn the lecture videos we emphasized the idea of matrix multiplication in terms of columns-times-rows and the related idea of breaking a matrix into a sum of rank-1 matrices. Now, let's see how this gives a different way of looking at the LU decomposition.The idea is that we would like to split $A$ into a rank-1 piece that picks out the first row and first column, plus a rank-1 piece that picks out the next row and column, and so on:$$\begin{aligned}A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 10 \end{bmatrix} &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & \_ & \_ \\ 7 & \_ & \_ \end{bmatrix}   + \begin{bmatrix} 0 & 0 & 0 \\ 0 & \_ & \_ \\ 0 & \_ & \_ \end{bmatrix}    + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \_ \end{bmatrix}\end{aligned}$$We can fill in all the blanks here by insisting that each term is rank-1 and that we recover $A$. Doing so, we get$$\begin{aligned}A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 10 \end{bmatrix} &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & \_ & \_ \\ 7 & \_ & \_ \end{bmatrix}   + \begin{bmatrix} 0 & 0 & 0 \\ 0 & \_ & \_ \\ 0 & \_ & \_ \end{bmatrix}    + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \_ \end{bmatrix}\\   &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & 8 & 12 \\ 7 & 14 & 21 \end{bmatrix}   + \begin{bmatrix} 0 & 0 & 0 \\ 0 & \_ & \_ \\ 0 & \_ & \_ \end{bmatrix}    + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \_ \end{bmatrix} \quad \text{(rank-1)}\\   &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & 8 & 12 \\ 7 & 14 & 21 \end{bmatrix}   + \begin{bmatrix} 0 & 0 & 0 \\ 0 & -3 & -6 \\ 0 & -6 & \_ \end{bmatrix}    + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \_ \end{bmatrix} \quad \text{(=$A$)}\\   &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & 8 & 12 \\ 7 & 14 & 21 \end{bmatrix}   + \begin{bmatrix} 0 & 0 & 0 \\ 0 & -3 & -6 \\ 0 & -6 & -12 \end{bmatrix}    + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \_ \end{bmatrix} \quad \text{(rank-1)}\\   &= \begin{bmatrix} 1 & 2 & 3 \\ 4 & 8 & 12 \\ 7 & 14 & 21 \end{bmatrix}   + \begin{bmatrix} 0 & 0 & 0 \\ 0 & -3 & -6 \\ 0 & -6 & -12 \end{bmatrix}    + \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad \text{(=$A$)} \\   &= \begin{bmatrix} 1 \\ 4 \\ 7 \end{bmatrix} \begin{bmatrix} 1 & 2 & 3 \end{bmatrix}    + \begin{bmatrix} 0 \\ 1 \\ 2 \end{bmatrix} \begin{bmatrix} 0 & -3 & -6 \end{bmatrix}    + \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix} \begin{bmatrix} 0 & 0 & 1 \end{bmatrix} \\    &= l_1 u_1{}^T + l_2 u_2{}^T + l_3 u_3{}^T = LU\end{aligned}$$<jupyter_code>l1  = L[:,0:1]
u1T = U[0:1]
l2  = L[:,1:2]
u2T = U[1:2]
l3  = L[:,2:3]
u3T = U[2:3]
l1@u1T
l2@u2T
l3@u3T
l1@u1T + l2@u2T + l3@u3T<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup