# 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]]

In [ ]: