If you have one (or a small number of) right-hand sides, I would try to make LSQR work. It can find a minimum norm solution even if A is rank-deficient, and you can use preconditioning.
Unfortunately, in my case, the generalized inverse of AA' is the preconditioner for the system, which is why I need the factorization of A'. Essentially, I take this factorization and then run it through my own iterative method. When I run tests in MATLAB, SPQR scales fine for matrices of at least a few hundred thousand rows and columns. For larger, it would be nice to essentially have an incomplete Q-less QR factorization, which I don't think exists, but should be an extension of the incomplete Choleski work.
But, yes, LSQR or more fitting LSMR solves a similar problem, but they're the iterative solver and I need the preconditioner, which I'm using the factorization for.
Also, if your problem is a good fit for a method like this, it could be impetus to add it to PETSc. https://epubs.siam.org/doi/pdf/10.1137/120866580