# Gaussian elimination with elimination matrices

In :
import numpy as np
import numpy.linalg as la

In :
n = 3

np.random.seed(15)
A = np.round(5*np.random.randn(n, n))

A

Out:
array([[-2.,  2., -1.],
[-3.,  1., -9.],
[-5., -5., -2.]])


U is the copy of A that we'll modify:

In :
U = A.copy()


Now eliminate U[1,0]:

In :
M1 = np.eye(n)
M1[1,0] = -U[1,0]/U[0,0]
M1

Out:
array([[ 1. ,  0. ,  0. ],
[-1.5,  1. ,  0. ],
[ 0. ,  0. ,  1. ]])

In :
U = M1.dot(U)
U

Out:
array([[-2. ,  2. , -1. ],
[ 0. , -2. , -7.5],
[-5. , -5. , -2. ]])


Now eliminate U[2,0]:

In :
M2 = np.eye(n)
M2[2,0] = -U[2,0]/U[0,0]

In :
U = np.dot(M2, U)
U

Out:
array([[ -2. ,   2. ,  -1. ],
[  0. ,  -2. ,  -7.5],
[  0. , -10. ,   0.5]])


Now eliminate U[2,1]:

In :
M3 = np.eye(n)
M3[2,1] = -U[2,1]/U[1,1]

In :
U = M3.dot(U)
U

Out:
array([[ -2. ,   2. ,  -1. ],
[  0. ,  -2. ,  -7.5],
[  0. ,   0. ,  38. ]])


Try inverting one of the Ms:

In :
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:

In :
U2 = M3.dot(M2.dot(M1.dot(A)))
U2

Out:
array([[ -2. ,   2. ,  -1. ],
[  0. ,  -2. ,  -7.5],
[  0. ,   0. ,  38. ]])

In :
U

Out:
array([[ -2. ,   2. ,  -1. ],
[  0. ,  -2. ,  -7.5],
[  0. ,   0. ,  38. ]])


Now define L:

In :
L = la.inv(M1).dot(la.inv(M2).dot(la.inv(M3)))
L

Out:
array([[ 1. ,  0. ,  0. ],
[ 1.5,  1. ,  0. ],
[ 2.5,  5. ,  1. ]])


Observations? (Shape? Diagonal values?)

In :
np.dot(L, U) - A

Out:
array([[ 0.,  0.,  0.],
[ 0.,  0.,  0.],
[ 0.,  0.,  0.]])