[SciPy-User] preconditioned conjugate gradient

[Pauli Virtanen]

Tue, 23 Feb 2010 11:14:53 -0800, Jake VanderPlas wrote:
> I'm looking for a method to solve a sparse linear equation A*x=b, where
> A is a NxN symmetric scipy.sparse.LinearOperator object, and b is a 1D
> numpy vector.  The obvious choice would be something like
> scipy.sparse.linalg.cg.  The problem is, the condition number of A is
> very large - on order of 10^26.  From a search through relevant
> literature, I know that matlab's preconditioned conjugate gradient (pcg)
> routine works well for the type of problem I'm dealing with.  Is there
> any similar routine in scipy?

Supply a preconditioner M for the conjugate gradients:

	scipy.sparse.linalg.cg(A, b, M=M)

Matlab's `pcg` is completely identical to this, AFAIK -- also with it you 
have to supply a suitable preconditioner yourself.

For a suitable preconditioner, you may want to look at PyAmg [1] or try 
something simpler such as Jacobi or SOR. Right now, Scipy does not have 
interfaces to incomplete LU or Cholesky, though.

.. [1] http://code.google.com/p/pyamg/

-- 
Pauli Virtanen