LU and upper echelon form¶

In [1]:
#keep
import numpy as np
import scipy.linalg as la


Part I: The Problem¶

Here's a matrix:

In [2]:
#keep
A = np.array([[0.,  0., -1., -1.,  0.,  0.,  0., -1.,  0.],
[0.,  0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],
[0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[-1.,  0.,  1.,  0.,  0.,  0., -1.,  0.,  0.],
[1., -1.,  0.,  1., -1.,  1.,  0.,  0., -1.],
[0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
[0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
[0.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[0.,  0.,  0.,  0.,  1.,  0.,  0.,  0.,  0.]])

In [3]:
#keep
P, L, U = la.lu(A)


Check that it is actually a factorization:

In [4]:
#keep
la.norm(A-P.dot(L).dot(U))

Out[4]:
0.0

Now look at U:

In [5]:
U

Out[5]:
array([[-1.,  0.,  1.,  0.,  0.,  0., -1.,  0.,  0.],
[ 0., -1.,  1.,  1., -1.,  1., -1.,  0., -1.],
[ 0.,  0., -1., -1.,  0.,  0.,  0., -1.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
[ 0.,  0.,  0.,  0., -1.,  1., -1., -1., -1.],
[ 0.,  0.,  0.,  0.,  0., -1.,  0.,  1.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  1.,  0.,  0.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.],
[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.]])
• What do you notice about the last two rows?
• Would this be allowed if U were upper echelon?

Part II: Getting to Echelon Form¶

In [6]:
#keep
def swap_rows(mat, i, j):
temp = mat[i].copy()
mat[i] = mat[j]
mat[j] = temp

In [7]:
#keep
def m_echelon(A):
m, n = A.shape
M = np.eye(m)
U = A.copy()

row = 0
for col in range(min(m, n)):
piv_row = row + np.argmax(np.abs(U[row:, col]))

if abs(U[piv_row, col]) == 0:
# column is exactly zero
continue

swap_rows(U, row, piv_row)
swap_rows(M, row, piv_row)

for el_row in range(row + 1, m):
fac = -U[el_row, col]/U[row, col]
U[el_row] += fac*U[row]
M[el_row] += fac*M[row]

row += 1

return M, U


Compute M and U, and check that $MA=U$:

In [8]:
#keep
M, U = m_echelon(A)

diff = M.dot(A)-U

print(la.norm(diff))

0.0


Let's see if $U$ is actually in echelon form:

In [9]:
#keep
U

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

And what does $M$ look like?

In [10]:
#keep
M

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

... not much structure here.

But we can still have something a little like LU:

In [11]:
#keep
A - la.inv(M).dot(U)

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


In [12]:
A = np.random.rand(5,5)

In [13]:
M, U = m_echelon(A)

In [14]:
M

Out[14]:
array([[ 0.        ,  0.        ,  0.        ,  0.        ,  1.        ],
[ 1.        ,  0.        ,  0.        ,  0.        , -0.01735781],
[-0.94411529,  0.        ,  1.        ,  0.        , -0.08147027],
[-0.15695783,  0.        ,  0.30258661,  1.        , -0.59104586],
[-0.66098952,  1.        ,  0.18899847,  0.23650075, -0.90859199]])
In [16]:
U.round(3)

Out[16]:
array([[ 0.777,  0.194,  0.063,  0.318,  0.569],
[ 0.   ,  0.752,  0.825,  0.592,  0.306],
[-0.   ,  0.   , -0.73 , -0.034,  0.535],
[-0.   ,  0.   , -0.   ,  0.65 ,  0.839],
[-0.   ,  0.   , -0.   ,  0.   ,  0.155]])
In [17]:
L = np.linalg.inv(M)

In [18]:
L

Out[18]:
array([[ 0.01735781,  1.        ,  0.        ,  0.        ,  0.        ],
[ 0.76814609,  0.51299527, -0.11743651, -0.23650075,  1.        ],
[ 0.09785804,  0.94411529,  1.        ,  0.        ,  0.        ],
[ 0.56415977, -0.12871882, -0.30258661,  1.        ,  0.        ],
[ 1.        ,  0.        ,  0.        ,  0.        ,  0.        ]])
In [19]:
L.round(2)

Out[19]:
array([[ 0.02,  1.  ,  0.  ,  0.  ,  0.  ],
[ 0.77,  0.51, -0.12, -0.24,  1.  ],
[ 0.1 ,  0.94,  1.  ,  0.  ,  0.  ],
[ 0.56, -0.13, -0.3 ,  1.  ,  0.  ],
[ 1.  ,  0.  ,  0.  ,  0.  ,  0.  ]])
In [ ]: