1

I want to do constrained optimisation using a vector of constraints using the scipy.optimize library. In particular, I am supplying a vector of 3d coordinates r0 of N points -- hence a matrix of size N x 3 -- as input to the function. The coordinates are Cartesian, and I wish to freeze out all y dependence. So that means that I need the second column of my N x 3 matrix to be held to a constant, y0 say. How do I go about defining such a list of constraints?

To be concrete, let's the consider the COBYLA algorithm (https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.fmin_cobyla.html#scipy.optimize.fmin_cobyla). I tried the following construction:

cons = []
for i in range(xyz0.shape[0]):
    def f(x):
        return x[i,1]-xyz0cyl[i,1]
    cons.append(f)
fmin_cobyla(energy, xyz0, cons, rhoend=1e-7)

and got the error: 

41 for i in range(xyz0.shape[0]):
42     def f(x):
---> 43         return x[i,1]-xyz0cyl[i,1]
 44     cons.append(f)
 45 

IndexError: too many indices for array

What is going on?

ap21
  • 2,372
  • 3
  • 16
  • 32

1 Answers1

0

Your approach is wrong in quite a number of ways.

First, minimize takes a sequence as constraint, so that your Nx3 array is first flattened before it is passed to constraint functions leaving you with an array of only one dimension. Therefore you can't index with a tuple except you reshape your array inside the constraint functions to the original Nx3; could be pretty expensive for large N:

return x.reshape(-1, 3)[i,1] - xyz0cyl[i,1]

Secondly, closures in Python are late binding; all of the constraints functions will use the last value of i after the for loop as completed. You'll only finding out later on after fixing the first bug that the optimisation does not go as expected. See How do lexical closures work? to learn more.

A better approach is to actually make the y-axis (i.e. 1st column) stationary in your energy function or simply passing a Nx2 matrix instead to fmin_cobyla.

Moses Koledoye
  • 77,341
  • 8
  • 133
  • 139
  • I will try freezing out the y-axis. However, if I had to use the constraints, could I write `def f(x, i = i):` instead of `def f(x):`? Would that get around the lexical closure problem? – ap21 Jun 02 '18 at 22:16
  • @ap21 Yes, that would fix it. – Moses Koledoye Jun 02 '18 at 22:19
  • And when you say make the y-axis stationary, do you mean set it to a constant vector on every evaluation of the energy function? Or something smarter? – ap21 Jun 03 '18 at 14:59