[SciPy-dev] Abstract vectors in optimization

[Ben FrantzDale]

I've started playing with SciPy for numerical optimization. So far, it looks
like a great set of useful algorithm implementations.

If I understand the architecture correctly, there is one change that could
make these libraries even more powerful, particularly for large-scale
scientific computing. The optimization algorithms seem to be designed to
work only with NumPy vectors. Mathematically, though, many of these
algorithms can operate on an arbitrary Hilbert space, so the algoirhtms can
be slightly rewritten to work with any objects that support addition, scalar
multiplication, and provide an inner product, and a few other things.

If these functions could be made more general in this way, I could see the
SciPy optimization tools being extremely useful for large-scale optimization
problems in which people already have the cost function, its gradient, and
the representation of a state vector implemented (e.g., in C++). It would be
easy for them to provide a Hilbert-space interface to Python.

For example, the nonlinear CG solver, fmin_cg, takes a cost function, and an
initial state, x0, and among other things, a number describing a norm to use
to test for convergence. At present, it requires the following functions
(among others) work:
   x0 = asarray(x0).flatten()
   len(x0)
   vecnorm(gfk, ord=norm) # where gfk is the derivative vector and norm is
passed in
   numpy.dot(gfk,gfk)
It was easy for me (new to Python, seasoned in C++) to rewrite this function
to work with a class providing a Hilbert-space interface, thereby removing
the dependence on any particular representation of the vector and gradient.
For the sake of testing, I wrote a wrapper class that wraps a NumPy array.
When rewriting fmin_cg to use that class, I don't need
asarray(x0).flatten(), len(x0) is provided as __len__(self), norm(self,
norm=2) is provided, and inner_product(self,other) is provided, along with
the arithmatic operations.

Has anyone considered this approach? It really seems like a small code
change that would both simplify the optimization code and make it
significantly more flexible. I've tried similar approaches for generic
high-performance optimization functions in C++, but I think Python's type
system and memory management makes it a better language to write this sort
of solver.


―Ben
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20090105/ee2d75a6/attachment.html>