# LU Factorization¶

In [2]:
import numpy as np
import numpy.linalg as la


## Part 1: One Column of LU¶

In [3]:
n = 3

np.random.seed(15)
A = np.round(5*np.random.randn(n, n))

A

Out[3]:
array([[-2.,  2., -1.],
[-3.,  1., -9.],
[-5., -5., -2.]])

Initialize L and U with zeros:

In [4]:
L = np.zeros((n,n))
U = np.zeros((n,n))


Set U to be the first row of A:

In [5]:
U[0,:] = A[0,:]
U

Out[5]:
array([[-2.,  2., -1.],
[ 0.,  0.,  0.],
[ 0.,  0.,  0.]])

Compute the first column of L:

In [6]:
L[:,0] = A[:,0]/U[0,0]
L

Out[6]:
array([[1. , 0. , 0. ],
[1.5, 0. , 0. ],
[2.5, 0. , 0. ]])

Compare what we have to A:

In [7]:
print(A)
print(L@U)

[[-2.  2. -1.]
[-3.  1. -9.]
[-5. -5. -2.]]
[[-2.   2.  -1. ]
[-3.   3.  -1.5]
[-5.   5.  -2.5]]


Perform the Schur complement update and store the result in A1:

In [8]:
A1 = A - L @ U
A1

Out[8]:
array([[  0. ,   0. ,   0. ],
[  0. ,  -2. ,  -7.5],
[  0. , -10. ,   0.5]])

Take the second row of U to be the second row of A1:

In [9]:
U[1,1:] = A1[1,1:]
U

Out[9]:
array([[-2. ,  2. , -1. ],
[ 0. , -2. , -7.5],
[ 0. ,  0. ,  0. ]])

We can now compute the next column of L:

In [10]:
L[1:,1] = A1[1:,1]/U[1,1]
L

Out[10]:
array([[1. , 0. , 0. ],
[1.5, 1. , 0. ],
[2.5, 5. , 0. ]])

And finally, compute the bottom right elements of L and U

In [11]:
U[2,2] = A1[2,2] - L[2,1]*U[1,2]
L[2,2] = 1.0

In [12]:
print(L)
print(U)

[[1.  0.  0. ]
[1.5 1.  0. ]
[2.5 5.  1. ]]
[[-2.   2.  -1. ]
[ 0.  -2.  -7.5]
[ 0.   0.  38. ]]

In [13]:
print(A)
print(L@U)

[[-2.  2. -1.]
[-3.  1. -9.]
[-5. -5. -2.]]
[[-2.  2. -1.]
[-3.  1. -9.]
[-5. -5. -2.]]


## Part 2: The Full Algorithm¶

Implement the general LU factorization algorithm

In [14]:
n = 4
A = np.random.random((n,n))
L = np.zeros((n,n))
U = np.zeros((n,n))
M = A.copy()

In [15]:
for i in range(n):
U[i,i:] = M[i,i:]
L[i:,i] = M[i:,i]/U[i,i]
M[i+1:,i+1:] -= np.outer(L[i+1:,i:i+1],U[i:i+1,i+1:])

In [16]:
print(L)
print(U)
print(A-L@U)

[[  1.           0.           0.           0.        ]
[  0.87952614   1.           0.           0.        ]
[  0.04296101  -8.67490157   1.           0.        ]
[  0.82877669 -11.69108892   1.80202512   1.        ]]
[[ 0.9176299   0.26414685  0.71777369  0.86571503]
[ 0.         -0.02177348 -0.46405769 -0.71471261]
[ 0.          0.         -3.05794766 -5.86446657]
[ 0.          0.          0.          2.43970137]]
[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
[ 0.00000000e+00  0.00000000e+00 -2.22044605e-16 -3.33066907e-16]
[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  2.22044605e-16]]


## Part 3: LU with Partial Pivoting¶

When a divisor (U_ii) becomes zero, while nonzeros remain in the trailing matrix, the LU factorization does not exist. Further, if U_ii is small, the magnitude of elements in L and U can blow up, increasing round-off error. Partial pivoting circumvents both problems.

In [31]:
n = 4
A = np.random.random((n,n))
L = np.zeros((n,n))
U = np.zeros((n,n))
P = np.eye(n,n)
M = A.copy()

def perm(i,j):
P = np.eye(n,n)
P[i,i] = 0
P[j,j] = 0
P[i,j] = 1
P[j,i] = 1
return P

for i in range(n):
imax = np.argmax(np.abs(M[i:,i]))
Pi = perm(i,i+imax)
P = Pi @ P
L = Pi @ L
M = Pi @ M
U[i,i:] = M[i,i:]
L[i:,i] = M[i:,i]/U[i,i]
M[i+1:,i+1:] -= np.outer(L[i+1:,i:i+1],U[i:i+1,i+1:])


Check residual error for the permuted A, i.e., PA-LU. What property do off-diagonal entries in L satisfy?

In [32]:
print(P@A - L @ U)

[[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
[ 0.00000000e+00  1.11022302e-16  0.00000000e+00  0.00000000e+00]
[ 0.00000000e+00  0.00000000e+00 -1.11022302e-16  0.00000000e+00]
[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]]

In [33]:
L

Out[33]:
array([[ 1.        ,  0.        ,  0.        ,  0.        ],
[ 0.93673604,  1.        ,  0.        ,  0.        ],
[ 0.91055425, -0.35441189,  1.        ,  0.        ],
[ 0.26450396, -0.02409714, -0.02133641,  1.        ]])
In [34]:
P

Out[34]:
array([[0., 1., 0., 0.],
[0., 0., 1., 0.],
[0., 0., 0., 1.],
[1., 0., 0., 0.]])
In [ ]: