# A Surrogate Objective Framework for Prediction+Optimization with Soft
Constraints

Kai Yan
Department of Computer Science
University of Illinois at Urbana-Champaign
Urbana, IL 61801
<EMAIL_ADDRESS>
&Jie Yan
Microsoft Research
Beijing, China
<EMAIL_ADDRESS>
&Chuan Luo
Microsoft Research
Beijing, china
<EMAIL_ADDRESS>
&Liting Chen∗
Microsoft Research
Beijing, China
<EMAIL_ADDRESS>
&Qingwei Lin†
Microsoft Research
Beijing, China
<EMAIL_ADDRESS>
&Dongmei Zhang
Microsoft Research
Beijing, China
<EMAIL_ADDRESS>Contributed during internship at Microsoft Research.
$\dagger$ Corresponding Author.

###### 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.

##  1 Introduction

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].

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\in\mathbb{R}^{n}$ and context
variables $C\in\mathbb{R}^{m\times n},d\in\mathbb{R}^{m}$. 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(\cdot,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(\cdot,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.

##  2 Preliminaries

###  2.1 Real-world Optimization Problems with Soft Constraints

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 under-
provisioning, 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\leq 0$ is allowed to violate as long
as the objective improves.

The general problem formulation is

$\small\max_{x}g(\theta,x)-\alpha^{T}\max(z,0),\ z=Cx-d,\ \text{s.t. }Ax\leq
b,Bx=c,\ x\geq 0$ (1)

where $g$ is a utility function, $x\in\mathbb{R}^{n}$ the decision variable,
and $\theta\in\mathbb{R}^{n}$ 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\in\mathbb{R}^{m_{1}\times
n}\geq 0$, $b\in\mathbb{R}^{m_{1}}\geq 0$, $B\in\mathbb{R}^{m_{2}\times n}\geq
0$ ,and $c\in\mathbb{R}^{m_{2}}\geq 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\in\mathbb{R}^{m_{3}}$, $C\in\mathbb{R}^{m_{3}\times
n}$, and $d\in\mathbb{R}^{m_{3}}$. 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,\theta)=\theta^{T}x$. The
problem formulation is

$\small\max_{x}\theta^{T}x-\alpha^{T}\max(Cx-d,0),\text{s.t. }Ax\leq
b,Bx=c,x\geq 0$ (2)

where $\alpha\geq 0$. Consider the application of logistics where $\theta$
represents goods’ prices, {$A$, $B$, $C$} are the capability of transportation
tools, and {$b$, $c$} are the thresholds. Obviously, $A,b,B,c\geq 0$ hold.

Quadratic programming with soft constraints, where
$g(\theta,x)=\theta^{T}x-x^{T}Qx$. One example is the classic minimum variance
portfolio problem [16] with semi-definite positive covariance matrix $Q$ and
expected return $\theta$ 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 $\alpha\geq 0$, we have:

$\displaystyle\max_{x}\theta^{T}x-x^{T}Qx-\alpha^{T}\max(Cx-d,0),\text{s.t.
}Bx=c,~{}x\geq 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 $\alpha_{1},\alpha_{2}>0$, we
have:

$\displaystyle\max_{x}-(\alpha_{1}^{T}\max(Cx-d,0)+\alpha_{2}^{T}\max(d-Cx,0)),\text{
s.t. }Bx=c,~{}x\geq 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$.

###  2.2 Prediction+Optimization

For compactness we write Eq.1 as $\max_{x\in\mathcal{X}}f(x,\theta)$, where
$f$ is the objective function and $\mathcal{X}$ is the feasible domain. The
solver for $f$ is to solve $x^{*}=\text{argmax}_{x\in\mathcal{X}}f(x,\theta)$.
With parameters $\theta$ known, Eq.2–LABEL:eq:prob3 can be solved by mature
solvers like Gurobi [6] and CPLEX [7].

In prediction+optimization, $\theta$ is unknown and needs to be predicted from
some observed features $\xi\in\Xi$. The prediction model is trained on a
dataset $D=\\{(\xi_{i},\theta_{i})\\}_{i=1}^{N}$. In this paper, we consider
the prediction model, namely $\Phi$, as a neural network parameterized with
$\psi$. In traditional supervised learning, $\Phi_{\psi}$ 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].

Figure 1: Computation graph of the decision-focused prediction methods.

Instead, decision-focused learning directly optimizes $\psi$ with respect to
the optimization objective $f$, that is,
$\max_{\psi}{\mathbf{E}}_{(\xi,\theta)\sim
D}[f(\hat{x^{*}}(\Phi_{\psi}(\xi)),\theta)]$. The full computational flow is
illustrated in Fig.1. In the gradient-based learning, update of $\psi$’s
gradient is $\frac{\partial
r}{\partial\psi}=\frac{\partial\hat{\theta}}{\partial\psi}~{}\frac{\partial\hat{x^{*}}}{\partial\hat{\theta}}~{}\frac{\partial
r}{\partial\hat{x^{*}}}$, where utility $r=f(\hat{x^{*}},\theta)$. The
Jacobian ${\partial\hat{\theta}}/{\partial\psi}$ is computed implicitly by
auto-differentiation of deep learning frameworks (e.g, PyTorch [17]), and
${\partial r}/{\partial\hat{x^{*}}}$ is analytical. The main challenge is to
compute ${\partial\hat{x^{*}}}/{\partial\hat{\theta}}$, which depends on
differentiating the $\operatorname*{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 $\hat{\theta}$ to the solution $\hat{x^{*}}$, and
then apply implicit function theorem to obtain
${\partial\hat{x^{*}}}/{\partial\hat{\theta}}$. 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.

##  3 Methodology

Our main idea is to derive a surrogate function for $f$ with a closed-form
solution such that the Jacobian $\frac{\partial x}{\partial\theta}$ 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
$\beta_{1},\beta_{2},\beta_{3}$, such that Eq.1 can be equivalently
transformed into an unconstrained problem:

$\small\max_{x}L(x)=\max_{x}g(x,\theta)-\alpha^{T}\max(Cx-d,0)-\beta_{1}^{T}\max(Ax-b,0)-\beta_{2}^{T}|Bx-c|-\beta_{3}^{T}\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
$\beta_{1},\beta_{2},\beta_{3}$; for this paper, we will assign each entry of
the three vectors the equal value (with a slight abuse of notation, we denote
$\beta_{1}=\beta_{2}=\beta_{3}=\beta$ 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(\cdot,0)$, such that the analytical form of
$\frac{\partial x}{\partial{\theta}}$ 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 $\frac{\partial
x}{\partial\theta}$, we develop our end-to-end learning algorithm of
prediction+optimization whose detailed procedure is described in Appendix C.

###  3.1 Softening the Hard Constraints

For any hard constraints $w=Ax-b\leq 0$, we denote its equivalent soft
constraints as $H(w)=\beta^{T}\max(w,0)$. $H(w)$ should satisfy two
conditions: 1) for $w\leq 0$ (feasible $x$), $H(w)=0$; 2) for $w\geq 0$
(infeasible $x$), $H(w)$ is larger than the utility gain(i.e., improvement of
the objective value) $R=f(x,\theta)-\max_{x_{1}:Ax_{1}\leq b}f(x_{1},\theta)$
by violating $Ax-b\leq 0$. Intuitively, the second condition requires a
sufficiently large-valued $\beta>0$ to ensure that the optimization on the
unconstrained surrogate objective never violates the original $Ax\leq 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.

(a) 2D acute angle

(b) 3D acute angle

(c) Non-degraded cone

(d) Degraded cone

Figure 2: (a) and (b) are $2$ and $3$-dimensional "acute angles"; (c) and (d)
shows two corresponding activated cones for given acute angles. The green area
is the feasible region, $x$ is the red point and $x_{0}$ is the yellow point;
the red $\theta$ is the derivative of an objective $g(x,\theta)=\theta^{T}x$.

The main effort of this subsection is to analyze and bound the effect caused
by such acute angles. Given a convex hull
$\mathcal{C}=\\{z\in\mathbb{R}^{n}|Az\leq b\\}$ ($A\geq 0$ is not required
here) and any point $x\not\in\mathcal{C}$, let $x_{0}\in\mathcal{C}$ be the
nearest point of $\mathcal{C}$ to $x$, and $A^{\prime}x\geq b^{\prime}$
represent all active constraints at $x$, then all such active constraints must
pass through $x_{0}$.111Otherwise, there must exist an active constraint $i$
and a small neighbourhood of $x_{0}$, say $B(x_{0},\epsilon)$, such that
$\forall z\in B(x_{0},\epsilon)$, $A_{i}z-b_{i}\neq 0$, which implies that
either constraint $i$ is inactive ($<0$) or $x_{0}\not\in\mathcal{C}$ ($>0$).
The domain $\mathcal{K}=\\{z\in\mathbb{R}^{n}|A^{\prime}z\geq b^{\prime}\\}$
is a cone or degraded cone where the tip of the cone is a subspace of $R^{n}$.
For the rest of the paper, we will call $\mathcal{K}$ activated cone, as shown
in Fig. 2. Note that for any degraded activated cone, $x-x_{0}$ 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^{\prime}x\leq b^{\prime}$ with the penalty
$\beta^{T}(A^{\prime}x-b^{\prime})$, _i.e.,_ ensure
$\beta^{T}(A^{\prime}x-b^{\prime})\geq R$ hold for any $x$. For the
compactness of symbols and simplicity, we will assume that the main objective
is $g(x,\theta)=\theta^{T}x$ 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,\theta)-\max_{x_{1}\in\mathcal{C}}f(x_{1},\theta)$ be the utility gain,
then $R\leq f(x,\theta)-f(x_{0},\theta)\leq\frac{E}{\cos
p_{0}}\sum_{i=1}^{n}(A^{\prime}_{i}x-b^{\prime}_{i})$, where $x$ is an
infeasible point, $A^{\prime}x\leq b^{\prime}$ the active constraints at $x$,
$p_{0}=\angle({A^{\prime}_{i}}^{*},\theta^{*})$ where ${A^{\prime}_{i}}^{*}$
and $\theta^{*}$ are the optimal solution of
$\max_{\theta}\min_{A^{\prime}_{i}}cos\angle(A^{\prime}_{i},\theta)$ (i.e.,
the maximin angle between $\theta$ and any hyperplane of the activated cone
$\mathcal{K}=\\{z\in\mathbb{R}^{n}|z-x_{0}\in cone(A^{\prime})\\}$), and
$x_{0}$ the projection of $x$ to the tip of cone
$\mathcal{C}=\\{z\in\mathbb{R}^{n}|A^{\prime}z\leq b^{\prime}\\}$. $E$ is the
upper bound of $||\theta||_{2}$. $\angle(\cdot,\cdot)$ denotes the angle of
two vectors.

Thus, $\frac{E}{\cos p_{0}}\mathbf{1}$ is a feasible choice of $\beta$, and it
suffices by finding the lower bound for $\cos p_{0}$. For the rest of the
section, we find the lower bound of $\cos p_{0}$ by exploiting the assumed
properties of $A^{\prime}$, _e.g.,_ $A^{\prime}\geq 0$; we give the
explanation of the full proofs for all theoretical results in Appendix A.

####  3.1.1 Conversion of Inequality Constraints $Ax\leq b$

Let us first consider the constraints $w=Ax-b\leq 0$, where $A\geq 0$, $b\geq
0$. It is easy to prove that given $A\geq 0$ and $b\geq 0$, the distance of a
point to the convex hull $Ax\leq b$ is bounded. More rigorously, we have the
following theorem, which guarantees the feasibility of softening the hard
constraints of inequalities $Ax\leq b$:

######  Theorem 2.

Assume the optimization objective $\theta^{T}x$ with constraints $Ax\leq b$,
where $A\geq 0$, and $b\geq 0$. Then, the utility gain $R$ obtained from
violating $Ax-b\leq 0$ has an upper bound of $O(\sum_{i}\max(w_{i},0)E)$,
where $w=A^{\prime}x-b^{\prime}$, and $A^{\prime}x\leq b^{\prime}$ is the
active constraints.

####  3.1.2 Conversion of Inequality Constraints $x\geq 0$

With inequality constraints $Ax\leq b$ converted, we now further enforce
$x\geq 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\geq 0$ into soft constraints by proving the
following theorem:

######  Theorem 3.

