import numpy as np
import numpy.linalg as la
n = 3
np.random.seed(15)
A = np.round(5*np.random.randn(n, n))
A
array([[-2., 2., -1.], [-3., 1., -9.], [-5., -5., -2.]])
Initialize L
and U
with zeros:
L = np.zeros((n,n))
U = np.zeros((n,n))
Set U
to be the first row of A
:
U[0,:] = A[0,:]
U
array([[-2., 2., -1.], [ 0., 0., 0.], [ 0., 0., 0.]])
Compute the first column of L
:
L[:,0] = A[:,0]/U[0,0]
L
array([[ 1. , 0. , 0. ], [ 1.5, 0. , 0. ], [ 2.5, 0. , 0. ]])
Compare what we have to A
:
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
:
A1 = A - L @ U
A1
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
:
U[1,1:] = A1[1,1:]
U
array([[-2. , 2. , -1. ], [ 0. , -2. , -7.5], [ 0. , 0. , 0. ]])
We can now compute the next column of L
:
L[1:,1] = A1[1:,1]/U[1,1]
L
array([[ 1. , 0. , 0. ], [ 1.5, 1. , 0. ], [ 2.5, 5. , 0. ]])
And finally, compute the bottom right elements of L
and U
U[2,2] = A1[2,2] - L[2,1]*U[1,2]
L[2,2] = 1.0
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. ]]
print(A)
print(L@U)
[[-2. 2. -1.] [-3. 1. -9.] [-5. -5. -2.]] [[-2. 2. -1.] [-3. 1. -9.] [-5. -5. -2.]]
Implement the general LU factorization algorithm
n = 4
A = np.random.random((n,n))
L = np.zeros((n,n))
U = np.zeros((n,n))
M = A.copy()
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:])
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]]