Check if binary equations are linearly independent in MATLAB

Question

I have a matrix A, with the coefficients of a binary equation system.

In order to solve it, I need to know if all of those equations are linearly independent. However, rank(A) doesn't work here, because MATLAB doesn't know they're binary equations.

Here's an example:

>> A
A =

   1   0   1   1
   0   1   1   1
   1   1   0   0

>> rank(A)
ans =  3

It's important to mention that it's a $Ax = 0 \pmod{2}$ type system.

As you can see, if you add rows 2 and 3, the result is the first one. This means the three of them aren't linearly independent, while rank(A) is saying they are.

How can I solve this?

Answer 1

rank(A) of Matlab works under the assuption that the matrix has elements in the field of complex numbers, so this won't work for you. But, the rank detection of a binary system is really easy (only algebra no numerical noise considerations). Just apply Gaussian algorithm with rank detection (i.e., in general with row and column permutations).

In your example the first non-trivial elimination step is adding the first equation to the last. In this way you get: \begin{align*} A_2=\pmatrix{1&0&1&1\\0&1&1&1\\0&1&1&1} \end{align*} The next elimination step consists of adding the second row to the just created modified third row. In this way you get just a zero row: \begin{align*} A_2=\pmatrix{1&0&1&1\\0&1&1&1\\0&0&0&0} \end{align*} You 've got two non-zero rows of the staircase system. Therefore, $\operatorname{rank}(A)=2$.

By applying the transformation steps also to the right-hand side you can use the algorithm to solve your system if it is solvable.

There is code for Gauss' algorithm at Mathworks. For completing the answer there follows a code for Octave. Maybe it runs in Matlab too:

## Usage: [x,r,LU,pr,pc,res] = gfGauss(A,y,p,pr,pc)
## Solve system Ax=y in the Galois Field GF(p).
## @param[in] A system matrix
## @param[in] y right-hand side
## @param[in] p size of field (must be prim)
## @param[in] pr optional row permutation; when specified A is already decomposed
## @param[in] pc optional column permutation
## @param[in] r optional rank; only required when A is already decomposed
## @return x solution
## @return r rank of the system
## @return LU LU-factorization of A
## @return pr row permutation vector
## @return pc column permutation vector
## @return res residuum The system cannot be solved if this is non-zero.
function [x,r,LU,pr,pc,res] = gfGauss(A,y,p,pr,pc,r)
  if ~exist("y","var")
    y = [];
  end
  if ~exist("p","var")
    p = 2;
  end
  n = min(size(A));
  nc = size(A,2);
  ny = size(y,2);
  if exist("pr","var") ## System already decomposed
    ## Solve L*z = y (z stored in y)
    LU = A;
    if ~exist("r","var")
      r = n;
    end
    x = zeros(nc,ny);
    y = y(pr,:);
    for iDiag=1:r
      y(iDiag,:) = y(iDiag,:) - LU(iDiag,1:(iDiag-1))*y(1:(iDiag-1),:);
    end
    x = [y(1:r,:);zeros(nc-r,ny)];
  else ## decomposition required
    ## Decompose L*U = A with permuation
    ## and simultaneously solve L*z = y (z stored in y)
    LU = mod([A,y],p);
    r = n;
    pr = 1:size(A,1);
    pc = 1:size(A,2);
    for iDiag=1:n
      [ir,ic] = find(LU(iDiag:end,iDiag:nc)); # Efficiency: search could stop at first non-zero element
      if length(ir)==0
    r = iDiag-1;
    break;
      else
    ir = ir(1)+iDiag-1;
    ic = ic(1)+iDiag-1;
      end;
      ## row/column permutation
      pr([ir,iDiag]) = pr([iDiag,ir]);
      LU([ir,iDiag],:) = LU([iDiag,ir],:);
      pc([ic,iDiag]) = pc([iDiag,ic]);
      LU(:,[ic,iDiag]) = LU(:,[iDiag,ic]);
      ## elimination step
      pivInv = mod(LU(iDiag,iDiag)^(p-2),p); ## multiplicative inverse of the pivot in GF(p); only works for pime p
      LU((iDiag+1):end,iDiag) = mod(pivInv.*LU((iDiag+1):end,iDiag),p);
      LU((iDiag+1):end,(iDiag+1):end) = mod(LU((iDiag+1):end,(iDiag+1):end)-LU((iDiag+1):end,iDiag)*LU(iDiag,(iDiag+1):end),p);
    end
    x = [LU(1:r,(nc+1):end);zeros(nc-r,ny)];
  end
  if ny > 0
    ## back-substitution
    for iDiag=r:-1:1
      pivInv = mod(LU(iDiag,iDiag)^(p-2),p);
      x(iDiag) = mod(pivInv.*(x(iDiag)-LU(iDiag,(iDiag+1):r)*x((iDiag+1):r)),p);
    end
    x(pc,:) = x;
  else
    x = [];
  end
  LU = LU(:,1:nc);
  if isargout(6)
    res = mod(A*x-y,p);
  end
endfunction

The following code segment shows the application of the algorithm to the problem in the question for deducing the rank of $A$:

A = [1   0   1   1
     0   1   1   1
     1   1   0   0];
[x,r,LU,pr,pc]=gfGauss(A)

If the size p of the Galois field is not included in the argument list it defaults to GF(2). That means the computation is based on dual numbers. If the right-hand side is omitted an empty vector x=[] is returned. Nevertheless, the LU-decomposition and the rank-detection are carried out.

The output of the program is:

>>     [x,r,LU,pr,pc]=gfGauss(A)
x = [](0x0)
r =  2
LU =

   1   0   1   1
   0   1   1   1
   1   1   0   0

pr =

   1   2   3

pc =

   1   2   3   4

The rank information r=2 that is included in the output is most interesting in this example with respect to the question.

Answer 2

use rank(gf(A)).

see MathWorks support page: Difference between rank(S) and rank(gf(S)) when S is a matrix