I'll have a look at Catamari and thanks for the link. Maybe you'll have a better idea, but essentially I need a generalized inverse of AA' where A has more columns than rows (short and fat.) Often, A becomes underdetermined enough where AA' no longer has full-rank, but I need a generalized inverse nonetheless. If A' was full rank, then the R in the QR factorization of A' is upper triangular. If A' is not full rank, but we can permute the columns, so that the R in the QR factorization of A' has the form [RR S] where RR is upper triangular and S is rectangular, we can still find the generalized inverse. As far as I know, the permutation that ensures this form requires a rank-revealing QR factorization.
For dense matrices, I believe GEQP3 in LAPACK pivots so that the diagonal elements of R are decreasing, so we can just threshold and figure out when to cut things off. For sparse, the only code I've tried that's done this properly is SPQR with its rank-revealing features.
In truth, there may be a better way to do this, so I might as well ask: Is there a good way to find the generalized inverse of AA' where A is rank-deficient as well as short and fat?
As far as where they come from, it's related to finding minimum norm solutions to Ax=b even when A is rank-deficient. In my case, I know the solution exists for a given b, even though the solution may not exist in general.
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.
For dense matrices, I believe GEQP3 in LAPACK pivots so that the diagonal elements of R are decreasing, so we can just threshold and figure out when to cut things off. For sparse, the only code I've tried that's done this properly is SPQR with its rank-revealing features.
In truth, there may be a better way to do this, so I might as well ask: Is there a good way to find the generalized inverse of AA' where A is rank-deficient as well as short and fat?
As far as where they come from, it's related to finding minimum norm solutions to Ax=b even when A is rank-deficient. In my case, I know the solution exists for a given b, even though the solution may not exist in general.