0
$\begingroup$

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?

$\endgroup$

2 Answers 2

1
$\begingroup$

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.

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

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 leading_ones = Vec::with_capacity(cols);
    let mut row = 0;
    let mut col = 0;
    while row < rows && col < cols
    {
        if (reduced[(row, col)] - convert(1.0)).abs() < epsilon
        {
            leading_ones.push(col);
            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 =
            if counter < leading_ones.len() - 1 { leading_ones[counter + 1] }
            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);

        for n in 0..leading_ones.len()
        {
            basis[(leading_ones[n], b)] = -reduced[(n, free_vars[b])];
        }
    }

    basis
}
```
$\endgroup$
5
  • $\begingroup$ 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. $\endgroup$ Commented Aug 18, 2023 at 21:12
  • $\begingroup$ @JohnHughes I agree, but doing it that way requires more computations than if you can work direactly on a sorted RREF which may have additional rows which are all 0. Not being a mathematician but a CS I prefer the ugliness of keeping the 0 rows around. It;s less elegant but less work for the rref algorithm to do. $\endgroup$
    – Makogan
    Commented Aug 18, 2023 at 21:49
  • $\begingroup$ I used to be a mathematician and now am a prof. of CS. While computing the RREF, you can keep track of the number of nonzero rows (which will be in an upper block of the matrix) at essentially O(1) extra cost. Or you can mess around with them later. I suppose it's a matter of which you prefer. You might say "but they'll only sort of be zero-rows," but by the time you start worrying about roundoff error, you're far beyond what OP was asking, I expect. $\endgroup$ Commented Aug 19, 2023 at 3:05
  • $\begingroup$ @JohnHughes Note I am the op :p The reason I implemented it this way is just ease of implementation while making sure it worked for both kinds of square matrices because that's what I need for my given problem. It is still highly based on your answer I just tweaked what I found necessary to make it work in the computational context I am working on. $\endgroup$
    – Makogan
    Commented Aug 19, 2023 at 3:24
  • $\begingroup$ Oops--apologies for not noticing you're the OP. TBH, I figured that anyone who knew stuff about roundoff and rank "estimation" would have known how to do the stuff I described, although now that I think about it, I don't recall ever seeing it written down anywhere else. I just figured it out back when I taught linear algebra classes an assumed that "everyone else knew." :) I'm sure I'm not the first, of course. $\endgroup$ Commented Aug 19, 2023 at 12:18

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .