In [1]:

from __future__ import division

import numpy as np
import numpy.linalg as la

In [2]:

A = np.array([
        [3, 4, 5],
        [2, 3, 1],
        [9, 2, 7]
        ], dtype=np.float64)

In [3]:

A = np.array([
        [1e-20, 1],
        [1, 1],
        ])

In [4]:

n = A.shape[0]

elim_matrices = []

W = A

for i in xrange(n-1):
    M = np.eye(n)
    M[i+1:, i] = -W[i+1:, i]/W[i,i]
    
    elim_matrices.append(M)
    
    W = np.dot(M, W)
    print "--------------------------------------"
    print "column %d" % i
    print "--------------------------------------"
    print "elimination mat:"
    print M
    print "left over:"
    print W

--------------------------------------
column 0
--------------------------------------
elimination mat:
[[  1.00000000e+00   0.00000000e+00]
 [ -1.00000000e+20   1.00000000e+00]]
left over:
[[  1.00000000e-20   1.00000000e+00]
 [  0.00000000e+00  -1.00000000e+20]]

In [5]:

L = reduce(np.dot, [la.inv(M) for M in elim_matrices])
U = W

print "L"
print L
print "U"
print U

L
[[  1.00000000e+00  -0.00000000e+00]
 [  1.00000000e+20   1.00000000e+00]]
U
[[  1.00000000e-20   1.00000000e+00]
 [  0.00000000e+00  -1.00000000e+20]]

In [6]:

A_rebuilt = np.dot(L, U)
A_rebuilt

Out[6]:

array([[  1.00000000e-20,   1.00000000e+00],
       [  1.00000000e+00,   0.00000000e+00]])

In [7]:

la.norm(A_rebuilt - A)

Out[7]:

1.0