import sys
########################################
# Change the string in the line below! #
########################################
sys.path.append("/Users/youngmin/Documents/CS/Introduction to Computational Science - F21/CS111-2021-fall/Python")
import os
import time
import math
import numpy as np
import numpy.linalg as npla
import scipy
from scipy import linalg as spla
import scipy.sparse
import scipy.sparse.linalg
from scipy import integrate
import networkx as nx
import cs111
##########################################################
# If this import for matplotlib doesn't work, try saying #
# conda install -c conda-forge ipympl #
# at a shell prompt on your computer #
##########################################################
import matplotlib
%matplotlib ipympl
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import axes3d
np.set_printoptions(precision = 4)
Although matrix inversion is not used for solving $A\mathbf{x} = \mathbf{b}$, there are a few applications where explicit knowledge of $A^{-1}$ is desirable.
We can use the $LU$ factors of a nonsingular $n$-by-$n$ matrix $A$ to compute $A^{-1}$ efficiently. The strategy is to solve the matrix equation $AX = I$. Recall that $AA^{-1} = I$ implies $A[A^{-1}]_{*j} = \mathbf{e}_j$, where $[A^{-1}]_{*j}$ is the $j^{th}$ column of $A^{-1}$, and $\mathbf{e}_j$ is the $j^{th}$ column of the identity matrix. Thus, $[A^{-1}]_{*j}$ is the solution of a system $A\mathbf{x}_j = \mathbf{e}_j$, for $j = 0, 1, ..., n-1$. Each of these $n$ systems has the same coefficient matrix, so, once the $LU$ factors for $A$ are known, each system $A\mathbf{x}_j = LU\mathbf{x}_j = \mathbf{e}_j$ can be solved by the standard two-step process:
Set $\mathbf{y}_j = U\mathbf{x}_j$, and solve $L\mathbf{y}_j = \mathbf{e}_j$ for $\mathbf{y}_j$ by forward substitution.
Solve $U\mathbf{x}_j = \mathbf{y}_j$ for $\mathbf{x}_j = [A^{-1}]_{*j}$ by back substitution.
The method has at least two advantages: it's efficient, and any code written to solve $A\mathbf{x} = \mathbf{b}$ can also be used to compute $A^{-1}$.
Write a function computeInverse( A ) that uses cs111.LUfactor, cs111.Lsolve, and cs111.Usolve to compute the inverse of the input matrix A by following the above algorithm. Here, however, you must take into account a permutation too. You may assume that A will always be square. Note that you should not call cs111.LUfactor more than once within your code.
Test your function on the following matrix:
$$A = \begin{pmatrix} 1 & 1 & 2 \\ 1 & 1 & 3 \\ 2 & 3 & 4 \\ \end{pmatrix} $$by computing $AA^{-1}$ and also by comparing your result to the inverse matrix calculated with numpy's builtin function np.linalg.inv( A ). This is the same matrix we saw in lecture, and we know that the two last rows must permute so that $LU$ factorization succeeds.
def computeInverse( A ):
L, U, p = cs111.LUfactor( A ) # Obtain the LU factors and permutation vector p.
I = np.eye( len(A) ) # Identity matrix with same size as A.
PI = I[p,:] # Identify matrix with permuted rows.
# Now, we want to solve the system(s) AX = I <=> PAX = PI <=> LUX = PI
X = np.zeros( A.shape )
for i in range( len(A) ): # Solving for the columns of X.
rhs = PI[:,i]
y = cs111.Lsolve( L, rhs )
x = cs111.Usolve( U, y )
X[:,i] = x
return X
# Let's test our inverse-computing function.
A = np.array( [[1,1,2], [1,1,3], [2,3,4]] )
X = computeInverse( A )
print( "A:\n", A, sep="" )
print( "\nX:\n", X, sep="" )
print( "\nA @ X:\n", A @ X, sep="" )
print( "\nnpla.inv(A):\n", npla.inv( A ), sep="" )
A: [[1 1 2] [1 1 3] [2 3 4]] X: [[ 5. -2. -1.] [-2. -0. 1.] [-1. 1. -0.]] A @ X: [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] npla.inv(A): [[ 5. -2. -1.] [-2. -0. 1.] [-1. 1. -0.]]
An $n$-by-$n$ matrix $A$ has an $LU$ factorization if and only if it is invertible.
Proof:
(No partial pivoting for simplicity here) First, we prove that if $LU$ factorization succeeds, then $A$ is invertible. To do so, recall that an elementary matrix of the form $M = I + \alpha\mathbf{e}_i \mathbf{e}_j^T$ adds a multiple ($\alpha$) of row $j$ to row $i$ when it premultiplies another matrix. For example, if $\alpha = -1$ and
$$ B = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 8 & 9 & 10 \end{pmatrix}, $$then
$$ M = I + (-1)\mathbf{e}_2\mathbf{e}_0^T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ -1 & 0 & 1 \end{pmatrix} $$and
$$ MB = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 7 & 7 \end{pmatrix}. $$Further, any elementary matrix of the above type (e.g., $M$) is always invertible and $M^{-1} = I - \alpha\mathbf{e}_i\mathbf{e}_j^T$. For example,
$$ M^{-1} = I - (-1)\mathbf{e}_2\mathbf{e}_0^T = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{pmatrix} \quad \textrm{and} \quad M^{-1}M = MM^{-1} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}. $$Thus, we can use elementary matrices like $M$ above to perform Gaussian elimination and transform $A$ into an upper triangular matrix $U$. If no zero pivot is found along the way (i.e., all diagonal elements in $U$ are nonzero), then
$$ M_m M_{m-1} \cdots M_\ell \cdots M_1 M_0 A = U $$Since any $M_\ell$ is invertible, we can write
$$ A = M_0^{-1} M_1^{-1} \cdots M_\ell^{-1} \cdots M_{m-1}^{-1} M_m^{-1}U = LU, $$which proves that $A$ is non singular because $A^{-1} = (LU)^{-1} = U^{-1}L^{-1} = U^{-1} M_m M_{m-1} \cdots M_\ell \cdots M_1 M_0$ exists.
Note: For the explanation involving permutations and partial pivoting, refer to https://piazza.com/class/kt26q57z3e27p3?cid=14.
To prove the converse (if $A$ is invertible, then its $LU$ factorization succeeds), recall that $A$ is invertible if $\textrm{rank}(A) = n$. This implies that all the rows in $A$ span $\mathbb{R}^n$ and no linear combination of $(n-1)$ rows can produce any other row. Therefore, the operations encoded in the elementary matrices $M$ above can never lead us to a zero row (and zero pivot, at least in exact arithmetic). This ensures that $M_m M_{m-1} \cdots M_\ell \cdots M_1 M_0 A = U$, and Gaussian elimination successfuly factors $A$ into $L$ and $U$ (with partial pivoting).
# Here's an example of an elementary matrix of the type above.
M = np.array([[1,0,0],[0,1,0],[-1,0,1]])
print( "M:\n", M, sep="" )
print( "\ninv(M):\n", npla.inv( M ), sep="" )
M: [[ 1 0 0] [ 0 1 0] [-1 0 1]] inv(M): [[1. 0. 0.] [0. 1. 0.] [1. 0. 1.]]
# Here's an example of a noninvertible matrix where LU=PA fails.
A = np.array([[1,4,7], [2,5,8], [3,6,9]])
L, U, p = cs111.LUfactor( A )
print( "L:\n", L, sep="" )
print( "\nU:\n", U, sep="" )
print( "\np:", p )
--------------------------------------------------------------------------- AssertionError Traceback (most recent call last) <ipython-input-5-7cf2afe67cc8> in <module> 1 # Here's an example of a noninvertible matrix where LU=PA fails. 2 A = np.array([[1,4,7], [2,5,8], [3,6,9]]) ----> 3 L, U, p = cs111.LUfactor( A ) 4 5 print( "L:\n", L, sep="" ) ~/Documents/CS/Introduction to Computational Science - F21/CS111-2021-fall/Python/cs111/LU.py in LUfactor(A, pivoting) 141 if pivoting: 142 piv_row = piv_col + np.argmax(np.abs(LU[piv_col:, piv_col])) --> 143 assert LU[piv_row, piv_col] != 0., "can't find nonzero pivot, matrix is singular" 144 LU[[piv_col, piv_row], :] = LU[[piv_row, piv_col], :] 145 p[[piv_col, piv_row]] = p[[piv_row, piv_col]] AssertionError: can't find nonzero pivot, matrix is singular
If Gaussian elimination (or $LU$ factorization with partial pivoting) fails, then $A$ is singular.
Proof: Suppose $A$ is an $n$-by-$n$ matrix. $LU$ factorization fails if in some step the algorithm produces a $\mathbf{0}$ row. In that case, $\textrm{rank}(A) < n$ or its nullspace has at least one nonzero basis vector. Similarly, if we are solving a system of linear equations, this means $A\mathbf{x} = \mathbf{b}$ has either no solution or infinitely many solutions depending on $\mathbf{b}$. Therefore, $A^{-1}$ does not exist.