MIPaaL: Mixed Integer Program as a Layer

TL;DR

This paper introduces a novel method for decision-focused learning that enables differentiation through Mixed Integer Linear Programs (MIPs) using a cutting planes approach, improving decision quality in complex optimization problems.

Contribution

It extends decision-focused learning to general MIPs by developing a differentiation technique via cutting planes, supporting arbitrary linear constraints over discrete and continuous variables.

Findings

Outperforms standard two-phase approaches in real-world tasks

Effective differentiation through MIPs improves decision quality

Supports broad class of problems with complex constraints

Abstract

Machine learning components commonly appear in larger decision-making pipelines; however, the model training process typically focuses only on a loss that measures accuracy between predicted values and ground truth values. Decision-focused learning explicitly integrates the downstream decision problem when training the predictive model, in order to optimize the quality of decisions induced by the predictions. It has been successfully applied to several limited combinatorial problem classes, such as those that can be expressed as linear programs (LP), and submodular optimization. However, these previous applications have uniformly focused on problems from specific classes with simple constraints. Here, we enable decision-focused learning for the broad class of problems that can be encoded as a Mixed Integer Linear Program (MIP), hence supporting arbitrary linear constraints over discrete and continuous variables. We show how to differentiate through a MIP by employing a cutting planes solution approach, which is an exact algorithm that iteratively adds constraints to a continuous relaxation of the problem until an integral solution is found. We evaluate our new end-to-end approach on several real world domains and show that it outperforms the standard two phase approaches that treat prediction and prescription separately, as well as a baseline approach of simply applying decision-focused learning to the LP relaxation of the MIP.

Tables

Table 1. Table 1 : Decision quality comparison of portfolio optimization (monthly percentage increase), and bipartite matching (number of pairs successfully matched). MIPaaL gives 2x monthly returns on SP500 and 8x on DAX compared to TwoStage, and successfully matches 40.3% more edges.

	SP500	DAX	Matching
MIPaaL-Exact	2.79 $\pm$ 0.17	5.70 $\pm$ 0.68	4.80 $\pm$ 0.71
MIPaaL-1000	2.60 $\pm$ 0.16	4.39 $\pm$ 0.66	3.45 $\pm$ 0.71
MIPaaL-100	1.25 $\pm$ 0.14	0.35 $\pm$ 0.63	2.57 $\pm$ 0.54
RootLP	1.97 $\pm$ 0.17	-1.97 $\pm$ 0.69	3.17 $\pm$ 0.60
TwoStage	1.19 $\pm$ 0.15	0.70 $\pm$ 1.46	3.42 $\pm$ 0.78

Table 2. Table 2 : Decision quality per instance win / loss percentages, with * indicating decision performance of winner is statistically significant. MIPaaL reliably improves or ties in these settings.

		MIPaaL-1000	MIPaaL-100	RootLP	TwoStage
	MIPaaL-Exact	57.2 / 42.3 *	87.2 / 12.8 *	82.8 / 17.2 *	90.6 / 9.4 *
SP500	MIPaaL-1000		84.4 / 15.6 *	80.0 / 20.0 *	85.0 / 15.0 *
	MIPaaL-100			27.8 / 72.2 *	52.8 / 47.2
	RootLP				69.4 / 30.6 *
	MIPaaL-Exact	66.7 / 33.3 *	86.6 / 13.4 *	88.9 / 11.1 *	84.4 / 15.6 *
DAX	MIPaaL-1000		81.1 / 18.9 *	85.6 / 14.4 *	73.9 / 26.1 *
	MIPaaL-100			69.4 / 30.6 *	43.3 / 56.7
	RootLP				28.3 / 71.7 *
	MIPaaL-Exact	11.7 / 3.9 *	15.6 / 2.8 *	15.0 / 3.9 *	12.2 / 5.6 *
Diverse Matching	MIPaaL-1000		11.7 / 3.9	9.4 / 8.3	10.0 / 7.2
	MIPaaL-100			3.9 / 11.1	5.0 / 8.9
	RootLP				7.2 / 8.3

Table 3. Table 3 : ML performance. TwoStage generalizes on metrics used for training (MSE, CE), whereas MIPaaL improves deployment loss and AUC. Edge prediction is hard as evidenced by AUCs near 0.5.

	SP500	DAX	Matching
	MSE	MSE	AUC	CE
MIPaaL-Exact	0.215 $\pm$ 0.043	0.126 $\pm$ 0.017	0.535 $\pm$ 0.004	0.658 $\pm$ 0.009
MIPaaL-1000	0.117 $\pm$ 0.020	0.349 $\pm$ 0.010	0.506 $\pm$ 0.007	0.614 $\pm$ 0.010
MIPaaL-100	0.983 $\pm$ 0.089	0.989 $\pm$ 0.060	0.503 $\pm$ 0.004	0.543 $\pm$ 0.013
RootLP	0.705 $\pm$ 0.178	1.055 $\pm$ 0.137	0.513 $\pm$ 0.001	0.493 $\pm$ 0.007
TwoStage	0.086 $\pm$ 0.017	0.022 $\pm$ 0.066	0.514 $\pm$ 0.005	0.392 $\pm$ 0.004

Table 4. Table 4 : Transfer learning across data distributions: trained on SP- 30 a superscript 30 𝑎 \text{30}^{a}

	Decision quality		MSE
	SP-${\text{30}}^b$	DAX	SP-${\text{30}}^b$	DAX
MIPaaL	2.02 $\pm$ 0.48	2.77 $\pm$ 0.40	4.81 $\pm$ 8.59	4.59 $\pm$ 8.80
RootLP	1.81 $\pm$ 0.44	1.74 $\pm$ 0.43	5.14 $\pm$ 1.02	5.39 $\pm$ 1.04
TwoStage	0.71 $\pm$ 0.04	0.82 $\pm$ 0.54	0.079 $\pm$ 0.052	0.065 $\pm$ 0.032

Table 5. Table 5 : Transfer learning across problem sizes: trained on SP- 30 a superscript 30 𝑎 \text{30}^{a} . Decision quality, monthly rate of return, on different problem sizes. MIPaaL ensures that the predictions perform well in multiple deployment settings whereas the TwoStage approach is unable to extract relevant information for larger scale problems.

	SP-50	SP-100	SP-200	SP500
MIPaaL	1.93 $\pm$ 0.13	2.27 $\pm$ 0.11	2.17 $\pm$ 0.48	2.26 $\pm$ 0.37
RootLP	1.50 $\pm$ 0.09	1.58 $\pm$ 0.08	1.82 $\pm$ 0.41	1.90 $\pm$ 0.29
TwoStage	1.58 $\pm$ 0.13	1.22 $\pm$ 0.09	1.50 $\pm$ 0.58	1.11 $\pm$ 0.35

Table 6. Table 6 : Transfer learning MSE results on different problem sizes: trained on SP- 30 a superscript 30 𝑎 \text{30}^{a} . The TwoStage approach clearly does better in approximating the distribution of objective coefficients even when the data distribution changes. However, the improvement in MSE does not ensure good deployment quality as show in 5

	SP-50	SP-100	SP-200	SP500
MIPaaL	5.421 $\pm$ 3.163	5.422 $\pm$ 2.368	5.254 $\pm$ 1.834	5.431 $\pm$ 1.670
RootLP	4.725 $\pm$ 3.169	4.876 $\pm$ 2.579	4.807 $\pm$ 1.906	4.834 $\pm$ 1.556
TwoStage	0.078 $\pm$ 0.019	0.074 $\pm$ 0.012	0.078 $\pm$ 0.011	0.076 $\pm$ 0.011

Equations

$$x^{*}(\hat{c};A,b,\mathcal{I})=\begin{array}[]{ll@{}ll}\text{argmin}_{x}&\hat{c}^{T}x&\\ \text{subject to}&Ax&\leq b\\ &x_{i}\in\mathbb{Z}&\,\forall i\in\mathcal{I}\end{array}$$

$$\begin{array}[]{ll@{}ll}\text{minimize}_{x}&\hat{c}^{T}x&\\ \text{subject to}&Ax&\leq b\\ &Sx&\leq t\end{array}$$

$$\hat{c}$$

$$\hat{x}$$

$$loss(\hat{c},c)$$

$$\begin{array}[]{rr@{}ll}\underset{w_{f},x,z_{1},z_{2},z_{3},f,y_{\text{tickets}},y_{\text{names}}}{\text{maximize}}&\sum_{i=1}^{n}\alpha_{i}w(i)&\\ \text{subject to}&\sum_{i=1}^{n}w_{f}(i)&=1\\ &\sum_{s\in S}x(s)&\leq\operatorname*{B_{\text{sector}}}\\ &\sum_{i=1}^{n}\operatorname*{y_{\text{names}}}(i)&\leq\operatorname*{B_{\text{name}}}\\ &\sum_{i=1}^{n}y_{\text{tickets}}(i)&\leq\operatorname*{B_{\text{ticket}}}\\ &w_{f}(i)-w_{t}(i)&\leq y(i)&\operatorname*{\forall i=1,\dots,n}\\ &-(w_{f}(i)-w_{t}(i))&\leq y(i)&\operatorname*{\forall i=1,\dots,n}\\ &w_{f}(i)-w_{0}(i)&\leq f(i)&\operatorname*{\forall i=1,\dots,n}\\ &-(w_{f}(i)-w_{0}(i))&\leq f(i)&\operatorname*{\forall i=1,\dots,n}\\ &\sum_{i=1}^{n}M_{s}(i)(w_{f}(i)-w_{t}(i))&\leq x(s)&\forall s\in S\\ &-\sum_{i=1}^{n}M_{s}(i)(w_{f}(i)-w_{t}(i))&\leq x(s)&\forall s\in S\\ &w_{f}(i)&\leq\operatorname*{y_{\text{names}}}(i)&\operatorname*{\forall i=1,\dots,n}\\ &f(i)&\leq y_{\text{tickets}}(i)&\operatorname*{\forall i=1,\dots,n}\\ &f(i)&=z_{1}(i)+z_{2}(i)+z_{3}(i)&\operatorname*{\forall i=1,\dots,n}\\ &0\leq z_{1}(i)&\leq 0.1\operatorname*{\text{vol}(i)}(i)&\operatorname*{\forall i=1,\dots,n}\\ &0\leq z_{2}(i)&\leq 0.2\operatorname*{\text{vol}(i)}(i)&\operatorname*{\forall i=1,\dots,n}\\ &0\leq z_{3}(i)&\leq 0.2\operatorname*{\text{vol}(i)}(i)&\operatorname*{\forall i=1,\dots,n}\\ &w_{f}(i),f(i),y(i)&\geq 0&\operatorname*{\forall i=1,\dots,n}\\ &z_{1}(i),z_{2}(i),z_{3}(i)&\geq 0&\operatorname*{\forall i=1,\dots,n}\\ &\operatorname*{y_{\text{names}}}(i),y_{\text{tickets}}(i)&\in\{0,1\}&\forall i=1,\dots,n\\ &x(s)&\geq 0&\forall s\in S\end{array}$$

$$\begin{array}[]{lr@{}ll}\underset{x}{\text{maximize}}&\sum_{i,j}c_{i,j}x_{i,j}&\\ \text{subject to}&\sum_{j}x_{i,j}&\leq 1&\forall i\in N_{1}\\ &\sum_{i}x_{i,j}&\leq 1&\forall j\in N_{2}\\ &\sum_{i,j}m_{i,j}x_{i,j}&\geq p\sum_{i,j}x_{i,j}\\ &\sum_{i,j}(1-m_{i,j})x_{i,j}&\geq q\sum_{i,j}x_{i,j}\\ &x_{i,j}&\in{0,1}&\forall i\in N_{1}\,,j\in N_{2}\end{array}$$

