# Constructing the basis vectors for the null space of a matrix

The algorithm to find a basis for the nullspace of a matrix $$A$$ is relatively simple for a human to do by hand on small matrices.

You find the RREF and then each row will have a leading 1 or be 0. If it has a leading one you express the constrained variable as $$x_i = \sum a_j s_j$$ where $$s_j$$ are your free variables.

You do this for each row and then you can construct a vector with as many rows as there are columns in $$A$$ expressing each equation, similar to this.

That's fine if you can do symbolic manipulation, but I need to get this basis numerically, i.e. I am telling a computer to do this and I cannot use a symbolic analyzer like sympy.

If I already have the RREF, how do I go about extracting the basis vectorss for the nullspace without invoking symbolic analysis?

Let's look at that example closely: the row-reduced matrix is $$\begin{pmatrix} 1 & -1 & 0 & 2 \\ 0 & 0 & 1 & -1 \end{pmatrix}$$ right? More generally, after you row-reduce, you'll need to get rid of any rows consisting of all zeroes at the bottom of your row-reduced matrix for the rest of this to work.

You can look at each row and find its leading "1", with pseudocode something like this:

for i = 1 to nrows
location(i) = 0  // says which column contains the "1" in the ith row
for j = 1 to ncols
if R(i, j) == 1
location(i) = j
break out of j-loop


When you're done, the array "location" will contain nrows column-indexes in increasing order, the indexes of the "non free" variables. In this example, we'll end up with location(1) = 1, location(2) = 3.

Now let's produce a list of "free"-variable indexes:

col = 1
counter = 1 // which entry in "location" are we looking at
for i = 1 to ncols-nrows
if location(counter) == col
counter = counter + 1
col = col + 1
else
free(i) = col
col = col + 1


When we're done with this, we'll have free(1) = 2, free(2) = 4.

Now comes the slightly tricky part. Let's suppose we set the free variables, $$x_2$$ and $$x_4$$ to $$x_2 = 1, x_4 = 0$$. (More generally: we'd set one of the free variables to $$1$$ and the others all to zero, but let's look at this concrete case as an example.) When we solve for $$x_1$$ and $$x_3$$, we get $$x_1 = 1$$ and $$x_3 = 0$$. Notice that these are in fact the entries in column $$2$$ of the row-reduced matrix, but negated! (Why column 2? Because we set $$x_2$$ to be 1.)

If we set $$x_2 = 0$$ and $$x_4 = 1$$, we get that $$x_1$$ and $$x_3$$ are $$-2$$ and $$+1$$, i.e., the entries in column $$4$$ of the row-reduced matrix, negated. There's a pattern here!

We want to produce $$ncols - nrows$$ basis vectors (where nrows is the number of rows in the reduced matrix after removing the all-zero rows!). We'll put each of these in a column of a matrix $$B$$. So $$B(1,1) ... B(1, ncol)$$ will be our first basis vector, and $$B(2, 1), ..., B(2, ncol)$$ will be our second, and so on. Here's some pseudocode:

B = matrix of zeros, with ncols rows, and ncols-nrows columns, one
for each basis vector
nBasis = ncols - nrows

for b = 1 to nBasis
// first fill in the zeroes-and-ones for the non-free var spots
B(free(b), b) = 1 // other entries are all zeroes already!

for n = 1 to nrows // number of entries in "location" matrix
B(location(n), b) = -R(n, free(b))



And that's pretty much it!

Here's that code in Matlab, which uses 1-based indexing (like most math books), but tends to do things with a matrix-at-a-time approach rather than using loops. Note that matlab has a built-in function for computing the row-reduces echelon form.

function B = nullBasis(A)
R = rref(A);
nrows = size(R, 1);
ncols = size(R, 2);

location = zeros(nrows, 1);
% First pseudocode fragment
% this loop is not very idiomatic matlab, alas.
for i = 1:nrows
location(i) = find(R(i, :), 1, 'first');
end
% second pseudocode fragment
free = 1:ncols;  % write all column indices in a sequence
free(location) = []; % and delete from the sequence all those that are non-free indices

% Third pseudocode fragment
nBasis = ncols - nrows;
B = zeros(ncols, nBasis);
B(free, :) = eye(nBasis);
B(location, :) = -R(:, free);


And here's an instance of it in operation:

>> A = [1 -1 -1 3; 2 -2 0 4]

A =

1    -1    -1     3
2    -2     0     4

>> nullBasis(A)

ans =

1    -2
1     0
0     1
0     1


By the way, most of this (in Matlab) is pointless, because not only is there a builtin rref function, there's also a builtin null, which returns an orthonormal basis for the nullspace of a matrix. More important, though, is that applying this code to large matrices is likely to lead to wrong results because of roundoff errors. It's probably fine for most 10 x 10 matrices you'll encounter at random, but bad for 1000 x 1000 matrices.

• I think Now let's produce a list of "free"-variable indexes: May have a bug, the variable col is never updated. Commented Aug 17, 2023 at 21:38
• Thanks...it's why I hate writing pseudocode. I believe I've fixed it, but who can be sure without testing? :) Commented Aug 18, 2023 at 21:14
• beware of this code I have only proven it correct, not run it. Commented Aug 18, 2023 at 21:47

John Hughes answer works but has some limitations, for example it only works for matrices where the rows are larger than the columns.

I tried generalizing the algorithm to work regardless of the matrix dimension, I made this rust code:

/// Compute a basis for the null space of a matrix.
pub fn null_basis<T: RealField, D: Dim, E: Dim>(
matrix: &OMatrix<T, D, E>,
) -> Matrix<T, na::Dyn, na::Dyn, na::VecStorage<T, na::Dyn, na::Dyn>>
where
DefaultAllocator: Allocator<T, D> + Allocator<T, E>,
DefaultAllocator: Allocator<T, D, E> + Allocator<T, D, E>,
DefaultAllocator: Allocator<T, na::base::dimension::Const<1>, D> + Allocator<T, D, E>,
DefaultAllocator: Allocator<T, na::base::dimension::Const<1>, E> + Allocator<T, D, E>,
DefaultAllocator: Allocator<T, <D as DimMin<E>>::Output, E>,
DefaultAllocator: Allocator<T, <D as DimMin<E>>::Output>,
DefaultAllocator: Allocator<T, D, <D as DimMin<E>>::Output>,
DefaultAllocator:
Allocator<T, <<D as DimMin<E>>::Output as DimSub<Const<1>>>::Output>,
T: Scalar + ClosedMul + Copy,
D: DimMin<E>,
<D as DimMin<E>>::Output: DimSub<Const<1>>,
{
let reduced = rref(matrix);

let epsilon = na::convert(f64::EPSILON * 10.0);
let rank = reduced.rank(epsilon);
let cols = reduced.ncols();
let rows = reduced.nrows();

// Find all leading 1's in the rref.
let mut row = 0;
let mut col = 0;
while row < rows && col < cols
{
if (reduced[(row, col)] - convert(1.0)).abs() < epsilon
{
row += 1;
}

col += 1;
}

// Find the columns corresponding to the free variables.
let mut free_vars = Vec::new();
let mut end = 0;
let mut counter = 0;
#[rustfmt::skip] {
while end < cols
{

end =
else { cols };

for c in leading_ones[counter] + 1..end
{
free_vars.push(c);
}

counter += 1;
}};

// Construct basis vectors initialized to 0.
let zero: T = convert(0.0);
let mut basis = DMatrix::from_element(cols, cols - rank, zero);

// Compute the entries in each basis vector.
for b in 0..cols - rank
{
basis[(free_vars[b], b)] = na::convert(1.0);

{
$$$$

• You're slightly mistaken in asserting that my solution has a limitation: ncols and nrows` are the number of columns and rows in the row-reduced matrix with all completely-zero rows removed. In this situation, the number of rows will never be greater than the number of columns, because row-rank and column-rank are the same. The ORIGINAL matrix may have many additional (linearly dependent) rows as well, of course. Commented Aug 18, 2023 at 21:12