Precept 5: LU

In this precept, we will investigate LU factorization. Recall the versions of LU factorization discussed in class:

Task 1

Implement LU factorization without pivoting. In particular, write a function

$$ L, U = \text{lu\_nopivot}(A) $$

that takes an $n \times n$ matrix $A$ and outputs $n \times n$ matrices $L$ and $U$. To save time, only implement LU factorization without pivoting. You can use scipy's $P, L, U = \text{scipy.linalg.lu}(A)$ for the version with partial pivoting in the following. Scipy factorizes $A = PLU$ instead of the more traditional form $ PA = LU \Rightarrow A = P^{-1} LU$

Task 2

Test LU factorization (both with and without pivoting) on $250 \times 250$ matrices $A$ with condition numbers $$ \kappa(A) = 10,10^6 $$ . Check the relative error in the Frobenius norm $$ \text{relerr} = \|LU - A \|_F / \|A\|_F $$ or $$ \text{relerr} = \| LU - PA\|_F / \|A\|_F $$ Repeat each of these tests several times.

import numpy as np
from scipy.linalg import lu as lu

def lu_nopivot(A):
    m = A.shape[0]
    L = np.eye(m)
    U = A.copy()
    ## add your code here 
    return L, U


def rel_err_LU(L,U,A):
    return np.linalg.norm( L @ U - A) / np.linalg.norm(A)

def rel_err_pivoted_LU(L,U,P, A):
    return np.linalg.norm( P@ L @ U - A) / np.linalg.norm(A)

# Change n 
n = 250
# Create random well conditioned matrix
U, _,Vt = np.linalg.svd(np.random.randn(n,n), full_matrices = False)

# Set singular values
S = np.diag(np.linspace(1,10,n))
A_well_cond = U@S@Vt
print('condition number of A: ', np.linalg.cond(A_well_cond))


## Test your lu and built-in pivoted LU

print('relative error without pivot:', rel_err_LU(L,U,A_well_cond))

print('relative error with pivot:', rel_err_pivoted_LU(L,U,P, A_well_cond))


# Set singular values
S = np.diag(10**np.linspace(0,6,n))
A_bad_cond = U@S@Vt
print('condition number of A: ', np.linalg.cond(A_bad_cond))

## Test your lu and built-in pivoted LU
print('relative error without pivot:', rel_err_LU(L,U,A_bad_cond))

print('relative error with pivot:', rel_err_pivoted_LU(L,U,P, A_bad_cond))

Task 3

For certain examples, partial pivoting has no effect and the values in $U$ can become massive. Consider the following matrix $$ A = \left( \begin{array}{rrrrrr} 1 & & & & 1 \\ -1 & 1 & & & 1 \\ -1 & -1& 1& & 1 \\ -1 & -1& -1& 1& 1 \\ -1 & -1& -1& -1& 1 \\ \end{array} \right) $$ the code to form A is given to you below. Use your $\text{lu\_nopivot}$ code on this matrix for $n = 25,50,100,200$ and check the relative error as before $$ \text{relerr} =\| LU - A\|_F / \|A\|_F $$ Additionally, visualize the matrices $U$ and $L$ using plt.imshow(..., interpolation='none') In general, when working with matrices, it is often useful to visualize them using imxhow if they are too big to print on the screen (say bigger than $6 \times 6$).

import matplotlib.pyplot as plt
def make_pathological_A(n):
    A = -np.ones((n,n))
    # Get the lower triangular part of A (zero are 0s)
    A = np.tril(A)
    # Fill diagonal with ones
    np.fill_diagonal(A, np.ones(n))

    A[:,-1] = np.ones(n)
    return A
n = 25
A = make_pathological_A(n)
# Plot A.
plt.imshow(A, interpolation='none')
plt.colorbar()

# Check relative error, plot L and U, increase n

Task 4

Implement a solver for lower and upper triangular systems of linear equations, and then use LU to solve linear systems. Initially, test your solver by creating a $6 \times 6$ linear system with condition number $\kappa(A) = 2$, preforming an $L U$ decomposition using LU with partial pivoting and solving the linear systems. After this code works, work up to testing your solver starting with a $250 \times 250$ linear system with condition number $10^6$.

def back_sub(U, b):
    # Solve Ux = b
    return x


def forward_sub(L, y):
    # Solve Lx = y
    return x



def solve_by_LU(A, b):
    # Compute A = LU, then solve the linear system Ax = b by forward sub and back sub
    return x



n = 6
# Create random well conditioned matrix
U, _,Vt = np.linalg.svd(np.random.randn(n,n), full_matrices = False)
# Set singular values
S = np.diag(np.linspace(1,2,n))
A_well_cond = U@S@Vt
print('condition number of A: ', np.linalg.cond(A_well_cond))

x_true = np.random.randn(n)
b = A_well_cond @ x_true
x = solve_by_LU(A_well_cond, b)
print('error:', np.linalg.norm(x - x_true) / np.linalg.norm(x))