Full text

1 Introduction

We propose a method of training predictive models to directly optimize the quality of decisions that are made based on the model’s predictions. We are particularly interested in decision-making problems that take the form of Mixed Integer Programs (MIPs) because they arise in settings as diverse as electrical grid load control [32], RNA string prediction [40], and many other industrial applications [34]. MIPs naturally arise in so many settings largely due to their flexibility, computational complexity (ability to capture NP-hard problems), and interpretability. In many practical situations it is often necessary to predict some component (e.g., the objective) of the MIP based on historical data, such as estimated demand [35], price forecasts [14], or patient readmission rate [10]. However, these predictive models are often trained without regard for the downstream optimization problem. A standard loss function may place emphasis on correctly predicting outliers or small values because they greatly impact the loss function. However, in practice it may be the case that decisions will never be made on outliers, and so model capacity may be going towards adequately predicting values for a region of the feature space that is readily disregarded in the decision problem. Furthermore, practitioners may want to enforce that the outputs of the predictions meet semantically meaningful objectives such as ensuring predictions result in making fair decisions downstream [5, 43, 46].

Machine learning components commonly appear in larger decision-making pipelines; however, the model training process typically focuses only on a loss that measures accuracy between predicted values and ground truth values. Decision-focused learning, introduced for a financial criterion in [6] and extended to the more general quadratic programs in [2] and linear programs in [47], explicitly integrates the downstream decision problem when training the predictive model, in order to optimize the quality of decisions induced by the predictions. In the commonly used gradient-based predictive models, the central challenge is in passing gradients back to give the predictive model an indication of how it should shift its weights in order to improve decision quality of the resulting optimal solution. The discrete and discontinuous solution space that makes the MIP so widely applicable also prevents us from easily differentiating through it, as has been done for embedding continuous optimization problems in neural networks [2, 47, 13]. Our approach to computing gradients relies on the fact that we can algorithmically generate an exact continuous surrogate for the original discrete optimization problem. We employ previous work on cutting plane approaches, which exactly solve a MIP by iteratively solving a continuous relaxation and cutting off the discovered integer infeasible solution until an integer feasible solution is found [17]. The final continuous, and convex, optimization problem can then be used for backpropagation by differentiating the KKT conditions [21] of the continuous surrogate, as has been suggested for convex problems [2]. While pure cutting plane approaches are often slower than alternate branch-and-bound MIP solvers in practice [12], we note that our approach only needs the cutting plane methodology for backpropagation during training. Indeed, at test time, we can make predictions and find optimal decisions based on those predictions using any state-of-the-art MIP solver, ensuring the running time is exactly the same were any other training method used. Due to the computational complexity of computing the cutting planes for backpropagation, we also analyze a hybrid approach that stops cut generation after a fixed number of cuts have been generated, trading off the exactness of the differentiable solver with improved training runtime. An extreme version of this is when we disregard integrality requirements entirely and just use the LP relaxation of the problem for training. Finally, we evaluate our decision-focused approaches against the baseline of a two-stage approach, that simply relies on training using a relevant classification or regression loss function.

We demonstrate the effectiveness of our approach on two different real-world MIP problem domains concerning investment portfolio optimization and bipartite matching with diversity constraints, showing significant improvements in solution quality over the baseline.

2 Problem description

We consider problems that combine learning with an optimization problem that can be modeled as a mixed integer program (MIP). Specifically, each instance of the problem is a triple $(\phi,c,D)$, where $\phi$ is a feature vector, $c$ is the objective coefficient vector of a MIP, and $D$ represents additional known data that plays a role in the downstream optimization. In a MIP, $D$ will include the left hand and right hand constraint coefficients in each train instance $A,b$. If $c$ were known a priori, we could simply solve the MIP; however, we consider the setting where $c$ is unknown and must be estimated from $\phi$. We assume that we observe training instances $\{(\phi_{1},c_{1},D_{1}),\ldots,(\phi_{m},c_{m},D_{m})\}$ drawn from some distribution. We will use this data to train a predictive model $f_{\theta}$ (where $\theta$ denotes the internal parameters of the model) which outputs an estimate $f_{\theta}(\phi)=\hat{c}$ on a test-time instance. The standard two stage approach in this setting is to train the machine learning model $f_{\theta}$ that minimizes a loss comparing predicted values $\hat{c}$ and ground truth values $c$. Our objective is to find model parameters $\theta$ which directly maximize the expected quality of MIP solution for $\hat{c}$, evaluated with respect to the (unknown) ground truth objective $c$. We formalize this problem and our proposed approach below.

3 MIPaaL: Encoding MIP in a Neural Network

We formulate the MIP as a differentiable layer in a neural network which takes objective coefficients $\hat{c}$ as input and outputs the optimal MIP solution. Formally, we consider the optimal solution $x^{*}(\hat{c};A,b,\mathcal{I})$ of the MIP as a function of the input coefficients $\hat{c}$ given linear constraints on the feasible region $Ax\leq b$ and the set of integral variables $\mathcal{I}$. We write a functional form of the layer as:

[TABLE]

We can perform a forward pass given input objective coefficients, which are potentially outputs of a neural network, and feasibility parameters of the MIP using any solver.

Standard practice in this setting is to first train a model $f_{\theta}$ to predict the coefficients based on embeddings $\phi_{i}$ of the different predicted components such that on average the model predictions $\hat{c}=f_{\theta}(\phi)$ are not far away from the ground truth objective coefficients $c$. Then, decisions are made during deployment based on the predicted values by finding the optimal solution with respect to the predicted values $x^{*}(\hat{c};A,b,\mathcal{I})$ based on 1.

