超大规模二进制矩阵模2求逆的性能优化咨询

我现在需要处理超大规模的二进制矩阵模2求逆的问题，矩阵尺寸甚至要达到100,000×100,000这么大。我试过sympy、numpy、numba这些库，但它们大多是基于浮点数运算的，不会在每一步都应用模2操作——比如numpy连一个1000×1000的随机可逆二进制矩阵，计算行列式都会返回INF。

我自己写了一个函数，逻辑上是没问题的：它不仅能求逆，还能在矩阵不可逆的时候得到它的秩，这点很实用。但问题是，当矩阵到10,000×10,000这个规模时，速度就慢得不行了。

我的代码如下：

import numpy as np

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:]

想请教各位大佬，有没有什么办法能把这个函数的速度提上去，尽可能快地处理超大矩阵？

备注：内容来源于stack exchange，提问作者Bobby Ocean