# coding: utf-8 # # 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)) # Now look at `U`: # In[10]: U # * 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)) # Let's see if \$U\$ is actually in echelon form: # In[25]: U # And what does \$M\$ look like? # In[27]: M # ... not much structure here. # -------------- # But we can still have *something* a little like LU: # In[34]: A - la.inv(M).dot(U) # In[ ]: