[Math] Solving for Moore Penrose pseudo inverse

data analysismatricesna.numerical-analysisst.statistics

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.

Best Answer

For black-box linear algebra (GMRES and the like) you don't need "sparse", you only need "can compute products quickly". If you check the docs for your sparse solver, I'm sure there'll be a version where you can provide directly the function $v\mapsto Av$ (and sometimes $w\mapsto A^Tw$ is needed, too). In your case, you can compute $(A+\frac{ee^T}{n})v$ quickly by exploting sparsity.

And yes, if you wish to use iterative solvers you'll have to solve the system more or less from scratch for every new right-hand side linearly independent from the previous ones. Some savings are possible, but getting them is still a hot research topic.