bitgauss is a library for doing linear algebra over the 2-element finite field. It consists of a Rust crate with Python bindings. To use the Rust crate, add bitgauss as a dependency to Cargo.toml, or to use the Python bindings, run:

pip install bitgauss

It provides a BitMatrix class, which is a bitpacked 2D matrix that implements fast linear algebraic operations using fast, vectorized bitwise operations. The main features are:

• getting and setting individual matrix elements (as bools)
• fast row operations and dot product using bitwise operations
• fast in-place and out-of-place matrix transpose using a recursive block method
• horizontal and vertical concatenation of matrices
• matrix multiplication
• Gaussian elimination and related methods (e.g. rank, inverse, and nullspace)

The goal of this library is to keep things small and simple, but if you need some core functionality that is missing, feel free to send a pull request!

Here is a simple example of using the BitMatrix class from Python:

from bitgauss import BitMatrix

# Construct a 300x400 matrix whose entries are given by the bool-valued function
m1 = BitMatrix.build(300, 400, lambda i, j: (i + j) % 2 == 0)

# Construct a random 80x300 matrix with an optional random seed
m2 = BitMatrix.random(80, 300, seed=1)

# Construct a random invertible 300x300 matrix
m3 = BitMatrix.random_invertible(300, seed=1)

m4 = m2 * m3           # Matrix multiplication
m3_inv = m3.inverse()  # Returns the inverse
m1_t = m1.transposed() # Returns transpose
m1.transpose_inplace() # Transpose inplace (padding if necessary)
m1.gauss()             # Transform to row-echelon form
m1.gauss(full=True)    # Transform to reduced row-echelon form
ns = m1.nullspace()    # Returns a spanning set for the nullspace

...or from Rust:

use bitmatrix::BitMatrix;
use std::ops::Mul;
use rand::{Rng, SeedableRng, rngs::SmallRng};

let mut rng = SmallRng::seed_from_u64(1);

// Construct a 300x400 matrix whose entries are given by the bool-valued function
let mut m1 = BitMatrix::build(300, 400, |i, j| (i + j) % 2 == 0);

// Construct a random 80x300 matrix using the given random number generator
let m2 = BitMatrix::random(&mut rng, 80, 300);

// Construct a random invertible 300x300 matrix
let m3 = BitMatrix::random_invertible(&mut rng, 300);

let m4 = &m2 * &m3;         // Matrix multiplication
let m3_inv = m3.inverse();  // Returns the inverse
let m1_t = m1.transposed(); // Returns transpose
m1.transpose_inplace();     // Transpose inplace (padding if necessary)
m1.gauss(false);            // Transform to row-echelon form
m1.gauss(true);             // Transform to reduced row-echelon form
let ns = m1.nullspace();    // Returns a spanning set for the nullspace

For more info, have a look at the Python documentation, the Rust documentation, or the Jupyter notebook getting_started.ipynb.

Graphic form

This library also contains methods for converting matrices to graphic form using the Bixby-Wagner algorithm. For a given matrix, this will try to find another matrix with the same rowspace whose columns have Hamming weight at most 2. This is a useful transformation e.g. for error correcting codes. Here's an example:

use bitmatrix::BitMatrix;

let m = BitMatrix::from_int_vec(&vec![
    vec![1, 1, 0, 1, 1],
    vec![0, 1, 1, 1, 0],
    vec![1, 1, 1, 0, 0],
]);

let n = m.graphic_form().unwrap();
println!("graphic form:\n{}", n);
// graphic form:
//  1  0  1  0  1 
//  0  1  0  0  1 
//  0  1  1  1  0

This method returns None if the matrix cannot be put into graphic form. Alternatively, the method graphic_form_partial will try to put the matrix into graphic form and return a pair consisting of a partial solution and a list of hyperedge columns, i.e. those whose Hamming weight is still above 2.

use bitmatrix::BitMatrix;

let m = BitMatrix::from_int_vec(&vec![
    vec![1, 1, 0, 1, 1, 0, 0],
    vec![0, 1, 1, 1, 0, 1, 0],
    vec![1, 1, 1, 0, 0, 0, 1],
]);

let (n, hyper) = m.graphic_form_partial();
println!("partial graphic form:\n{}", n);
println!("hyperedge columns: {:?}", hyper);
// partial graphic form:
//  1  1  0  1  1  0  0 
//  1  0  1  0  1  1  0 
//  0  1  0  0  1  1  1 
// 
// hyperedge columns: [4]

Performance

The BitMatrix class is substantially faster than a pure Python implementation of a binary matrix without bitpacking. For example, in getting_started.ipynb there is a calculation that is more than 100X faster using bitgauss:

m = BitMatrix.random(1000, 2000, seed=42) # bitgauss
m1 = Mat2(m.to_int_list()) # pure python

%timeit m.copy().gauss()
# 82.2 ms ± 1.7 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

%timeit m1.copy().gauss()
# 10.1 s ± 172 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

The difference is probably less pronounced with numpy or other fast compiled linear algebra backends. If you'd like to contribute some benchmarks vs. those, feel free.

A note on SIMD

bitgauss is designed to take advantage of the fact that modern hardware and compilers are able to process many bit-entries at once using SIMD. While this package doesn't explicitly use SIMD (via e.g. the experimental std::simd module in nightly Rust), it relies on the fact that LLVM, which the Rust compiler uses under the hood, is already very good at using SIMD automatically when needed and available. From my experiments and benchmarking, this seems to almost always do better than introducing SIMD instructions by hand. Please keep this in mind before you spend a lot of time writing a PR that adds explicit SIMD instructions. Basically, don't do what I did. :)