A Surrogate Objective Framework for Prediction+Optimization with Soft Constraints
Year: 2021

Abstract:
Prediction+optimization is a common real-world paradigm where we have to predict problem parameters before solving the optimization problem. However, the criteria by which the prediction model is trained are often inconsistent with the goal of the downstream optimization problem. Recently, decision-focused prediction approaches, such as SPO+ and direct optimization, have been proposed to fill this gap. However, they cannot directly handle the soft constraints with the max operator required in many real-world objectives. This paper proposes a novel analytically differentiable surrogate objective framework for real-world linear and semi-definite negative quadratic programming problems with soft linear and non-negative hard constraints. This framework gives the theoretical bounds on constraints’ multipliers, and derives the closed-form solution with respect to predictive parameters and thus gradients for any variable in the problem. We evaluate our method in three applications extended with soft constraints: synthetic linear programming, portfolio optimization, and resource provisioning, demonstrating that our method outperforms traditional two-staged methods and other decision-focused approaches.

heading: 1 Introduction
text: Mathematical optimization (a.k.a. mathematical programming), e.g., linear and quadratic programming, has been widely applied in decision-making processes, such as resource scheduling [1], goods production planning [2], portfolio optimization [3], and power scheduling [4]. In practice, problem parameters (e.g., goods demands, and equity returns) are often contextual and predicted by models with observed features (e.g., history time series). With the popularity of machine learning techniques and increasing available data, prediction+optimization has become a normal paradigm [5]. ∗Contributed during internship at Microsoft Research. † Corresponding Author.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Prediction becomes critical to the performance of the full prediction+optimization workflow since modern optimization solvers (e.g., Gurobi [6] and CPLEX [7]) can already efficiently find optimal solutions for most large scale optimization problems. Traditionally, prediction is treated separately as a general supervised learning problem and learned through minimizing a generic loss function (e.g., mean squared error for regression). However, studies have shown that minimization of the fitting errors does not necessarily lead to better final decision performance [8, 9, 10, 11].
Recently, a lot of efforts on using optimization objective to guide the learning of prediction models have been made, which are decision-focused when training prediction models instead of using traditional prediction metrics, e.g. mean squared error losses. For linear objectives, the ‘Smart Predict then Optimize’ ([10]) proposes the SPO+ loss function to measure the prediction errors against optimization objectives, while direct optimization ([12]) updates the prediction model’s parameters by perturbation. For quadratic objectives, OptNet [13, 14] implements the optimization as an implicit layer whose gradients can be computed by differentiating the KKT conditions and then back propagates to prediction neural network. CVXPY [15] uses similar technologies with OptNet but extends to more general cases of convex optimization. However, all the above state-of-the-art approaches do not contain the soft constraints in their objectives. In this paper, We consider linear ‘soft-constraints’, a penalty in the form of max(z, 0), where z = Cx− d is a projection of decision variable x ∈ Rn and context variables C ∈ Rm×n, d ∈ Rm. Such soft constraints are often required in practice: for example, they could be the waste of provisioned resources over demands or the extra tax paid when violating regulatory rules. Unfortunately, the max(·, 0) operator is not differentiable and thus cannot be directly handled by these existing approaches. To differentiate soft constraints is a primary motivation of this paper.
In this paper, we derive a surrogate objective framework for a broad set of real-world linear and quadratic programming problems with linear soft constraints and implement decision-focused differentiable predictions, with the assumption of non-negativity of hard constraint parameters. The framework consists of three steps: 1) rewriting all hard constraints into piece-wise linear soft constraints with a bounded penalty multiplier; 2) using a differentiable element-wise surrogate to substitute the piece-wise objective, and solving the original function numerically to decide which segment the optimal point is on; 3) analytically solving the local surrogate and obtaining the gradient; the gradient is identical to that of the piecewise surrogate since the surrogate is convex/concave such that the optimal point is unique.
Our main contributions are summarized as follows. First, we propose a differentiable surrogate objective function that incorporates both soft and hard constraints. As the foundation of our methodology, in Section 3, we prove that, with reasonable assumptions generally satisfied in the real world, for linear and semi-definite negative quadratic programming problems, the constraints can be transformed into soft constraints; then we propose an analytically differentiable surrogate function for the soft constraints max(·, 0). Second, we present the derived analytical and closed-form solutions for three representative optimization problems extended with soft constraints in Section 4 – linear programming with soft constraints, quadratic programming with soft constraints, and asymmetric soft constraint minimization. Unlike KKT-based differentiation methods, our method makes the calculation of gradients straightforward for predicting context parameters in any part of the problem. Finally, we apply with theoretical derivations and evaluate our approach in three scenarios, including synthetic linear programming, portfolio optimization, and resource provisioning in Section 4, empirically demonstrate that our method outperforms two-stage and other predict+optimization approaches.

