Quick and dirty minimization

Sometimes you just want to minimize a function. You're in the middle of your coding and you suddenly realize that you need to know the minimum of your 7 dimensional function, let's call it

def dirty_function(x[0], x[1], x[2]....):
  return some_complicated_value

So you hit the books. Problem is, when you start reading up on minimization algorithms, you'll find a whole smörgåsbord of complicated looking algorithms. You get all confused and you resign yourself to learning how to use a full-blown numerics library in order to minimize your dirty little function. There goes the rest of the week.

A couple of years ago, I had a conformation-analysis problem where I wanted to minimize some distances on a little peptide. I ended wasting a week reading up on Numerical Recipes in C. Their code was doubly terrible. It was terrible because their code was Fortran disguised in C, and it was doubly terrible because in their shitty license, they own the rights to their code and any modifications that you make.

It doesn't have to be like that. I've found that if you want something rough-as-guts, it's actually very simple to do a minimization in just a few lines of code.

There are roughly two types of minimization algorithms. Those that require a derivative, and those that don't. Since you have a dirty_function, you probably don't have an analytical derivative of dirty_function, so that rules out the algorithms that require a derivative (the best being Powell minimization). You have basically two options: downhill-simplex minimization and sequential-univariate minimization.

Although downhill-simplex is more robust, sequential-univariate minimization (even though its name is uglier) is way easier to implement.

The idea to seq-uni (as we'll call it from now on) is to pretend that your function is parabolic in all 7 dimensions. Let's say you are at

v = (x[0], x[1], x[2], ...)

So, for each iteration, you first pick two tiny steps in the first dimension, dimension 0:

x_old = x[0]
x_step1 = x_old + step
x_step2 = x_old + 2*step

Then you evaluate your dirty_function at:

y = dirty_function(x_old, x[1], x[2]....)
y1 = dirty_function(x_step1, x[1], x[2]...)
y2 = dirty_function(x_step2, x[1], x[2]...)

Now we pretend that (x_old y), (x_step1, y1), (x_step2, y2) are points on a parabola and calculate the x_min of the parabola.

def get_parabola_x_min(x1, y1, x2, y2, x3, y3):
  # assuming y = A*x^2 + B*x + C 
  term1 = (y3-y2) / ((x3-x2)*(x3-x1))
  term2 = (y1-y2) / ((x1-x2)*(x3-x1))
  A = term1 - term2
  B = (y3 - y2 + A*(x2*x2-x3*x3)) / (x2-x3)
  x_min = 0.5 * B / A
  return x_min;

x_min = get_parabola_x_min(x_old, y, 
                  x_step1, y1, x_step2, y2)

x_min is now your new point for value of your vector in dimension 0.

v_new = (x_min, x[1], x[2]....)

Now cycle through the rest of your dimensions for your first iteration. At the end of the cycle, we've found the minimum of the fake parabolic minimums for each dimension.

Repeat the iteration until the x_min for each dimension is smaller than a certain tolerance (I've found that the tiny step you take has to be smaller than the tolerance):

error = abs(x_min - x_old) 
if error < tolerance:
  finish

If the error is small, it means that approximating a new parabola from the starting point brings us back to the same point. As the approximation is good only near the real minimum, we should be there.

So there you have it. A fast and dirty minimization. But beware: seq-uni can converge slowly with some functions, but this is balanced by the fact that you don't make a lot of calls to dirty_function. And sometimes, seq-univ cycles around the minimum without terminating. Oh well, that's why it's a dirty method.