Just noting that this here is quite inefficient: It would be more efficient to compute the row echelon form for as many columns as we need and extracting the information on the linearly independent columns from the row echelon form.

Keeping track of the row transformations would allow to add additional columns as needed.

I do not know whether such a variant of row echelon form is available in SageMath.

Furthermore, in verify_linear_combination, we could find out what column we need to add to get the next invertible submatrix.

Perhaps add some comment into the code, we do not want to delay this PR by rewriting everything.

Originally posted by @cheuberg in #35682 (comment)