heading: 2 Preliminaries

heading: 2.1 Real-world Optimization Problems with Soft Constraints
text: Our target is to solve the broad set of real-world mathematical optimization problems extended with a max(z, 0) term in their objectives where z depends on decision variables and predicted context parameters. In practice, max(z, 0) is very common; for example, it may model overhead of underprovisioning, over-provisioning of goods, and penalty of soft regulation violations in investment portfolios. We call the above max(z, 0) term in an objective as soft constraints, where z ≤ 0 is allowed to violate as long as the objective improves.
The general problem formulation is
max x
g(θ, x)− αT max(z, 0), z = Cx− d, s.t. Ax ≤ b,Bx = c, x ≥ 0 (1)
where g is a utility function, x ∈ Rn the decision variable, and θ ∈ Rn the predicted parameters. Based on observations on a broad set of practical problem settings, we impose two assumptions on the formulation, which serves as the basis of following derivations in this paper. First, we assume A ∈ Rm1×n ≥ 0, b ∈ Rm1 ≥ 0, B ∈ Rm2×n ≥ 0 ,and c ∈ Rm2 ≥ 0 hold. This is because for problems with constraint on weights, quantities or their thresholds, these parameters are naturally non-negative. Second, we assume the linearity of soft constraints, that is z = Cx − d, where z ∈ Rm3 , C ∈ Rm3×n, and d ∈ Rm3 . This form of soft constraints make sense in wide application situations when describing the penalty of goods under-provisioning or over-provisioning, a vanish in marginal profits, or running out of current materials.
Now we look into three representative instances of Eq.1, extracted from real-world applications.
Linear programming with soft constraints, where g(x, θ) = θTx. The problem formulation is
max x
θTx− αT max(Cx− d, 0), s.t. Ax ≤ b,Bx = c, x ≥ 0 (2)
where α ≥ 0. Consider the application of logistics where θ represents goods’ prices, {A, B, C} are the capability of transportation tools, and {b, c} are the thresholds. Obviously, A, b,B, c ≥ 0 hold. Quadratic programming with soft constraints, where g(θ, x) = θTx − xTQx. One example is the classic minimum variance portfolio problem [16] with semi-definite positive covariance matrix Q and expected return θ to be predicted. We extend it with soft constraints which, for example, may represent regulations on portions of equities for some fund types. Formally with α ≥ 0, we have:
max x
θTx− xTQx− αT max(Cx− d, 0), s.t. Bx = c, x ≥ 0. (3)
Optimization of asymmetric soft constraints. This set of optimization problems have the objective to match some expected quantities by penalizing excess and deficiency with probably different weights. Such formulation represents widespread resource provisioning problems, e.g., power[4] and cloud resources[1], where we minimize the cost of under-provisioning and over-provisioning against demands. Formally with α1, α2 > 0, we have:
max x −(αT1 max(Cx− d, 0) + αT2 max(d− Cx, 0)), s.t. Bx = c, x ≥ 0. (4)
In this paper, we consider a challenging task where C is a matrix to be predicted with known constants d. In reality, the Cx− d term may represent the "wasted" part when satisfying the actual need of d.

