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 [ ]: