Computing the Nullspace

In [41]:
import numpy as np
import numpy.linalg as la
In [46]:
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]
In [47]:
from scipy.linalg import lu
In [48]:
PT, L, U = lu(A.T)
P = PT.T
In [49]:
la.norm(
    P.dot(A.T) - L.dot(U))
Out[49]:
1.1322097734007351e-15
In [50]:
np.set_printoptions(precision=3)
U
Out[50]:
array([[  1.027e+00,   6.761e-01,  -2.323e-01,   1.202e+00,  -1.347e-01],
       [  0.000e+00,  -3.498e+00,  -1.468e+00,   1.749e+00,  -7.342e+00],
       [  0.000e+00,   0.000e+00,  -2.140e+00,   4.441e-16,  -1.070e+01],
       [  0.000e+00,   0.000e+00,   0.000e+00,   5.761e-16,   0.000e+00],
       [  0.000e+00,   0.000e+00,   0.000e+00,   0.000e+00,   8.882e-16]])

Now define NUT as vectors spanning the nullspace of $N(U^T)$.

In [55]:
NUT = np.eye(n)[:, 3:]

Check that it's actually a nullspace:

In [56]:
U.T.dot(NUT)
Out[56]:
array([[  0.000e+00,   0.000e+00],
       [  0.000e+00,   0.000e+00],
       [  0.000e+00,   0.000e+00],
       [  5.761e-16,   0.000e+00],
       [  0.000e+00,   8.882e-16]])

Now define NA as some vectors spanning $N(A)$:

In [58]:
NA = P.T.dot(la.solve(L.T, NUT))

And check:

In [59]:
A.dot(NA)
Out[59]:
array([[ -1.110e-16,   0.000e+00],
       [  0.000e+00,   0.000e+00],
       [  0.000e+00,   1.110e-16],
       [  1.110e-16,   5.551e-17],
       [  0.000e+00,   4.441e-16]])