[SciPy-user] Least squares for sparse matrices?

[Andy Fraser]

Thank you for the quick response.  I spent some time thinking and
looking things up, and I think I understand your suggestions (with the
exception of preconditioning.)  I can afford to do a sparse LU.  But I
think that your suggestion does not provide a control equivalent to
the rcond argument to numpy.linalg.lstsq.

>>>>> "NB" == Nathan Bell <wnbell at gmail.com> writes:

    NB> 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.
    >> 

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

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

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

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

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


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

-- 
Andy Fraser				ISR-2	(MS:B244)
afraser at lanl.gov			Los Alamos National Laboratory
505 665 9448				Los Alamos, NM 87545