Gaussian elimination with elimination matrices

In [56]: 
import numpy as np
import numpy.linalg as la
In [57]: 
n = 3

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

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

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

In [58]: 
U = A.copy()

Now eliminate U[1,0]:

In [59]: 
M1 = np.eye(n)
M1[1,0] = -U[1,0]/U[0,0]
M1
Out[59]:
array([[ 1. ,  0. ,  0. ],
       [-1.5,  1. ,  0. ],
       [ 0. ,  0. ,  1. ]])
In [60]: 
U = M1.dot(U)
U
Out[60]:
array([[-2. ,  2. , -1. ],
       [ 0. , -2. , -7.5],
       [-5. , -5. , -2. ]])

Now eliminate U[2,0]:

In [61]: 
M2 = np.eye(n)
M2[2,0] = -U[2,0]/U[0,0]
In [62]: 
U = np.dot(M2, U)
U
Out[62]:
array([[ -2. ,   2. ,  -1. ],
       [  0. ,  -2. ,  -7.5],
       [  0. , -10. ,   0.5]])

Now eliminate U[2,1]:

In [63]: 
M3 = np.eye(n)
M3[2,1] = -U[2,1]/U[1,1]
In [64]: 
U = M3.dot(U)
U
Out[64]:
array([[ -2. ,   2. ,  -1. ],
       [  0. ,  -2. ,  -7.5],
       [  0. ,   0. ,  38. ]])

Try inverting one of the Ms:

In [65]: 
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 [66]: 
U2 = M3.dot(M2.dot(M1.dot(A)))
U2
Out[66]:
array([[ -2. ,   2. ,  -1. ],
       [  0. ,  -2. ,  -7.5],
       [  0. ,   0. ,  38. ]])
In [67]: 
U
Out[67]:
array([[ -2. ,   2. ,  -1. ],
       [  0. ,  -2. ,  -7.5],
       [  0. ,   0. ,  38. ]])

Now define L:

In [68]: 
L = la.inv(M1).dot(la.inv(M2).dot(la.inv(M3)))
L
Out[68]:
array([[ 1. ,  0. ,  0. ],
       [ 1.5,  1. ,  0. ],
       [ 2.5,  5. ,  1. ]])

Observations? (Shape? Diagonal values?)

In [69]: 
np.dot(L, U) - A
Out[69]:
array([[ 0.,  0.,  0.],
       [ 0.,  0.,  0.],
       [ 0.,  0.,  0.]])