Search code examples
scipylexical-scope

scipy - why isn't COBYLA respecting constraint?


I'm using COBYLA to do a cost minimization on a linear objective function with constraints. I'm implementing lower and upper bounds by including a constraint for each.

import numpy as np
import scipy.optimize

def linear_cost(factor_prices):
    def cost_fn(x):
        return np.dot(factor_prices, x)
    return cost_fn


def cobb_douglas(factor_elasticities):
    def tech_fn(x):
        return np.product(np.power(x, factor_elasticities), axis=1)
    return tech_fn

def mincost(targets, cost_fn, tech_fn, bounds):

    n = len(bounds)
    m = len(targets)

    x0 = np.ones(n)  # Do not use np.zeros.

    cons = []

    for factor in range(n):
        lower, upper = bounds[factor]
        l = {'type': 'ineq',
             'fun': lambda x: x[factor] - lower}
        u = {'type': 'ineq',
             'fun': lambda x: upper - x[factor]}
        cons.append(l)
        cons.append(u)

    for output in range(m):
        t = {'type': 'ineq',
             'fun': lambda x: tech_fn(x)[output] - targets[output]}
        cons.append(t)

    res = scipy.optimize.minimize(cost_fn, x0,
                                  constraints=cons,
                                  method='COBYLA')

    return res

COBYLA doesn't respect the upper or lower-bound constraints, but it does respect the technology constraint.

>>> p = np.array([5., 20.])
>>> cost_fn = linear_cost(p)

>>> fe = np.array([[0.5, 0.5]])
>>> tech_fn = cobb_douglas(fe)

>>> bounds = [[0.0, 15.0], [0.0, float('inf')]]

>>> mincost(np.array([12.0]), cost_fn, tech_fn, bounds)
       x: array([ 24.00010147,   5.99997463])
 message: 'Optimization terminated successfully.'
   maxcv: 1.9607782064667845e-10
    nfev: 75
  status: 1
 success: True
     fun: 239.99999999822359

Why wouldn't COBYLA respect the first factor constraint (i.e. upper-bound @ 15)?


Solution

  • COBYLA is in fact respecting all the bounds you give.

    The problem lies in the construction of the cons list. Namely, binding of variables in lambda and other inner-scoped functions in Python (and Javascript) is lexical, and does not behave in the way you assume: http://eev.ee/blog/2011/04/24/gotcha-python-scoping-closures/ After the loop is finished, the variables lower and upper have values 0 and inf, and variable factor has value 1, and these values are what is then used by all of the lambda functions.

    One workaround is to explicitly bind the specific values of variables to dummy keyword arguments:

    for factor in range(n):
        lower, upper = bounds[factor]
        l = {'type': 'ineq',
             'fun': lambda x, a=lower, i=factor: x[i] - a}
        u = {'type': 'ineq',
             'fun': lambda x, b=upper, i=factor: b - x[i]}
        cons.append(l)
        cons.append(u)
    
    for output in range(m):
        t = {'type': 'ineq',
             'fun': lambda x, i=output: tech_fn(x)[i] - targets[i]}
        cons.append(t)
    

    A second way is to add a factory function generating the lambdas.