Friday, 8 June 2007

Minimizing an objective function using Python

Thanks largely to physicists, Python has very good support for efficient scientific computing. The following code shows how to use the brute-force optimization function of scipy to minimize the value of some objective function with 4 parameters. Since it is a grid-based method, it's likely that you may have to rerun the optimization with a smaller parameter space.

import numpy
import scipy.optimize

def my_objective_fn(params):
print params,
value_of_objectivefn = calc_objectivefn_value(params)
print value_of_objectivefn
return value_of_objectivefn

if __name__=="__main__":
# Method 1 (same number of steps for each parameter)
myranges = ((21,49), (0,0.1), (95,123), (1.6, 2.0))
scipy.optimize.brute(my_objective_fn, myranges, Ns=5)

# Method 2 (can have different number of steps for
# different parameters)

slice_obj = numpy.s_[21:56:7, 0:0.1:0.02,
95:130:7, 1.6:2.1:0.1]
scipy.optimize.brute(my_objective_fn, slice_obj)

A number of cleverer optimization functions are also available in scipy, including constrained optimization and stochastic methods. For more information, install scipy (on Windows or Linux), and type the following at the Python prompt:

>>> import scipy.optimize
>>> help(scipy.optimize)
Optimization Tools

A collection of general-purpose optimization routines.

fmin -- Nelder-Mead Simplex algorithm
(uses only function calls)
fmin_powell -- Powell's (modified) level set method
(uses only function calls)


Felix said...

this is not really related to your post but to your category python:

I feel like I should know a different programming new language as a future theoretical chemist. Visual basic is all I know now but it doesn't seem to be the best choice

do you have a recommendation for what I should take a look at? python? C, fortran?


baoilleach said...

Languages differ in speed of execution, speed to write a program, how they make you think, features (e.g. object orientation), platforms supported, domain of applicability, idioms and conventions.

I think it is safe to say that VB is not a good choice simply on the basis that it is not available for Linux. So what should you do...

If you are going to be doing hard-core numerical algorithms, you'll be using either C, C++ or FORTRAN; the choice will usually depend on whether you are interfacing to some legacy code and on what your supervisor knows. Note that C++ is object-oriented, unlike the others.

For day-to-day work, you need to learn a scripting language. It's also useful for prototyping an algorithm. I recommend Python, which has libraries for efficient numerical computation, and encourages clean and elegant code. Other choices are Perl and Ruby. These are all free and available cross-platform. The number of third-party libraries available for a language like Python is unbelievable. With something like pyVTK you can visualise 3D objects and spin them around. matplotlib brings matlib-like plotting to python. And so on...

Also for day-to-day work, it's useful to learn common bash commands for text manipulation (I use these on Windows too by installing Cygwin), e.g. cut, paste, tail, head, awk, grep, find, sed.

Note also that Python (and other scripting languages) can interface to C, C++, and FORTRAN (although this is more difficult) libraries. In fact, this is how scipy links to the BLAS or ATLAS numerical libraries.

And...don't get caught in the language wars. A lot of hot air is wasted on C vs. FORTRAN, and Perl vs. Python. You will have to make your own decision.

Felix said...

thanks for the information. python seems pretty cool according to what you said and what I read on wikipedia. I'll check it out.

if I ever happen to do some serious programming then I'll just have to get used to a new syntax again

baoilleach said...

As an addition to the original article...

The Rosenbrock function in (x,y) is often used to test optimization routines. It has a minimum at (1,1):
f(x,y) = (1-x)2 + 100(y-x2)2

You can optimize functions in multiple variables in scipy if you use tuples for x0.
import scipy.optimize
lambda testfn (x,y): (1-x)**2 + 100*(y-x**2)**2
print scipy.optimize.fmin_bfgs(testfn, (2,2))

Dmitrey said...

see also scikits.openopt - new Python-written numerical optimization module from scipy developers