In _remove_redundancy_sparse, spsolve is repeatedly called within a loop.
pi = scipy.sparse.linalg.spsolve(B.transpose(), e).reshape(-1, 1)
Since B does not change, it would be a lot faster to have an LU factorization outside of the loop