Forward propagation in this setting is straightforward; however, the highly nonconvex and discrete structure of the MIP, which enable it’s flexibility, hinders straightforward gradient computation for the backward pass.

We propose an approach to compute the backward pass which relies on finding a continuous optimization problem which is equivalent to the original optimization problem at the optimal integral solution. In particular, we use a pure cutting plane approach which is guaranteed to give the optimal solution (at the expense of potentially generating exponentially many cutting planes along the way) [9, 48].

The cutting plane approach iteratively solves the linear programming relaxation of the current problem. If the found solution is integral then the algorithm terminates since the found solution is both feasible to the original MIP and optimal for a relaxation of the original problem. Otherwise, a cut is generated which removes the discovered fractional solution and leaves all feasible integral solutions. Since the individual cuts do not remove any integral solutions, the final LP retains all integral solutions. In the extreme case, assuming the feasible region is bounded, we could describe the convex hull of the integral solutions yielding a linear program equivalent to the original MIP, thus ensuring that all potentially optimal integral solutions lie on extreme points of the feasible region. The simplex algorithm, a practically efficient LP solver, finds the optimal extreme point. In practice, finding the convex hull is intractable but we can obtain cutting planes that yield an exact integer solution via Gomory cuts or other globally valid cuts [17, 3].

We can then consider that we have generated cuts $Sx\leq t$ and write out the following equivalent linear program:

[TABLE]

Given this continuous optimization problem, we can now find the gradient of the optimal solution with respect to the input parameters by differentiating through the KKT conditions which are necessary and sufficient for the optimal solution of 2 (and hence also of the original MIP). This is done via the quadratic smoothing approach proposed in [47] for linear programs based on differentiating the KKT conditions of a quadratic program as shown in [2].

4 Decision-Focused Learning with MIPaaL

Training to minimize the error between predicted and ground truth values $\hat{c}$ and $c$ can yield decisions that have poor performance in practice in that the standard error metrics like mean squared error or crossentropy may not be directly aligned with the decision quality of the solutions obtained. To remedy this, we can train the predictive component to perform well by defining the loss to be the solution quality of the solution given the predicted coefficients. In other words, we can use the MIPaaL formulation to train the model to directly minimize the deployment objective. Using a neural network parametrized by $\theta$, $f_{\theta}$ to predict objective coefficients based on embeddings $\phi$ of the decision variables, we compute one forward pass as shown in 3 utilizing the cutting plane solver to generate the LP corresponding to the original MIP.

[TABLE]

Since the dot product in 3c and prediction of the neural network in 3a are differentiable functions of their inputs, the parameters of the neural network predictor $\theta$ can be trained via backpropagation using KKT conditions of the computed surrogate LP found in 2, relying on a small quadratic regularization term for the LP proposed in [47] to enforce strong convexity and perform backpropagation through 3b.

5 Empirical Evaluation

We instantiate MIPaaL for a range of tasks which require predicting and optimization with the overall goal of improving the objectives upon deployment of the recommendations. Specifically, we run experiments on combinatorial portfolio optimization and diverse bipartite matching. The portfolio optimization setting accounts for various combinatorial constraints enforcing small cardinality of the assets in the portfolio and limiting rebalancing transactions to maximize the overall predicted return of the portfolio. Diverse bipartite matching enforces diversity in the types of pairs that are recommended to maximize an overall predicted utility of the suggested pairing. In each setting, the predictive problem is nontrivial since the features do not contain much signal for the predictive task. However, we demonstrate that in these settings, a predictive model can perform well when given the specific task of ensuring that the predictions yield high-quality decisions.

5.1 Decision-Making Settings

Combinatorial Portfolio Optimization arises in several financial applications, with convex variants motivating early work in developing loss functions directly tied to decision quality [6]. However, the previously considered Markowitz portfolio optimization problem [31] is unable to capture discrete decisions which are desired or necessary to meet combinatorial deployment requirements as done in [8], where the authors use MIP to maximize percentage increase while limiting fixed costs incurred on trades, and limiting position changes for individual industries. We use a slightly modified version of the MIP from [8] which we specify fully in the appendix. The resulting formulation leverages the real-world modeling flexibility that MIPs afford: including big-M constraints, piecewise linear function modeling, and logical implication.

In the combined prediction and optimization problem, the next period’s percent increase in different assets (i.e., the objective) are unknown and must be learned from historical data. We gather features from the Quandl WIKI and FSE datasets [38] (which aggregate indicators for each company). We evaluate on the SP500, a collection of the 505 largest companies representing the American market, and the DAX, a set of 30 large German companies. We split the data temporally into train, test, and validation sets. More details about the data collection may be found in the appendix.

Diverse Bipartite Matching

Bipartite matching is used in many applications, where the benefit of a particular matching is often estimated from data. Without additional constraints, bipartite matching can be formulated as an LP; this formulation has previously been used for decision-focused learning [47] since the LP relaxation is exact. In practice though, matching problems are often subject to additional constraints. For instance, finding fair and high-quality housing allocations or kidney matching with additional contingency plans [5, 15] require additional decision variables and constraints on the solutions which make the integer problem NP-Hard and thus not solvable using a polynomially-sized LP, unless P=NP.

We use the problem of bipartite matching with added diversity constraints, which enforce a maximum and minimum bound on the percent of edges selected with a specified property. We use the experimental setup of [47], who did not include diversity constraints. Specifically, the matching instances are constructed from the CORA citation network [41] by partitioning the graph into 27 instances on disjoint sets of nodes (split into train, test and validation). Diversity constraints are added to enforce that at least some percent of edges selected are between papers in the same field and some percent are between papers in different fields.

