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 [59]:
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 [61]:
A1 = A - L @ U
A1
Out[61]:
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 [62]:
U[1,1:] = A1[1,1:]
U
Out[62]:
array([[-2. ,  2. , -1. ],
       [ 0. , -2. , -7.5],
       [ 0. ,  0. ,  0. ]])

We can now compute the next column of L:

In [64]:
L[1:,1] = A1[1:,1]/U[1,1]
L
Out[64]:
array([[ 1. ,  0. ,  0. ],
       [ 1.5,  1. ,  0. ],
       [ 2.5,  5. ,  0. ]])

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

In [65]:
U[2,2] = A1[2,2] - L[2,1]*U[1,2]
L[2,2] = 1.0
In [66]:
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 [67]:
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 [7]:
n = 4
A = np.random.random((n,n)) 
L = np.zeros((n,n)) 
U = np.zeros((n,n)) 
M = A.copy()
In [108]:
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 [110]:
print(L)
print(U)
print(A-L@U)
[[ 1.          0.          0.          0.        ]
 [ 0.48283957  1.          0.          0.        ]
 [ 1.0179142   0.08898352  1.          0.        ]
 [ 0.08538441  0.70860694 -0.76449667  1.        ]]
[[ 0.81336415  0.0571982   0.97858297  0.90071912]
 [ 0.          0.86381599  0.00137876  0.32479681]
 [ 0.          0.         -0.66417473 -0.90299951]
 [ 0.          0.          0.         -0.03547261]]
[[  0.00000000e+00   0.00000000e+00   0.00000000e+00   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   0.00000000e+00   0.00000000e+00   1.11022302e-16]]