[SciPy-user] Least squares for sparse matrices?

[Nathan Bell]

On Fri, Nov 21, 2008 at 1:14 PM, Andy Fraser <afraser at lanl.gov> wrote:
> I am trying to understand the behavior of matrices that approximate
> Radon transforms and their inverses.  So far, I've used commands like
> the following to experiment with various values of SVDcond:
>
> V_1,resids,rank,s = numpy.linalg.lstsq(M_Radon,U_0,rcond=SVDcond)
>
> Is there a sparse version of lstsq that will let me exploit the
> sparseness of M_Radon?
>
> It would be even better if there were something like a sparse SVD that
> would let me pre-calculate an SVD decomposition and then later use it to
> invert several U_0 vectors.
>

I think the best we can currently offer you is an iterative solver
(e.g. cg()) on the normal equations.  Something like (untested)

from scipy.sparse.linalg import LinearOperator, cg
from scipy.sparse import csr_matrix

M_Radon = csr_matrix(M_Radon)
M,N = M_Radon.shape

def matvec(x):
    return M_Radon.T * (M_Radon * x)

A = LinearOperator( (N,N), matvec)
x,info = cg(A, U_0)


If you can afford to do a sparse LU decomposition of M_Radon, then you
can either use that directly, or as a preconditioner to cg.


-- 
Nathan Bell wnbell at gmail.com
http://graphics.cs.uiuc.edu/~wnbell/