5.2 Baselines

**Two-Stage: ** We compare against the standard predict-then-optimize approach which treats prediction and optimization components separately. The predictive component is trained to minimize a standard loss between predicted objective coefficients and the ground truth (e.g., mean squared error or cross-entropy). Afterwards, we solve the MIP to optimality using the predicted coefficients.

**RootLP: ** Next, we compare against an alternate decision-focused learning method from [47], which uses only the naive LP relaxation of the MIP (disregarding integrality constraints). The predictive model is trained using the LP relaxation, and at test time we solve the true MIP using the predicted objective coefficients to obtain an integral decision. This tests the impact of our cutting plane method, which allows us to fully account for combinatorial constraints in MIPaaL’s gradients.

**MIPaaL - k cuts: ** Given that the cut generation process is time consuming and must be done for each forward pass, we examine the tradeoff of decision quality when limiting the generated cuts to a fixed number $k$, and stopping cut generation after the first $k$ cuts have been generated. This essentially attempts to solve the original generally NP-Hard problem in polynomial time since each LP is solvable in polynomial time and we generate a constant number of cuts. We experiment with two settings ($k=100$ and $k=1000$) to determine how the exactness of the cut generation process impacts the decision quality at test-time. Note that in the case where $k=0$ cuts are generated, this method is equivalent to the RootLP approach where no cuts are generated.

5.3 Metrics

We evaluate the deployment quality of MIPaaL’s outputs as well as the similarity of the predicted MIP inputs to the ground truth objective coefficients. Where applicable, we provide $95\%$ confidence half-widths to indicate our confidence in the mean evaluation. For comparative metrics, we use a one-sided paired t-test with significance level of 0.05 to indicate whether we can reject the hypothesis that observed improvement is due to random chance.

Decision quality: The decision quality of a given model’s outputs is determined by the objective value after deployment of the model. A trained ML model predicts the objective coefficients of the particular decision problem, and CPLEX, a fast commercial optimization solver, gives an optimal solution based on those predictions. We evaluate the proposed solution based on the objective value with respect to realized ground truth objective coefficients. For portfolio optimization, the decision quality corresponds to the percentage change in the portfolio value from one time period to the next so the $2.79\%$ that MIPaaL-Exact achieves in 1 indicates that with this method we increase the portfolio value by $2.79\%$ every month on average. In bipartite matching, our goal is to maximize the number of successful matches and so the solution quality, as in previous work [47], corresponds to the average number of successful matches in the given instance. Top entries are bolded along with any method whose confidence interval overlaps with that of the best entry.

ML performance: We show the predictive performance for the given models in 3. For portfolio optimization (which is a regression problem), we report MSE, while for bipartite matching (a classification problem), we report Cross-Entropy (CE) and AUC. We also visualize each model’s predictions in 1 to better understand performance improvement.

5.4 Experimental Results

**Setup: **Experiments are run on a cluster of five identical 32-core machines with Intel 2.1 GHz processors and 264 GB of memory. We use the C API of IBM ILOG CPLEX 12.6.1 to generate and record cutting planes during training. Neural networks are trained and tested with the PyTorch Python API [36] and evaluated with the CPLEX’s Python API. We average over 5 training and testing iterations per problem setting with different seeds to evaluate the given approaches. This results in 180 portfolio optimization instances, and 135 diverse bipartite matching instances.

**Predictive architectures: **The neural networks for the portfolio optimization problems consist of two fully-connected layers with 100 nodes each, represented as a fully-connected layer, followed by batch normalization [20], LeakyRelu as a nonlinearity [30], and dropout with probability 0.5 as proposed and suggested in [42]. As portfolio optimization is built on a regression problem, we add a linear layer on the output of the predictive component. The architecture for matching is similar although it has only one layer and uses a sigmoid activation function on the output of the predictive component since the task is edge classification. These architectures were selected from grid search based on the task-based validation loss of the standard two stage model. The grid search was done considering models with between one and three layers, ReLU, LeakyReLU, and Sigmoid activation functions, and either 100 or 200 nodes per layer. We additionally, ran experiments disregarding dropout and batch normalization and found that without either the networks were unable to generalize. The models are trained with the Adam optimizer [26] with learning rate 0.01 and l2 normalization coefficient of 0.01 which were selected based on TwoStage decision quality.

Results: As shown in 1, the exact MIPaaL gives uniformly higher decision quality than baseline methods, with more than 2x the average return of the TwoStage or RootLP methods for portfolio optimization and 40.3% more successful matches for biparite matching. These results drive home the importance of integrating the full combinatorial problem into training, as enabled by MIPaaL. Comparing between versions of MIPaaL, we find that stopping cut generation at $k=1000$ results in competitive performance in an aggregate sense with the exact MIPaaL layer. Interestingly, stopping cut generation too early at $k=100$ underperforms RootLP on SP500 and Matching.

To further tease out the benefit of using variants of MIPaaL on a case-by-case basis, we compare decision quality in terms of win/loss rates in 2. This table shows the percent of instances the given approach on the left won/lost to the approach on the top, matching instances based on optimization problem and seed. Additionally, we compare the means of these values using a one-sided t-test with significance level 0.05. Asterisks indicate statistically significant results. Note that in portfolio optimization the allocated weights are somewhat continuous and so getting a tie is very rare. However, in the matching problem, we assign only integral pairs and get integral rewards if the match is successful or not so the objectives here are discrete and thus ties are more common.

As is apparent in 2, the MIPaaL method tends to outperform the other methods on all problem settings on an individual instances basis. Also useful to note is that the RootLP method beats the MIPaaL-100 approach in both the portfolio optimization settings, potentially indicating that limiting cut generation too much can be detrimental to the performance of the overall system. Additionally, in the diverse matching setting, the standard two stage approach does comparably to the $k$-cut decision-focused methods which stop cut generation early. This suggests that while the diverse matching problem is similar to the original bipartite matching LP, adding a few diversity constraints breaking the integrality property of the LP relaxation can remove gains made from decision-focused learning unless the exact MIPaaL layer is used.

Looking at the machine learning metrics in 3, we notice that the final machine learning performance of the decision-focused methods vary widely. In particular, the testing mean squared error for the portfolio optimization problems is quite high compared to the two stage approach. This mismatch between the MSE and decision quality exemplifies the need for training with the task in mind in that even though the MIPaaL model has worse MSE than TwoStage in both portfolio optimization settings, it results in much higher-return decisions.

1 gives a better understanding of how this occurs. Examining 1 and 1, we notice that even though the model learned with MIPaaL does poorly in general, it focuses in on identifying points that yield very high return while disregarding elements of mediocre return, leaving them in the center. The two-stage approach seen in 1 and 1 is focused more heavily on predicting many of the values somewhat well, disregarding that incorrectly predicting an asset to have high return may mean that it incorrectly ends up in the portfolio. Looking comparatively at the predictions of MIPaaL in 1 and TwoStage in 1, we can see that TwoStage predicts a large column of assets to have high returns when in practice they have relatively low or even negative returns. On the other hand, MIPaaL predicts relatively high returns correctly while mostly disregarding assets with lower returns. This is even more apparent looking at 1 and 1 as MIPaaL correctly predicts values for very profitable assets while TwoStage predicts relatively low or negative returns for most assets.

In terms of the matching problem, we see that even though TwoStage has better Cross-Entropy loss at test time, as it was trained with that specific classification loss in mind, it lacks in AUC which both corresponds to the findings in previous work [47], and indicates that the predictions learned by MIPaaL may sometimes also be accurate in a traditional sense.

Transfer learning: We evaluate MIPaaL, RootLP, and TwoStage on two transfer learning tasks for portfolio optimization. Here, models are trained on 30 assets randomly drawn from SP500 (SP-$\text{30}^{a}$), with data from January 2005 - December 2010. We then evaluate each model on data from December 2013 to November 2016 based on 1) SP-$\text{30}^{b}$, a set of 30 randomly drawn assets from the SP500, disjoint from SP-$\text{30}^{a}$, and 2) the DAX, a separate index comprising 30 companies from a different country. The transfer learning results in 4 demonstrate that MIPaaL not only performs well across time periods, but generalizes to unseen assets as well as unseen countries. On SP-$\text{30}^{b}$ MIPaaL gives more than double the improvement in the average rate of return over the standard TwoStage approach, and a 59% improvement over RootLP, indicating MIPaaL’s good generalization performance. Furthermore, while the transferred MIPaaL applied to DAX doesn’t beat MIPaaL trained on the same task, shown in 1, it improves performance over the RootLP and TwoStage in both settings, showing that MIPaaL learns a useful model for portfolio optimization as a whole rather than approximating the training data distribution. Lastly, MIPaaL’s performance improvement over TwoStage occurs despite a higher MSE, highlighting the importance of the decision-focused loss function.

Transfer Learning: Problem Size Generalization We evaluate the abilit of MIPaaL to perform well when predictions are required in differently sized settings. In particular, we evaluate when a model is trained to perform well on the portfolio optimization with 30 assets, but used in deployment on a larger number of assets. SP-50, SP-100, and SP-200 are instances drawn on 3 sets of assets which are all disjoint from both SP-$\text{30}^{a}$ and SP-$\text{30}^{b}$, with 50, 100, and 200 assets each. Again, we only evaluate on testing time periods to prevent data leakage and present results on the monthly rates of return in 5 and MSE values in 6.

As shown in 5, decision-focused approaches perform well when the predictions are used to find the optimal MIP solution. In this case the performance improvement is not as drastic as when the MIP is trained directly for the task in question. However, these results demonstrate that MIPaaL is able to efficiently extract information relevant to the general task of portfolio optimization, regardless of the set or number of assets in question.

6 Related Work

The interaction of machine learning and optimization has recently provided a new set of tools for efficiently solving a wide variety of problems. Recent work has framed traditionally heuristic components of optimization algorithms as machine learning tasks such as learning in branch and bound for MIPs [1, 25, 24, 29, 19, 39, 4, 7]. These approaches work inside an exact solver when the objective is known and not predicted from data. Deep reinforcement learning methods have also been used to generate high-quality data-driven solutions to hard problems such as graph optimization [23], vehicle routing [44, 27, 33], and realtime patrol planning in security games [49]. These problems also take as input different known features such as the true objective coefficients and are mainly used to help quickly generate high-quality solutions with a known input, which would be computationally intractable or take an unknown amount of time for an exact solver.

Several approaches have been proposed which embed optimization components in neural networks. This includes specific problems such as Markowitz portfolio optimization [6] or physically feasible state transitions [13], as well as larger classes which exhibit properties like convexity or submodularity. Problem classes used in end-to-end training include quadratic programs [2], linear programs [47], and zero-sum games [28, 37]. In addition, a solution is proposed for solving submodular optimization problems in [47]. In [16], the authors create a surrogate loss function and connect it to a decision-focused regret bound. Additionally, in [18] the authors instantiate end-to-end learning with linera models predicting components of a quadratic optimization function. Lastly, [45] incorporates a differentiable semidefinite programming problem as a relaxation for MAXSAT instances. To our knowledge, MIPaaL is the first approach for imbuing neural networks with the highly flexible Mixed-Integer Program, a widely-used class of potentially inapproximable NP-Hard optimization problems, while providing exact feedback on the decision quality.

7 Conclusion