heading: 2.2 Prediction+Optimization
text: For compactness we write Eq.1 as maxx∈X f(x, θ), where f is the objective function and X is the feasible domain. The solver for f is to solve x∗ = argmaxx∈X f(x, θ). With parameters θ known, Eq.2–4 can be solved by mature solvers like Gurobi [6] and CPLEX [7].
In prediction+optimization, θ is unknown and needs to be predicted from some observed features ξ ∈ Ξ. The prediction model is trained on a dataset D = {(ξi, θi)}Ni=1. In this paper, we consider the prediction model, namely Φ, as a neural network parameterized with ψ. In traditional supervised learning, Φψ is learned by minimizing a generic loss, e.g., L1 (Mean Absolute Error) or L2 (Mean Squared Error), which measures the expected distance between predictions and real values. However, such loss minimization is often inconsistent with the optimization objective f , especially when the prediction model is biased [11].
Instead, decision-focused learning directly optimizes ψ with respect to the optimization objective f , that is, maxψ E(ξ,θ)∼D[f(x̂∗(Φψ(ξ)), θ)]. The full computational flow is illustrated in Fig.1. In
the gradient-based learning, update of ψ’s gradient is ∂r∂ψ = ∂θ̂ ∂ψ
∂x̂∗
∂θ̂ ∂r ∂x̂∗ , where utility r = f(x̂∗, θ).
The Jacobian ∂θ̂/∂ψ is computed implicitly by auto-differentiation of deep learning frameworks (e.g, PyTorch [17]), and ∂r/∂x̂∗ is analytical. The main challenge is to compute ∂x̂∗/∂θ̂, which depends on differentiating the argmax operation. One recent approach is to rewrite the objective to be convex (by adding a quadratic regularizer if necessary), build and differentiate the optimality conditions (e.g.,
KKT conditions) [18] which map θ̂ to the solution x̂∗, and then apply implicit function theorem to obtain ∂x̂∗/∂θ̂. Alternatively, in this paper we propose a novel approach that rewrites the problem as an unconstrained problem with soft constraints and derives analytical solutions, based on our observations on the real-world problem structures and coefficient properties.

heading: 3 Methodology
text: Our main idea is to derive a surrogate function for f with a closed-form solution such that the Jacobian ∂x∂θ is analytical, making the computation of gradient straightforward. Unlike other recent work [13, 14, 15, 3], our method does not need to solve KKT optimality condition system. Instead, by adding reasonable costs for infeasibility, we convert the constrained problem into an unconstrained one. With the assumption of concavity, we prove that there exist constant vectors β1, β2, β3, such that Eq.1 can be equivalently transformed into an unconstrained problem:
max x L(x) = max x
g(x, θ)−αT max(Cx−d, 0)−βT1 max(Ax−b, 0)−βT2 |Bx−c|−βT3 max(−x, 0) (5)
The structure of this section is as follows. Section 3.1 proves that the three types of hard constraints can be softened by deriving bounds of β1, β2, β3; for this paper, we will assign each entry of the three vectors the equal value (with a slight abuse of notation, we denote β1 = β2 = β3 = β for the proofs of bounds; we align with the worst bound applicable to the problem formulation.) Section 3.2 proposes a novel surrogate function of max(·, 0), such that the analytical form of ∂x∂θ can be easily derived via techniques of implicit differentiation [18] and matrix differential calculus [19] on equations derived by convexity [20]. Based on such derived ∂x∂θ , we develop our end-to-end learning algorithm of prediction+optimization whose detailed procedure is described in Appendix C.

