# LU and upper echelon form

In [5]:
import numpy as np
import scipy.linalg as la


## Part I: The Problem

Here's a matrix:

In [6]:
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 [12]:
P, L, U = la.lu(A)


Check that it is actually a factorization:

In [13]:
la.norm(A-P.dot(L).dot(U))

Out[13]:
0.0


Now look at U:

In [10]:
U

Out[10]:
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 [19]:
def swap_rows(mat, i, j):
temp = mat[i].copy()
mat[i] = mat[j]
mat[j] = temp

In [36]:
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 [29]:
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 [25]:
U

Out[25]:
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 [27]:
M

Out[27]:
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 [34]:
A - la.inv(M).dot(U)

Out[34]:
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 []: