Basis complement of image of sparse matrix

I have a map of vector space $$R:T\to M$$, and I'm interested in the quotient $$E=M/R(T)$$. In particular I want to construct a section $$s:E\to M$$. Theoretically this is not very difficult to solve, but I need a fast implementation. A simple theoretical solution is using Gram-Schmidt to build an orthogonal projection $$p:M\to R(T)$$ and taking $$s=(\mathrm{Id}_M-p)^\top$$, but this is computationally extremely inefficient.

Right now my approach is as follows. I take the matrix R and then augment by putting the identity matrix of M on the right, and then row reducing. $$\left(\begin{array}{c|c}R & \mathrm{Id}_M\end{array}\right) \to \left(\begin{array}{c|c}b^\top&*\\\hline0&s^\top\end{array}\right)$$ Here $$b$$ is a basis of the transpose $$R^\top(M)$$, and $$s$$ is the section I wanted.

In my particular case both $$T$$ and $$M$$ can have high dimension (on the order of 10,000 but perhaps much larger), but $$R$$ is very sparse (the number of non-zero entries is roughly twice the dimension of $$T$$). Doing row reduction with dense matrices of this size is very slow.

In python (scipy) there seems to be an efficient implementation of finding the LU decomposition of a square, non-singular, sparse matrix. This seems to be very useful, but doesn't help in my approach because the matrix is rectangular. One could add rows of zeroes to make it square, but then the resulting matrix is singular.

My question: is there a modification of my approach, or a different approach entirely, where I can make use of the sparsity of $$R$$ to speed up computation? If it would work with built-in scipy functions that would be a big bonus, but I don't mind implementing an algorithm.

A little more context, in case it is relevant. I have a differential complex $$d_i:M_i\to M_{i+1}$$, and each $$M_i$$ has a subspace $$T_i$$. The differential descends to the subspace $$E_i=M_i/T_i$$. I want to compute the Betti numbers $$\dim H^i(E)$$. If $$s_i:E_i\to M_i$$ is a section, then we can compute this if we know the dimensions of the kernels/cokernels (i.e. the rank) of $$\widetilde{d}_i =s_{i+1}^\top\circ d\circ s_i$$. Note that $$\widetilde d_i$$ is not the differential of $$E_i$$, but its (co)kernel is the same.

Finding a basis for the complement $$R(T)^\perp$$ is equivalent to finding a basis of the kernel $$\ker R^\top$$ (recall that $$\ker R^\top\oplus \mathrm{im}\,R=M$$). Hence the problem reduces to finding the kernel of a sparse rectangular matrix.
In fact, there it is not even necessary to compute a basis of $$\ker R^\top$$. Using the identity $$\mathrm{im}\,A^\top\oplus \ker A=V$$ for $$A:V\to W$$ repeatedly, together with the observation that $$\dim \mathrm{im}\,R_i^\top+\dim \mathrm{im}\,s^\top_i=\dim M_i$$, one can derive a formula for the rank of $$\widetilde d_i$$ in terms of $$\dim T_i$$, $$\dim M_i$$ and the ranks of $$R_{i+1}^\top d_i R_i^\top$$, $$R_{i+1}^\top d_i$$, $$d_i R_i^\top$$, and $$d_i$$. This still leaves the problem of computing ranks of large sparse rectangular matrices, which in fact doesn't seem to be much easier than computing the basis of the kernel.