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