5
$\begingroup$

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.

Improve this question
$\endgroup$

2 Answers 2

Reset to default
6
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ The fundamental point here is that you don't want to actually compute the Moore-Penrose pseudoinverse, since that matrix will almost certainly be fully dense. Rather, you should consider using an iterative method to compute a minimum norm least squares solution. You could do this by solving a system involving $(A+\frac{ee^{T}}{n})$, or you might want to consider a more straight forward approach. For example, LSQR can be used with a damping factor to minimize a weighted sum of $\| Ax - -b\|^2$ and $\| x \|^2$. $\endgroup$
    – Brian Borchers
    Commented Aug 11, 2011 at 4:19
  • 2
    $\begingroup$ indeed - I coded Federico's suggestion and it works well :) I used conjugate gradient since A+ee'/n is SPD (diagonal dominant with positive diagonal). $\endgroup$
    – WhitAngl
    Commented Aug 11, 2011 at 11:07
3
$\begingroup$

Are you looking for a fast practical method, or a fast theoretical method? If the former, there are very fast solvers based on sparse Cholesky or sparse LDL decomposition (both of which can be reused for many $b$). You should check out Tim Davis' beautiful book called something like "Sparse direct solvers".

There are very few theoretical bounds, since the running time depends heavily on the sparsity pattern.

Improve this answer
$\endgroup$
3
  • 1
    $\begingroup$ To my knowledge, Cholesky doesn't work for rank deficient matrices (otherwise, there is a 0 pivot at some point). $\endgroup$
    – WhitAngl
    Commented Aug 10, 2011 at 21:57
  • 2
    $\begingroup$ Yes, but LDL does, I am pretty sure. $\endgroup$
    – Igor Rivin
    Commented Aug 11, 2011 at 0:06
  • 1
    $\begingroup$ It does, but only after you build-in symmetric pivoting. Even then, it's easy to choke $LDL^\top$; try $\begin{pmatrix}0&1\\1&0\end{pmatrix}$. $\endgroup$
    – J. M. isn't a mathematician
    Commented Dec 9, 2011 at 4:22

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .