Week 4 Self check questions and solutions

Question 1:

Write a Python function to implement the LU Decomposition with partial pivoting. You can take the following skeleton as starting point.

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

def partial_pivot_lu(A):
    """Implementation of a pivoted LU decomposition."""
    
    import numpy as np
    
    # Return an error if matrix is not square
    if not A.shape[0]==A.shape[1]:
        raise ValueError("Input matrix must be square")
        
    n = A.shape[0] # The number of rows/columns of A
    p = np.arange(n) # The permutation vector
    
    L = np.zeros((n,n),dtype='float64') # Reserve space for L
    U = np.zeros((n,n),dtype='float64') # Reserve space for U
    U[:] = A # Copy A into U as we do not want to modify A
    np.fill_diagonal(L,1) # fill the diagonal of L with 1
    
    for i in range(n-1):
        # IMPLEMENT THE ROW EXCHANGE FOR PARTIAL PIVOTING HERE.
        # STORE THE PERMUTATION IN THE p VECTOR.
        for j in range(i+1,n):
            L[j,i] = U[j,i]/U[i,i]
            U[j,i:] = U[j,i:]-L[j,i]*U[i,i:]
            U[j,i] = 0
    # The vector p stores the forward permutations.
    # The following code takes an identity matrix and swaps the rows accordingly.
    P = np.eye(n)[p,:] 
    return (P,L,U)

n = 10

A = np.random.randn(n, n)
b = np.random.randn(n)

# Example solution with Scipy
# Note that Scipy computes A = PLU and not PA=LU.
P, L, U = lu(A)
x_compare = solve_triangular(
    U, solve_triangular(L, P.T @ b, lower=True))
# Print the residual
print(f"Residual with Scipy LU: {np.linalg.norm(b - A @ x_compare)}")

P2, L2, U2 = partial_pivot_lu(A)
x_pivot = solve_triangular(
    U2, solve_triangular(L2, P2 @ b, lower=True))
print(f"Residual with partial pivot LU: {np.linalg.norm(b - A @ x_pivot)}")

Residual with Scipy LU: 3.0356551825401023e-15
Residual with partial pivot LU: 8.041205793198706e-13

Solution:

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

def partial_pivot_lu(A):
    """Implementation of a pivoted LU decomposition."""
    
    import numpy as np
    
    # Return an error if matrix is not square
    if not A.shape[0]==A.shape[1]:
        raise ValueError("Input matrix must be square")
        
    n = A.shape[0] # The number of rows/columns of A
    p = np.arange(n) # The permutation vector
    
    L = np.zeros((n,n),dtype='float64') # Reserve space for L
    U = np.zeros((n,n),dtype='float64') # Reserve space for U
    U[:] = A # Copy A into U as we do not want to modify A
    np.fill_diagonal(L,1) # fill the diagonal of L with 1
    
    for i in range(n-1):
        # The outer iteration
        # Permute the rows to ensure the element with largest magnitude is on top
        max_index = i+np.argmax(np.abs(U[i:,i]))
        U[[i,max_index],:] = U[[max_index,i],:]
        L[[i,max_index],:i] = L[[max_index,i],:i]
        p[[i,max_index]] = p[[max_index,i]]
        for j in range(i+1,n):
            L[j,i] = U[j,i]/U[i,i]
            U[j,i:] = U[j,i:]-L[j,i]*U[i,i:]
            U[j,i] = 0
    # The vector p stores the forward permutations.
    # The following code takes an identity matrix and swaps the rows accordingly.
    P = np.eye(n)[p,:] 
    return (P,L,U)

n = 10

A = np.random.randn(n, n)
b = np.random.randn(n)

# Example solution with Scipy
# Note that Scipy computes A = PLU and not PA = LU
P, L, U = lu(A)
x_compare = solve_triangular(
    U, solve_triangular(L, P.T @ b, lower=True))
# Print the residual
print(f"Residual with Scipy LU: {np.linalg.norm(b - A @ x_compare)}")

P2, L2, U2 = partial_pivot_lu(A)
x_pivot = solve_triangular(
    U2, solve_triangular(L2, P2 @ b, lower=True))
print(f"Residual with partial pivot LU: {np.linalg.norm(b - A @ x_pivot)}")

Residual with Scipy LU: 1.1123893155135927e-15
Residual with partial pivot LU: 2.7294119545946647e-15