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.]])
U
is the copy of A
that we'll modify:
U = A.copy()
Now eliminate U[1,0]
:
M1 = np.eye(n)
M1[1,0] = -U[1,0]/U[0,0]
M1
array([[ 1. , 0. , 0. ], [-1.5, 1. , 0. ], [ 0. , 0. , 1. ]])
U = M1.dot(U)
U
array([[-2. , 2. , -1. ], [ 0. , -2. , -7.5], [-5. , -5. , -2. ]])
Now eliminate U[2,0]
:
M2 = np.eye(n)
M2[2,0] = -U[2,0]/U[0,0]
U = np.dot(M2, U)
U
array([[ -2. , 2. , -1. ], [ 0. , -2. , -7.5], [ 0. , -10. , 0.5]])
Now eliminate U[2,1]
:
M3 = np.eye(n)
M3[2,1] = -U[2,1]/U[1,1]
U = M3.dot(U)
U
array([[ -2. , 2. , -1. ], [ 0. , -2. , -7.5], [ 0. , 0. , 38. ]])
Try inverting one of the M
s:
print(M2)
print(la.inv(M2))
[[ 1. 0. 0. ] [ 0. 1. 0. ] [-2.5 0. 1. ]] [[ 1. -0. -0. ] [ 0. 1. 0. ] [ 2.5 0. 1. ]]
So we've built M3*M2*M1*A=U
. Test:
U2 = M3.dot(M2.dot(M1.dot(A)))
U2
array([[ -2. , 2. , -1. ], [ 0. , -2. , -7.5], [ 0. , 0. , 38. ]])
U
array([[ -2. , 2. , -1. ], [ 0. , -2. , -7.5], [ 0. , 0. , 38. ]])
Now define L
:
L = la.inv(M1).dot(la.inv(M2).dot(la.inv(M3)))
L
array([[ 1. , 0. , 0. ], [ 1.5, 1. , 0. ], [ 2.5, 5. , 1. ]])
Observations? (Shape? Diagonal values?)
np.dot(L, U) - A
array([[ 0., 0., 0.], [ 0., 0., 0.], [ 0., 0., 0.]])