import numpy as np
import numpy.linalg as la
n = 5
np.random.seed(25)
A = np.random.randn(n, n)
# Decrease the rank
A[4] = A[0] + 5 * A[2]
A[1] = 3 * A[0] -2 * A[3]
from scipy.linalg import lu
PT, L, U = lu(A.T)
P = PT.T
la.norm(
P.dot(A.T) - L.dot(U))
np.set_printoptions(precision=3)
U
Now define NUT
as vectors spanning the nullspace of $N(U^T)$.
NUT = np.eye(n)[:, 3:]
Check that it's actually a nullspace:
U.T.dot(NUT)
Now define NA
as some vectors spanning $N(A)$:
NA = P.T.dot(la.solve(L.T, NUT))
And check:
A.dot(NA)