[SciPy-user] Least squares for sparse matrices?
Nathan Bell
wnbell@gmail....
Fri Nov 21 14:17:37 CST 2008
On Fri, Nov 21, 2008 at 1:14 PM, Andy Fraser <afraser@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@gmail.com
http://graphics.cs.uiuc.edu/~wnbell/
More information about the SciPy-user
mailing list