Precept 6: LU and Woodbury¶
Q1: Rank-1 perturbations¶
The Woodbury matrix identity is:
$$(A + uv^T)^{-1} = A^{-1} - \frac{1}{1 + v^T A^{-1} u }A^{-1}u v^T A^{-1}$$ where $A \in \mathbb{R}^{n \times n} $ , $u, v \in \mathbb{R}^n$, and both $A$ and $A + uv^T$ are invertible. (Note that $v^T A^{-1} u $ is a scalar). The Woodbury identity is commonly used in situations where the inverse of $A$ is easy to compute, and we wish to invert a rank-1 perturbation of A.
Prove the Woodbury matrix identity. Hint: $C = B^{-1} \Leftrightarrow BC = I$.
Consider a matrix $B \in \mathbb{R}^{n\times n}$ that is circulant except in the last column. For example, a $5 \times 5$ matrix of this form is:
$$B = \begin{bmatrix} c_0 & c_4 & c_3 & c_2 & l_0 \\ c_1 & c_0 & c_4 & c_3 & l_1 \\ c_2 & c_1 & c_0 & c_4 & l_2 \\ c_3 & c_2 & c_1 & c_0 & l_3 \\ c_4 & c_3 & c_2 & c_1 & l_4 \end{bmatrix} $$
Use the Woodbury matrix identity to write a solver that computes the solution $x$ of $Bx =d$ in $\mathcal{O}(n \log n)$ operations.
import numpy as np
import scipy.linalg
def almost_circulant_solve(c, l, d ):
# c is a vector of size n, it is the first column of B.
# v is a vector of size n, it is the last column of B.
# d is a vector of size n, it is the right hand side vector Bx = d
n = c.size
## ADD YOUR CODE HERE
x = np.zeros(n)
# Should return just a vector
return x
# ALL THE CODE BELOW IS PROVIDED FOR YOU. NO NEED TO CHANGE IT
# Test for accuracy (there is no test for scaling provided but make sure the complexity is correct!)
n = 5
c = np.random.randn(5)
v = np.random.randn(5)
# Form a circulant matrix with first column c
B = scipy.linalg.circulant(c)
# Change last column
B[:,-1] = v
x = np.random.randn(n)
d = B @ x
x_sol = almost_circulant_solve(c, v, d )
print('Q3 error:', np.linalg.norm(x - x_sol) / np.linalg.norm(x))
Programming task¶
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))