# Python solvers for mixed-integer nonlinear constrained optimization

I want to minimize a black box function $f(x)$, which takes a 8$\times$3 matrix of non-negative integers as input. Each row specifies a variable, whereas each column specifies a certain time period so that $x_{ij}$ is the $i$th variable in the $j$th time period.

$f(x)$ is the activation function for a convolutional filter in the 2. convolutional layer of a convolutional neural network. It is therefore expected to be nonlinear and nonconvex. The input is subject to constraints as seen below, where $\text{sgn}(x)$ is the sign function as defined at Wikipedia.

\begin{aligned} \text{Objective:} \hspace{1cm} & \text{minimize} \hspace{0.2cm} f(x)\\ \text{Constraints:} \hspace{1cm} & \sum_{i=3}^{6}\sum_{j=1}^{3} x_{ij} = 10 & \\ & x_{3j} + x_{5j} \geq x_{1j} &\forall j = 1,2,3\\ & x_{4j} + x_{6j} \geq x_{2j} &\forall j = 1,2,3\\ & x_{7j} \leq 15 &\forall j = 1,2,3\\ & x_{8j} \leq 15 &\forall j = 1,2,3\\ & \text{sgn}(x_{7j}) = \text{sgn}(x_{3j}) &\forall j = 1,2,3 \\ & \text{sgn}(x_{8j}) = \text{sgn}(x_{4j}) &\forall j = 1,2,3 \\ & x_{ij} \in \mathbb{N_0} &\forall i = 1,2,\cdots,8 \hspace{0.2cm} \forall j = 1,2,3 \end{aligned}

Can anyone recommend any Python packages that would be able to solve this problem? Any commercial software with an interface to Python and a free academic license/evaluation period would also be great.

EDIT: It should be noted that the optimization does not have to find a global minimum (although that is, of course, preferred).

• I might not have explained the problem well enough (see the updated question). It is nonlinear and I am certain that the optimal solution will have D's unequal to 0. – pir Jun 7 '15 at 8:01
• See the updated question. – pir Jun 7 '15 at 11:52
• Thanks for finding time to look at the problem! I've looked more into CVXPY, but it cannot seem to find whether it can handle black box optimization? – pir Jun 7 '15 at 12:03
• @felbo: I'm guess that your nonlinearity consists of the restrictions on signs of certain pairs of unknowns. If you are asking how to solve this problem with off-the-shelf software packages, one approach would be breaking each of those pair-sign constraints into pairs of inequalities, e.g. $x_{7,j} \gt 0$ and $x_{3,j} \gt 0$ OR $x_{7,j} \lt 0$ and $x_{3,j} \lt 0$ (trusting you to fill in the details about how signs for zero values should be treated). But perhaps you are not as much interested in this one problem as in a broader survey of integer optimization software. – hardmath Jun 7 '15 at 14:00
• @felbo: I see. Given the difficult (nonlinear, nonconvex) behavior of your function $f(x)$, I would focus on effciently generating all feasible points $(x_{i,j})$ in your domain and evaluating the function at those points. Off-the-shelf software would not be likely to do anything more intelligent than such exhaustive search. – hardmath Jun 7 '15 at 14:19

I want to minimize a black box function $f(x)$, which takes a 8×3 matrix of non-negative integers as input.

In general, this sort of problem should be solved with a derivative-free MINLP (or if the function is linear, an MILP) solver. A very cursory glance at the literature suggests that the algorithm DFL, presented in a JOTA paper (also see preprint) could work; there might be other algorithms that also solve derivative-free MINLPs. The authors' implementation is in Fortran 90, so you would need to write some wrappers.

Generally speaking, though, you need some solver that can solve:

• black-box/derivative-free problems (I assume here that derivatives are not available, otherwise it would not be "black box")
• that are also mixed-integer

Since your problem contains no continuous decision variables, exhaustive sampling, as proposed by @hardmath, is another option that is probably easier to implement if you'd rather not write Python wrappers to a Fortran package (I wouldn't blame you).

However, IPOPT would not work, because it requires at least gradient information, and solves continuous problems, not mixed-integer problems. Even ignoring that estimating both gradient and Hessian information with finite-difference-type approaches is likely to yield poor outcomes due to numerical errors in these approximations, IPOPT will only solve the continuous relaxation of the problem. IPOPT (like any other continuous optimization solver) would have to be augmented with branch-and-cut or branch-and-bound-type methods, which is a great deal of work. If you're going to forge ahead and use a solver that assumes you have gradient information (even though you don't), you would be much better served using an MINLP solver (e.g., BARON, Couenne, Bonmin) that would at least give you integer-feasible solutions.