heading: 3.1 Softening the Hard Constraints
text: For any hard constraints w = Ax − b ≤ 0, we denote its equivalent soft constraints as H(w) = βT max(w, 0). H(w) should satisfy two conditions: 1) for w ≤ 0 (feasible x), H(w) = 0; 2) for w ≥ 0 (infeasible x), H(w) is larger than the utility gain(i.e., improvement of the objective value) R = f(x, θ) − maxx1:Ax1≤b f(x1, θ) by violating Ax − b ≤ 0. Intuitively, the second condition requires a sufficiently large-valued β > 0 to ensure that the optimization on the unconstrained surrogate objective never violates the original Ax ≤ b; to make this possible, we assume that the l2-norm of the derivative of the objective f before conversion is bounded by constant E. The difficulty of requirement 2) is that the distance of a point to the convex hull l is not bounded by the sum of distances between the point and each hyper-plane in general cases, so the utility gain obtained from violating constraints is unbounded. Fig. 2-(a) shows such an example which features the small angle between hyper-planes of the activated cone. We will refer such kind of ’unbounding’ as "acute angles" below.
The main effort of this subsection is to analyze and bound the effect caused by such acute angles. Given a convex hull C = {z ∈ Rn|Az ≤ b} (A ≥ 0 is not required here) and any point x 6∈ C, let x0 ∈ C be the nearest point of C to x, and A′x ≥ b′ represent all active constraints at x, then all such active constraints must pass through x0.2 The domain K = {z ∈ Rn|A′z ≥ b′} is a cone or degraded
2Otherwise, there must exist an active constraint i and a small neighbourhood of x0, say B(x0, ), such that ∀z ∈ B(x0, ), Aiz − bi 6= 0, which implies that either constraint i is inactive (< 0) or x0 6∈ C (> 0).
cone where the tip of the cone is a subspace of Rn. For the rest of the paper, we will call K activated cone, as shown in Fig. 2. Note that for any degraded activated cone, x− x0 is always perpendicular to the tip subspace; therefore, we may project the cone onto the complementary space of the tip and get the same multiplier bound on the projected activated cone with lower dimensions.
Ideally, we aim to eliminate the utility gain R obtained from violating A′x ≤ b′ with the penalty βT (A′x− b′), i.e., ensure βT (A′x− b′) ≥ R hold for any x. For the compactness of symbols and simplicity, we will assume that the main objective is g(x, θ) = θTx in this section; however note that our proofs apply with the existence of soft constraints and quadratic terms in g, which is discussed at the beginning of Appendix A. With such assumption, we now give a crucial lemma, which is the core of our proof for most of our theorem:
Lemma 1. (Bounding the utility gain) Let R = f(x, θ) − maxx1∈C f(x1, θ) be the utility gain, then R ≤ f(x, θ) − f(x0, θ) ≤ Ecos p0 ∑n i=1(A ′ ix − b′i), where x is an infeasible point, A′x ≤ b′ the active constraints at x, p0 = ∠(A′i ∗ , θ∗) where A′i
∗ and θ∗ are the optimal solution of maxθ minA′i cos∠(A ′ i, θ) (i.e., the maximin angle between θ and any hyperplane of the activated cone K = {z ∈ Rn|z − x0 ∈ cone(A′)}), and x0 the projection of x to the tip of cone C = {z ∈ Rn|A′z ≤ b′}. E is the upper bound of ||θ||2. ∠(·, ·) denotes the angle of two vectors.
Thus, Ecos p01 is a feasible choice of β, and it suffices by finding the lower bound for cos p0. For the rest of the section, we find the lower bound of cos p0 by exploiting the assumed properties of A′, e.g., A′ ≥ 0; we give the explanation of the full proofs for all theoretical results in Appendix A.

heading: 3.1.1 Conversion of Inequality Constraints Ax ≤ b
text: Let us first consider the constraints w = Ax − b ≤ 0, where A ≥ 0, b ≥ 0. It is easy to prove that given A ≥ 0 and b ≥ 0, the distance of a point to the convex hull Ax ≤ b is bounded. More rigorously, we have the following theorem, which guarantees the feasibility of softening the hard constraints of inequalities Ax ≤ b: Theorem 2. Assume the optimization objective θTx with constraints Ax ≤ b, where A ≥ 0, and b ≥ 0. Then, the utility gain R obtained from violating Ax − b ≤ 0 has an upper bound of O( ∑ i max(wi, 0)E), where w = A ′x− b′, and A′x ≤ b′ is the active constraints.

heading: 3.1.2 Conversion of Inequality Constraints x ≥ 0
text: With inequality constraints Ax ≤ b converted, we now further enforce x ≥ 0 into soft constraints. It seems that the constraint of this type may form a very small angle to those in w = Ax− b. However, as −x is aligned with the axes, we can augment x ≥ 0 into soft constraints by proving the following theorem:
Theorem 3. When there is at least one entry of x ≥ 0 in the activated cone, the utility gain R from deviating the feasible region is bounded by O(n 1.5E ∑ imax(wi,0)
sin p ), where p is the smallest angle between axes and other constraint hyper-planes and w is the union of Ax− b and −x.
Hence, we can set β = O(n 1.5E sin p ). Specially, for binary constraints we may set β = O(n 1.5E):
Corollary 4. For binary constraints where the entries of A are either 0 or 1, the utility gain R of violating x ≥ 0 constraint is bounded by O(n1.5E ∑ i max(wi, 0)), where wi = Aix− bi or −x.
which gives a better bound for a set of unweighted item selection problem (e.g. select at most k items from a particular subset).

