CODEWITHSUNDEEP
pythonmathoptimizationscipymathematical-optimization
wa4557
wa4557

Reputation: 1045

scipy minimize with constraints

Let's suppose I have a matrix

arr = array([[0.8, 0.2],[-0.1, 0.14]])

with a target function

def matr_t(t):
    return array([[t[0], 0],[t[2]+complex(0,1)*t[3], t[1]]]

def target(t):
    arr2 = matr_t(t)
    ret = 0
    for i, v1 in enumerate(arr):
          for j, v2 in enumerate(v1):
               ret += abs(arr[i][j]-arr2[i][j])**2
    return ret

now I want to minimize this target function under the assumption that the t[i] are real numbers, and something like t[0]+t[1]=1.

Upvotes: 44

Views: 110336

Answers (2)

cottontail
cottontail

Reputation: 23351

Instead of writing a custom constraint function, you can construct a scipy.optimize.LinearConstraint object and pass it as the constraint. Its construction asks for upper and lower bounds; also, the vector of independent variables has to have the same length as the variable length passed to the objective function, so the constraint such as t[0] + t[1] = 1 should be reformulated as follows (because t is length 4 as can be seen from its manipulation in matr_t()):

constraint

Also minimize optimizes over the real space, so the restriction of t[i] being real is already embedded into the algorithm. Then the complete code becomes:

import numpy as np
from scipy.optimize import minimize, LinearConstraint

def matr_t(t):
    return np.array([[t[0], 0],[t[2]+complex(0,1)*t[3], t[1]]])

def target(t):
    arr2 = matr_t(t)
    ret = 0
    for i, v1 in enumerate(arr):
          for j, v2 in enumerate(v1):
               ret += abs(arr[i][j]-arr2[i][j])**2
    return ret

arr = np.array([[0.8, 0.2],[-0.1, 0.14]])
linear_constraint = LinearConstraint([[1, 1, 0, 0]], [1], [1])
result = minimize(target, x0=[0.5, 0.5, 0, 0], constraints=[linear_constraint])

x_opt = result.x            # array([ 0.83,  0.17, -0.1 , 0.])
minimum = result.fun        # 0.0418

For a more involved example, let's use a common problem in economics, Cobb-Douglas utility maximization as an illustrative example. This is actually a constrained maximization problem but because minimize is a minimization function, it has to be coerced into a minimization problem (just negate the objective function). Also in order to pass the constraints as a scipy.optimize.LinearConstraint object, we have to write them to have lower and upper bounds. So the optimization problem is as follows:

problem

In this function, there are two variables x and y; the rest are all hyperparameters that have to be externally supplied (alpha, beta, px, py and B). Among them only alpha and beta are parameters of the objective function, so they must be passed to it through args= argument of minimize (alphas in the example below).

The optimization problem solves for x and y values where the objective function attains its minimum value given the constraint. They must be passed as a single object (variables in the function below) to the objective function. As mentioned before, we must pass an educated guess for these variables (x and y of the objective function) in order for the algorithm to converge.

from scipy.optimize import minimize, LinearConstraint

def obj_func(variables, hyperparams):
    x, y = variables
    alpha, beta = hyperparams
    return - x**alpha * y**beta

B, px, py, alphas = 30, 2, 5, [1/3, 2/3]
linear_constraint = LinearConstraint([[px, py], [1, 0], [0, 1]], [-np.inf, 0, 0], [B, np.inf, np.inf])
result = minimize(obj_func, x0=[1, 1], args=alphas, constraints=[linear_constraint])

x_opt, y_opt = result.x     # 4.9996, 4.0001
optimum = result.fun        # -4.3088

result.x are the minimizers and result.fun is the local minimum.

Cobb-Douglas has a closed form solution where for the example input, the correct solution is (x_opt, y_opt) = (5, 4). The result of minimize is not quite equal to that but since minimize is an iterative algorithm, this is as close as it got before it stopped.

Upvotes: 5

askewchan
askewchan

Reputation: 46578

This constraint

t[0] + t[1] = 1

would be an equality (type='eq') constraint, where you make a function that must equal zero:

def con(t):
    return t[0] + t[1] - 1

Then you make a dict of your constraint (list of dicts if more than one):

cons = {'type':'eq', 'fun': con}

I've never tried it, but I believe that to keep t real, you could use:

con_real(t):
    return np.sum(np.iscomplex(t))

And make your cons include both constraints:

cons = [{'type':'eq', 'fun': con},
        {'type':'eq', 'fun': con_real}]

Then you feed cons into minimize as:

scipy.optimize.minimize(func, x0, constraints=cons)

Upvotes: 73

Related Questions