matlab's rref function modified to operate in gf(2)

g2rref.m

% This is a modified version of matlab's building rref which calculates
% row-reduced echelon form in gf(2).  Useful for linear codes.
% Tolerance was removed because yolo, and because all values
% should only be 0 or 1.  @benathon

function [A] = g2rref(A)
%G2RREF   Reduced row echelon form in gf(2).
%   R = RREF(A) produces the reduced row echelon form of A in gf(2).
%
%   Class support for input A:
%      float: with values 0 or 1
%   Copyright 1984-2005 The MathWorks, Inc. 
%   $Revision: 5.9.4.3 $  $Date: 2006/01/18 21:58:54 $

[m,n] = size(A);

% Loop over the entire matrix.
i = 1;
j = 1;

while (i <= m) && (j <= n)
   % Find value and index of largest element in the remainder of column j.
   k = find(A(i:m,j),1) + i - 1;

   % Swap i-th and k-th rows.
   A([i k],j:n) = A([k i],j:n);
   
   % Save the right hand side of the pivot row
   aijn = A(i,j:n);
   
   % Column we're looking at
   col = A(1:m,j);
   
   % Never Xor the pivot row against itself
   col(i) = 0;
   
   % This builds an matrix of bits to flip
   flip = col*aijn;
   
   % Xor the right hand side of the pivot row with all the other rows
   A(1:m,j:n) = xor( A(1:m,j:n), flip );

   i = i + 1;
   j = j + 1;
end

Thank you very much for this function!

I slightly modified the function to return the matrix M of row operations corresponding to reducing A to the reduced row echelon form. I also return a matrix N such that if it is applied to the right of the rref matrix A, then the first m columns have the form [ I_k, 0; 0, 0]. This is useful in some cases for matrix decomposition. I use it for decomposing a binary symplectic matrix into a product of elementary symplectic transformations.

Please update your function from my fork if you like these additions.

It's very useful!

Hi and thanks for the code.
I think I found a bug in the code, for example for the given matrix:
1 1 1 1 0 0 1 0
0 1 0 0 1 1 0 0
1 0 0 1 0 1 1 1
0 1 0 1 0 0 1 0
0 1 1 0 0 0 1 0
0 1 0 1 0 1 1 0
At the 5-th iteration, the following line fails (no such index is found)
% Find value and index of largest element in the remainder of column j.
k = find(A(i:m,j),1) + i - 1;

If we are interested in linear codes only, I think column swapping can solve this. Your thoughts?

Thanks!

Thank you very much,
I have been able to port this to numba python link to gist
Literally i have been wrestling with big binary matrices inversion for almost two days, you saved my day !!!
I had a pure python implementation which was not compatible with numba , this can be compiled with numba and now is blazing fast.

Hi,

Thanks for the code. I happened to try this code and found a bug of it. Like the original rref.m in MATLAB, I think in the beginning we need a condition "if isempty(find(A(i:m,j),1)))'' to make some situations work. Moreover, it would be great if the pivots can be one of the outputs of the program (if we add "if.., else..." part there the pivots can be saved). Thanks.

Hi,
Thank you for providing this efficient well working code. I was dealing with this problem for a while and your code helped me a lot. However, I need the result to be in the form I have sent you the picture for and the result that I get using your code is not like that. can you help me with that? thanks