We propose MIPaaL, a principled method for incorporating Mixed Integer Programs as a differentiable layer in neural networks. We approach the task of differentiating with respect to this flexible, discrete, and potentially inapproximable problem by algorithmically generating an equivalent continuous optimization problem via cutting planes. We instantiate our proposed approach for decision-focused learning wherein a predictive model is trained with a loss function that directly corresponds to the quality of the decisions made based on the predictions. MIPaaL is evaluated on two settings of portfolio optimization and one setting of bipartite matching with additional diversity constraints, which contain many modeling techniques widely used in combinatorial optimization that make the problem more complex but also more realistic. We demonstrate empirically that MIPaaL is able to outperform the standard approach of decoupling the prediction and decision components, as well as an approach of just using a continuous relaxation of the original combinatorial optimization problem. To better understand the impact of the cutting plane technique, we explore hybrid strategies that stop the cutting plane generation early. Ultimately, we find empirically that our approach can give high-quality solutions in the investigated settings.

Appendix A Portfolio Optimization Specifics

A.1 Data specification

We simulate our approach on this problem setting we use historical price and volume data downloaded from the Quandl WIKI dataset [38]. Our goal is to generate monthly portfolios of stocks in a given market index which targets a portfolio of these stocks weighted based on market capitalization, assuming that we start out with a portfolio that is weighted by market capitalization from the previous month. We train our model on data collected from January 2005 to December 2010, validate on data collected from January 2011 to November 2013 and test our model based on data from December 2013 to November 2016 resulting in 72 time periods used for training, 35 for validation and 36 for testing. We split the data temporally to ensure that data from the future isn’t used to inform predictions in the current time period. The 11 features used are made up of historical price averages and rates of return over different time horizons, which are commonly used indicators for simple trading strategies. Furthermore, to evaluate the generality of the approach we evaluate on two market indices: the SP500 and DAX. The SP500 is an index of 505 large companies in the united states which are widely representative of the US trading market. The DAX is composed of 30 large German companies and is much smaller than the SP500, being somewhat less representative of its respective market as a whole, but still representative of the central publicly traded companies in terms of trading volume and amount invested in those companies.

We use the following 11 indicators commonly used in simple momentum-based and mean-reversion trading strategies [11]:

percentage increase of price from previous 5 time periods

percentage increase of price from previous year

percentage increase of price from previous quarter

mean of price percentage increase over previous year

variance of price percentage increase over previous year

mean of price percentage increase over previous quarter

variance of price percentage increase over previous quarter

We randomly sample 60 indices from the SP500 for SP-$\text{30}^{a}$ and SP-$\text{30}^{b}$, these indices are:

SP-$\text{30}^{a}$ = LRCX, PHM, WY, SPG, EMN, CME, AEP, F, CAG, FISV, WBA, XOM , NVDA, ETN, MDT, FL, HBAN, FFIV, BLK, IPG, EXPD, IRM, PH , DLTR, COST, NBL, INCY, CSX]