When there is at least one entry of $x\geq 0$ in the activated cone, the
utility gain $R$ from deviating the feasible region is bounded by
$O(\frac{n^{1.5}E\sum_{i}max(w_{i},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 $\beta=O(\frac{n^{1.5}E}{\sin p})$. Specially, for binary
constraints we may set $\beta=O(n^{1.5}E)$:

######  Corollary 4.

For binary constraints where the entries of $A$ are either $0$ or $1$, the
utility gain $R$ of violating $x\geq 0$ constraint is bounded by
$O(n^{1.5}E\sum_{i}\max(w_{i},0))$, where $w_{i}=A_{i}x-b_{i}$ 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).

####  3.1.3 Conversion of Equality Constraints $Bx=c$

Finally, we convert $Bx=c$ into soft constraints. This is particularly
difficult, as $Bx=c$ implies both $Bx\leq c$ and $-Bx\leq-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
$I^{n\times n}$.

######  Theorem 5.

If there is only one equality constraint $Bx=c$ (i.e., $B$ degrades as a row
vector, such like $\sum_{i}x_{i}=1$) and special inequality constraints $x\geq
0$, $Ix\leq b$, then the utility gain $R$ from violating constraints is
bounded by $O(\frac{n^{1.5}E\sum_{i}\max(w_{i},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(\frac{n^{1.5}E}{\sin p})$ for $\beta$.

Finally, we give bounds for general $Ax\leq b$, $Bx=c$ with $A,b,B,c\geq 0$ as
below:

######  Theorem 6.

Given constraints $Ax\leq b$, $x\geq 0$, and $Bx=c$, where $A,B,b,c\geq 0$,
the utility gain $R$ obtained from violating constraints is bounded by
$O(\sqrt{n}\lambda_{max}\sum_{i}\max(w_{i},0))$, where $\lambda_{max}$ is the
upper bound for eigenvalues of $P^{T}P$ ($P:x\rightarrow 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\leq b$, $x\geq
0$, $Bx\leq c$ and $-Bx\leq-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
$\binom{n+m_{1}+m_{2}}{n}$ different cases of $P$, and $\lambda_{max}$ is the
upper bound of eigenvalues of $P^{T}P$ over all possible cases of $P$. Note
that there are no direct bounds on $\lambda_{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,
$\lambda_{max}=O(n^{2})$ follows. Therefore, empirically we can use a bound
$O(\frac{n^{2.5}E}{\sin p})$ for $\beta$.222For real optimization problems, we
can sample $A^{\prime}$ from $A$ and estimate the largest eigenvalue of
$P^{T}P$. 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:

$\displaystyle L(x)=g(x,\theta)-\gamma^{T}\cdot
max(C^{\prime}x-d^{\prime}),\text{ where }~{}\gamma={\begin{bmatrix}\alpha\\\
O(\frac{n^{2.5}E}{\sin{p}})\mathbf{1}\\\
\end{bmatrix}},~{}C^{\prime}={\scriptsize\begin{bmatrix}C\\\ A\\\ -B\\\ B\\\
-I\end{bmatrix}},~{}d^{\prime}={\scriptsize\begin{bmatrix}d\\\ b\\\ -c\\\ c\\\
0\end{bmatrix}}$ (6)

###  3.2 The Unconstrained Soft Constraints

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^{\prime}x-d^{\prime}$ 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:

$\small S(z)=\begin{cases}0\quad&\text{if}\,z<-\frac{1}{4K}\\\
K(z+\frac{1}{4K})^{2}\quad&\text{if}\,-\frac{1}{4K}\leq z\leq\frac{1}{4K}\\\
z\quad&\text{if}\,z\geq\frac{1}{4K}\\\ \end{cases}$ (7)

Let $M$ and $U$ be diagonal matrices as the indicator of $S(z)$, which are
$\small M_{i,i}=2K[-\frac{1}{4K}\leq z_{i}\leq\frac{1}{4K}]$ and
$U_{i,i}=[\frac{1}{4K}<z_{i}]$, where $[\cdot]$ 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, $\frac{\partial
f}{\partial\hat{\theta}}|_{\hat{x}}$ might be too steep. Then consider the
optimal point for the unconstrained optimization problem maximizing
$\theta^{T}x-\gamma^{T}\max(C^{\prime}x-d^{\prime},0)$, by differentiaing on
both sides, we can obtain:

$\small\theta=C^{\prime
T}M\text{diag}(\gamma)(C^{\prime}x-d^{\prime})+C^{\prime
T}(\frac{1}{4K}M+U)\gamma$ (8)

This equation reveals the closed-form solution of $x$ with respect to
$\theta$, and thus the calculation of $\frac{\partial x}{\partial\theta}$
becomes straightforward:

$\frac{\partial x}{\partial\theta}=(C^{\prime
T}M\text{diag}(\gamma)C^{\prime})^{-1}$ (9)

$C^{\prime T}M\text{diag}(\gamma)C^{\prime}$ 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$ in
$\mathcal{R}^{n}$. 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.

##  4 Applications and Experiments

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,333The 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. 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\sqrt{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 $\sqrt{n}$, where such "reasonable constant" is
the product of a constant factor (e.g. $5-10$) and a roughly estimated upper
bound of $||\theta||_{2}$ (which corresponds to $E$ in our bounds) with
specific problem settings; then alternately increase the constant and time an
extra $\sqrt{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.

###  4.1 Synthetic Linear Programming

Problem setup. The prediction dataset $\\{(\xi_{i},\theta_{i})\\}_{i=1}^{N}$
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.

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$.

|  | Regret (the lower, the better)
---|---|---
$N$ | Problem Size | L1 | L2 | SPO+ [10] | DF [9] | Ours
100 | (40, 40, 0) | 2.454$\pm$0.232 | 2.493$\pm$0.295 | 2.506$\pm$0.294 | 2.478$\pm$0.425 | 2.258$\pm$0.311
| (40, 40, 20) | 2.626$\pm$0.307 | 2.664$\pm$0.303 | 2.667$\pm$0.281 | 2.536$\pm$0.376 | 2.350$\pm$0.263
| (80, 80, 0) | 5.736$\pm$0.291 | 5.831$\pm$0.361 | 5.711$\pm$0.309 | 5.756$\pm$0.317 | 5.200$\pm$0.506
| (80, 80, 40) | 4.786$\pm$0.403 | 4.786$\pm$0.596 | 4.939$\pm$0.382 | 4.902$\pm$0.537 | 4.570$\pm$0.390
1000 | (40, 40, 0) | 1.463$\pm$0.143 | 1.447$\pm$0.155 | 1.454$\pm$0.148 | 1.434$\pm$0.268 | 1.346$\pm$0.144
| (40, 40, 20) | 1.626$\pm$0.141 | 1.613$\pm$0.110 | 1.618$\pm$0.103 | 1.529$\pm$0.151 | 1.506$\pm$0.102
| (80, 80, 0) | 3.768$\pm$0.132 | 3.718$\pm$0.117 | 3.573$\pm$0.113 | 3.532$\pm$0.102 | 3.431$\pm$0.100
| (80, 80, 40) | 2.982$\pm$0.176 | 2.913$\pm$0.172 | 2.879$\pm$0.148 | 3.351$\pm$0.212 | 2.781$\pm$0.165
5000 | (40, 40, 0) | 1.077$\pm$0.105 | 1.080$\pm$0.109 | 1.090$\pm$0.105 | 1.078$\pm$0.092 | 1.037$\pm$0.100
| (40, 40, 20) | 1.283$\pm$0.070 | 1.277$\pm$0.077 | 1.298$\pm$0.077 | 1.291$\pm$0.091 | 1.220$\pm$0.071
| (80, 80, 0) | 2.959$\pm$0.086 | 2.943$\pm$0.091 | 2.926$\pm$0.079 | 2.869$\pm$0.085 | 2.845$\pm$0.064
| (80, 80, 40) | 2.239$\pm$0.122 | 2.224$\pm$0.106 | 2.234$\pm$0.122 | 2.748$\pm$0.165 | 2.172$\pm$0.098
Table 1: Performance comparison (regret mean with std. deviation) for the
synthetic linear programming problem. $N$ is the size of the training dataset,
and problem size is a triplet (# of decision variables’ dimension, # of hard
constraints, # of soft constraints).

###  4.2 Portfolio Optimization

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
$\alpha^{T}\max(Cx-d,0)$, $\alpha=\frac{15}{n}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 $\sum_{i}x_{i}=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 $\sum_{i}x_{i}=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.

| Regret measured in % (the lower, the better)
---|---
#Equities | L1 | L2 | DF [9] | ours($K=100$)
50 | 4.426$\pm$0.386 | 4.472$\pm$0.385 | 4.016$\pm$0.389 | 3.662$\pm$0.238
100 | 4.262$\pm$0.231 | 4.320$\pm$0.229 | 3.500$\pm$0.252 | 3.214$\pm$0.138
150 | 3.878$\pm$0.281 | 3.950$\pm$0.287 | 3.419$\pm$0.281 | 3.109$\pm$0.162
200 | 3.755$\pm$0.236 | 3.822$\pm$0.273 | 3.406$\pm$0.287 | 3.152$\pm$0.183
250 | 3.721$\pm$0.205 | 3.751$\pm$0.212 | 3.335$\pm$0.175 | 3.212$\pm$0.135
Table 2: Performance comparison (regret mean with std. deviation) for
portfolio optimization.

###  4.3 Resource Provisioning

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\in\mathbb{R}^{24\times 8}$, the loads of $24$ hours in $8$ regions. The
decision vairable $x$ is $8$-dimensional, and
$d=0.5\times\mathbf{1}+0.1N(0,1)$. We test five sets of
$(\alpha_{1},\alpha_{2})$, with $\alpha_{1}/\alpha_{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\times 24\times
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
$\alpha_{1}/\alpha_{2}$. Our method is better than other methods, except for
$\alpha_{1}/\alpha_{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 $\alpha_{1}/\alpha_{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 $\alpha_{1}/\alpha_{2}=0.01$ is caused by a pair of outliers.

| Regret (the lower, the better)
---|---
$\alpha_{1}/\alpha_{2}$ | L1 | L2 | Weighted L1 | Ours($K=0.05$)
100 | 105.061$\pm$21.954 | 93.193$\pm$29.815 | 79.014$\pm$32.069 | 20.829$\pm$8.289
10 | 13.061$\pm$2.713 | 13.275$\pm$6.208 | 7.743$\pm$1.305 | 2.746$\pm$1.296
1 | 4.267$\pm$0.618 | 5.136$\pm$0.722 | 4.267$\pm$0.618 | 5.839$\pm$0.512
0.1 | 10.846$\pm$1.606 | 13.619$\pm$2.195 | 16.462$\pm$2.093 | 10.240$\pm$1.248
0.01 | 99.145$\pm$21.159 | 118.112$\pm$29.957 | 230.825$\pm$91.184 | 94.341$\pm$29.821
Table 3: Performance comparison (regret mean with std. deviation) for resource
provisioning.

##  5 Related Work

Differentiating $\bm{\operatorname*{argmin}}$/$\bm{\operatorname*{argmax}}$
through optimality conditions. For convex optimization problems, the KKT
conditions map coefficients to the set of solutions, and thus can be
differentiated for $\operatorname*{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 $\operatorname*{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.

##  6 Conclusion

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.

## References

  * Luo et al. [2020] Chuan Luo, Bo Qiao, Xin Chen, Pu Zhao, Randolph Yao, Hongyu Zhang, Wei Wu, Andrew Zhou, and Qingwei Lin.  Intelligent virtual machine provisioning in cloud computing.  In _Proceedings of IJCAI 2020_ , pages 1495–1502, 2020.
  * Chansombat et al. [2018] Sirikarn Chansombat, Pupong Pongcharoen, and Christian Hicks.  A mixed-integer linear programming model for integrated production and preventive maintenance scheduling in the capital goods industry.  _International Journal of Production Research_ , 57:1–22, 04 2018.  doi: 10.1080/00207543.2018.1459923.
  * Wang et al. [2020] Kai Wang, Bryan Wilder, Andrew Perrault, and Milind Tambe.  Automatically learning compact quality-aware surrogates for optimization problems.  In _NeurIPS 2020_ , Vancouver, Canada, 2020.
  * Demirovi et al. [2020] Emir Demirovi, Peter Stuckey, Tias Guns, James Bailey, Christopher Leckie, Kotagiri Ramamohanarao, and Jeffrey Chan.  Dynamic programming for predict+optimise.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 34:1444–1451, 04 2020.  doi: 10.1609/aaai.v34i02.5502.
  * Ganeshapillai et al. [2013] Gartheeban Ganeshapillai, John Guttag, and Andrew Lo.  Learning connections in financial time series.  In Sanjoy Dasgupta and David McAllester, editors, _Proceedings of the 30th International Conference on Machine Learning_ , volume 28 of _Proceedings of Machine Learning Research_ , pages 109–117, Atlanta, Georgia, USA, 17–19 Jun 2013. PMLR.
  * Gurobi Optimization [2021] LLC Gurobi Optimization.  Gurobi optimizer reference manual, 2021.  URL http://www.gurobi.com.
  * Cplex [2020] IBM ILOG Cplex.  Release notes for cplex 20.1.0.  2020\.  URL https://www.ibm.com/docs/en/icos/20.1.0?topic=2010-release-notes-cplex.
  * Bengio [1997] Yoshua Bengio.  Using a financial training criterion rather than a prediction criterion.  _Int. J. Neural Syst._ , 8(4):433–443, 1997.  doi: 10.1142/S0129065797000422.  URL https://doi.org/10.1142/S0129065797000422.
  * Wilder et al. [2019] Bryan Wilder, Bistra Dilkina, and Milind Tambe.  Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization.  In _The Thirty-Third AAAI Conference on Artificial Intelligence_ , pages 1658–1665. AAAI Press, 2019.
  * Elmachtoub and Grigas [2017] Adam Elmachtoub and Paul Grigas.  Smart "predict, then optimize".  _Management Science_ , 10 2017.  doi: 10.1287/mnsc.2020.3922.
  * Donti et al. [2017] Priya L. Donti, Brandon Amos, and J. Zico Kolter.  Task-based end-to-end model learning in stochastic optimization.  In _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , NIPS’17, page 5490–5500, Red Hook, NY, USA, 2017. Curran Associates Inc.  ISBN 9781510860964.
  * Song et al. [2016] Yang Song, Alexander Schwing, Richard, and Raquel Urtasun.  Training deep neural networks via direct loss minimization.  In _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pages 2169–2177, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * Amos and Kolter [2017] Brandon Amos and J. Zico Kolter.  OptNet: Differentiable optimization as a layer in neural networks.  In Doina Precup and Yee Whye Teh, editors, _Proceedings of the 34th International Conference on Machine Learning_ , volume 70 of _Proceedings of Machine Learning Research_ , pages 136–145, International Convention Centre, Sydney, Australia, 06–11 Aug 2017. PMLR.  URL http://proceedings.mlr.press/v70/amos17a.html.
  * Amos [2019] Brandon Amos.  _Differentiable Optimization-Based Modeling for Machine Learning_.  PhD thesis, Carnegie Mellon University, May 2019.
  * Agrawal et al. [2019] Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J. Zico Kolter.  Differentiable convex optimization layers.  In _Advances in Neural Information Processing Systems_ , volume 32, pages 9562–9574. Curran Associates, Inc., 2019.
  * Markowitz [1959] Harry M. Markowitz.  _Portfolio Selection: Efficient Diversification of Investments_.  Yale University Press, 1959.  ISBN 9780300013726.  URL http://www.jstor.org/stable/j.ctt1bh4c8h.
  * Paszke et al. [2019] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Kopf, Edward Yang, Zachary DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32_ , pages 8024–8035. Curran Associates, Inc., 2019.
  * Griewank and Walther [2008] A. Griewank and A. Walther.  _Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation, Second Edition_.  Other Titles in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM, 3600 Market Street, Floor 6, Philadelphia, PA 19104), 2008.  ISBN 9780898717761.
  * Magnus [2019] Jan R. Magnus Jan R. Magnus.  _Matrix Differential Calculus with Applications in Statistics and Econometrics, Third Edition_.  Wiley Series in Probability and Statistics. Wiley, 2019.  ISBN 9781119541202.
  * Boyd and Vandenberghe [2004] Stephen Boyd and Lieven Vandenberghe.  _Convex Optimization_.  Cambridge University Press, 2004.
  * Peters et al. [2017] Jonas Peters, Dominik Janzing, and Bernhard Schölkopf.  _Elements of causal inference_.  The MIT Press, 2017.
  * Nair and Hinton [2010] Vinod Nair and Geoffrey E. Hinton.  Rectified linear units improve restricted boltzmann machines.  In _Proceedings of the 27th International Conference on International Conference on Machine Learning_ , ICML’10, page 807–814, Madison, WI, USA, 2010. Omnipress.  ISBN 9781605589077.
  * Duchi et al. [2011] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  _Journal of Machine Learning Research_ , 12:2121–2159, 07 2011.
  * sp5 [2021] Sp500 data set.  2021\.  URL https://docs.quandl.com/.
  * erc [2021] Ercot data set.  2021\.  URL https://github.com/kmcelwee/mediumBlog/tree/master/load_forecast/data.
  * McElwee [2020] Kevin McElwee.  Predict daily electric consumption with neural networks, 2020.  URL https://www.brownanalytics.com/load-forecasting.
  * Ferber et al. [2020] Aaron Ferber, Bryan Wilder, Bistra Dilkina, and Milind Tambe.  Mipaal: Mixed integer program as a layer.  In _AAAI Conference on Artificial Intelligence_ , 2020.
  * Mandi and Guns [2020] Jayanta Mandi and Tias Guns.  Interior point solving for lp-based prediction+optimisation.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems_ , volume 33, pages 7272–7282. Curran Associates, Inc., 2020.
  * Chen et al. [2020] Di Chen, Yada Zhu, Xiaodong Cui, and Carla Gomes.  Task-based learning via task-oriented prediction network with applications in finance.  In Christian Bessiere, editor, _Proceedings of the Twenty-Ninth International Joint Conference on Artificial Intelligence, IJCAI-20_ , pages 4476–4482. International Joint Conferences on Artificial Intelligence Organization, 7 2020.  doi: 10.24963/ijcai.2020/617.  URL https://doi.org/10.24963/ijcai.2020/617.  Special Track on AI in FinTech.
  * Mensch and Blondel [2018] Arthur Mensch and Mathieu Blondel.  Differentiable dynamic programming for structured prediction and attention.  In _35th International Conference on Machine Learning_ , volume 80, pages 3459–3468, 2018.
  * Djolonga and Krause [2017] Josip Djolonga and Andreas Krause.  Differentiable learning of submodular models.  In _Advances in Neural Information Processing Systems_ , volume 30, pages 1013–1023. Curran Associates, Inc., 2017.
  * Rezende et al. [2014] Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra.  Stochastic backpropagation and approximate inference in deep generative models.  In _Proceedings of the 31st International Conference on International Conference on Machine Learning - Volume 32_ , ICML’14, page II–1278–II–1286. JMLR.org, 2014.
  * Berthet et al. [2020] Quentin Berthet, Mathieu Blondel, Olivier Teboul, Marco Cuturi, Jean-Philippe Vert, and Francis Bach.  Learning with differentiable perturbed optimizers.  _arXiv preprint arXiv:2002.08676_ , 2020.
  * Bertsekas [1997] Dimitri P Bertsekas.  Nonlinear programming.  _Journal of the Operational Research Society_ , 48(3):334–334, 1997.

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [Yes]

    3. (c)

Did you discuss any potential negative societal impacts of your work? [Yes]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [Yes]

    2. (b)

Did you include complete proofs of all theoretical results? [Yes]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [Yes]

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [Yes]

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [Yes]

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [Yes]

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [Yes]

    2. (b)

Did you mention the license of the assets? [Yes]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [Yes]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [Yes]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [Yes]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A] We did not use crowdsourcing, nor did we
conduct research with human subjects.

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]

##  Appendix A Mathematical Proofs

Throughout the following proofs, for the convenience of description, we denote
two symbols:

  * •

the symbol $\angle(a,b)$, which means the angle of two vectors
$a,b\in\mathbb{R}^{n}$;

  * •

$cone(A)$, which means the conic combination of the row vectors of matrix $A$;
$cone(A,B)$ stands for the conic combination of the union of the row vectors
of matrix $A,B$.

Throughout this section, we will assume the same variables as in Eq.1-5 and
Section 3 by default, and further make the following assumptions or
simplification for ease of description.

  * •

Assume the main objective $g(x,\theta)=\theta^{T}x$, where $||\theta||_{2}\leq
E$. In reality, elements of the predictive vector $\theta$ have a known range
and thus $l_{2}$-norm of $\theta$ is bounded by a constant, namely
$E\in\mathbb{R}$.

  * •

We omit the soft constraints in original objective. Since all proofs in this
section consider the worst direction of $\frac{\partial L}{\partial x}$ in
Eq.5, for soft constraints which are linear, their derivatives are constant
and can be integrated into $E$ for any bounded penalty multiplier $\alpha$.

  * •

Similarly, we omit quadratic term $Q$. At first glance, it seems that $Q$ will
bring unbounded derivative of decision variable $x$. However, note that we are
finding the maximum value of function, and the matrix $Q$ is semi-definite
positive, we have $\frac{\partial-x^{T}Qx}{\partial x_{i}}<0$ for any $i$.
Thus, denote the maximum point of $-x^{T}Qx$ to be $x_{0}$, we can find a
finite radius $r>0$, such that any point $x\in\mathcal{R}^{n}$ out of the
circle with $x_{0}$ as the center and $r$ as radius is impossible to be the
global optimal point; intuitively, the norm of the derivative of $Q$ becomes
so large as to "draw the optimal point back" from too far, and thus we can
ignore the points outside the circle. Since the circle is a bounded set, the
norm of derivative $E$ is also bounded inside the circle.

  * •

When we consider the effect of violating constraints, we only consider non-
degraded activated cone, _i.e.,_ , for an activated cone in a $n$-dimensional
Euclidean space, there are at least $n$ non-redundant active constraints.
Otherwise, as we stated in the main paper, we can project the activated cone
$\mathcal{K}=\\{z\in\mathbb{R}^{n}|z-x_{0}\in cone(A^{\prime})\\}$ onto the
complementary space of the cone tip $x_{0}$ (_i.e.,_ the solution space of
$A^{\prime}$) and proceed on a low-dimensional subspace.

  * •

As for constraint itself, we assume that for any row vector of $A$ in $Ax\leq
b$, and any row vector of $B$ in $Bx=c$, we have $||A_{i}||_{2}=1$ and
$||B_{i}||_{2}=1$ for any row vector $A_{i},B_{i}$, $i\in\\{1,2,...,m\\}$
where $m$ is the number of constraints. If not, we can first normalize the row
vectors of $A$ and $B$ while adequately scaling $b,c$, then apply our proofs.

Thus it is straightforward to apply our proofs to the general form of Eq.5.

###  A.1 Lemma 1

######  Lemma 1.

(Bounding the utility gain) Let
$R=f(x,\theta)-\max_{x_{1}\in\mathcal{C}}f(x_{1},\theta)$ be the utility gain,
then $R\leq f(x,\theta)-f(x_{0},\theta)\leq\frac{E}{\cos
p_{0}}\sum_{i=1}^{n}(A^{\prime}_{i}x-b^{\prime}_{i})$, where $x$ is an
infeasible point, $A^{\prime}x\leq b^{\prime}$ the active constraints at $x$,
$p_{0}=\angle({A^{\prime}_{i}}^{*},\theta^{*})$ where ${A^{\prime}_{i}}^{*}$
and $\theta^{*}$ are the optimal solution of
$\max_{\theta}\min_{A^{\prime}_{i}}cos\angle(A^{\prime}_{i},\theta)$ (i.e.,
the maximin angle between $\theta$ and hyperplanes of the activated cone
$\mathcal{K}=\\{z\in\mathbb{R}^{n}|z-x_{0}\in cone(A^{\prime})\\}$), and
$x_{0}$ the projection of $x$ to the tip of cone
$\mathcal{C}=\\{z\in\mathbb{R}^{n}|A^{\prime}z\leq b^{\prime}\\}$. $E$ is the
upper bound of $||\theta||_{2}$.

###### Proof.

The first inequality is trivial as $x_{0}\in\mathcal{C}$ and directly followed
from the definition of utility gain $R$.

For the second inequality, by the law of sines, we have
$||x-x_{0}||_{2}=\frac{d_{i}}{\sin\tau_{i}}$, where $d_{i}$ is the distance of
$x$ to the $i$-th hyper-plane in the activated cone $cone(A^{\prime})$ with
$A^{\prime}_{i}$ (_i.e.,_ the $i$-th row of $A^{\prime}$) being the normal
vector of the $i$-th hyper-plane in $A^{\prime}x\leq b^{\prime}$, and
$\tau_{i}$ is the angle between $i$-th hyper-plane and $\theta$. Then the
utility gain $R$ from violating the constraints follows:

$\displaystyle R$ $\displaystyle\leq\theta^{T}(x-x_{0})$ (10)
$\displaystyle\leq||\theta||_{2}||x-x_{0}||_{2}$
$\displaystyle=||\theta||_{2}\frac{d_{i}}{\text{sin}~{}\tau_{i}}\ (\text{By
the law of sines; $i$ is arbitrary})$
$\displaystyle=||\theta||_{2}\frac{d_{j}}{\text{sin}~{}\tau_{j}}\ (\text{By
selecting
$j=\operatorname*{argmin}_{i}{\text{sin}~{}\tau_{i}}=\operatorname*{argmax}_{i}{cos\angle(A^{\prime}_{i},x-x_{0})}$})$
$\displaystyle\leq E\frac{\sum_{i=1}^{n}d_{i}}{\text{sin}~{}\tau_{j}}\text{ (E
is the upper bound of }||\theta||_{2}\text{)}$
$\displaystyle=E\frac{\sum_{i=1}^{n}(A^{\prime}_{i}x-b^{\prime}_{i})}{\text{sin}~{}\tau_{j}}(\text{$d_{i}=\frac{|A^{\prime}_{i}x-b^{\prime}_{i}|}{||A^{\prime}_{i}||_{2}}=A^{\prime}_{i}x-b^{\prime}_{i}$;
note the fifth assumption.})$

Thus, setting all entries of penalty multiplier $\beta$ in the hard constraint
conversion to
$O(\frac{E}{sin~{}\tau_{j}})=O(\frac{E}{cos\angle(A^{\prime}_{j},x-x_{0})})$
will give us an upper bound of the minimum feasible $\beta$. ∎

The rest of the work is to find the lower bound for
$\min_{\theta}\max_{i}cos\angle(A^{\prime}_{i},x-x_{0})$ for any given
activated cone $cone(A^{\prime})$ and $x_{0}$; as we have assumed
$||A^{\prime}_{i}||=1\ \forall i\in\\{1,2,...,n\\}$, the objective
$\min_{\theta}\max_{i}cos\angle(A^{\prime}_{i},x-x_{0})$ can be equivalently
written as444With an abuse of notation, the $\theta$ in Eq.11 represents
$x-x_{0}$, as the worst situation for deriving upper bound of $\beta$ is that
$x-x_{0}$ has the same direction with the objective $\theta^{T}x$’s derivative
$\theta$.

$F=\min_{\theta}\max_{i}A_{i}^{\prime}\theta$ (11)

Deriving the lower bound of $F$ is the core part of proving Theorem 3 and 6.

###  A.2 Lemma 7

In order to prove Theorem 2, we give a crucial lemma.

######  Lemma 7.

Consider a set of hyper-planes with normal vectors
$A^{\prime}_{1},A^{\prime}_{2},...,A^{\prime}_{n}\in\mathbb{R}^{n}\geq 0$,
$||A^{\prime}_{i}||_{2}=1$ for $i\in\\{1,2,...,n\\}$ which forms a cone. Let
$d_{i}$ be the distance of a point $x$ to the hyper-plane $A^{\prime}_{i}$,
where $x$ is in the cone of $A^{\prime}_{i}$ (_i.e.,_
$x=k_{1}A^{\prime}_{1}+...+k_{n}A^{\prime}_{n}$, $k_{i}\geq 0~{}\forall
i\in\\{1,2,...,n\\}$), then we have $||x||_{2}\leq\sum_{i=1}^{n}d_{i}$.

###### Proof.

Without loss of generality, let
$x=k_{1}A^{\prime}_{1}+k_{2}A^{\prime}_{2}+...+k_{n}A^{\prime}_{n}$, and
$||x||_{2}=1$ (we can scale $x$ if $||x||_{2}\neq 1$). The distance
$d_{i}=\frac{{A^{\prime}_{i}}^{T}x}{||A^{\prime}_{i}||_{2}}={A^{\prime}_{i}}^{T}x$.
Therefore, we have

$\displaystyle\sum_{i=1}^{n}d_{i}$
$\displaystyle=\sum_{i=1}^{n}{A^{\prime}_{i}}^{T}(k_{1}A^{\prime}_{1}+...+k_{n}A^{\prime}_{n})$
(12)
$\displaystyle=\sum_{j=1}^{n}k_{j}({A^{\prime}_{j}}^{T}A^{\prime}_{j}+\sum_{i=1,i\neq
j}^{n}{A^{\prime}_{i}}^{T}A^{\prime}_{j})$
$\displaystyle\geq\sum_{j=1}^{n}k_{j}||A^{\prime}_{j}||_{2}\text{
($A^{\prime}_{i}\geq 0$, therefore ${A^{\prime}_{i}}^{T}A^{\prime}_{j}\geq
0$)}$ $\displaystyle=k_{1}+k_{2}+...+k_{n}$

We next prove that $k_{1}+k_{2}+...+k_{n}\geq 1$. As $||x||_{2}=1$, we have

$\displaystyle||x||_{2}$ $\displaystyle=1$ (13) $\displaystyle
k_{1}{A^{\prime}_{1}}^{T}x+k_{2}{A^{\prime}_{2}}^{T}x+...+k_{n}{A^{\prime}_{n}}^{T}x$
$\displaystyle=1$

As ${A^{\prime}_{i}}^{T}x\in[0,1]$, $k_{1}+...+k_{n}\geq 1$ must hold.
Therefore, $\sum_{i=1}^{n}d_{i}\geq k_{1}+...+k_{n}\geq 1=||x||_{2}$. ∎

###  A.3 Theorem 2

######  Theorem 2.

Assume the optimization objective $\theta^{T}x$ with constraints $Ax\leq b$,
where $A\geq 0$, and $b\geq 0$. Then, the utility gain $R$ obtained from
violating $Ax-b\leq 0$ has an upper bound of $O(\sum_{i}\max(w_{i},0)E)$,
where $w=A^{\prime}x-b^{\prime}$, and $A^{\prime}x\leq b^{\prime}$ is the
active constraints.

###### Proof.

Let $d_{i}=(d_{i,1},d_{i,2},...,d_{i,n})^{T}$ be the vector of distances to
the $i$-th hyper-plane of $A^{\prime}$ with normal vector $A^{\prime}_{i}$
from $x$. By definition of $x$ and $x_{0}$, the distance vector of the point
to the convex hull $\mathcal{C}$ is $x-x_{0}$. When $A\geq 0$ and $b\geq 0$,
obviously Lemma 7 holds. Consider the active constraints
$w=A^{\prime}x-b^{\prime}\geq 0$ which is the subset of active constraints in
$Ax\leq b$. The utility gain $R$ from violating the rules
$w=A^{\prime}x-b^{\prime}$ satisfies:

$\displaystyle R$ $\displaystyle\leq\theta^{T}(x-x_{0})$ (14)
$\displaystyle\leq||\theta||_{2}||(x-x_{0})||_{2}$
$\displaystyle\leq||\theta||_{2}\sum_{i}||d_{i}||_{2}\text{ (Lemma 7)}$
$\displaystyle\leq E\sum_{i}||d_{i}||_{2}$
$\displaystyle=E\sum_{i}\max(w_{i},0)\text{
($d_{i}=\frac{|A^{\prime}_{i}x-b^{\prime}_{i}|}{||A^{\prime}_{i}||_{2}}=A^{\prime}_{i}x-b^{\prime}_{i}$;
note the fifth assumption.)}$

This holds for any $x$. Therefore, we turn the hard constraint $Ax\leq b$ into
a soft constraint with penalty multiplier $\beta=E$. ∎

###  A.4 Theorem 3

######  Theorem 3.

Assume at least one constraint in $x\geq 0$ is active, then the utility gain
$R$ by deviating from the feasible region is bounded by
$O(\frac{n^{1.5}E\sum_{i}max(w_{i},0)}{\sin{p_{0}}})$, where
$p_{0}=\min_{i,j}\angle(A_{i},e_{j})$ (_i.e.,_ the smallest corner between
axes and other constraint hyper-planes), and $w=\begin{bmatrix}A^{\prime}\\\
-I^{\prime}\end{bmatrix}x-\begin{bmatrix}b^{\prime}\\\ 0\end{bmatrix}$, where
$A^{\prime}x\leq b^{\prime}$ and $-I^{\prime}x\leq 0$ are active constraints
in $Ax\leq b$ and $x\geq 0$ respectively.

###### Proof.

As Lemma 1 implies, the key point of proving this theorem is to prove that
$\theta$ (the derivative of the objective), is always close to some rows of
$A^{\prime}$ and $-I^{\prime}$.

Consider the activated cone $\mathcal{K}=\\{z\in\mathbb{R}^{n}|z-x_{0}\in
cone(A^{\prime},-I^{\prime})\\}$; all hyper-planes of $A^{\prime}x=b^{\prime}$
and violated entries of $x\geq 0$ must pass through $\mathcal{K}$’s tip
$x_{0}$. Figure 3 gives such an illustrative example. For the rest of the
proof, it is enough to assume that the activated cone $\mathcal{K}$ is a non-
degraded cone, since otherwise, as the main paper states, we can project the
activated cone onto the supplementary space of its tip. By the projection, we
actually reduce the problem to the same one with lower dimensions of $x$, and
can apply the following proof in the same way. Let $q$ be the number of rows
in $I^{\prime}$, i.e., the number of active constraints in $x\geq 0$; and $r$
be the number of rows in $A^{\prime}$, i.e., the number of active constraints
in $Ax\leq b$. Since $\mathcal{K}$ is non-degraded, we have $r+q\geq n$.

Without loss of generality, we transform the problem equivalently by
rearranging the order of dimensions of $x$ so that $I^{\prime}$ is the first
$q$ rows of $I_{n\times n}$. Let’s see an example after this transformation.
As shown in Figure 3, the three solid black vectors are $(-1,0,0)$,
$(0,-1,0)$, and $A^{\prime}_{1}\geq 0$, where the first and second entries of
$A^{\prime}_{1}$ are strictly greater than $0$ (otherwise the cone is
degraded). Similarly, for the $n$-dimensional non-degraded activated cone, we
consider the maximin angle between $\theta$ and all active inequality
constraints $A^{\prime}$ and $-I^{\prime}$ (which is equivalent to finding the
minimax cosine value of angles between $\theta$ and all active inequality
constraints).

Figure 3: A $3$-dimensional illustration of the proof for Theorem 3. The green
polytope is the feasible region; the black solid vectors are the related
vectors (i.e. rows of $A^{\prime}$) and the red line is the utility vector
$\theta$. Apparently, $(0,0,1)^{T}\theta\geq 0$, otherwise the intersection
point will move down to origin. The worse case appears when the red $\theta$
is on the circumcenter of the triangle formed by the three solid unit vectors.

Then, guided by Equation 11 in Lemma 1, the upper bound of distance between
$\theta$ and hyper-planes with normal vector $A^{\prime}_{j}$, is the solution
of the following optimization problem w.r.t. $\theta\in\mathbb{R}^{n}$ and
$A^{\prime}_{j}$ (since the cone is non-degraded, we assume that
$A^{\prime}_{j}$ are linearly independent):

$\displaystyle F$ $\displaystyle=\min_{\theta}\max_{j\in J=\\{1,..,r\\},k\in
K=\\{1,..,q\\}}\\{A^{\prime}_{j}\theta,-e_{k}^{T}\theta\\},$ (15)
$\displaystyle s.t.~{}||\theta||_{2}$
$\displaystyle=1,\theta\in\mathcal{K},r+q\geq n.$

where $e_{k}$ is the unit vector where the $k$-th entry is $1$ and others are
$0$. Note that according to our presumptions, $A^{\prime}_{j}\geq 0$,
$||A^{\prime}_{j}||_{2}=1$ for $j\in\\{1,2,...,r\\}$ have already satisfied.
The constraints $\theta\in\mathcal{K}$ come from the following consideration:
if with given $\theta$ we update $x$ in an iterative manner such as gradient
ascending, $x$ would eventually leave current activated cone where
$\theta_{i}<0~{}~{}\text{for}~{}~{}i\in\\{q+1,...,n\\}$ (see Figure 3 for an
illustration). Since for any dimension index $i\in\\{q+1,...,n\\}$, all
entries of $A^{\prime}_{j}$ and $-e_{k}$ ($j\in J,k\in K$) are non-negative,
we know that $\theta_{i}\geq 0,~{}\forall i\in\\{q+1,...,n\\}$. To further
relax the objective $F$ and derive the lower bound, we will assume $r+q=n$ in
the rest of the paper; if $r\geq n-q$, we can simply ignore all
$A^{\prime}_{j}$ in the max operator with $j>n-q$.

Given the signs of $\theta$, we next derive a lower bound of $F$ with respect
to any given $A^{\prime}_{j}$ and set of $k$. Now we relax $F$ by setting all
entries of $A^{\prime}_{j}$ to $0$, except for the entries with index
$\\{1,2,...,q\\}$ and $s_{j}$ to get $A^{\prime\prime}_{j}$, where
$\\{s_{j}\\}(j\in\\{q+1,...,n\\})$ is a permutation of $\\{q+1,q+2...,n\\}$.
We have:

$\displaystyle F$
$\displaystyle\geq\min_{\theta}\max_{j\in\\{1,..,r\\},k\in\\{1,..,q\\}}\\{{A^{\prime\prime}}_{j}^{T}\theta,-e_{k}^{T}\theta\\}\text{
(Relaxing within the $\max$ operator; note the sign of $\theta$)}$ (16)
$\displaystyle\geq\min_{\theta}\max_{j\in\\{1,..,r\\},k\in\\{1,..,q\\}}\\{\sum_{i=1}^{q}A^{\prime}_{j,i}\theta_{i}+A^{\prime}_{j,s_{j}}\theta_{s_{j}},-\theta_{k}\\}$
$\displaystyle\geq\min_{\theta}\max_{s_{j}\in\\{q+1,n\\},k\in\\{1,..,q\\}}\\{\sum_{i=1}^{q}\theta_{i}+\alpha_{0}\theta_{s_{j}},-\theta_{k}\\}$

where $||\theta||_{2}=1$, $\alpha_{0}$ is the smallest non-zero entry among
$A^{\prime}_{j,s_{j}}$ (and note that $A\geq 0$ and $\theta_{s_{j}}\geq 0$,
which means $\alpha_{0}\theta_{s_{j}}\geq 0$). The inequalities hold for the
relaxation within the $\max$ operator. Note that the last line of inequality
assumes $\theta_{i}^{*}\leq 0$ for $i\in\\{1,2,...,q\\}$ for the optimal
$\theta^{*}$; otherwise, we can scale further, ignore the index $i$ of the
optimal point $\theta^{*}$ where $\theta^{*}_{i}>0$ and proceed with $n-1$
dimensions for the following two facts:

  1. 1.

$\sum_{i=1}^{q}\theta_{i}+\alpha_{0}\theta_{s_{j}}$ term becomes smaller after
ignoring such dimension;

  2. 2.

$\theta\in\mathcal{K}$, which means for any $\theta$ and any vector $y$ of the
cone we have $\theta^{T}y\geq 0$, and the last line still corresponds to a
cone. Therefore the optimal value of last line is no less than $0$; the
removal of $-\theta_{k}$ term given $\theta_{k}>0$ will not affect the optimal
value.

Then, according to the property of minimax,
$\sum_{i=1}^{q}\theta_{i}+\alpha_{0}\theta_{s_{j}}$ should be equal to
$-\theta_{k}$; Otherwise, if
$\sum_{i=1}^{q}\theta_{i}+\alpha_{0}\theta_{s_{j}}$ is larger than
$-\theta_{k}$, we can adjust the value of $\sum_{i=1}^{q}|\theta_{i}|$ and
$\sum_{j=q+1}^{n}|\theta_{j}|=\sum_{j=q+1}^{n}\theta_{j}$ by a small amount
such that the result is better and $||\theta||_{2}=1$ is still satisfied,
which can be repeated until no optimization is possible, and vice versa. With
$\sum_{i=1}^{q}|\theta_{i}|$ and $\sum_{j=q+1}^{n}\theta_{j}$ fixed, the
entries of $\theta$ are equal to each other within each group by symmetry.

Therefore, the solution $\theta$ should be in the form of
($\frac{c_{1}}{\sqrt{n}},\frac{c_{1}}{\sqrt{n}},...,\frac{c_{2}}{\sqrt{n}},\frac{c_{2}}{\sqrt{n}}$)
where $c_{1}<0,c_{2}>0$, and we have the following set of equations:

$\displaystyle qc^{2}_{1}+(n-q)c^{2}_{2}$ $\displaystyle=n$ (17)
$\displaystyle qc_{1}+\alpha_{0}c_{2}$ $\displaystyle=-c_{1}$

With equations in Eq. 17, we get $c_{2}=-(1+\alpha_{1}q)c_{1}/\alpha_{0}$,
$c^{2}_{1}(q+(n-q)\frac{(1+q)^{2}}{\alpha^{2}_{0}})=n$. As $q\leq n$, and as
$||\alpha_{j}||=1$, $\alpha_{0}$ can be further relaxed to the sine of the
smallest angle $p_{0}$ between axes $x\geq 0$ and the other inequality
constraints, and thus we have the lower bound
$F\geq\frac{c_{1}}{\sqrt{n}}=O(\sin p/n^{1.5})$ for the optimization problem
listed in Equation 16, which is also the cosine lower bound to the nearest
normal vector $A^{\prime}_{j}$ and $-e_{k}$ of the activated cone. This is the
denominator of the desired result; the final result follows as we apply such
bound to Lemma 1. ∎

###  A.5 Corollary 4

######  Corollary 4.

For binary constraints where the entries of $A$ (before normalization) are
either $0$ or $1$, the utility gain $R$ of violating $x\geq 0$ constraint is
bounded by $O(n^{1.5}E\sum_{i}\max(w_{i},0))$, where $w$ is the same as
Theorem 3.

###### Proof.

The proof of Corollary 4 is almost the same with Theorem 3, except that if the
constraint $A$ is cardinal, then each entry of $A$ is either $0$ or $1$. Thus,
as all the non-zero entry are the same, we shall replace the third line in
Equation 16 with

$\min_{\theta}\max_{j}\\{\frac{\sum_{i=1}^{q}\theta_{i}}{\sqrt{h_{j}}}+\frac{\theta_{s_{j}}}{\sqrt{h_{j}}},-\theta_{k}\\}$
(18)

where $h_{j}$ is the number of non-zero entry of $A^{\prime}_{j}$, each entry
being $\frac{1}{\sqrt{h_{j}}}$ as normalized $||A^{\prime}_{j}||_{2}=1$. By
substituting $h_{j}$ with $n$555Note in Theorem 3 we have already mentioned
the non-positivity of $\theta_{i}$ where $i\in\\{1,2,...,q\\}$, thus bigger
$h_{j}$ brings smaller objective., this optimization problem can be further
relaxed to

$\min_{\theta}\max_{j}\\{\frac{\sum_{i=1}^{q}\theta_{i}}{\sqrt{n}}+\frac{\theta_{s_{j}}}{\sqrt{n}},-\theta_{k}\\}$
(19)

and we can change Equation 17 to

$\displaystyle qc^{2}_{1}+(n-q)c^{2}_{2}$ $\displaystyle=n$ (20)
$\displaystyle qc_{1}+c_{2}$ $\displaystyle=-\sqrt{n}c_{1}$

Therefore, similar to the proof of Theorem 3, we get
$\frac{c_{1}}{\sqrt{n}}=\sqrt{\frac{1}{q+(n-q)(\sqrt{n}+q)^{2}}}$, which leads
to the bound $O(n^{1.5}E\sum_{i}\max(w_{i},0))$, removing the $\sin p$ in the
denominator. ∎

###  A.6 Theorem 5

######  Theorem 5.

If there is only one equality constraint $B^{T}x=c$ (e.g. $\sum_{i}x_{i}=1$)
and special inequality constraints $x\geq 0$, $Ix\leq b$, then the utility
gain $R$ from violating constraints is bounded by
$O(\frac{n^{1.5}E\sum_{i}\max(w_{i},0)}{\sin{p}})$, where $p$ is the same with
theorem 3, $w$ is the union of active $B^{T}x-c$ and $-x$.666See Theorem 3 for
the meaning of union.

###### Proof.

We first prove the situation where the inequality constraint is only $x\geq
0$. If we only have one equality constraint, then it can be seen as two
separate inequality constraints, which are $B^{T}x\leq c$ and $-B^{T}x\leq-c$;
for any non-degraded activated cone, as $B^{T}x<c$ and $B^{T}x>c$ cannot hold
simutaneously, there is at most one normal vector of constraint in the
activated cone.

If $B^{T}x\leq c$ is active (_i.e.,_ this constraint is violated, now we have
$B^{T}x>c$), then with $B\geq 0$, the case is exactly the same with Theorem 3.
Otherwise, if $-B^{T}x\leq-c$ (_i.e.,_ $B^{T}x\geq c$) is active (_i.e.,_
$B^{T}x\geq c$ is violated, now $B^{T}x<c$), then the normal vector of the
current active constraint $-B^{T}x\leq-c$ is $-B\leq 0$. Note that all other
normal vectors of active constraints are $-e_{k}\leq 0$ (which is the same
situation as that in Equation 15 in the proof of Theorem 3, except that
$A^{\prime}_{j}$ is substituted with $-B^{T}\leq 0$); this indicates that the
condition of Lemma 7 is satisfied under such scenario. Therefore, we can apply
the proof of Theorem 2 in this case and get a better bound than that in the
first scenario.

Then, the full theorem is proved as follows: consider any activated cone
$\mathcal{K}$. For any $i\in\\{1,2,...,n\\}$, $x_{i}\leq b_{i}$ and $x_{i}\geq
0$ cannot be active simultaneously. If the former is activated, we replace
$x_{i}$ with $b_{i}-x_{i}$; if the latter (or neither) is active, we remain
$x_{i}$ as normal. Then, for this activated cone, the scenario is exactly the
same with the situation where the inequality constraints are only $x\geq 0$.
Similarly, if the equality constraint is cardinal, with the same proof of
Corollary 4, we can remove the $\sin p$ in the denominator.

∎

###  A.7 Theorem 6

######  Theorem 6.

Given constraints $Ax\leq b$, $x\geq 0$, and $Bx=c$, where $A,B,b,c\geq 0$,
the utility gain $R$ obtained from violating constraints is bounded by
$O(\sqrt{n}\lambda_{max}\sum_{i}\max(w_{i},0))$, where $\lambda_{max}$ is the
upper bound for eigenvalues of $P^{T}P$ ($\mathcal{P}:x\rightarrow 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\leq
b$, $x\geq 0$, $Bx\leq c$ and $-Bx\leq-c$.

###### Proof.

For any non-degraded activated cone $\mathcal{K}$, we have the following
optimization problem (it is worth noting that reformulating our original
problem to this optimization one is inspired by Equation 11 in Lemma 1):

$\displaystyle F=\min_{\theta}\max_{j\in J,k\in K}\\{\max_{i\in
D,n_{i}\in\mathbb{R}^{n}}$
$\displaystyle\\{n_{i}^{T}\theta\\},A^{\prime}_{j}\theta,-e_{k}^{T}\theta\\}$
(21) $\displaystyle s.t.||\theta||_{2}=1,\theta\in\mathcal{K}$ $\displaystyle
d+r+q\geq n,||n_{i}||_{2}=1$ $\displaystyle~{}\forall i\in
D=\\{1,2,...,d\\},<B_{1},...,B_{d}>=<n_{1},...,n_{d}>$

where $A^{\prime}_{j}$ is the normal vector of the $j$-th inequality
constraints in $Ax\leq b$, $n_{i}$ is the $i$-th normal vector for the
subspace of $Bx=c$ (and thus represents the same subspace that is represented
by the rows of $B$). According to the assumption, we scale $Ax\leq b$ and
$Bx=c$ so that $||n_{i}||=||A^{\prime}_{j}||=1$ for any $i,j$. Note that
different from the optimization problem stated in Theorem 3, in this problem
$q$ can be $0$, which means that there can be no entry of $x\geq 0$ active. To
get the lower bound of $F$, we may relax the function by keeping the normal
vectors $\\{n_{i}\\}$ in their original direction $B$ to $F\geq
G=\min_{\theta}\max_{i\in D,j\in J,k\in
K}\\{B^{\prime}_{i}\theta,A^{\prime}_{j}\theta,-e_{k}^{T}\theta\\}$ where the
$i$-th row of $B$ satisfies $B^{\prime}_{i}=B_{i}$ or $B^{\prime}_{i}=-B_{i}$
(the sign of $B^{\prime}_{i}$ is decided by the direction of $\theta$), and
thus $B^{\prime}_{i}\geq 0$ or $B^{\prime}_{i}\leq 0$ (_i.e.,_ any pair of
elements in $B^{\prime}_{i}$ will not have opposite signs) for any $i$.
Moreover, $||B^{\prime}_{i}||_{2}=1$. We denote the positive $B^{\prime}_{i}$
as $B^{p}_{i_{1}}$, $i_{1}\in D_{1}$ and the negative $B^{\prime}_{i}$ as
$B^{n}_{i_{2}}$, $i_{2}\in D_{2}$. Then similar to Theorem 3,

$G=\min_{\theta}\max_{i_{1}\in D_{1},i_{2}\in D_{2},j\in J,k\in
K}\\{B^{n}_{i_{1}}\theta,B^{p}_{i_{2}}\theta,A^{\prime}_{j}\theta,-e_{k}^{T}\theta\\}$
(22)

where $D_{1}\bigcap D_{2}=\emptyset$, $D_{1}\bigcup D_{2}=D$. we aggregate
$A^{\prime}$, $B^{p}$ and $-e_{k}$, $B^{n}$ into $\alpha$, and get

$G=\min_{\theta}\max_{j\in S}\\{{\alpha}_{j}^{T}\theta\\}$ (23)

where $|S|\geq n$, $\alpha_{j}\geq 0$ or $\leq 0$ for any $j\in S$. Therefore,
we can further relax $G$ by selecting $n$ linearly independent vectors that
are the closest to the minimum product777The existence of such set of vectors
comes from the non-degradation of the activated cone. and ignore the others
for the $\max$ operator, and for the rest of the proof we may assume that
$|S|=n$.

We will next consider a linear transformation $\mathcal{P}:x\rightarrow Px$
from $\mathbb{R}^{n}$ to $\mathbb{R}^{n}$ that transforms $\\{\alpha_{j}\\}$
into an orthogonal normal basis. Then we have

$\displaystyle\lambda_{max}G$
$\displaystyle=\min_{\theta}\max_{j}\lambda_{max}{\alpha}_{j}^{T}\theta$ (24)
$\displaystyle\geq\min_{\theta}\max_{i}{\alpha}_{j}P^{T}P\theta$
$\displaystyle=\frac{1}{\sqrt{n}}\sum_{k=1}^{n}{\alpha}_{j}P^{T}P{\alpha}_{k}$
$\displaystyle=\frac{1}{\sqrt{n}}{\alpha}_{j}P^{T}P{\alpha}_{j}$
$\displaystyle=\frac{1}{\sqrt{n}}$

where $\lambda_{max}$ is the maximum eigenvalue of $P^{T}P$. The row vectors
of $A^{\prime\prime}$ are $\\{\alpha^{\prime}_{i}\\}$. Therefore, our problem
$G$ has a lower bound of $\frac{1}{\sqrt{n}\lambda_{max}}$ where
$\lambda_{max}$ is the maximum eigenvalue for $P^{T}P$; $A^{\prime\prime}$
consists of $n$ of the vectors in $A^{\prime}$ and has the largest maximum
eigenvalue for $P^{T}P$. ∎

Though we do not derive the bound of $\lambda_{max}$ with respect to the
number of dimension $n$ and the angle $p_{0}$ between hyper-planes in the
activation cone, we empirically evaluate the behavior of $\lambda_{max}$ with
respect to $n$ on randomly generated data. We first generate a normal vector
$n\in R^{n}$, with each entry generated independently at random from the
distribution $D$; then we generate $n$ vectors
$\\{\alpha_{1},...,\alpha_{n}\\}$ in $R^{n}$ with their entries either all not
smaller than $0$ or all not greater than $0$; the orthant is chosen with
probability $0.5$. Each entry is independently generated from $D$ and shifted
by a constant to enforce the signs. We ensure that $\forall i$,
$\alpha_{i}^{T}n\geq 0$ by discarding the vectors that do not satisfy such
constraint. $D$ can be uniform distribution $U(0,1)$, the distribution for
absolute value of normal distribution $|N(0,1)|$, or beta distribution
$B(2,2)$. The $\alpha_{i}$ is then normalized, and we record the largest
eigenvalue of $R^{T}R$, where $A=\begin{bmatrix}\alpha_{1}^{T}\\\
\alpha_{2}^{T}\\\ ...\\\ \alpha_{n}^{T}\end{bmatrix}$, $A=RQ$ is the RQ
decomposition of $A$. We repeat $1000$ times for each $n$ and record the mean
and maximum eigenvalue for $R^{T}R$. Below is the result of our evaluation; it
shows that the maximum eigenvalue $\lambda_{max}$ is approximately $O(n^{2})$.
See Figure 4,5, and 6 for illustration.

Figure 4: Empirical estimation of the max maximum eigenvalue over $1000$
trials with $D$ being uniform. The $x$ axis is $n$, and the $y$ axis is
$\lambda_{max}$. The $2$-degree polynomial fitting curve is
$0.1282x^{2}+2.392x-32.89$. Figure 5: Empirical estimation of the max maximum
eigenvalue over $1000$ trials with $D$ being Gaussian. The $x$ axis is $n$,
and the $y$ axis is $\lambda_{max}$. The $2$-degree polynomial fitting curve
is $0.2693x^{2}+13.95x-716.8$. Figure 6: The max maximum eigenvalue over
$1000$ trials with $D$ being $Beta(2,2)$. The $x$ axis is $n$, and the $y$
axis is $\lambda_{max}$. The $2$-degree polynomial fitting curve is
$0.1325x^{2}+3.403x-69.95$.

##  Appendix B Choices of Surrogate $\max$ Functions

The final goal of surrogate $max$ function is to relax the objective with the
term $\alpha^{T}\max(z=Cx-d,0)$, making it differentiable over
$\mathbb{R}^{n}$. Such objective is a piecewise function with respect to $z$
with two segments: one is constant $0$ with derivative $0$, the other is
linear with a constant derivative.

At first glance, sigmoidal surrogates are seemingly the most straightforward
candidate for modeling the derivative of such a piecewise function. For a
sigmoidal approximation of soft constraints, $S(z)$ should satisfy the
following four conditions, among which the first two are compulsory, and the
third and fourth can be slightly altered (e.g. by setting $\epsilon_{2}=0$ and
remove the fourth condition).

  1. 1.

When $z\rightarrow\infty$, $S^{\prime}(z)\rightarrow
M^{-}=(1+\epsilon_{1})^{-}$. $\epsilon_{1}>0$ should be a small amount, and it
serves as a perturbation since otherwise $g^{\prime}(z)$ would be always
greater than $0$. However, $\epsilon_{1}$ should not be too small to let the
optimal point be too far away from the original constraint (for the gradient
will reach $\alpha$ too late).

  2. 2.

$S$ must be differentiable, and must have a closed-form inverse function.

  3. 3.

When $z\rightarrow-\infty$, $S^{\prime}(z)\rightarrow\epsilon_{2}^{+}<0$,
where $\epsilon_{2}$ is constant and very close to $0$. The $S^{\prime}(z)$
should cause as small impact as possible when $z<0$ (when there is no waste),
and meanwhile it should (very slightly) encourage $z$ to grow to the limit
instead of discouraging.

  4. 4.

When $z=0$, $S^{\prime}(z)=0$, which means the penalty function $S$ has no
influence on the clipping border. This condition is important for linear
programming; it can be removed for quadratic programming.

To satisfy the conditions above, we need a closed-form sigmoidal function with
closed-form inverse function to be $S^{\prime}(z)$, the derivative of the
surrogate. Such function can be Sigmoid function
($S^{\prime}(z)=\frac{1}{1+e^{-z}}$), Tanh/arctan function, or fractional
function ($S^{\prime}(z)=0.5(\frac{z}{\sqrt{z^{2}+1}}+1)$, or equivalently
$S(z)=0.5\sqrt{z^{2}+1}+0.5z$). With the form of $S(z)$ confirmed, our
function should look like888To satisfy the requirements above, the function
needs scaling and translation; we omit them for the simplicity of formulae by
simply set $\alpha=1$, $\epsilon_{1}=\epsilon_{2}=0$.
$L(x)=g(x,\theta)-\alpha^{T}S(z)$, with the derivative $\frac{\partial
L}{\partial x}=\frac{\partial g}{\partial x}-\frac{\partial z}{\partial
x}S^{\prime}(z)\alpha$.

Take $g(x)=\theta^{T}x$, $S(z)=0.5\sqrt{1+z^{2}}+0.5z$ as an example. The
optimal point satisfies

$\theta=C^{T}\text{diag}(\frac{0.5z}{(1+z^{2})^{0.5}}+0.5)\alpha\\\ $ (25)

Let $y_{i}=(\frac{0.5z}{(1+z_{i}^{2})^{0.5}}+0.5)\alpha_{i}$,
$w_{i}=(\frac{0.5x}{(1+x_{i}^{2})^{0.5}}+0.5)$, then $\theta=C^{T}y$. $S(z)$
can be substituted by any function with sigmoidal function (e.g.
$arctan,tanh$, and $sigmoid$) as its derivative. However, $C$ is impossible to
be invertible, for $C$ includes an identity matrix representing $x\geq 0$ and
must have strictly more rows than columns if $C$ has any other constraint.
Therefore, we cannot solve $z$ by first solving $y$ exactly in a linear
system; we need to solve the equation of fractional function, or estimate $y$
with pseudo-inverse. Hence, the closed-form solution with sigmoidal surrogate
is at least non-attractive.

Figure 7: S(z), the proposed surrogate function of $max(z,0)$

To derive a closed-form solution $x$ for the equation $0=\frac{\partial
L}{\partial x}=\frac{\partial g}{\partial x}-\frac{\partial z}{\partial
x}S^{\prime}(z)\alpha$, $S^{\prime}(z)$ should be in a simple form (e.g.
polynomial form with degrees lower or equal to $4$ with respect to $x$).
Unfortunately, no basic function except the sigmoid function can satisfy the
conditions listed above; the degree limit makes higher-order Taylor expansion
infeasible. Thus, our only choice left is to find a differentiable piecewise
function with linear parts on both sides; we shall first solve the optimal
point of the surrogate numerically, and extend the segment where the optimal
point locates to $R^{n}$; the gradient $\frac{\partial x}{\partial\theta}$ is
correct due to the uniqueness of optimality (which requires a convex
surrogate).

Therefore, we choose the function $S(z)$ to be

$S(z)=\begin{cases}0\quad&\text{if}\,z<-\frac{1}{4K}\\\
K(z+\frac{1}{4K})^{2}\quad&\text{if}\,-\frac{1}{4K}\leq z\leq\frac{1}{4K}\\\
z\quad&\text{if}\,z\geq\frac{1}{4K}\\\ \end{cases}$ (26)

as illustrated in Figure 7, where $K$ is a constant. The surrogate is chosen
for its simplicity. Not only does such choice satisfy our conditions with
minimum matrix computation, but moreover, $S^{\prime}(z)$ is linear, which
means the equation can be solved in a linear system. The limitation of such
method is that it introduces a hyper-parameter $K$ to tune: when $K$ is large,
if the predicted optimal point is far from ground-truth, the gradient
$\frac{\partial f(x,\theta_{\text{real}})}{\partial x}$ will be too large,
making the training process instable. On the other hand, if $K$ is small, the
optimal point may deviate from the feasible region too far, as the hard
constraints are not fully enforced within the region $z\in[-1/4K,4K]$.
Fortunately, our empirical evaluation shows that the effect of $K$ is smooth,
and can be tuned by a grid search.

##  Appendix C Learning and Inference Algorithm

In general, the learning and inference procedures based on our method have the
same prediction+optimization workflow as the computational graph we gave in
Fig.1. In this section, we propose the detailed algorithm with emphasis on the
specific steps in our method.

Alg. 1 is the learning algorithm of stochastic gradient descent in mini-batch.
For each sample $(\xi_{i},\theta_{i})$, it predicts the context variables as
$\hat{\theta_{i}}$ (Statement 3) and then solve the program $f$ to get the
optimal decision variable as $\hat{x_{i}^{*}}$ (Statement 4). Note that here
$f$ is solved by solvers such as Gurobi and CPLEX that are capable to get the
deterministic optimal solutions. With $\hat{x_{i}^{*}}$, the piece of $S(z)$
is determined, and thus the surrogate objective $\bar{f}$. From Statement 5 to
7, we decide the form of surrogate $\bar{f}$ for $\hat{x_{i}^{*}}$ by first
calculating $\hat{z_{i}^{*}}$, then use $\hat{z_{i}^{*}}$ to decide $M_{i}$
and $U_{i}$ in the surrogate. $[\cdot]$ is an element-wise indicator function.
Specially, if $C$ or $d$ is the predicted parameter, $\hat{z_{i}^{*}}$ should
be calculated by the predicted value of $C$ or $d$; moreover, in statement 8
we should calculate $\frac{\partial r}{\partial x}$ on the estimated segment,
but with ground-truth value of $C$ or $d$. In Statement 8, we compute the
gradient by back-propagation. Specifically, we directly compute
$\frac{\partial x}{\partial\theta}$ and $\frac{\partial r}{\partial x}$ by
their analytical forms, while $\frac{\partial\theta}{\partial\psi}$ by auto-
grad mechanism of end-to-end learning software such as PyTorch [17], which we
use throughout our experiments. Finally, in Statement 9, we update the
parameters of $\psi$ with accumulated gradients in this mini-batch.

Input : dataset $D=\\{(\xi_{i},\theta_{i})\\}_{i=1}^{N}$

Input : optimization settings $\\{A,b,C,d,\alpha\\}$

Input : derived hyper parameters $\\{K\\}$

Input : learning rate $a$ and batch size $s$

Output : the predictive model $\Phi_{\psi}$ with parameters $\psi$

begin

1 Sample mini-batch $D_{k}\sim D$;

2 foreach  _$(\xi_{i},\theta_{i})\in D_{k}$_ do

// Estimate $\theta_{i}$ with the model $\Phi_{\psi}$.

3 $\hat{\theta_{i}}\leftarrow\Phi_{\psi}(\xi_{i})$;

// Solve the original $f$ with $\hat{\theta}_{i}$.

4 $\hat{x_{i}^{*}}\leftarrow{\textrm{arg}}\max_{x}{f(x,\hat{\theta_{i}})}$;

// decide surrogate $\bar{f}$ for $\hat{x_{i}^{*}}$ by calculating $M_{i}$ and
$U_{i}$

5$\hat{z_{i}^{*}}\leftarrow C\hat{x_{i}^{*}}-d$;

6 $M_{i}\leftarrow diag([-1/4K\leq\hat{z_{i}^{*}}\leq 1/4K])$;

7 $U_{i}\leftarrow diag([\hat{z_{i}^{*}}\geq 1/4K])$;

// Compute $\frac{\partial\hat{r_{i}^{*}}}{\partial\psi}$, where
$\bar{r_{i}^{*}}=\bar{f}(\hat{x_{i}^{*}},\theta_{i})$

8
$\frac{\partial\bar{r_{i}^{*}}}{\partial\psi}\leftarrow\frac{\partial\hat{\theta_{i}}}{\partial\psi}\times\frac{\partial\hat{x_{i}^{*}}}{\partial\hat{\theta_{i}}}\times\frac{\partial\bar{r_{i}^{*}}}{\partial\hat{x_{i}^{*}}}\mid_{\theta_{i},~{}\hat{\theta_{i}},~{}\hat{x_{i}^{*}}}$

// Update $\psi$ with ascent gradients

9
$\psi\leftarrow\psi+a\times\frac{1}{s}\times\sum\limits_{i=1}^{s}\frac{\partial\bar{r_{i}^{*}}}{\partial\psi}$

Algorithm 1 SGD Learning with the surrogate $\bar{f}$

Alg. 2 is the program with context variable prediction for inference; it is
the same with the traditional predict-then-optimize paradigm. Note that the
surrogate is no longer needed in this phase.

Input : dataset $D=\\{\xi_{i}\\}_{i=1}^{N}$

Input : optimization settings $\\{A,b,C,d,\alpha\\}$

Input : derived hyper parameters $\\{K\\}$

Input : learned predictive model $\Phi_{\psi}$

Output : Solutions {$(\hat{x_{i}^{*}},\hat{r_{i}^{*}})$}

begin

1 foreach  _$\xi_{i}\in D$_ do

// Estimate $\theta_{i}$ with the model $\Phi_{\psi}$.

2 $\hat{\theta_{i}}\leftarrow\Phi_{\psi}(\xi_{i})$;

// Solve $f$ with the estimated $\hat{\theta}_{i}$.

3 $\hat{x_{i}^{*}}\leftarrow{\textrm{arg}}\max_{x}{f(x,\hat{\theta_{i}})}$;

// Evaluate $f$ with the real $\theta_{i}$.

4 $\hat{r}_{i}^{*}\leftarrow f(\hat{x_{i}^{*}},\theta_{i})$;

Algorithm 2 optimization with predictive variables

##  Appendix D Derivation of Gradients

In this section, we will show the derivation process of the gradient with
respect to the three problems stated in the main paper.

###  D.1 Linear Programming with Soft Constraints

The problem formulation is:

$\max_{x}\theta^{T}x-\alpha^{T}\max(Cx-d,0)\text{, s.t. }Ax\leq b$ (27)

where $\alpha\in\mathbb{R}^{n}\geq 0,A\in\mathbb{R}^{m_{1}\times n}\geq
0,b\in\mathbb{R}^{m_{1}}\geq 0$, and $\theta\in\mathbb{R}^{n}$ is to be
predicted. In this formulation and the respective experiment (synthetic linear
programming), we assume that there is no equality constraint.

Let $C^{\prime}=\begin{bmatrix}C\\\ A\\\
-I\end{bmatrix},d^{\prime}=\begin{bmatrix}d\\\ b\\\
0\end{bmatrix},\gamma=\begin{bmatrix}\alpha\\\ O(\frac{n^{1.5}E}{\sin p})\\\
O(\frac{n^{1.5}E}{\sin{p}})\end{bmatrix}$, where $E$ is the upper bound of
$||\theta||_{2}$, $p$ is the minimum angle between hyper-planes of $Ax\leq b$
and the axes. Then, we can write the surrogate as

$\displaystyle\begin{split}\gamma^{T}S(z)&=\gamma^{T}(0.5M(z+\frac{1}{4K})^{2}+Uz)\\\
\end{split}$ (28)

Calculating the derivative of $\theta^{T}x-\gamma^{T}S(z)$ at the optimal
point, we get

$\displaystyle\begin{split}\theta&=C^{\prime
T}(M(\text{diag}(z)+\frac{1}{4K}I)+U)\gamma=C^{\prime
T}M\text{diag}(\gamma)(C^{\prime}x-d^{\prime})+C^{\prime
T}(\frac{1}{4K}M+U)\gamma\end{split}$ (29)

This equation gives us the analytical solution of $x$:

$x=(C^{\prime T}M\text{diag}(\gamma)C^{\prime})^{-1}(\theta+C^{\prime
T}M\text{diag}(\gamma)d^{\prime}-C^{\prime T}(\frac{1}{4K}M+U)\gamma)$ (30)

Based on such derivative solution of $x$, we differentiate with respect to
$\theta$ on both sides:

$\displaystyle\begin{split}\frac{\partial x}{\partial\theta}&=(C^{\prime
T}M\text{diag}(\gamma)C^{\prime})^{-1}\end{split}$ (31)

On the other hand, given real parameter $\theta_{real}$, we calculate the
derivative of $f(x,\theta_{\text{real}})$:

$\displaystyle\begin{split}\frac{\partial f(x,\theta_{\text{real}})}{\partial
x}&=\theta_{\text{real}}-C^{\prime
T}M\text{diag}(\gamma)(C^{\prime}x-d^{\prime})-C^{\prime
T}(\frac{1}{4K}M+U)\gamma\end{split}$ (32)

###  D.2 Portfolio Selection with Soft Constraints

We consider minimum variance portfolio ([16]) which maximizes the return while
minimizes risks of variance. The problem formulation is:

$\displaystyle\begin{split}&\max_{x}\theta^{T}x-x^{T}Qx-\alpha^{T}\max(Cx-d,0)\\\
\text{s.t.}&~{}x^{T}\mathbf{1}=1,~{}x\geq 0,~{}Q\geq 0,~{}\alpha\geq 0\\\
\end{split}$ (33)

where $x\in\mathbb{R}^{n}$ is the decision variable vector – equity weights,
$\theta\in\mathbb{R}^{n}$ is the equity returns, the semi-definite positive
$Q\in\mathbb{R}^{n\times n}$ is the covariance matrix of returns
$\theta\in\mathbb{R}^{n}$. Equivalently, we rewrite the constraints of $x$ to
fit our surrogate framework, as $Ax\leq b$ where
$A=\begin{bmatrix}-\mathbf{1}_{1\times n}\\\ \mathbf{1}_{1\times n}\\\
-I\end{bmatrix},b=\begin{bmatrix}-1\\\ 1\\\ 0\end{bmatrix}$. Let
$C^{\prime}=\begin{bmatrix}C\\\ A\end{bmatrix},d^{\prime}=\begin{bmatrix}d\\\
b\end{bmatrix},\gamma=\begin{bmatrix}\alpha\\\
O(n^{1.5}E)\textbf{1}\end{bmatrix}$, and we have surrogate
$\gamma^{T}S(z)=\gamma^{T}(0.5M(Z+\frac{1}{4K})^{2}+Uz)$. Then, we may derive
the optimal solution $x$ and the gradient of
$f(x,\theta_{\text{real}},Q_{\text{real}})$ with real data with respect to $x$
as

$\displaystyle x$ $\displaystyle=(2Q+C^{\prime
T}M\text{diag}(\gamma)C^{\prime})^{-1}(\theta+C^{\prime
T}M\text{diag}(\gamma)d^{\prime}-C^{\prime T}(U+\frac{M}{4K})\gamma)$ (34)
$\displaystyle\frac{\partial
f(x,\theta_{\text{real}},Q_{\text{real}})}{\partial x}$
$\displaystyle=\theta_{\text{real}}-2Q_{\text{real}}x-C^{\prime
T}M\text{diag}(\gamma)(C^{\prime}x-d^{\prime})-C^{\prime
T}(\frac{1}{4K}M+U)\gamma$

Differentiating on both sides of the analytical solution of $x$, we get:

$\displaystyle\frac{\partial x}{\partial\theta}$ $\displaystyle=(2Q+C^{\prime
T}M\text{diag}(\gamma)C^{\prime})^{-1}$ (35) $\displaystyle\frac{\partial
x}{\partial Q}$ $\displaystyle=(\theta+C^{\prime
T}M\text{diag}(\gamma)d^{\prime}-C^{\prime
T}(U+\frac{M}{4K})\gamma)\frac{\partial(2Q+C^{\prime
T}M\text{diag}(\gamma)C^{\prime})^{-1}}{\partial Q}$

To derive a simplified norm of $\frac{\partial x}{\partial Q}$, let
$R=(2Q+C^{\prime T}M\text{diag}(\gamma)C^{\prime})^{-1}$, $S=C^{\prime
T}M\text{diag}(\gamma)C^{\prime}$, $\beta=\theta+C^{\prime
T}M\text{diag}(\gamma)d^{\prime}-C^{\prime T}(U+\frac{M}{4K})\gamma$. With
such notations, we can simplify the previous results to

$x=(2Q+S)^{-1}\beta=R\beta,\ \frac{\partial x}{\partial\theta}=R^{T}=R,\
\frac{\partial f(x,\theta_{\text{real}},Q_{\text{real}})}{\partial
x}=\beta_{\text{real}}-R_{\text{real}}^{-1}x$ (36)

and the derivative $\frac{\partial x_{i}}{\partial Q_{j,k}}$ and
$\frac{\partial f(x,\theta_{\text{real}},Q_{\text{real}})}{\partial Q_{j,k}}$
can be derived as follows:

$\displaystyle\frac{\partial x_{i}}{\partial Q_{j,k}}$
$\displaystyle=\sum_{x}\sum_{y}\frac{\partial R_{x,y}}{\partial
Q_{j,k}}\frac{\partial x_{i}}{\partial R_{x,y}}$ (37)
$\displaystyle=\sum_{y}\frac{\partial R_{i,y}}{\partial Q_{j,k}}\beta_{y}$
$\displaystyle=\sum_{y}\beta_{y}\sum_{p}\sum_{q}\frac{\partial(2Q+S)_{p,q}}{Q_{j,k}}\frac{\partial
R_{i,y}}{\partial(2Q+S)_{p,q}}$
$\displaystyle=\sum_{y}2\beta_{y}\frac{\partial
R_{i,y}}{\partial(2Q+S)_{j,k}}$
$\displaystyle=-\sum_{y}2\beta_{y}R_{i,j}R_{k,y}$ $\displaystyle\frac{\partial
f(x,\theta_{\text{real}},Q_{\text{real}})}{\partial Q_{j,k}}$
$\displaystyle=-\sum_{i}(\sum_{y}2\beta_{y}R_{i,j}R_{k,y})(\beta_{\text{real},i}-(R_{\text{real}}^{-1}x)_{i})$
(38)
$\displaystyle=2\sum_{i}((R_{\text{real}}^{-1}x)_{i}-\beta_{\text{real},i})R_{i,j}\sum_{y}\beta_{y}R_{k,y}$

Let
$p_{j}=\sum_{i}{R^{T}_{j,i}((R^{-1}_{\text{real}}x)_{i}-\beta_{\text{real},i})},t_{k}=\sum_{y}\beta_{y}R_{k,y}$
($p=R^{T}(R_{\text{real}}^{-1}x-\beta_{\text{real}})$, $t=x$), and finally

$\frac{\partial f(x,\theta_{\text{real}},Q_{\text{real}})}{\partial
Q_{j,k}}=2p_{j}x_{k},\ \frac{\partial
f(x,\theta_{\text{real}},Q_{\text{real}})}{\partial Q}=2px^{T}$ (39)

###  D.3 Resource Provisioning

The problem formulation of resource provisioning is

$\min_{x}\alpha_{1}^{T}\max(Cx-d,0)+\alpha_{2}^{T}\max(d-Cx,0),\
s.t.~{}x^{T}\mathbf{1}=1,x\geq 0$ (40)

Let $C^{\prime}\begin{bmatrix}C\\\ -C\\\ -I\\\ \mathbf{1_{1\times n}}\\\
\mathbf{-1_{1\times n}}\end{bmatrix}$, $d^{\prime}=\begin{bmatrix}d\\\ -d\\\
0\\\ 1\\\ -1\end{bmatrix}$, $\gamma=\begin{bmatrix}\alpha_{1}\\\ \alpha_{2}\\\
O(n^{1.5}E)\mathbf{1}\\\ O(n^{1.5}E)\\\ O(n^{1.5}E)\end{bmatrix}$,
$P=(C^{\prime T}M\text{diag}(\gamma)C^{\prime})^{-1}$, $\beta=(C^{\prime
T}M\text{diag}(\gamma)d^{\prime}-C^{\prime T}(\frac{M}{4K}+U)\gamma)$,
$\eta=M\text{diag}(\gamma)d^{\prime}-(\frac{M}{4K}+U)\gamma$. Then for the
derivative $\frac{\partial f(x,C^{\prime}_{\text{real}})}{\partial x}$ and the
analytical solution of $x$, we have

$\displaystyle\frac{\partial f(x,C^{\prime}_{\text{real}})}{\partial x}$
$\displaystyle=-{C^{\prime}}_{\text{real}}^{T}M_{\text{real}}\text{diag}(\gamma)(C^{\prime}_{\text{real}}x-d^{\prime})-{C^{\prime}}_{\text{real}}^{T}(\frac{1}{4K}M_{\text{real}}+U_{\text{real}})\gamma$
(41) $\displaystyle x$ $\displaystyle=(C^{\prime
T}M\text{diag}(\gamma)C^{\prime})^{-1}(C^{\prime
T}M\text{diag}(\gamma)d^{\prime}-C^{\prime T}(\frac{M}{4K}+U)\gamma)=P\beta$

According to the analytical solution of optimal point $x$, the derivative
$\frac{\partial x_{i}}{\partial C^{\prime}_{k,l}}$ is

$\frac{\partial x_{i}}{\partial C^{\prime}_{k,l}}=\sum_{j}\frac{\partial
P_{i,j}}{\partial
C^{\prime}_{k,l}}\beta_{j}+\sum_{j}P_{i,j}\frac{\partial\beta_{j}}{\partial
C^{\prime}_{k,l}}$ (42)

For the first term of the derivative above, we have:

$\displaystyle\frac{\partial P_{i,j}}{\partial C^{\prime}_{k,l}}$
$\displaystyle=\frac{\partial(C^{\prime
T}M\text{diag}(\gamma)C^{\prime})_{i,j}^{-1}}{C^{\prime}_{k,l}}$ (43)
$\displaystyle=\sum_{p}\sum_{q}\frac{\partial(C^{\prime
T}M\text{diag}(\gamma)C^{\prime})_{p,q}}{\partial
C^{\prime}_{k,l}}\frac{P_{i,j}}{P^{-1}_{p,q}}$
$\displaystyle=\sum_{p}\sum_{q}\frac{\partial\sum_{x}\sum_{y}C^{\prime}_{x,p}M\text{diag}(\gamma)_{x,y}C^{\prime}_{y,q}}{\partial
C^{\prime}_{k,l}}(-P_{i,p}P_{q,j})$
$\displaystyle=\sum_{q}\sum_{y}[p==l]M\text{diag}(\gamma)_{k,y}C^{\prime}_{y,q}(-P_{i,p}P_{q,j})+\sum_{p}\sum_{x}[q==l]M\text{diag}(\gamma)_{x,k}C^{\prime
T}_{p,x}(-P_{i,p}P_{q,j})$
$\displaystyle=\sum_{q}-P_{i,l}(M\text{diag}(\gamma))_{k,*}C^{\prime}_{*,q}P_{q,j}+\sum_{p}-P_{l,j}C^{\prime
T}_{p,*}(M\text{diag}(\gamma))_{*,k}P_{i,p}$
$\displaystyle=-(P_{i,l}(M\text{diag}(\gamma))_{k,*}C^{\prime}P_{*,j}+P_{l,j}P_{i,*}C^{\prime
T}(M\text{diag}(\gamma))_{*,k})$
$\displaystyle=-(P_{i,l}(M\text{diag}(\gamma)C^{\prime}P)_{k,j}+P_{l,j}(PC^{\prime
T}M(\text{diag}(\gamma)))_{i,k})$

Therefore, the simplified result for the first term of the derivative in
Equation 42 is:

$\sum_{j}\frac{\partial P_{i,j}}{\partial
C^{\prime}_{k,l}}\beta_{j}=-((M\text{diag}(\gamma)C^{\prime}P\beta)_{k}P_{i,l}+(P\beta)_{l}(PC^{\prime
T}M\text{diag}(\gamma))_{i,k})$ (44)

For the second term, we have:

$\frac{\partial\beta_{j}}{\partial C^{\prime}_{k,l}}=\frac{\partial(C^{\prime
T}\eta)_{j}}{C^{\prime}_{k,l}}=\frac{\partial\sum_{p}C^{\prime}_{p,j}\eta_{p}}{\partial
C^{\prime}_{k,l}}=\eta_{k}[j==l]$ (45)

The simplified second term is thus

$\sum_{j}P_{i,j}\frac{\partial\beta_{j}}{\partial
C^{\prime}_{k,l}}=P_{i,l}\eta_{k}$ (46)

Finally, the derivative can be written as

$\frac{\partial x_{i}}{\partial
C^{\prime}_{k,l}}=-((M\text{diag}(\gamma)C^{\prime}P\beta)_{k}P_{i,l}+(P\beta)_{l}(PC^{\prime
T}M\text{diag}(\gamma))_{i,k})+P_{i,l}\eta_{k}$ (47)

##  Appendix E Benchmark Details : Dataset and Problem Settings

Our code is public in the repo:
https://github.com/PredOptwithSoftConstraint/PredOptwithSoftConstraint.

###  E.1 Synthetic Linear Programming

####  E.1.1 Prediction Dataset

We generate the synthetic dataset $\\{\xi_{i},\theta_{i}\\}_{i=1}^{N}$ under a
general structural causal model ([21]). In the original form, it is like:

$\displaystyle\begin{split}z\sim N(0,\Sigma)\\\ \xi=g(z)+\epsilon_{1}\\\
\theta=h(z)+\epsilon_{2}\\\ \end{split}$ (48)

where $z$ is the latent variable, $\xi$ observed features of $z$, and $\theta$
the result variable caused by $z$. According to physical knowledge, $h$ can be
a process of linear, quadratic, or bi-linear form. However, it is difficult to
get an explicit form of reasonable $g$. Instead, $g^{-1}$ can be well-
represented by deep neural networks.

Thus, alternatively, we use the following generative model:

$\displaystyle\begin{split}\xi^{*}\sim N(0,\Sigma)\\\ z=m(\xi^{*})\\\
\theta=h(z)+\epsilon_{2}\\\ \xi=\xi^{*}+\epsilon_{1}\end{split}$ (49)

where $m$ behaves as $g^{-1}$. In our experiment settings, $\Sigma=I+QQ^{T}$,
where each element of $Q$ is generated randomly at uniform from $(0,1)$. We
set $m(x)=\sin(2\pi xB)$, where $B$ is a matrix whose elements are generated
randomly at uniform in $\\{0,1\\}$, and $\sin$ is applied element-wisely. We
implement $h(z)$ as a MLP with two hidden layers, and the output is normalized
to $(0,1]$ for each dimension through different data points. Finally, we add a
noise of $0.01\epsilon_{x}$ to $x$, where $\epsilon_{x}\sim N(0,1)$; and
$0.01\epsilon_{\theta}$ to $\theta$, where $\epsilon_{\theta}$ follows a
truncated normal distribution which truncates a normal distribution $N(0,1)$
to $[0,1.5]$.

The dataset is split into training, validation and test sets with the
proportions $50\%$, $25\%$, $25\%$ in respect. The batch size is set to $10$
for $N=100$, $50$ for $N=1000$, and $125$ for $N=5000$.

####  E.1.2 Problem Settings

We generate hard constraint $Ax\leq b$, and soft constraint $Cx\leq d$. Each
element of $A$ or $C$ is first generated randomly at uniform within $(0,1)$,
then set to $0$ with probability of $0.5$. We generate $b,d$ as
$b=0.5A\mathbf{1}$ and $d=0.25C\mathbf{1}$. The soft constraint coefficient
$\alpha$ is generated randomly at uniform from $(0,0.2)$ for each dimension.

###  E.2 Portfolio Optimization

####  E.2.1 Dataset

The prediction dataset is daily price data of SP500 from 2004 to 2017
downloaded by Quandl API [24] with the same settings in [3]. Most settings are
aligned with those in [3], including dataset configuration, prediction model,
learning rate (initial $0.01$ with scheduler), optimizer (Adam), gradient clip
($0.01$), the number of training epochs ($20$), and the problem instance size
(the number of equities being $\\{50,100,150,200,250\\}$).

####  E.2.2 Problem Settings

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
$\alpha^{T}\max(Cx-d,0)$, $\alpha=\frac{15}{n}v$, where each element of $v$ is
generated randomly at uniform from $(0,1)$; the elements of matrix $C$ are
generated independently from $\\{0,1\\}$, where the probability of $0$ is
$0.9$ and $1$ is $0.1$. $K$ is set as $100$.

###  E.3 Resource Provisioning

####  E.3.1 Dataset

The ERCOT energy dataset [25] contains hourly data of energy output from 2013
to 2018 with $52535$ data points. We use the first $70\%$ as the training set,
the middle $10\%$ as the validation set, and the last $20\%$ as the test set.
We normalize the dataset by dividing the labels by $10^{4}$, which makes the
typical label becomes value around $(0.1,1)$. We train our model with a batch
size of $256$; each epoch contains $144$ batches. We aim to predict the matrix
$C\in\mathbb{R}^{24\times 8}$, where $24$ represents the following $24$ hours
and $8$ represents the $8$ regions, which are $\\{\text{COAST, EAST, FWEST,
NCENT, NORTH, SCENT, SOUTH, WEST}\\}$. The data is drawn from the dataset of
the corresponding region. the decision variable $x$ is $8$-dimensional, and
$d=0.5\mathbf{1}+0.1N(0,1)$.

####  E.3.2 Problem Settings

We test five sets of $(\alpha_{1},\alpha_{2})$, which are
$(50\times\mathbf{1},0.5\times\mathbf{1})$,
$(5\times\mathbf{1},0.5\times\mathbf{1})$, $(\mathbf{1},\mathbf{1})$,
$(0.5\times\mathbf{1},5\times\mathbf{1})$, and
$(0.5\times\mathbf{1},50\times\mathbf{1})$, against two-stage with $L1$-loss,
$L2$-loss, and weighted $L1$-loss. We use AdaGrad as optimizer with learning
rate $0.01$, and clip the gradient with norm $0.01$. The feature is a
$(8\times 24\times 77)$-dimensional vector for each matrix $C$; we adopted
McElWee’s Blog [26] for the feature generation and the model. Weighted L1-loss
has the following objective:

$\alpha_{2}^{T}\max(Cx-d,0)+\alpha_{1}^{T}\max(d-Cx,0)$ (50)

Note that $\alpha_{2}$ and $\alpha_{1}$ are exchanged in the objective.
Intuitively, this is because an under-estimation of entries of $C$ will cause
a larger solution of $x$, which in turn makes $C_{\text{real}}x-d$ larger at
test time, and vice versa. Another thing worth noting is that SPO+ cannot be
applied to the prediction of the matrix $C$ or vector $d$, for it is designed
for the scenario where the predicted parameters are in the objective.

##  Appendix F Supplementary Experiment Results

###  F.1 Linear Programming with Soft Constraints

The effect of $K$. Table 4 provides the mean and standard deviation of regrets
under varying values of $K$. we can see that our method performs better than
all other methods with most settings of $K$.

Empirically, to find an optimal $K$ for a given problem and its experimental
settings, a grid search with roughly adaptive steps suffices. For example,in
our experiments, we used a proposal where neighbouring coefficient has $5x$
difference (e.g. $\\{0.2,1,5,25,125\\}$) works well. Quadratic objective, as
in the second experiment, usually requires larger $K$ than linear objective as
in the first and the third experiment.

|  | Regret |
---|---|---|---
$N$ | Problem Size | ours($K=0.2$) | $K=1.0$ | $K=5.0$ | $K=25.0$ | $K=125.0$
100 | (40, 40, 0) | 2.423$\pm$0.305 | 2.378$\pm$0.293 | 2.265$\pm$0.238 | 2.301$\pm$0.325 | 2.258$\pm$0.311
| (40, 40, 20) | 2.530$\pm$0.275 | 2.574$\pm$0.254 | 2.384$\pm$0.277 | 2.321$\pm$0.272 | 2.350$\pm$0.263
| (80, 80, 0) | 5.200$\pm$0.506 | 5.529$\pm$0.307 | 5.610$\pm$0.351 | 5.563$\pm$0.265 | 5.502$\pm$0.353
| (80, 80, 40) | 4.570$\pm$0.390 | 4.747$\pm$0.470 | 4.758$\pm$0.531 | 4.780$\pm$0.450 | 4.796$\pm$0.624
1000 | (40, 40, 0) | 1.379$\pm$0.134 | 1.387$\pm$0.130 | 1.365$\pm$0.147 | 1.346$\pm$0.144 | 1.352$\pm$0.140
| (40, 40, 20) | 1.587$\pm$0.112 | 1.553$\pm$0.114 | 1.523$\pm$0.109 | 1.507$\pm$0.102 | 1.506$\pm$0.102
| (80, 80, 0) | 3.617$\pm$0.125 | 3.431$\pm$0.100 | 3.540$\pm$0.105 | 3.513$\pm$0.097 | 3.507$\pm$0.117
| (80, 80, 40) | 2.781$\pm$0.165 | 2.817$\pm$0.167 | 2.819$\pm$0.155 | 2.820$\pm$0.157 | 2.826$\pm$0.197
5000 | (40, 40, 0) | 1.041$\pm$0.094 | 1.065$\pm$0.098 | 1.045$\pm$0.099 | 1.038$\pm$0.102 | 1.037$\pm$0.100
| (40, 40, 20) | 1.278$\pm$0.072 | 1.248$\pm$0.076 | 1.223$\pm$0.075 | 1.220$\pm$0.071 | 1.221$\pm$0.074
| (80, 80, 0) | 3.074$\pm$0.098 | 2.845$\pm$0.064 | 2.884$\pm$0.100 | 2.889$\pm$0.083 | 2.864$\pm$0.086
| (80, 80, 40) | 2.217$\pm$0.123 | 2.222$\pm$0.111 | 2.172$\pm$0.098 | 2.174$\pm$0.113 | 2.177$\pm$0.105
Table 4: The mean and standard deviation of the regret of all $K$ in our
experiment.

Prediction error. Table 5 and Table 6 are the mean and standard deviation of
prediction MSE for all methods. We can find that two-stage methods have an
advantage over end-to-end methods on prediction error. however, as
demonstrated in our experiments, such advantage does not necessarily lead to
the advantage on the regret performance. Surprisingly, although DF[9] works
generally on par with other baselines in regret performance, it generates a
very large MSE error compared with our method; closer examination suggests
that there are some outlier cases where the prediction loss is magnitudes
higher than others. Also, the prediction error of DF gets much higher when no
soft constraint exists, which is probably because the optimal solution remains
the same when all parameters to predict are scaled by a constant factor.

|  | MSE of $\theta$
---|---|---
$N$ | Problem Size | L1 | L2 | SPO+ | DF | ours ($K=0.2$)
100 | (40, 40, 0) | 0.062$\pm$0.006 | 0.063$\pm$0.008 | 0.079$\pm$0.150 | 14.592$\pm$9.145 | 0.094$\pm$0.039
| (40, 40, 20) | 0.061$\pm$0.005 | 0.062$\pm$0.006 | 0.071$\pm$0.006 | 6.414$\pm$4.451 | 0.061$\pm$0.007
| (80, 80, 0) | 0.058$\pm$0.004 | 0.059$\pm$0.004 | 0.066$\pm$0.011 | 28.768$\pm$13.90 | 0.247$\pm$0.009
| (80, 80, 40) | 0.059$\pm$0.004 | 0.059$\pm$0.004 | 0.085$\pm$0.008 | 8.98$\pm$ 9.99 | 0.061$\pm$0.005
1000 | (40, 40, 0) | 0.031$\pm$0.001 | 0.031$\pm$0.001 | 0.034$\pm$0.002 | 34.931$\pm$14.309 | 0.036$\pm$0.003
| (40, 40, 20) | 0.031$\pm$0.001 | 0.031$\pm$0.001 | 0.033$\pm$0.001 | 9.255$\pm$8.497 | 0.033$\pm$0.001
| (80, 80, 0) | 0.030$\pm$0.001 | 0.029$\pm$0.001 | 0.032$\pm$0.002 | 90.590$\pm$35.298 | 0.222$\pm$0.009
| (80, 80, 40) | 0.030$\pm$0.001 | 0.029$\pm$0.001 | 0.033$\pm$0.001 | 1.060$\pm$2.395 | 0.033$\pm$0.002
5000 | (40, 40, 0) | 0.022$\pm$0.000 | 0.022$\pm$0.000 | 0.025$\pm$0.001 | 72.083$\pm$46.685 | 0.025$\pm$0.001
| (40, 40, 20) | 0.022$\pm$0.000 | 0.022$\pm$0.000 | 0.238$\pm$0.001 | 235.789$\pm$842.939 | 0.243$\pm$0.001
| (80, 80, 0) | 0.021$\pm$0.000 | 0.021$\pm$0.000 | 0.023$\pm$0.001 | 147.144$\pm$51.101 | 0.196$\pm$0.011
| (80, 80, 40) | 0.020$\pm$0.000 | 0.020$\pm$0.000 | 0.022$\pm$0.001 | 0.372$\pm$0.311 | 0.024$\pm$0.001
Table 5: The mean and standard deviation of the MSE loss of the predicted variable $\theta$ in the synthetic experiment on the test set. The two-stage methods have better performance on the accuracy of prediction. |  | MSE of $\theta$
---|---|---
$N$ | Problem Size | ours($K=1.0$) | ours($K=5.0$) | ours($K=25.0$) | ours($K=125.0$)
100 | (40, 40, 0) | 0.068$\pm$0.008 | 0.091$\pm$0.037 | 0.188$\pm$0.089 | 0.204$\pm$0.079
| (40, 40, 20) | 0.094$\pm$0.031 | 0.204$\pm$0.074 | 0.244$\pm$0.105 | 0.343$\pm$0.200
| (80, 80, 0) | 0.129$\pm$0.049 | 0.257$\pm$0.182 | 0.641$\pm$0.476 | 0.816$\pm$0.351
| (80, 80, 40) | 0.094$\pm$0.013 | 0.151$\pm$0.035 | 0.172$\pm$0.091 | 0.177$\pm$0.042
1000 | (40, 40, 0) | 0.036$\pm$0.002 | 0.056$\pm$0.010 | 0.143$\pm$0.048 | 0.195$\pm$0.095
| (40, 40, 20) | 0.060$\pm$0.010 | 0.127$\pm$0.051 | 0.172$\pm$0.071 | 0.193$\pm$0.078
| (80, 80, 0) | 0.114$\pm$0.070 | 0.220$\pm$0.097 | 1.124$\pm$0.341 | 1.407$\pm$0.439
| (80, 80, 40) | 0.039$\pm$0.002 | 0.060$\pm$0.006 | 0.078$\pm$0.010 | 0.080$\pm$0.008
5000 | (40, 40, 0) | 0.026$\pm$0.001 | 0.033$\pm$0.004 | 0.076$\pm$0.025 | 0.089$\pm$0.037
| (40, 40, 20) | 0.036$\pm$0.004 | 0.073$\pm$0.019 | 0.110$\pm$0.038 | 0.113$\pm$0.026
| (80, 80, 0) | 0.152$\pm$0.079 | 0.097$\pm$0.040 | 0.703$\pm$0.177 | 0.895$\pm$0.205
| (80, 80, 40) | 0.026$\pm$0.001 | 0.032$\pm$0.003 | 0.038$\pm$0.004 | 0.040$\pm$0.007
Table 6: The mean and standard deviation of the MSE loss of the predicted
variable $\theta$ in the synthetic experiment for different values of $K$.

Detailed behaviour of SPO+. As mentioned in the main paper, SPO+ quickly
becomes overfited in our experiments, shown by Table 7 as empirical evidence.
In our experiments, early stopping starts at epoch $8$, and maximum number of
run epochs is $40$; the results show that SPO+ stops much earlier than other
methods under most settings. Figure 8 plots the timely comparison of test
performance during training, where SPO+ get overfited quickly, although it
performs better in the earlier period of training.

Figure 8: The average regret on test set of $K=25.0$, training set size $5000$ and problem size $(80,80,0)$ with respect to the number of epochs (methods other than ours does not reach $40$ in this figure, for all runs are early-stopped before epoch $40$). While SPO+ is better than two-stage in the best performance, it overfits rather quickly. |  | Average Episode
---|---|---
$N$ | Problem Size | L1 | L2 | SPO+ | DF | ours ($K=0.2$)
100 | (40, 40, 0) | 20.46$\pm$10.20 | 12.53$\pm$3.00 | 19.6$\pm$8.58 | 29.53$\pm$10.24 | 20.53$\pm$10.39
| (40, 40, 20) | 21.47$\pm$9.58 | 14.13$\pm$4.47 | 17.8$\pm$9.98 | 25.46$\pm$11.41 | 22.47$\pm$10.29
| (80, 80, 0) | 22.6$\pm$5.78 | 22.6$\pm$8.52 | 20.07$\pm$7.23 | 30.87$\pm$8.83 | 24.53$\pm$8.02
| (80, 80, 40) | 25.2$\pm$8.38 | 24.73$\pm$8.88 | 19.33$\pm$764 | 24$\pm$9.06 | 28.33$\pm$8.89
1000 | (40, 40, 0) | 14.67$\pm$2.61 | 13.2$\pm$1.74 | 10$\pm$2.04 | 27.93$\pm$7.42 | 13.6$\pm$2.53
| (40, 40, 20) | 13.2$\pm$2.37 | 13.4$\pm$2.50 | 10.33$\pm$1.40 | 24.6$\pm$8.14 | 10.67$\pm$1.23
| (80, 80, 0) | 19.07$\pm$2.25 | 18.07$\pm$1.83 | 11.87$\pm$2.00 | 25.2$\pm$6.41 | 24.27$\pm$11.02
| (80, 80, 40) | 19.93$\pm$2.66 | 19$\pm$1.31 | 12.33$\pm$1.95 | 25.4$\pm$9.22 | 17.73$\pm$4.08
5000 | (40, 40, 0) | 16.47$\pm$2.59 | 16.27$\pm$1.33 | 11.87$\pm$2.10 | 23.33$\pm$7.34 | 11.8$\pm$2.60
| (40, 40, 20) | 17.13$\pm$2.59 | 15.27$\pm$1.79 | 12.8$\pm$2.88 | 21.67$\pm$7.95 | 13.4$\pm$2.77
| (80, 80, 0) | 20.67$\pm$3.20 | 20.8$\pm$2.18 | 13.7$\pm$2.16 | 29.33$\pm$7.22 | 30.07$\pm$9.48
| (80, 80, 40) | 20.47$\pm$2.64 | 21$\pm$2.88 | 12.2$\pm$1.61 | 17.6$\pm$6.84 | 18.53$\pm$5.26
Table 7: The mean and standard deviation of the number of epochs run in each
instance (capped at $40$). Note that SPO+ is significantly more prone to
overfitting in our experiment settings.

Statistical significance tests. Table 9 shows the results of significance test
(one-tailed paired $t$-test), with the assumption that the means of two
distributions are the same. The results show that our method is significantly
better than other methods under almost all settings.

|  |  $t$-value
---|---|---
$N$ | Problem Size | L1 | L2 | SPO+ | DF
100 | (40, 40, 0) | 3.884 | 4.615 | 5.194 | 1.614
| (40, 40, 20) | 5.173 | 5.248 | 6.836 | 1.569
| (80, 80, 0) | 6.024 | 5.661 | 4.523 | 3.610
| (80, 80, 40) | 2.308 | 1.724 | 3.622 | 1.936
1000 | (40, 40, 0) | 8.438 | 5.271 | 4.860 | 1.115
| (40, 40, 20) | 6.547 | 6.212 | 8.562 | 2.620
| (80, 80, 0) | 12.620 | 10.925 | 7.979 | 2.751
| (80, 80, 40) | 7.116 | 4.465 | 6.223 | 8.230
5000 | (40, 40, 0) | 7.130 | 7.808 | 7.546 | 1.175
| (40, 40, 20) | 10.410 | 10.238 | 6.555 | 2.382
| (80, 80, 0) | 8.361 | 7.502 | 6.874 | 0.869
| (80, 80, 40) | 6.550 | 4.388 | 5.066 | 11.620
Table 8: The $t$-value of the one-tailed paired $t$-test between all other methods and our methods under optimal $K$. |  |  $p$-value
---|---|---
$N$ | Problem Size | L1 | L2 | SPO+ | DF
100 | (40, 40, 0) | 8.272 $\times 10^{-4}$ | 2.006$\times 10^{-4}$ | 6.809 $\times 10^{-5}$ | 0.058
| (40, 40, 20) | 7.076$\times 10^{-5}$ | 6.162$\times 10^{-5}$ | 4.054$\times 10^{-6}$ | 0.064
| (80, 80, 0) | 1.561$\times 10^{-5}$ | 2.936$\times 10^{-5}$ | 2.388$\times 10^{-4}$ | 5$\times 10^{-4}$
| (80, 80, 40) | 0.018 | 0.053 | 0.001 | 0.032
1000 | (40, 40, 0) | 3.662$\times 10^{-7}$ | 5.917$\times 10^{-5}$ | 1.263 $\times 10^{-4}$ | 0.137
| (40, 40, 20) | 6.489$\times 10^{-6}$ | 1.133$\times 10^{-5}$ | 3.078$\times 10^{-7}$ | 0.007
| (80, 80, 0) | 2.444$\times 10^{-9}$ | 1.544$\times 10^{-8}$ | 7.064$\times 10^{-7}$ | 0.005
| (80, 80, 40) | 2.600$\times 10^{-6}$ | 2.667$\times 10^{-4}$ | 1.114$\times 10^{-5}$ | 2.937$\times 10^{-9}$
5000 | (40, 40, 0) | 2.545$\times 10^{-6}$ | 9.074$\times 10^{-7}$ | 1.342$\times 10^{-6}$ | 0.125
| (40, 40, 20) | 2.834$\times 10^{-8}$ | 3.487$\times 10^{-8}$ | 6.406$\times 10^{-6}$ | 0.012
| (80, 80, 0) | 4.082$\times 10^{-7}$ | 1.435$\times 10^{-6}$ | 3.818$\times 10^{-6}$ | 0.196
| (80, 80, 40) | 6.460$\times 10^{-6}$ | 3.097$\times 10^{-4}$ | 8.606$\times 10^{-5}$ | 1.585$\times 10^{-12}$
Table 9: The $p$-value of the one-tailed paired $t$-test between all other
methods and our methods under optimal $K$.

###  F.2 Portfolio Optimization

Statistical significance tests. Table 10 shows the result of significance test
(one-tailed paired $t$-test), with the assumption that the means of two
distributions are the same. The results show that our method is significantly
better than other methods under every setting.

|  $t$-value |  $p$-value
---|---|---
#Equities | L1 | L2 | DF | L1 | L2 | DF
50 | 10.329 | 10.628 | 6.141 | 3.213$\times 10^{-8}$ | 2.186$\times 10^{-8}$ | 1.279$\times 10^{-5}$
100 | 21.414 | 23.995 | 6.200 | 2.125$\times 10^{-12}$ | 4.495$\times 10^{-13}$ | 1.157$\times 10^{-5}$
150 | 13.402 | 14.503 | 5.690 | 1.119$\times 10^{-9}$ | 3.971$\times 10^{-10}$ | 2.789$\times 10^{-5}$
200 | 13.617 | 13.104 | 4.329 | 9.086$\times 10^{-10}$ | 1.500$\times 10^{-9}$ | 3.465$\times 10^{-4}$
250 | 13.904 | 14.483 | 3.205 | 6.915$\times 10^{-10}$ | 4.044$\times 10^{-10}$ | 0.003
Table 10: The $t$-value and $p$-value of the one-tailed paired $t$-test
between all other methods and our methods.

###  F.3 Resource Provisioning

Statistical significance tests. Table 11 shows the result of significance test
(one-tailed paired $t$-test), with the assumption that the means of two
distributions are the same.

|  $t$-value |  $p$-value
---|---|---
$\alpha_{1}/\alpha_{2}$ | L1 | L2 | Weighted L1 | L1 | L2 | Weighted L1
100 | 13.133 | 9.357 | 6.429 | 2.914$\times 10^{-9}$ | 2.115$\times 10^{-7}$ | 1.576$\times 10^{-5}$
10 | 12.484 | 14.231 | 10.897 | 5.622$\times 10^{-9}$ | 1.493$\times 10^{-6}$ | 3.190 $\times 10^{-8}$
0.1 | 1.184 | 5.073 | 11.001 | 0.127 | 8.486 $\times 10^{-5}$ | 1.416 $\times 10^{-8}$
0.01 | 0.535 | 1.998 | 5.083 | 0.300 | 0.032 | 8.33$\times 10^{-5}$
Table 11: The $t$-value ($p$-value) of the one-tailed paired $t$-test between
all other methods and our methods.

##  Appendix G Computing Infrastructure

All experiments are conducted on Linux Ubuntu 18.04 bionic servers with 256G
memory and 1.2T disk space with no GPU, for GPU does not suit well with
Gurobi.999see Gurobi’s official support website:
https://support.gurobi.com/hc/en-us/articles/360012237852-Does-Gurobi-support-
GPUs- Each server has 32 CPUs, which are Intel Xeon Platinum 8272CL @ 2.60GHz.

For the first experiment, we use a few minutes to get one set of data
101010Running all methods simultaneously under one particular parameter
setting. with training set size $100$, and about $4-5$ hours to get one set of
data with training set size $5000$. For the second experiment, we use about an
hour to get one set of data with $N=50$, and $2-3$ days to get one set of data
with $N=250$. For the third experiment, we use around $6-7$ hours to obtain
one set of data. Though our method is slower to train than two-stage methods,
it is 2-3x faster to train than KKT-based decision-focused method.