• Why is it mixed-integer? All the parameters are integers. – pir Jun 18 '15 at 15:11
• Purely integer programs are a subset of mixed-integer programs. Since methods for solving integer programs include approaches like relaxing the integer variables to continuous variables, then solving using a branch-and-bound (or branch-and-cut, or other) scheme, from an implementation perspective, implementing an integer programming solver frequently means mixed-integer programs can also be solved. – Geoff Oxberry Jun 18 '15 at 22:15
• Can you recommend an MINLP solver that has an interface to something more accessible than Fortran and C? A solver with a Python interface would be amazing. – pir Mar 7 '16 at 13:45
• @pir: I know there are black-box MINLP solvers implemented in MATLAB, but I don't know if they've been released. If you're just looking for MINLP solvers in the standard case where you know your objective and constraints (and these are not black box, plus you can compute derivatives), you might try interfacing to them via modeling languages like GAMS, AMPL, JuMP (Julia), or Pyomo (Python). There should be a MATLAB interface as well, I just don't remember what it is off the top of my head. – Geoff Oxberry Mar 7 '16 at 20:51
• Thanks! I am able to compute the first-order derivatives, but the function would still be black box in the sense that it cannot be formulated in a concise mathematical way (output of a neural network). Do you have any recommendations for which solver to use of the ones available in Pyomo? (software.sandia.gov/downloads/pub/pyomo/…) – pir Mar 7 '16 at 22:19

The feasible region can be broken up into several "convex" integer lattice subdomains by setting out the specific combinations of sign-agreement pairs. [That is, each sign agreement between $x_{7,j}$ and $x_{3,j}$, resp. between $x_{8,j}$ and $x_{4,j}$ can occur in one of two ways: both positive or both zero. Given the three "time" indices $j=1,2,3$, there are $2^6 = 64$ "components" which are the integer lattice solutions of linear constraints.]

When the region described by linear constraints is too complicated or large for direct generation of all feasible points, a popular technique is to generate a random sample of the points, e.g. by a Monte Carlo Markov Chain (MCMC) algorithm.

There is a Python implementation of MCMC samplers as pymc by Christopher Fonnesbeck, Anand Patil and David Huard. Although targeted to Bayesian estimation, the sampling portion of the toolkit could be used to generate "random" feasible points from each subdomain described above.

The minimum would then be approximated by the least value of $f(x)$ attained over all the sampled points.

Here is a sketch of how one might generate all the feasible points for these particular constraints.

Note that the "odd" row and "even" row indices $i$ are linked only by the sum-to-ten constraint at top. We can rearrange the unknowns from the $8\times 3$ matrix $x$ into a $4\times 6$ matrix:

$$\begin{pmatrix} x_{1,1} & x_{1,2} & x_{1,3} & x_{2,1} & x_{2,2} & x_{2,3} \\ x_{3,1} & x_{3,2} & x_{3,3} & x_{4,1} & x_{4,2} & x_{4,3} \\ x_{5,1} & x_{5,2} & x_{5,3} & x_{6,1} & x_{6,2} & x_{6,3} \\ x_{7,1} & x_{7,2} & x_{7,3} & x_{8,1} & x_{8,2} & x_{8,3} \end{pmatrix}$$

Each column of non-negative integers $(a,b,c,d)^T$ in this form obeys the same constraints:

$$a \le b + c$$

$${b = 0} \Leftrightarrow {d = 0}$$

$$d \le 15$$

We can add the constraint $b+c \le 10$ deduced from the combined middle two row sums in the new form being exactly ten, and we get a finite set of possibilities $S \subset \mathbb{N}_0^4$ from which all six columns are chosen.

Furthermore the sum-to-ten constraint for the two middle rows allows us to assign parts that add to $10$ for the six columns, thereby generating all the feasible points for this problem. There are $35$ partitions of $10$ in not more than six parts, and $\binom{15}{5} = 3003$ weak compositions of $10$ in exactly six summands, so this appears to be a tractable task computationally.