SP-$\text{30}^{b}$ = [DVA, DE, BAC, KLAC, ADBE, FIS, IT, KR, FMC, HOG, SHW, RE , ETR, BK, ACN, NWL, ESS, VMC, C, EW, IR, SWKS, SNPS, ARE , SCHW, WEC, IVZ, SLB

Indices for the SP-50, SP-100, and SP-200 are:

SP-50 = VLO, EIX, RSG, PKG, CMI, AMZN, CHRW, ATVI, ADS, RHI, DVN , WMB, ILMN, LOW, SIVB, T, NTRS, HCP, ALXN, CPB, MRO, MU, CB , HES, APA, VFC, CHD, CMCSA, ALK, PBCT, BBY, MNST, UDR, CTXS , AZO, XRX, SBAC, M, MTD, COP, UTX, KO, MCO, TGT, CELG, HSIC , WYNN, YUM, ECL, ABT

SP-100 = RTN, K, ROST, WAB, MMM, WMT, JPM, MMC, AMD, BA, NEM, JBHT , STI, ITW, NKTR, ETFC, PNW, HSY, EL, MLM, UPS, FRT, ES, ROL , CMS, MAC, PVH, TSN, ANSS, SJM, DISH, LUV, MOS, TTWO, SEE , PKI, BAX, AMG, EOG, AMGN, CCI, EA, AGN, HP, HAL, AXP, SRE , INTU, AOS, APH, ROK, ABMD, CINF, UNM, ZION, HRL, ADSK, MSFT , JNJ, EBAY, ADI, EXC, CVS, PEG, CVX, PNR, PAYX, CPRT, DHR , JCI, SYK, HST, MO, INTC, GD, MCHP, PFG, PG, QCOM, STX, HAS , WM, CL, CDNS, BF_B, HIG, REGN, COG, RL, GPS, APD, APC, GPC , TSCO, JKHY, FAST, NOV, TMK, AVY, PGR

SP-200 = SYMC, WDC, CSCO, STZ, NUE, MAA, ARNC, DGX, EQIX, PWR, ED , LLY, AIV, ISRG, FITB, EQR, BXP, ORLY, MHK, JEC, UNP, TXT , IDXX, WHR, OMC, KEY, NI, FLIR, SNA, MCK, PXD, MCD, XEC, LEG , KSU, TSS, IP, AVB, RF, GPN, AFL, UNH, DTE, NEE, ZBH, NDAQ , DRE, MSI, VRSN, A, XEL, MAS, LKQ, FDX, LMT, WFC, NRG, RHT , COF, RMD, SLG, TROW, ANTM, BRK_B, PSA, RJF, BSX, AJG, FLR , BLL, AES, ORCL, MAR, HD, GILD, RCL, TIF, PRU, SWK, TMO, NOC , NKE, IFF, DHI, STT, AMT, TJX, PRGO, GS, NFLX, GLW, PFE, DOV , OKE, PLD, VZ, LNC, PNC, CTSH, JNPR, LEN, ABC, CI, CCL, LH , FCX, FLS, PCAR, HRS, EMR, GRMN, AKAM, BBT, PEP, OXY, KSS, MS , URI, AEE, NSC, CAH, VTR, O, USB, SYY, FE, TXN, BDX, BWA , CMA, COO, NTAP, SBUX, GWW, CNP, HOLX, CNC, D, TAP, MET, AAP , TFX, IBM, XLNX, LLL, AMAT, HUM, AON, DIS, KMB, BEN, GIS, ADM , MXIM, PPL, HPQ, ALL, VRTX, L, KIM, DOW, HON, DRI, AAPL, XRAY , MKC, MRK, HFC, FOX, CTAS, REG, HRB, KMX, SO, MYL, CLX, CERN , ALB, BIIB, LNT, BMY, MAT, DUK, MTB, AME, LB, GE, VNO, ROP , EFX, AIG, JWN, CAT, UHS, WAT, ADP, ATO, VAR, FOXA, MDLZ

A.2 Optimization model

We specify the optimization model used for the portfolio optimization task.

For our applications, we set the ticket budget $\operatorname*{B_{\text{ticket}}}$, and $\operatorname*{B_{\text{name}}}$ to be half of the number securities considered, 200 for SP500, and 15 for DAX and SP-$30$. We set the sector budget $\operatorname*{B_{\text{sector}}}$ to be 0.1, in that for any given sector we can change at most 10% of the portfolio weight into or out of that sector. Overall we found that this gave us a balance of reasonable problems that weren’t trivial to solve but still had a non-empty feasible region.

Given: Set of $n$ assets $i=1,\dots,n$, set of sectors $S$ represented as a partitioning of the $n$ assets

Ticket (trading) limit $\operatorname*{B_{\text{ticket}}}$

Name (unique item changes) limit $\operatorname*{B_{\text{name}}}$

Sector deviation limit $\operatorname*{B_{\text{sector}}}$

For each asset $i=1,\dots,n$:

Expected returns $\alpha_{i}$

Trading volume $\operatorname*{\text{vol}(i)}$

Indicator whether asset $i$ belongs to sector $s$

$M_{s}(i)\equiv i\in s$

Initial portfolio asset weight $w_{0}(i)$

Target portfolio asset weight $w_{t}(i)$

Goal: find final portfolio weights $w_{f}(i)$ which is close to the target portfolio $w_{t}$(i) in terms of weight, doesn’t incur too much cost from $w_{0}(i)$ to execute, respects budget constraints in terms of sector exposure, cardinality, and trading limits, and ensuring all weights add up to 1

Decision variables: final weight in asset $i$ $w_{f}(i)\in[0,1]\forall i=1,\dots,n$

auxiliary ticket-counting variables $z_{1}(i),z_{2}(i),z_{3}(i)$

auxiliary indicator of whether asset $i$ is used $y_{\text{names}}(i)$

auxiliary indicator of whether asset $i$ changes weight in the portfolio $y_{\text{tickets}}$

auxiliary variable $f(i)$ representing absolute value of weight change from original, keeping track of tickets bought

auxiliary variable $y(i)$ representing absolute value of deviation from target

auxiliary variable $x(s)$ representing absolute value of sector weight change

auxiliary variables $z_{1}(i),z_{2}(i),z_{3}(i)$ to put weights in three different compartments of a piecewise linear function related to volume.

Formulation:

[TABLE]

Overall this problem has $\lvert S\rvert+6n$ continuous decision variables, $\lvert 2n\rvert$ binary decision variables, and $10n+\lvert 2S\rvert+4$ constraints.

Appendix B Bipartite matching

B.1 Data specification

We run experiments on a version of the bipartite matching problem used in [47] which requires additional diversity constraints on the proposed matching. The matching problem is done on graphs sampled from the CORA citation network dataset [41] which consist of nodes and edges representing publications and citations respectively. In addition, we use information present about a given paper’s field of study. The full network we consider consists of 2708 nodes which are partitioned into 27 matching instances using metis [22]. Each instance is a bipartite matching problem on a complete bipartite graph with 50 nodes on each side of the graph. The 100 nodes in each instance are divided to maximize the number of edges between the two partitions. To ensure the diversity of our recommendations, we require a minimum $p=25\%$ percentage of the suggested pairings belong to distinct fields of study and similarly impose a constraint that a minimum percentage $q=25\%$ of suggested citations belong to the same field of study. The predicted values in this case are the edge weights which correspond to whether one paper actually cites another. The node representations in this case are binary feature vectors which correspond to whether a given word from a lexicon of 1433 words appears in the paper. Edge values are predicted based on the concatenation of the node features on the edge endpoints.

B.2 Optimization model

Given: two sets of nodes representing publications $N_{1},N_{2}$, weights on pairs of nodes $c_{i,j}\forall i\in N_{1},j\in N_{2}$ corresponding to how likely it is that one paper will cite another, same field indicator $m_{i,j}=\{1\text{ if i and j are in the same field},0\text{ otherwise}\}$

Goal: find matching of nodes such that 1) each node is matched at most once, 2) at least $p\in[0,1]$ proportion of selected edges connect nodes of the same field ($m_{i,j}=1$) 3) at least $q\in[0,1]$ proportion of selected edges connect nodes of different fields ($m_{i,j}=0$)

Decision Variables: $x_{i,j}\in{0,1}\forall(i,j)\in N_{1}\times N_{2}$ corresponding to whether we use edge $(i,j)$ in the matching or not.

Formulation:

[TABLE]

Overall this problem has $\lvert N_{1}\rvert\times\lvert N_{2}\rvert$ decision variables, and $\lvert N_{1}\rvert+\lvert N_{2}\rvert+2$ constraints. In our setting, we had $p=q=0.25$ to validate our experiments as this setting resulted in the problem not simply being solved too easily while yielding feasible regions.