heading: 3.1.3 Conversion of Equality Constraints Bx = c
text: Finally, we convert Bx = c into soft constraints. This is particularly difficult, as Bx = c implies both Bx ≤ c and −Bx ≤ −c, which will almost always cause acute angles. Let’s first consider a special case where there is only one equality constraints and A is an element matrix In×n. Theorem 5. If there is only one equality constraint Bx = c (i.e., B degrades as a row vector, such like ∑ i xi = 1) and special inequality constraints x ≥ 0, Ix ≤ b, then the utility gain R from violating constraints is bounded by O(n 1.5E ∑ i max(wi,0)
sin p ), where p is the same with theorem 3, w is the union of Bx− c and −x.
Intuitively, when there is only one equality constraint, the equality constraint can be viewed as an unequal one, for at most one side of the constraint can be in a non-degraded activated cone. Thus, we can directly apply the proof of Theorem 2 and 3, deriving the same bound O(n
1.5E sin p ) for β.
Finally, we give bounds for general Ax ≤ b, Bx = c with A, b,B, c ≥ 0 as below: Theorem 6. Given constraints Ax ≤ b, x ≥ 0, and Bx = c, where A,B, b, c ≥ 0, the utility gain R obtained from violating constraints is bounded by O( √ nλmax ∑ i max(wi, 0)), where λmax is the upper bound for eigenvalues of PTP (P : x→ Px is an orthogonal transformation for an n-sized subset of normalized row vectors in A,B and −I), and w is the union of all active constraints from Ax ≤ b, x ≥ 0, Bx ≤ c and −Bx ≤ −c.
In this theorem, P is generated by taking an arbitrary n-sized subset from the union of row vectors in A,B and −I , orthogonizing the subset, and using the orthogonal transformation matrix as P ; there are ( n+m1+m2
n
) different cases of P , and λmax is the upper bound of eigenvalues of PTP
over all possible cases of P . Note that there are no direct bounds on λmax with respect to n and the angles between hyper-planes. However, empirical results (see Appendix A for details) show that for a diverse set of synthetic data distributions, λmax = O(n2) follows. Therefore, empirically we can use a bound O(n
2.5E sin p ) for β. 3 So far, we have proven that all hard constraints can be transformed into soft constraints with bounded multipliers. For compactness, Eq.5 is rewritten in a unified form:
L(x) = g(x, θ)− γT ·max(C′x− d′), where γ = [
α
O(n 2.5E sin p )1
] , C′ =  C A −B B −I , d′ =  d b −c c 0  (6)

heading: 3.2 The Unconstrained Soft Constraints
text: As Eq.5 is non-differentiable for the max operator, we need to seek a relaxing surrogate for differentiation. The most apparent choice of the derivative of such surrogate S(z) for z = C ′x − d′ is sigmoidal functions; however, it is difficult to derive a closed-form solution for such functions, since sigmoidal functions yield z in the denominator, and z cannot be directly solved because C is not invertible (referring to Appendix B for detailed reasons). Therefore, we have to consider a piecewise roundabout where we can first numerically solve the optimal point to determine which segment the optimal point is on, and then expand the segment to the whole space. To make this feasible, two assumptions must be made: 1) this function must be differentiable, and 2) the optimal point must be unique; to ensure this, the surrogate should be a convex/concave piece-wise function. The second property is for maintaining the optimal point upon segment expansion. Fortunately, there is one simple surrogate function satisfying our requirement:
S(z) =  0 if z < − 1 4K K(z + 1 4K )2 if − 1 4K ≤ z ≤ 1 4K
z if z ≥ 1 4K
(7)
3For real optimization problems, we can sample A′ from A and estimate the largest eigenvalue of PTP .
Let M and U be diagonal matrices as the indicator of S(z), which are Mi,i = 2K[− 14K ≤ zi ≤ 1 4K ] and Ui,i = [ 14K < zi], where [·] is an indicator function. K > 0 is a hyper-parameter that needs to be balanced. Larger K makes the function closer to original; however, if K is too large, then the training process would be destabilized, because when the prediction error is large at the start of the training process, ∂f
∂θ̂ |x̂ might be too steep. Then consider the optimal point for the unconstrained optimization
problem maximizing θTx− γT max(C ′x− d′, 0), by differentiaing on both sides, we can obtain:
θ = C′TMdiag(γ)(C′x− d′) + C′T ( 1 4K M + U)γ (8)
This equation reveals the closed-form solution of x with respect to θ, and thus the calculation of ∂x∂θ becomes straightforward:
∂x ∂θ = (C ′TMdiag(γ)C ′)−1 (9)
C ′TMdiag(γ)C ′ is invertible as long as at least n soft constraints are on the quadratic segment (i.e., active), which is the necessary condition to fix a particular optimal point n inRn. With such solution, we can train our prediction model with stochastic gradient descent (SGD). For each batch of problem instances, we first solve optimization problems numerically using solvers like Gurobi to get the matrix M and U , and then calculate gradients with the analytical solution. The parameters of the prediction model are updated by such gradients. The sketch of our algorithm is outlined in Appendix C.

heading: 4 Applications and Experiments
text: We apply and evaluate our approach on the three problems described in Section 2, i.e.,linear programming, quadratic programming, and asymmetric soft constraint minimization. These problems are closely related to three applications respectively: synthetic linear programming, portfolio optimization, and resource provisioning, which are constructed using synthetic or real-world datasets. The detailed derivation of gradients for each application can be found in Appendix D.
In our experiments, the performance is measured in regret, which is the difference between the objective value when solving optimization over predicted parameters and the objective value when solving optimization over actual parameters. For each experiment, we choose two generic two-stage methods with L1-loss and L2-loss, as well as decision-focused methods for comparison baselines. We choose both SPO+[10] and DF proposed by Wilder et al. [9] for synthetic linear programming and DF only for portfolio optimization,4 as the former is specially designed for for linear objective. For resource provisioning, we use a problem-specific weighted L1 loss, as both SPO+ and DF are not designed for gradients with respect to variables in the soft constraints. All reported results for each method are obtained by averaging on 15 independent runs with different random seeds.
As real-world data is more lenient than extreme cases, in practice we use a much lower empirical bound than the upper bound proved in section 3.1., e.g., constants of around 20 and 5 √ n where n is the number of dimensions of decision variables. One rule of thumb is to start from a reasonable constant or a constant times √ n, where such "reasonable constant" is the product of a constant factor (e.g. 5− 10) and a roughly estimated upper bound of ||θ||2 (which corresponds to E in our bounds) with specific problem settings; then alternately increase the constant and time an extra √ n while resetting the constant until the program stops diverging, and the hard constraints are satisfied. In our experiments, we hardly observe situations where such process goes for two or more steps.

heading: 4.1 Synthetic Linear Programming
text: Problem setup. The prediction dataset {(ξi, θi)}Ni=1 is generated by a general structural causal model ([21]), ensuring it is representative and consistent with physical process in nature. The programming parameters are generated for various scales in numbers of decision variables, hard constraints, and soft constraints. Full details are given in Appendix E.
4The method [3] and its variant with dimension reduction [9] have same performance in terms of regret on this problem, thus we choose the former for comparison convenience.
Experimental setup. All five methods use the same prediction model – a fully connected neural network of two hidden layers with 128 neurons for each and ReLU [22] for activation. We use AdaGrad [23] as the optimizer, with learning rate 0.01 and gradient clipped at 1e− 4. We train each method for 40 epochs, and early stop when valid performance degrades for 4 consecutive epochs. Specially, to make DF[9] work right on the non-differentiable soft constraints, we first use black-box solvers to determine whether each soft constraint is active on the optimal point, and then optimize with its local expression (i.e. 2nd-order Taylor expansion at optimal point).
Performance analysis. Our experiments cover four programming problem scales with three prediction dataset sizes. Results are summarized in Table 1. In all cases, our method performs consistently better than two-stage methods, DF and SPO+. Even for the cases with only hard constraints (i.e., the third parameter of problem size is 0), our method still has significant advantage, demonstrating its effectiveness on handling hard constraints. Surprisingly, although the main objective is linear, SPO+ often performs even worse than two-stage methods. Detailed analysis (see the appendix) shows that SPO+ quickly reaches the test optimality and then over-fits. This may be due to that, unlike our method, SPO+ loss is not designed to align with the soft constraints. This unawareness of soft constraint is also why DF is performing worse than our method, as DF is working on an optimization landscape that is non-differentiable at the boundary of soft constraints, on which the optimal point usually lies. Besides, with the increment of the samples in train data, the performance of all methods is improved significantly and the performance gap among ours and two-stage methods becomes narrow, which implies that prediction of two-stage methods becomes better and with lower biases. Even so, our method has better sample efficiency than two-stage methods.
We also investigated effect of the hyper-parameter K in our surrogate max function, detailed in the appendix. Through our experiments, K’s effect to regret is not monotonic, and its optimal value varies for different problem sizes. Interestingly,K’s effect is approximately smooth. Thus, in practice, we use simple grid search to efficiently find the best setting of K.