Added: A program written to carry out these calculations found some 8 trillion feasible solutions. More precisely there are 8,054,848,375,116 black box function evaluations needed to exhaust the search space.

• Sounds interesting! Note that $x_{i,j} \in \mathbb{N_0}$ $\forall i = 1,2,..,8$ $\forall j = 1,2,3$, allowing only two ways of sign agreement. Based on my understanding there would then be $2^{2 \cdot 3} = 64$ components. – pir Jun 7 '15 at 19:33
• How do you know that each integer lattice subdomain is "convex"? – pir Jun 7 '15 at 19:51
• The mention of "convex" was perhaps too terse/cryptic. Any nonempty region of $\mathbb{R}^n$ defined by linear equalities and inequalities is necessarily convex. What I meant by "convex" is the nonmepty intersection of the integer lattice $\mathbb{Z}^n$ with a convex subset of $\mathbb{R}^n$. Unfortunately there seems to be no succinct name for this. Because these sets are disconnected if they contain more than one point, generally these are not truly convex sets (though their theory and the application of MCMC algorithms is tied to convexity). – hardmath Jun 8 '15 at 0:09
• The constraint is not across all $x_{ij}$, but only for $i = 3,\cdots,6$. – pir Jun 8 '15 at 16:24
• @hardmath: The terminology I'd use is that the continuous relaxation of the problem is convex. – Geoff Oxberry Jun 17 '15 at 1:07

I agree with all of the answers provided here but I wanted to supplement with a Python implementation. We developed the Python GEKKO package for solving similar problems. We're also working on machine learning functions that may be able to combine a convolutional neural network with this constrained mixed-integer problem as a single optimization. Here is a potential solution with Python GEKKO (>0.2rc4).

import numpy as np
from gekko import GEKKO
m = GEKKO()
ni = 8
nj = 3
x = [[m.Var(lb=0,integer=True) for j in range(nj)] for i in range(ni)]
s = 0
for i in range(ni):
for j in range(nj):
s += x[i][j]
m.Equation(s==10)
m.Equations([x[2][j]+x[4][j]>=x[0][j] for j in range(nj)])
m.Equations([x[3][j]+x[5][j]>=x[1][j] for j in range(nj)])
for j in range(nj):
x[6][j].upper=15
x[7][j].upper=15
m.Equations([(m.sign3(x[6][j])==m.sign3(x[2][j])) for j in range(nj)])
m.Equations([(m.sign3(x[7][j])==m.sign3(x[3][j])) for j in range(nj)])
m.Obj(sum([m.tanh(x[i][0]) for i in range(ni)]))
m.options.SOLVER=1
m.solve()
print(x)


When you switch to IPOPT with m.options.SOLVER=3, it gives a non-integer solution, as expected:

---------------------------------------------------
Solver         :  IPOPT (v3.12)
Solution time  :   0.109899999995832      sec
Objective      :   7.933869530630442E-009
Successful solution
---------------------------------------------------

[[0.0], [0.2838769309], [0.2838769309]]
[[0.0], [0.2838769309], [0.2838769309]]
[[3.9669348007e-09], [0.99734110398], [0.99734110398]]
[[3.96693473e-09], [0.99734110398], [0.99734110398]]
[[0.0], [0.52439022369], [0.52439022369]]
[[0.0], [0.52439022369], [0.52439022369]]
[[0.0], [0.69439174379], [0.69439174379]]
[[0.0], [0.69439174379], [0.69439174379]]


However, when you switch to the APOPT MINLP solver with m.options.SOLVER = 1, it gives an integer solution:

---------------------------------------------------
Solver         :  APOPT (v1.0)
Solution time  :   3.900000000430737E-002 sec
Objective      :   0.000000000000000E+000
Successful solution
---------------------------------------------------

[[0.0], [2.0], [1.0]]
[[0.0], [1.0], [1.0]]
[[0.0], [0.0], [0.0]]
[[0.0], [0.0], [0.0]]
[[0.0], [2.0], [1.0]]
[[0.0], [1.0], [1.0]]
[[0.0], [0.0], [0.0]]
[[0.0], [0.0], [0.0]]


As you can see, the solution is fast but it is likely because the problem size is small and the objective function is not correct. For small problems, an exhaustive search may be the right solution. However, if you scale-up to larger problems with more integer variables or include the convolutional neural network model, each NLP branch and bound sub-problem will take longer to solve or may not converge.