I am attempting to invert (mod 2) a non-sparse binary matrix. I am trying to do it with very very large matricies, like 100_000 by 100_000. I have tried libraries, such as sympy, numpy, numba. Most of these do this with floats and don't apply mod 2 at each step. In general, numpy reports the determinant of a random binary invertible array of size 1000 by 1000 as INF. Here is my attempted code, but this lags at 10_000 by 10_000. This function is also nice, since even if I am not invertible, I can see the rank. def binary_inv(A): n = len(A) A = np.hstack((A.astype(np.bool_),np.eye(n,dtype=np.bool_))) for i in range(n): j = np.argmax(A[i:,i]) + i if i!=j: A[[i,j],i:] = A[[j,i],i:] A[i,i] = False A[A[:,i],i+1:] ^= A[i:i+1,i+1:] return A[:,n:] Any thoughts on making this as fast as possible? Continue reading...