I have a system to solve, set up as : $$Ax = b$$ with a square rank deficient matrix $A$. The paper suggests to use a Moore Penrose pseudo inverse, which in my case can be computed using the traditional inverse : $$ A^+ = (A+\frac{ee^T}{n})^{-1} - \frac{ee^T}{n} $$ where $e$ is a vector containing only ones, and $n$ is the dimension of the matrix. This matrix comes from the solution of a Multidimensional Scaling problem using the SMACOF method (the Guttman transform).

However, in my case, my matrix $A$ is very sparse (and rank deficient) : what method can I use to efficiently solve the original system without making my matrix dense (as would be the case by using an SVD, by using the above formula for the pseudo inverse, by computing $A^TA$ or by QR factorization) ? $A$ is also symmetric, has a positive diagonal, and the other values are either -1 or 0, and such as the sum of each row (resp. column) is 0.

Preferably, since I'll need to solve for multiple right hand sides with this same matrix, I would like to avoid performing the resolution from scratch for each right hand side. I would also like to get exactly the same result as the one obtained with the Moore Penrose pseudoinverse.

Thanks.