[SciPy-user] conjugate gradients solver - operator adjoint

[Pauli Virtanen]

Tue, 21 Apr 2009 23:42:01 +0200, lists kirjoitti:
> in my program for solving the integral equations in the liquid state
> theory (based on scipy/numpy; check pyoz.vrbka.net if you are interested
> in details) i need to implement a solver of a linear system AX=B -
> generally, this shouldn't be a big problem. there are some issues,
> however...
[clip]
> having a look at the sparse.linalg.solve.cg, there's no way of
> specifying the adjoint

The conjugate gradient method assumes that the matrix is hermitian, hence 
no need for specifying adjoint separately.

The algorithm you refer to is probably something different. A reference 
to the paper would come handy. Anyway, as far as I see, none of the 
solvers in sparse.linalg.* make use of the adjoint, so the particular 
algorithm is probably not implemented in scipy.

> (which should be, if i understand the method
> correctly, necessary for getting the result, since it's impossible to
> invert something that big efficiently).

This is probably debatable. There are many iterative algorithms that can 
invert non-symmetric matrices not needing to know the adjoint.

You could try using one of

	scipy.sparse.linalg.cgs
	scipy.sparse.linalg.bicgstab
	scipy.sparse.linalg.gmres

to see if one of them works well with your problem.

Note that in general, large problems may require preconditioning for the 
solvers to work efficiently. If you have the sparse matrix at hand and 
don't want to think, possibilities include incomplete (LU)
decompositions (scipy.sparse.linalg.splu), algebraic multigrid (see 
http://code.google.com/p/pyamg/), etc... Especially PyAMG could be worth 
trying out, as it's fairly easy to use.

-- 
Pauli Virtanen