heading: 4.2 Portfolio Optimization
text: Problem and experimental setup. The prediction dataset is daily price data of SP500 from 2004 to 2017 downloaded by Quandl API [24] with the same settings in [3]. We use the same fix as that in linear programming experiment to make DF[9] work with non-differentiable soft constraints, which was also used in [3] for non-convex optimization applications. Most settings are aligned with those in [3], including dataset configuration, prediction model, learning rate (0.01), optimizer (Adam), gradient clip (0.01) and number of training epochs (20). We set the number of soft constraints to 0.4 times of n, where n is the number of candidate equities. For the soft constraint αT max(Cx− d, 0), α = 15n v, where each element of v is generated randomly at uniform from (0, 1); the elements of matrix C is generated independently from {0, 1}, where the probability of 0 is 0.9 and 1 is 0.1. K is set as 100.
Performance analysis. Table 2 summarizes the experimental results. In total, on all problem sizes (#equities), our method performs consistently and significantly better than both two-stage (L1 and L2) methods and the decision focused DF[9]. Among the three baselines, DF is significantly better than two-stage methods, similar to results in [3]. In fact, DF under this setting can be viewed as a variant of our method with infinite K and no conversion of softening ∑ i xi = 1. The comparison to DF also demonstrates the advantage of our method on processing such non-differentiable cases against the simple 2nd-order Taylor expansion. Besides, with the increment of the number of equities, regrets of all methods decease, which indicates that for the constraint ∑ i xi = 1, larger number of equities brings smaller entries of x on average (with the presence of Q, there are many non-zero entries of x), lowering the impact of prediction error for any single entry.

heading: 4.3 Resource Provisioning
text: Problem setup. We use ERCOT energy dataset [25], which contains hourly data of energy output from 2013 to 2018, 52535 data points in total. We use the last 20% samples for test. We aim to predict the matrix C ∈ R24×8, the loads of 24 hours in 8 regions. The decision vairable x is 8-dimensional, and d = 0.5× 1 + 0.1N(0, 1). We test five sets of (α1, α2), with α1/α2 ranging from 100 to 0.01. Experimental setup. We use AdaGrad with learning rate of 0.01, and clip the gradient with norm 0.01. For the prediction model, we use McElWee’s network [26] which was highly optimized for this task, with (8× 24× 77)-dimensional numerical features and embedding ones as input. Performance analysis. Table 3 shows the experimental results. The absolute value of regret differs largely across different ratios of α1/α2. Our method is better than other methods, except for α1/α2 = 1, where the desired objective is exactly L1-loss and thus L1 performs the best. Interestingly, compared to L1/L2, the Weighted L1 performs better when α1/α2 > 1, but relatively worse otherwise. This is probably due to the dataset’s inherent sample bias (e.g., asymmetric distribution and occasional peaks), which causes the systematic bias (usually underestimation) of the prediction model. This bias exacerbates, when superposed with weighted penalty multipliers which encourage the existing bias to fall on the wrong side. Besides, the large variance for weighted L1 at α1/α2 = 0.01 is caused by a pair of outliers.

heading: 5 Related Work
text: Differentiating argmin/argmax through optimality conditions. For convex optimization problems, the KKT conditions map coefficients to the set of solutions, and thus can be differentiated for argmin using implicit function theorem. Following this idea, existing work developed implicit layers of argmin in neural network, including OptNet [13] for quadratic programs (QP) problems and CVXPY [14] for more general convex optimization problems. Further with linear relaxation and
QP regularization, Wilder et al. derived an end-to-end framework for combinatorial programs [9], which accelerates the computation by leverage the low-rank properties of decision vectors [3], and is further extended to mixed integer linear programs in MIPaaL [27]. Besides, for the relaxed LP problems, instead of differentiating KKT conditions, IntOpt [28] proposes an interior point based approach which computes gradients by differentiating homogeneous self-dual formulation.
Optimizing surrogate loss functions. Elmachtoub and Grigas [10] proposed a convex surrogate loss function, namely SPO+, measuring the decision error induced by a prediction, which can handle polyhedral, convex and mixed integer programs with linear objectives. TOPNet [29] proposes a learned surrogate approach for exotic forms of decision loss functions, which however is hard to generalize to handle constrained programs.
Differentiating argmin is critical for gradient methods to optimize decision-focused prediction models. Many kinds of efforts, including regularization for specific problems (e.g., differentiable dynamic programming [30], differentiable submodular optimization [31]), reparameterization [32], perturbation [33] and direct optimization ([12]), are spent for optimization with discrete or continuous variables and still actively investigated.
As a comparison, our work proposes a surrogate loss function for constrained linear and quadratic programs extended with soft constraints, where the soft constraints were not considered principally in previous work. Also, unlike OptNet [13] and CVXPY [15], our method does not need to solve KKT conditions. Instead, by adding reasonable costs for infeasibility, we convert the constrained optimization problem into an unconstrained one while keeping the same solution set, and then derive the required Jacobian analytically. To some degree, we reinvent the exact function [34] in prediction+optimization.

heading: 6 Conclusion
text: In this paper, we have proposed a novel surrogate objective framework for prediction+optimization problems on linear and semi-definite negative quadratic objectives with linear soft constraints. The framework gives the theoretical bounds on constraints’ multipliers, and derives the closed-form solution as well as their gradients with respect to problem parameters required to predict by a model. We first convert the hard constraints into soft ones with reasonably large constant multipliers, making the problem unconstrained, and then optimize a piecewise surrogate locally. We apply and empirically validate our method on three applications against traditional two-stage methods and other end-to-end decision-focused methods. We believe our work is an important enhancement to the current prediction+optimization toolbox. There are two directions for the future work: one is to seek solutions which can deal with hard constraint parameters with negative entries, and the other is to develop the distributional prediction model rather than existing point estimation, to improve robustness of current prediction+optimization methods.

Review summaries:
This paper provides a method for addressing the "prediction + programming" problem in the case where the programming problem contains a $\max(z,0)$ term in the objective. To do so, the authors rewrite the optimization problem as an unconstrained optimization problem by converting all hard constraints of the optimization problem to soft constraints in the objective. They then pick the coefficients of those soft constraints in a way that is meant to prevent constraint violations, and then relax these soft constraints via a surrogate function. They then embed the resultant approximate optimization problem within a method where some $\theta$ is output by a neural network, the approximate optimization problem is solved for that value of $\theta$, and gradients are updated by differentiating end-to-end through the approximate optimization problem. The authors demonstrate their method on three synthetic settings, and report improved performance over two-stage methods and SPO+.

This paper presents an approach for "predict, then optimize" problems that is amenable to derivative-based learning methods. It involves converting hard constraints into soft constraints (i.e. implicitly enforced through objective penalization). The main theoretical contributions are bounds on the magnitude of objective change that can be attributed to infeasibility with respect to the original hard constraints.

The paper develops a surrogate function with closed-form solution for objective functions with "soft constraints," here defined as those with the term $\alpha^T \\max\\{0,z(x)\\}$ for affine $z$, component-wise $\\max$, and $\alpha > 0$. The goal is to apply such functions within the context of smart "predict-then-optimize" (SPO) by also relaxing the hard constraints with a sufficiently large penalty $\alpha$. To develop their surrogate, the authors evaluate bounds on the utility gained when violating the soft constraints using primarly a geometric perspective, specifically evaluating the angles or corners of infeasible points with respect to the convex hull of feasible solutions (and related properties). Numerical results compare their approach with the standard SPO+ loss function and a DF metric for portfolio optimization.

The paper proposes a new method to integrate soft constraints in the predict+programming paradigm that focuses on mathematical programming models with a prediction component that needs to be considered during optimization/solving. The new method works by reformulating the problem into piecewise linear constraints and then optimizing a surrogate for the original problem. A theoretical analysis and motivation is provided as well as an experimental evaluation that underlines the effectiveness of the method.