Conjugate gradients minimizer

[Carl Banks]

Andrew Walkingshaw wrote:
> Hi. Ages since I last posted to c.l.p, and this is a long shot, but here
> goes.
> 
> Does anyone have a conjugate-gradients minimization routine (failing
> that, even a line-search-using-gradients algorithm would be extremely
> useful):

You don't want to use steepest gradient.  If you get it, take the time
to modify it into conjugate gradient.

Just out of curiosity: how many variables are we talking about?


> a) in "pure" python (NumPy is *just about* OK for my purposes; C
> extensions aren't);
[snip license issues I don't care about]


Numeric is your friend for optimization algorithms.  Numeric does the
iteration though large arrays for you, which is useful because doing
it in Python isn't very fast.  Because of this, most optimization
algorithms would be terrible written in pure Python.  I highly suggest
upgrading the "just about" if you're going to do numerical work.

However, conjugate-gradient might be the exception, seeing that it
uses vectors only; no matrices.

In Python, you could use the map and reduce functions, with some
functions from the operator module, to get vector-like operations.
Unfortunately, I can't think of a good way to multiply a vector by a
constant (maybe there's a fast iterator that always returns the same
item that can be passed to map)?  Using a lambda is about the best
you can do there.

This is probably the best you can do in pure Python.  Here's an
example of the outer loop of conjugate gradient; maybe that'll get you
started:

    from __future__ import nested_scopes
    import operator

    g1 = gradient(x)
    s = map(operator.neg, g1)

    while not reached_termination(g1):
        x = linesearch(s)
        g2 = gradient(x)
        beta = (reduce(operator.add,map(operator.mul(g2,g2)) /
                reduce(operator.add,map(operator.mul(g1,g1)))
        s = map(lambda si,gi:si*beta-gi, s, g2)
        g1 = g2


Feel free to use this little code snippet as you wish; I put it into
the public domain.


> ... which I could either adapt or even just *read*? I've looked at
> Numerical Methods in C, and though it makes the algorithms pretty clear,
> it's so unPythonic it sets my teeth on edge. 

One thing you should keep in mind: in numerical algorithms, where
speed is at a premium, you shouldn't worry about things being
"Pythonic".  In fact, Numerical Algorithms code is probably not just
unPythonic, but also unCic.  In such code, it is often worth it to
make your code ugly just a save a few CPU cycles here and there.

The language of choice for a lot of people doing numerical work is
still Fortran.  That should tell you something about their priorities.
And having worked with some of them, I can tell you they do know what
they're doing.  (But God help them if they want to write regular
code.)


-- 
CARL BANKS