[SciPy-dev] The future of the scipy.sandbox and a reminder of upcoming doc-day

[Robert Cimrman]

Nathan Bell wrote:
> On Jan 24, 2008 2:31 AM, Robert Cimrman <cimrman3 at ntc.zcu.cz> wrote:
> 
>> Besides the sparse stuff, I would also like to have support for the
>> symmetric eigenvalue functions of lapack that are currently in the
>> symeig package, but are not in scipy.linalg. lobpcg needs them to work
>> correctly.
> 
> We should formalize your idea to make a standard "dummy" matrix soon
> also.  How about the name LinearOperator?  Our goal with this should
> be to eliminate all of the code that checks for matvec(), psolve(),
> etc.

Good name!

> I imagine a definition like so:
> 
> class LinearOperator:
>     def __init__(self, shape, matvec, rmatvec=None,psolve=None):
>         self.shape = shape
>         self.matvec = matvec
> 
>         if rmatvec is not None:
>             self.rmatvec

self.rmatvec = rmatvec

>         else:
>             def rmatvec(x):
>                  raise NotImplementedError('LinearOperator does not
> define the operation x*A")
> 
> 
> A question arises when dealing with psolve() (the preconditioner).  We
> could either continue checking for the existence of psolve() in each
> method, or we could make LinearOperator have a dummy routine psolve(x)
> -> x and then write the methods so that the preconditioner is always
> applied.
 >
> The downside to this approach is that unnecessary copies may be
> performed.  OTOH one could write the method to avoid such problems (at
> worst, by checking to psolve() as is currently done).  Ideas?

I am not familiar with internals of scipy.linalg, but looking at 
iterative.py, it seems to me that functions in it check for psolve and 
if it is not defined in A, they use the no-op method psolve(x) -> x. So 
there is a precedent :)

IMHO if one must use sparse matrices, a (possible) copy of a vector or 
two does not matter much, as the main memory+performance hogger is the 
matrix.

-> I would start with the dummy psolve. When it is done, profiling may 
indicate some better way.


A side-note:

currently, the iterative solvers are defined as, for example:
bicg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None)
- the preconditioner is passed in as an attribute of A. This is not too 
transparent, IMHO.

It might be better to use instead
bicg(A, b, x0=None, tol=1e-5, maxiter=None, xtype=None, callback=None, 
precond = None),
where both A, precond would be LinearOperator instances. The 
preconditioning would then be performed by precond.matvec/rmatvec.

just my 1.5c
r.