import numpy as np
import scipy.linalg as la
Here's a matrix:
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.]])
P, L, U = la.lu(A)
Check that it is actually a factorization:
la.norm(A-P.dot(L).dot(U))
Now look at U
:
U
U
were upper echelon?def swap_rows(mat, i, j):
temp = mat[i].copy()
mat[i] = mat[j]
mat[j] = temp
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$:
M, U = m_echelon(A)
diff = M.dot(A)-U
print(la.norm(diff))
Let's see if $U$ is actually in echelon form:
U
And what does $M$ look like?
M
... not much structure here.
But we can still have something a little like LU:
A - la.inv(M).dot(U)