Differentiating Through Integer Linear Programs with Quadratic Regularization and Davis-Yin Splitting
Authors: Daniel McKenzie
Howard Heaton
Samy Wu Fung
Venue: Accepted by TMLR ()

Abstract: In many applications, a combinatorial problem must be repeatedly solved with similar, but distinct parameters. Yet, the parameters $w$ are not directly observed; only contextual data $d$ that correlates with $w$ is available. It is tempting to use a neural network to predict $w$ given $d$. However, training such a model requires reconciling the discrete nature of combinatorial optimization with the gradient-based frameworks used to train neural networks. We study the case where the problem in question is an Integer Linear Program (ILP). We propose applying a three-operator splitting technique, also known as Davis-Yin splitting (DYS), to the quadratically regularized continuous relaxation of the ILP. We prove that the resulting scheme is compatible with the recently introduced Jacobian-free backpropagation (JFB). Our experiments on two representative ILPs: the shortest path problem and the knapsack problem, demonstrate that this combination---DYS on the forward pass, JFB on the backward pass---yields a scheme which scales more effectively to high-dimensional problems than existing schemes.

Published in Transactions on Machine Learning Research (07/2024)
Differentiating Through Integer Linear
Programs with
Quadratic Regularization and Davis- Yin Splitting
Daniel McKenzie
dmccenniee@minnes.el
Deportment of Applied Mathematics and Statistics
Colorado School of Mines
Samy Wu Fung
sww.fureg@mmnnes.coo edu
Deportment of Applied Mathematics and Statistics
Department of Computer Science
Colorado School of Mines
Howard Heaton
research @typal academy
Typal Academy
Reviewed on OpenReview: //tp:////pperevviee net/ forun? Ho.rrrlW1SSSt
Abstract
In many applications, a combinatorial problem must be repeatodly solved with similar, but distinct parameters. Yet, the parameters w not directly observed; only contextual data d that correlates with w is available. It is tempting to use neural network to predict given d. However, training such model requires reconciling the discrete nature of com- binatorial optimization with the gradient-based frameworks used to train neural networks. We study the case where the problem in question an Integer Linear Program (ILP). We propose applying a three-operator splitting technique, also known as Davis- Yin split- ting (DYS), to the quadratically regularized continuous relaxation of the ILP. We prove that the resulting scheme compatible with the recently introduced Jacobian-free back- propagation (JFB). Our experiments lem and the knapsack problem, demonstrate that this combination DYS on the forward pass, JFB on the backward pass-yields a scheme which scales more effectively to high- dimensional problems than existing schemes. All code associated with this paper is available at /ectps://githin com/mines coom/inees-optml///ee opt-ml/fpo-dys
representative IT.Ps: the shortest path prob-
1
Introduction
Many high-stakes decision problems in healthcare (Zhngg & Tang 2021 logistics and scheduling (Siihi & Eglese 2010 Kacem et al. 2021), and transportation (Wang & Tang 2021 can be viewed as a two step process. In the first step, one gathers data about the situation at hand. This data is used to assign a value (or cost) to outcomes arising from each possible action. The second step is to select the action yielding maximum value (alternatively, lowest cost). Mathematically, this can be framed as an optimization problem with a data-dependent cost function:
x x²) arg min f(;;d) xEX
(1)
In this work, we focus on the case where X' C Rn is a finite constraint set and f(x; = w(dd To is a linear function. This class of problems is quite rich, containing the shortest path, traveling salesperson, and sequence alignment problems, name few. Given w(d), solving equation may be straightforward (e.g. shortest path) or Na-dard (e.g. traveling salesperson problem (Karp) 1972)). However, our present interest is settings where the dependence of w(d) on d unknown. In such settings, it is intuitive to learn mapping

---

Published in Transactions on Machine Learning Research (07/2022 2024)
20
80
10
10
Figure The shortest path prediction problem (Pogančić al., 2019h. The goal is to find the shortest path (from top-left to bottom-right) through randomly generated terrain map from the Warcraft II tileset Guyomarch)) The contextual data d, shown in (a), is an image sub-divided into 8-by-8 squares, each representing a vertex in 12-by-12 grid graph. The cost of traversing each square, d.e. w(d), is shown in (b), with darker shading representing lower cost. The true shortest path is shown in (c).
(d) ze w(d) and then solve
(d) arg min wes (e) =EX
in lieu of w(d). The observed data d called the contert. As an illustrative running example, consider the shortest path prediction problem shown in Figure 1 which is studied in Berthet et al. (2020) and Pogančić et al. (011)).
At first glance, it may appear natural to tune the weights θ so as to minimize the difference between we(d) and w(d). However, this only possible if w(d) is available at training time. Even if this approach is feasible, it may not be advisable. If (d)
near-perfect predictor w(d) we can expect (d( x x (d). However, if there are even small differences between (e(d) and w(d) this can manifest in wildly different solutions (Bengio| 1997 Wilder al.) 2019p. Thus, we focus on methods where θ is tuned by minimizing the discrepancy between ze
(e(d) (d and x² (d). Such methods are instances of the decision-focused learning paradigm Mandi al. (2023), as the criterion fidelity of we(d) (the "prediction").
optimizing for is the quality of (d) (the "decision") not the
The key obstacle in using gradient-based algorithms to tune in such approaches "differentiating through" the solution (d) equation to obtain gradient with which to update O. Specifically, the combinatorial nature of X' may cause the solution (d) to remain unchanged for many small perturbations to +; yet, for some perturbations (d) may "jump" to different point in X. Hence, the gradient dre /dwe always either zero undefined. To compute informative gradient, we follow Wilder al. (2019) by relaxing equation 2 to quadratic program over the convex hull of X with added regularizer (see equation 5).
Contribution Drawing upon recent advances in convex optimization Ryu & Yin 2022 and implicit neural networks (Fung et al. 2022 McKenzie et al. 2021), we propose method designed specifically for "differentiating through" large-scale integer lincar programs. Our approach is fast, casy to implement using our provided code, and, unlike several prior works (e.g. see Pogančić et al. (2019); Berthet al. (2020)) ) trains completely on GPU. Numerical examples herein demonstrate our approach, run using only standard computing resources, easily scales to problems with tens of thousands of variables. Theoretically, we verify our approach computes informative gradient via refined analysis of Jacobian-free Backpropagation
(JFB) (Fung et al. 2022 In summary, we do the following.
Building upon McKenzie 2021 we use the three operator splitting technique of Davis & Yin (2017) to propose DYS-Net.
2

---

Published in Transactions on Machine Learning Research (07/2024)
We provide theoretical guarantees for differentiating through the fixed point of certain non- expansive, but not contractive, operator.
We numerically show DYS-Net easily handles combinatorial problems with tens of thousands of variables.
2 Preliminaries
LP Reformulation We focus on optimization problems of the form equation 1 where f(;; d) = w(d) Tx and is the integer or binary points of a polytope, which we may assume to expressed in standard form Ziegler
2012
X =n= or = where C= fx I Ax=b and x 2 0}.
(3)
In other words, equation is an Integer Linear Program (ILP). Similar to works by Elmachtoub & Grigas 2022;; Wilder et al. (2019) Mandi Guns (2020) and others replace the model equation with its continuous relaxation, and redefine
re (d) _ arg min (d) x
(4)
sEC
step towards making the computation of dre/dwe feasible.
Losses and Training Data We assume access to training data the tuple form (d, (d)). Such data may be obtained by observing experienced agents solve equation I for example taxi drivers in shortest path problem Piorkowski et al. (2009). We focus on the & loss:
(()) A EAED ()) where (0; = (1' zo(d)]
(5)
and D is the distribution of contextual data, although more sophisticated losses can be used with only superficial changes to our results. In principle we select the optimal weights by solving arg ming (()). In practice, the population risk is inaccessible, and so we minimize empirical risk instead (Vapnik 1999):
arg min
1 N
Σ le (O:d).
(6)
Argmin Differentiation For notational brevity, we temporarily suppress the dependence on The gradient of la with respect O, evaluated on single data point is
d de [|xe z²11 (re x²)T dre dwe dwe de
[()(
de
As discussed in Section 1 zg piecewise constant function of we. and this remains true for the LP relaxation equation 4 Consequently, for all we, either are 000 0 or Oze/Oe /80 is undefined -eefineedneeihher case yields an informative gradient. To remedy this, Wilder et al. 22019 Mandi & Guns (2020) propose adding a small amount of regularization to the objective function in equation A to make strongly convex. This ensures ze is a continuously differentiable function wa. We add a small quadratic regularizer, modulated by y 2 0, to obtain the objective
fo(x; A we(d) T x + /
(7)
and henceforth replace equation by
(d A arg arg min fe(æ; (1 d.d sEC
(8)
During training, aim solve equation and efficiently compute the derivative Dre/00.
3

---

Published in Transactions on Machine Learning Research (07/2024)
Pitfalls of Relaxation Finally, offer a note of caution when using the quadratically regularized problem equation S as proxy for the original integer linear program equation
1 Temporarily defining
(o(d) _ arg min (e(d) T Te (d)
arg arg we we(d) Tx (d) _ arg min fe(z; Y,))
#EX
sEC
sEC
we observe that two sources of error have arisen; one from replacing X' with C (i.e. the discrepancy between že and te), and one from replacing the linear objective function with its quadratically regularized counterpart (Ree the discrepancy between Fie and æa). When A is tolally unimodular, Fe is guaranteed to be integral,
in which case èe že. Moreover, Wilder et al. 2019 Theorem 2) shows that
0 S we(d) (Fe(d) - -aa(d) S max || - yl12 z.yEC
whence, in the proof Berthet al. 2200 Proposition 2.2), lim,- xe(d) ( (at least when C is compact). Thus, for totally unimodular ILPs (e.g. the shortest path problem), it is reasonable to expect that (d( obtained from equation S for well-tuned @ and small y, are of acceptable quality. For ILPs which are not totally unimodular no such theoretical guarantees exist. However, we observe that this relaxation works reasonably well, at least for the knapsack problem.
3 Related Works
We highlight recent works in decision-focused learning particularly for mixed integer programs and related areas.
Optimization Layers Introduced in Amos Kolter (017)) and studied further in Agrawal al.) c2 Bolte et al.] (2021); Blondel et al. (2022), an optimization layer is a modular component of a deep learning framework where forward propagation consists of solving a parametrized optimization prob- lem. Consequently, backpropagation entails differentiating the solution to this problem with respect to the parameters. Optimization layers are a promising technology as they are able to incorporate hard constraints into deep neural networks. Moreover, as their output is a (constrained) minimizer of the objective function, it is casier to interpret than the output a generic layer.
Decision-focused learning for LPs When an optimization layer is fitted to a data-dependent optimiza- tion problems of the form equation [ with the goal of maximizing the quality of the predicted solution, this is usually referred to deciston-focused learning. We narrow our focus to the case where the objective function in equation is linear Wilder et al. (2019) and Elmachtoub & Grigas (2022) are two of the earli- est works applying deep lcarning techniques to data-driven LPs, and delincate two fundamentally different approaches. Specifically, WWilder et al. (2019) proposes replacing the LP with a continuous and strongly convex relaxation, as described in Section 2 This approach extended to ILPs in Ferber et al. (2020), and non-quadratic regularizers Mandi & Guns (2020). On the other hand, Elmachtoub & Grigas (2022) avoid modifying the underlying optimization problem and instead propose new loss function: dubbed the Smart Prediet-then-Opfiee Opthimii- (SPO) loss, to be used instead of the &2 loss. We emphasize that the SPO loss requires access to the true cost vectors w(d), setting which we do not consider.
Several works define a continuously differentiable proxy for the solution to the unregularized LP equation which we rewrite here as
x² (w) = arg min w xx sEc
(9)
In Berthet et al. 2020), stochastic perturbation is considered:
x (1 = Bz arg min (w + eZ) zec
(10)
which is somewhat analogous to Nesterov-Spokoiny smoothing (Nesterov & Spokoiny 2017) in zeroth-order optimization. For some draws of Z the additively perturbed cost vector w + E2 may have negative entries,
For sake of notational simplicity, the dependence implicit.
4

---

Published in Transactions on Machine Learning Research (07 (07/2024)
which may cause solver (e.g. Dijkstra's algorithm) applied to the associated LP to fail. To avoid this,
Dalle et al. (2022) suggest using multiplicative perturbation of w instead. Pogančié et. al. 2019 define a piecewise-affine interpolant to (w)), where E is a suitable loss function.
We note a bifurcation in the literature; Wilder et al.] (2019); Ferber et al. 2020 Elmachtoub & Grigas
2022) assume access to training data of the form (d, w())), whereas Pogančić et al. (2019): Berthet el al.
g020 : Sahoo et al. (2022) use training data of the form (d, (d,a"(d)) The former sctting is frequently referred
as the
problem. Our focus on the latter setting. We refer the reader to Kotary
et al.
2021;; Mandi et al. 2023 for recent surveys of this area.
Deep Equilibrium Models Bai ct al. 2019 El Chaoui et al. (2021) propose the notion of deep equilib-
rium model (DEQ), also known as an implicit neural network. A DEQ is neural network for which forward propagation consists of (approximately) computing a fixed point of parametrized operator. We note that equation may be reformulated as a fixed point problem,
Find ze such that ze = Pe (ag - aV- V=
(11)
where Pe is the orthogonal projection onto C. Thus, DEQ techniques may be applied to constrained optimization layers (Chen et al., 2021 Blondel et al. 2022 McKenzie et al. 2021 However, the cost of
computing Pe can be prohibitive, see the discussion in Section
Learning-to-Opee Optimize (L20) Our work is related to the bearning-lo-optining (I.20) framework (Chen et al. 2022 Li & Malik 2017 Chen et al. 2018 Liu et al. 2023 where an optimization algorithm is learned and its outputs are used perform inferences. However, traditional L20 methods use fixed number of layers (i.e. unrolled algorithms). Our approach attempts to differentiate through the solution z4 and is therefore most closely aligned with works that employ implicit networks within the L20 framework {Amos) & Kolter (2017: Heaton & Fung, 2023 McKenzie et al., 2021 Gilton et al. 2021 Liu et al. 2022). We highlight that, unlike many L20 works, our goal is not to learn faster optimization method for fully
specified problems. Rather, we seek to solve partially specified problems given contextual data by combining learning and optimization techniques.
Computing the derivative of a minimizer with respect to parameters In all of the aforementioned works, the same fundamental problem encountered: Ore/00 must be computed where ze is the solution of (constrained) optimization problem with objective function parametrized by θ. The most common approach to doing so, proposed in Amos & Kolter (011)) and used in Ruthotto et al. (2018); Wilder et al.
2019 Mandi & Guns 2020 Ferber et al. 2020) starts with the KKT conditions for constrained optimality:
afe (re) + A x + e =0, 0, Az b =0 0, D(D) O,
Or
where Å and & 2 O are Lagrange multipliers associated to the optimal solution ze (erttskas 1997] and D(0) is matrix with & along its diagonal. Differentiating these equations with respect to Q and rearranging yields
A
dza de di de de de
0²
fa
(12)
D(D) D(xe)
The matrix and right hand side vector in equation 12 are computable, thus enabling one to solve for das de (as well as and de). We emphasize that both primal (i.e., æe) and dual (i.e., Å and D) variables at di
optimality are required. This constrains the choice of optimization schemes to those that track primal and dual variables, and prohibits the use of fast, primal-only schemes. Following (Amos & Kolter] 2017). most works employing this approach use a primal-dual interior point method, which compules acceptable approximations to Ie, Å, and at a cost of O (max(n, m) assuming x e R² and A E Raxn
Using the stochastic proxy equation 10 Berthet et al. 20200) derive a formula for da//dw which is also an
For C R., (a( A migmin...o E. || =l]-
5

---

Published in Transactions on Machine Learning Research (07 (07/2024)
expectation, and hence can be efficiently approximated using Monte Carlo methods. Pogančié et al. (0119)
show that the gradients of their interpolant are strikingly easy to compute, requiring just one additiona solve of equation 9 with perturbed cost w'. This approach is extended by Sahoo et al. (2022) which proposes to avoid computing /dw entirely, replacing it with the negative identity matrix. This is similar in spirit to our use of Jacobian-Iree Backpropagation (see Section 4
4 DYS-Net
We now introduce our proposed model, DYS-net. We use this term to refer to the model and the training procedure. Fixing an architecture for we, and an input d, DYS-net computes an approximation to re(d) way that is easy to backpropagate through:
DYs-net(d; Q) ~ re
2 arg min Se (e(r;,,))
(13)
#EC
The Forward Pass As we wish compute ze and 0re/00 for high dimensional settings (i.e. large n): eschew second-order methods (e.g. Newton's method) in favor of first-order methods such as projected gradient descent (PGD). With PGD, sequence {a*} of estimates of ze are computed so that ne lima- _ where
x+t = Pe (r aVaf(=" for all KEN.
This approach works for simple sets C for which there exists an explicit form of Pc, e.g. when C is the probability simplex Duchi et al. 2008 Condat 2016 Ti et al. 2023 However, for general polytopes no such form exists, thereby requiring second iterative procedure to compute Pe(x*). This projection must be computed at every iteration of every forward pass for every sample of every epoch, and this dominates the computational cost, see McKenzie et al. 2021 Appendix 4).
To avoid this expense, we draw on recent advances in the convex optimization literature, particularly Davis- Yin splitting or DYS (Davis & Yin 2017). DYS three operator splitting technique, and hence is able to decouple the polytope constraints into relatively simple constraints Davis & Yin (2017). This technique has been used elsewhere, e.g. Pedregosa & Gidel (2018); Yurtsever et al. (2021). ). Specifically, adapt the
architecture incorporating DYS proposed in McKenzie et al. (001)). To this end, we rewrite C as
= fx: Ar = b and x 2 0} = {x: Ax b) n fx: x x20) = =c
(14)
C
C
While Pe is hard to compute, both Po. and Pc2 can be computed cheaply via explicit formulas, once a singular value decomposition (SVD) computed for A. We verify this via the following lemma (included for completeness, the two results are already known).
Lemma 1. C1, are as in equation 14 and A is full-rank then:
=2) = - - b), where Al = VT and UEVT is the compact SVD of A.
2. (2) RcLU((2) A m0 2).
We remark that the SVD of A is computed offline, and is once-off cost. Alternatively, one could compute Pc, by using an iterative method (e.g. GMRES) to find the (least squares) solution to Av = Az 1---b b, whence (2) = z-u, but this would need to be done at each iteration of each forward pass for every sample of every epoch. Thus, it is unclear whether this procides computational savings over the once-off computation of the SVD. However, this could prove useful when only provided implicitly, i.e, via access to a matrix-vector product oracle.
The following theorem allows to approximate se using only Pe, and Pc,; not Pc-

---

Published in Transactions on Machine Learning Research (07/2024)
Theorem 2. Let C1. C2 be as in equation 14 and suppose fe(x; we(d) =111=112 for any neural
network we we(d). (d). For any E (0,2/7) define the sequence {=*} by:
= =T(2z*) for all =
(15)
where
x Pc, ((2 (2- (2) -
(16)
If t & Pc, (2 then the sequence { converges to re in equation 8 and || - -01/)))
A full proof is included in Appendix A here we provide proof sketch. Substituting the gradient
V v=Je(Ezn,d)) we +
for Fe(. d) into the formula for (-(-) given in McKenzic et al. [2021 Theorem 3.2) and rearranging yields
equation 16
Verifying the assumptions of McKenzie et al. (2021 Theorem 3.3) are satisfied is straight-
forward. For example, Vzfe(z; d) is 1/y-cocoercive by the Baillon- Haddla theorem (Baillon & Haddad 1977 Bauschke & Combettes 2009
The forward pass of DYS-net consists of iterating equation 15 until suitable stopping criterion met. For simplicity, in presenting Algorithm 1 we assume this to be maximum number of iterations is reached. One may also consider stopping criterion of ||x++1 z < tol.
The Backward Pass
From Theorem
we deduce the following fixed point condition:
re = Pc,( for some =o E (=: T-()))
(17)
As discussed McKenzie et 2021 instead of backpropagating through every iterate of the forward pass, may derive formula for dre/ de by appealing to the implicit function theorem and differentiating both sides of equation 17
dze aTe
aT dze
dze 8Te where Je(z) = I- aTe 0z
+
Te(ze)
(18)
de
00
Oz de
dO ae
We notice two immediate problems not addressod by McKenzie et al. 2001 p: (i) To not everywhere
differentiable with respect to 2, as Pc, is not; (ii) if Te were a contraction (i.e. Lipschitz in & with constant less than 1), then To would be invertible. However, this is not necessarily the case. Thus, it is not clear priori that equation 18 can be solved for dze/dO. Our key result (Theorem below) is to provide reasonable conditions under which Je(ze) is invertible.
Assuming these issues can be resolved, one may compute the gradient of the loss using the chain rule:
de dr de dr dz de = de dPc. To oTe
de dl dre de (dec, dze
dæ de
d0
This approach requires solving linear system with Te which becomes particularly expensive when is large. Instead, use Jacobion-free Backpropagation (JFB) (Fung et al., 2022), also known one-step differentiation (Bolte et al. 2024). in which the Jacobian Je is replaced with the identity matrix. This leads to an approximation of the true gradient dé/d@ using
0e dPc, ÜT
(e(xo)
(19)
dr dz 0e
(x.s)=((e.=e)
This update can be seen as a zeroth order Neumann expansion of J Fung et al. (2022) Bolte et al. (2024)
We show equation 10 valid descent direction by resolving the two problems highlighted above. Wo begin by rigorously deriving a formula for 0Ta/Oz Recall the following generalization of the Jacobian non-smooth operators due to Clarke (1983).

---

Published in Transactions on Machine Learning Research (07/2024)
Algorithm DYS-Net
1: Inputs:A and b defining C, hyperparamcteers a, K
2: Initialize: Attach neural network we(+). Compute SVD of A for Pc, formula
3: Forward Pass: Initialize zº. Given we(d) compute z',...., using equation 15 Return x²
Pe, (2*) 2 ze.
Backward Pass: Compute p=(x*) pe(ae) using equation 19 and return.
Definition For any locally Lipschitz F; 12 R" let Dy denote the set upon which F is differentiable. The Clarke Jacobian of F is the set-valued function defined as
= 0F(3)={ Con {lima-
if 2 DF if 2 DF
0F(
==
dF da
Where Con { denotes the conver hull of a set.
The Clarke Jacobian of Pc2 is easily computable, see Lemma
Define the (sel-valued) functions
if a >
and ā(a) {1
f a S
d(a) = mmxx(0,,))
if a < 0
[0,1] if a 0
Then
aPc. (z)
dz ReLU(x)
diag(c(Z)),
(20)
where e is applied element-wise. We shall now provide a consistent rule for selecting, at any ž, an element of PP), Tf 2 # 0 for all then OPc, is a singleton. one or more 2; = 0 then apc, is multi-valued, so we choose the element of OPc, with in the (i,i) position for every 2; = 0. We write dPc, for this chosen element, and note that
dPc.
diag(e(T)) e (2)
dz ===
This aligns with the default rule for assigning sub-gradient to ReLU used in the popular machine learning libraries TensorFlow AaasorFlow(Abad)) et al. 2016), PyTorch (Paszke et al.] 2019 and JAX (Bradbury et al.] 2018)) and has been observed to yicld networks which are more stable to train than other choices (Bertoin et al. 2021
Given the above convention, we can compute 0Te/O2, interpreted as an element of the Clarke Jacobian of Te chosen according to consistent rulc. Surprisingly, 0Te/Oz may be expressed using only orthogonal projections to hyperplanes. Throughout, we let e; E R² be the one-hot vector with in the i-th position and zeros elsewhere, and a. be the i-th row of A. For any subspace H. we denote the orthogonal subspace as H-- The following theorem, the proof of which is given in Appendix A characterizes 0Te/Oz.
Theorem 4. Suppose A is full-rank and Hi A Null(A), H2.z Span (c: :4 0). Then for all & E R.
öTe 0z
Pat Put + (1 an) Pn, PHiPHa..
(21)
==
To show JFB is applicable, it suffices to verify ||0e/001| when evaluated ze. Theorem eeaalles to show this inequality holds when ze satisfies commonly-used regularity condition, which formalize follows.
Definition (LICQ condition, specialized to case). Let ze denote the solution to equation A(ze) C (I. denole the set of active posilivily construints:
Let
A(re) & (i: [re] 0).
(22)
8

---

Published in Transactions on Machine Learning Research (07/2024)
The point ze satisfies the Linear Independence Constraint Qualification (LICQ) condition if the vectors
{ai, (a1 am) U fe:::e:e: A(ro))
(23)
are linearly independent.
Theorem 6. If the LICQ condition holds at ze, is full-rank and a (0,2/7), then ||0Ta/02||==se
The significance of Theorem d which is proved in Appendix A is outlined by the following corollary, stating that using pe instead of the true gradient dl/d@ still guaranteed to decrease {(0); least for small enough step-size.
Corollary 7. 1f Te is continuously differentiable with respect to θ at ze, the assumptions in Theorem hold and (dTe/00) (OTa/00) has condition number sufficiently small, then pe (=( a descent direction for l with respect to O.
Theorem alsso proves that Je(zg) in equation 18 invertible, thus also justifying the use of Jacobian-basee backprop as well as numerous other gradient approximation techniques (Liao ct al. 2018 Geng et al. 2021
Finally, we note that in practice (=(**) is used instead of pe(re)- Fung et. al. (2022 Corollary 0.1)
guarantees that p=(=*) will still be descent direction for K sufficiently large, see also Bolte et al. 2024 Cor. 1).
Synergy between forward and backward pass We highlight that the pseudogradient computed in the backward pass (i.c., pe) does not require any Lagrange multipliers, in contrast with the gradient computed using equation 12 This is why we may use Davis Yin splitting, as opposed a primal-dual interior point method, on the forward pass. This allows efficiency gains the forward pass, Davis Vinis-Yin splitting scales favorably with the problem dimension.
Implementation DYS-net is implemented an abstract PyTorch (Paskke et al. 2019 layer in the pro- vided code. To instantiate it, a user only need specify A and b (see equation 3) and choose architecture for we- At test time, it can be beneficial to solve the underlying combinatorial problem equation / exactly using a specialized algorithm, e.g. commercial integer programming software for the knapsack prediction problem. This can casily be done using DYS-net by instantiating the best-time-forvlte abstract method.
5 Numerical Experiments
We demonstrate the effectiveness of DYS-Net through four experiments. In all our experiments with DYS-net we used a 0.05 and K 100, with an early stopping condition of | - ze/2 S tol with tol 0.01. We compare DYS-net against three other methods: Perturbed Optimization (Berthet et all 2020 PertOpt-net;
approach using Blackbox Backpropagation (Vlastelica et al.] 2019 BB BB-net and
approach using cvxpylayers (Agrawal et al. 2019a) CVx-net. All code needed to reproduce our exper- iments is available at /ttts:///i tttgggtthuub.... www.coo.mmnnneeeee.. opt http://nnssooottm/oooo/-0----
5.1 Experiments using PyEPO
First, verify the efficioncy of DYS-net by performing benchmarking experiments within the PyEPO TTng
& Khalil 2022 framework.
Models We use the PyEPD implementations of Perturbed Optimization, respectively Blackbox Backprop- agation, as the core of PertOpt-nee respectively BB-net We use the same architecture for we (d) for all model and use an appropriate combinatorial solver at test time, via the PyEPO interface to Gurobi. Thus, the methods differ only in training procedure. The precise architectures used for (e) for each problem are described in Appendix and are also viewable in the accompanying code repository. Interestingly,
TTee architecture we shortest path prediction problems differs from that used for knapsack prediction problems. See the appendix for details.
9

---

Published in Transactions on Machine Learning Research (07/2024)
observed that adding drop-out during training to the output layer proved beneficial for the knapsack prediction problem. Without this, we find that we tends to output sparse approximation to supported on a feasible set of items, and so does not generalize well.
Training We train for a maximum of 100 epochs or 30 minutes, whichever comes first. We use a validation set for model selection as we observe that, for all methods, the best loss is seldom achieved at the final iteration. Exact hyperparameters are discussed in Appendix [ For each problem size, we perform three repeats with different randomly generated data. Results are displayed in Figure 2
Evaluation Given test data {(di, w(d) x² (d:)) and fully trained model we(d) we define the regret: r(d) ((:) (re(d;) - (d:)).
(24)
Note that regret is non-negalive, and is zero precisely when æe (d)) is a solution of equal quality to x (d:). Following Tang & Khalil 2022), we evaluate the quality of (e(d) using normalized regret, defined as
E r(d))
F
(25)
1i=1 [uo(d,) TT (d,))
5.1.1 Shortest Path Prediction
The shortest path between two vertices in a graph 9 = EV, can be found by:
Tx; x arg min wdd x; X' = {ar E {0, 11 E1 Ex = b),
(26)
ZEN
where E is the vertex-edge adjacency matrix, b encodes the initial and terminal vertices, and w(() E R is a vector encoding (d-dependent) edge lengths. Here, xt will encode the edges included in the optimal path. In this experiment we focus on the case where 9 is the & x grid graph. We use the PyEPO Tang & Khalil (2022) benchmarking software to generate training datasets of the form {(d, R x' (d:)) where N 1000 for each E (5, 10, 15, 20, 30). The d, are sampled from the five-dimensional multivariate Gaussian distribution with mean 0 and identity covariance. Note that the number of variables in equation 26 scales like k², not k.
5.1.2 Knapsack Prediction
In the (0-1, single) knapsack prediction problem, we are presented with a container (i.e. knapsack) of size and I items, of sizes 81,..., 87 and values wi wr)()),, The goal is to select the subset of maximum value that fits in the container, t.e. to solve:
xx ara max w(d) x; X E (0.13' sT
ZEX
Here, x² encodes the selected items. In the (0-1) k-knapsack prediction problem we imaginc the container having various notions of "size" (i.e. length, volume, weight limit) and hence a k-tuple of capacities E R* Correspondingly, the items each have a k-tuple of sizes $1,...,87 E Rt. We aim to select a subset of maximum value, amongst all subsets satisfying the capacity constraints:
x arg max w(d) REX
X = {a E (0,1)' Sr S c}
(27)
S = [81
ss] E R kx1
In Appendix B we discuss how to transform X' into the canonical form discussed in Section 2 We again use PyEPO to generate datasets of the form {(d,, xt Rs x² (d,))) where N 1000 for and varying from 50 750 inclusive in increments of 50. The d, from the five-dimensional multivariate Gaussian distribution with mean 0 and identity covariance.
10

---

Published in Transactions on Machine Learning Research (07 (024)
0.9
0.8
10a
PertOpt
CVX
0.4
PertOpt
10a
CVX
0.7
0.6
BB
DYS
0.3
BB
0.5
10a
BB
10a
0.4
PertOpt
0.2
0.3
CVX
PertOpt
10
BB
0.1
10a
0.1
DYS
CVX
10
20
30
20
400
600
800
200
400 600
800
Grid Size
Grid Size
Number of Items
Number of Items
(a) Normalized Regret
(b) Train Time (in seconds)
(c) Normalized Regret
(d) Train Time (in seconds)
Figure Results for the shortest path and knapsack prediction problems. Figures (a) and (b) show nor- malized regret and train time for the shortest path prediction problem, while Figures (c) and (d) show normalized regret and train time for the knapsack prediction problem. Note that the train time in figures (b) and (d) is the time till the model achieving best normalized regret the validation set is reached.
100%
DYS
80%
CVX
60%
PertOpt
40%
BB
20%
0%
20
40
60
80
100
Epochs
Figure 3: Accuracy (in percentage) of predicted paths on 5-by-5 grid during training.
5.1.3 Results
2 DYS-net trains the fastest among the four methods, and achieves the best or near-best normalized regret in all cases. We attribute the accuracy to the fact that is not a zeroth order method. Moreover, avoiding computation and backpropagation through the projection Pe allows DYS to have fast runtimes.
As is clear from Figure
5.2 Large-scale shortest path prediction
Our shortest path prediction experiment described in Section 5.1.1 is bottlenecked by the fact that the "academie" license of Gurobi does not allow for shortest, path prediction problems on grids larger than x0 To further explore the limits of DYS-net, we redo this experiment using a custom pyTorch implementation of Dijkstra's algorithm" as our base solver for equation 26
Data Generation We generate datasets D = ((d,, x² 1,000 for k-by-k grids where E (5, 10, 20, 30, 50, 100} and the d, are sampled uniformly at random from [o, 115, the true edge weights are computed as wr(d) = Wd for fixed W E R č|x5, and (d) computed given w(d) using the aforementioned pyTorch implementation of Dijkstra's algorithm. Thus, all entries of w(d) are non-negative, and are positive with overwhelming probability.
"adapted from Tensorflox code available peruuraaiiii///ggrrriaattttttttr
11

---

Published in Transactions on Machine Learning Research (07 (07/202) /224)
grid size
number of variables
neural network size
5-by-3
40
500
10-by-10
180
2040
20-by-20
760
8420
30-by-30
1740
19200
50-by-50
4900
53960
100-by-100
19800
217860
Table 1: Number of variables (í.c. number of odges) per grid size for the shortest path prediction problem of Section 5.2 Third column: number of parameters in (e(d) for DYS-Net, CVX-net and PertOpt-net. For BB BB-Net, we found latent dimension that is 20-times larger than the aforementioned three to be more effective.
10a
PertOpt
10a
BB
DYS
PertOpt
10
CVX
BB
10-1
BB
PertOpt
10a
10-1
CVX
DYS
CVX
DYS
10-3
10
10a
50
100
50
100
50 Grid Size
Grid Size
Grid Size
Test MSE Loss
(b) Training Time (min)
(c) Regret Values
Figure 4: Results for the shortest path prediction problem. a) Test MSE loss (left), b) training time in minutes (middle), and c) regret values (right) vs. gridsize for DYS-Net (proposed) and approaches using
cvxpylayers (Agrawal et al. (019b) labeled CVX; Perturbed Optimization (Berthet et al. 2020) labeled PertOpt; and Blackbox Backpropagation (Vlastelica et al. 2019) labeled BB. Note CVX unable to load or run problems with gridsize over 30. Dimensions of the variables can be found in Table
12

---

Published in Transactions on Machine Learning Research (07/2024)
Models and Training We test the same four approaches as in Sections 5.1.1 and 5.1.2 but unlike in these sections we do not use the PyEPO implementations of PertOpt-net or BB-net. We cannot, as the PyEPO implementations call Gurobi to solve equation 26 which, as mentioned above, cannot handle grids larger than 30 x 30. Instead, we use custom implementations based on the code accompanying the works introducing PertOpt (Berthet et al. 2020) and BB (Vlastelica et al.) 2019). We use the same architecture for (e(d) for DYS-net, PertOpt-net, and Cvx-net; a two layer fully connected neural network with leaky ReLU activation functions. For DYS-net and Cvx-net we do not modify the forward pass at test time. For BB-net we use larger network by making the latent dimension 20-times larger than that of the first three as we found this more effective. Network sizes can be seen in Table [
We tuned the hyperparameters for each approach to the best of our ability on the smallest problem (5-by-5 grid) and then used these hyperparameter values for all other graph sizes. See Figure 3 for the results of this training run. We train all approaches for 100 epochs total on each problem using the lz loss equation
1 While Cvx-net and PertOpt-ner net achieve low regret for small grids, DYS-net model achieves low regret for all grids. In addition to training faster, DYS-net can also be trained for much larger problems, e.g., 100-by-100 grids, as shown in Figure 4 We found that CVX-net could not handle grids larger than 30-by-30, f.e. problems with more than 1740 variables (see Table D. Importantly, PertOpt-net takes close to a week to train for the 100-by-100 problem, whereas DYS-net takes about a day (see right Figure A On the other hand, the training speed of BB-net is comparable to that of DYS-net, but does not lead to competitive accuracy as shown in Figure 4a). Interpreting the outputs of DYS-net and CVX-net as (unnormalized) probabilities over the grid, one can use a greedy algorithm to determine the most probable path from top-left to bottom-right. For small grids, e.g. 5-by-5, this path coincides exactly with the true path for most d (sce Fig. 3
Results The results are displayed in Figure
5.3 Warcraft shortest path prediction
Finally, as an illustrative example, we consider the Warcraft terrains dataset first studied in Vlastelica et al. (2019h. As shown in Figure 1 tile has a different cost to traverse, and the goal is to find the quickest path from the top-left corner to the bottom-right, corner.
a 96-by-96 RGB image, divided into a 12-by-12 grid of terrain tiles. Each
Models, Training, and Evaluation We build upon the code provided as part of the PyEPO package (Tang & Khalil, 2022). We use the same architecture for (e) for BB-net, PertOpt-net, and DYS-net- truncated Res, 1e (ee et al. 2016) as first proposed in Vlastelica et al. (2019). We train each network for 50 epochs, except for the bascline. The initial learning rate is 5 x 10 and it is decreased by a factor of 10 after 30 and 40 epochs respectively. PertOpt-net and DYS-net are trained to minimize the & loss equation 5 while BB-net is trained minimize the Hamming loss as described in Vlastelica et al. (0019) We evaluate the quality of (dl using normalized regret equation 25
Results The results for this experiment are shown in Table 2 Interestingly, BB-net and PertOpt are much more competitive in this experiment than in the experiments of Sections 5.1.1 5.1.2 and 5.2 We attribute this to the "discrete" nature of the true cost vector w(d) for the Warcraft problem entries of w(d) can only take on four, well-scparated values-aa as opposed to the more "contimuuus" nature of w(d) in the previous experiments. Thus, an algorithm that learns a rough approximation to the true cost vector will be more successful in the Warcraft problem than in the shortest path problems of Scctions 5.1.1 and 5.2 This difference is illustrated in Figure 5 We also note that model selection using a validation set is beneficial for all approaches, but particularly so for BB-net and DYs-net, which achieve their best-performing models in the first several epochs of training.
5rhii is to expected, as discussed Amos Kolter 2017
Agrawal et al. 2019a
13

---

Published in Transactions on Machine Learning Research (07 (07/202) 2024)
Figure 5: Left two figures: Sample cost matrices for shortest problem considered in Section 5.1.1 Right two figures: Sample cost matrices for the Warcraft shortest path prediction problem considered in Sec- tion 5.3 Note that Section 5.3 the node weighted shortest path prediction problem is considered, while in Section 5.1.1 the edge weighted variant is solved. For ease of comparison, in the left two figures we have reshaped the edge cost vector inlo a node cost. matrix.
Algorithm
Test Normalized Regret
Time (in hours)
BB-net
0.1204
0.11
PertOpt-net
0.1089
2.07
DYS-net
0.0889
0.09
CVX-net
0.0214
2.33
Table 2: Results for the 12-by-12 Warcraft shortest, path prediction problem. We select the model achieving best normalized regret on the validation set. The time displayed is the time till this best normalized regret is achieved. All results are averaged over three runs.
6 Conclusions
This work presents a new method for learning to solve ILPs using Davis-Yin splitting which we call DYS-net. We prove that the gradient approximation computed the backward pass of DYS-net indeed a descent di- rection, thus advancing the current understanding of implicit networks. Our experiments show that DYS-net is capable of scaling to truly large ILPs, and outperforms existing state-of-the-art methods, at least when training data is of the form (d, (d)). We note that in principle DYS-net may be applied given training data of the form (d, w(d); common paradigm for many trediit-then-pptimizz problems Mandi et al. (2023) However, preliminary experiments showed that SPO+ Ellachhoub & Grigas (2022) achieved a substantially
lower regret than DYS-net in this setting. We leave the extension of DYS-net to this data type to future work. Our experiments also reveal an interesting dichotomy between problems in which entries of w(d) may take on only handful of discrete values and problems which w(d) more "continuous". Future work could explore this dichotomy further, as well apply DYS-net to additional ILPs, for example the traveling salesman problem.
References
Martin Abadi, Paul Barham, Jianmin Chen, Zhifeng Chen, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Geoffrey Irving, Michael Isard, et al. Tensorflow: system for large-scale machine lcarning. In 12th USENIX symposium on operating systems design and implementation (OSDI 16), pp. 265-283, 2016.
A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and Z. Kolter. Differentiable convex optimization layers. In Advances in Newral Informalion Processing Systems, 2019a. 9a..
Akshay Agrawal, Brandon Amos, Shane Barratt, Stephen Boyd, Steven Diamond, and J Zico Kolter. Dif- ferentiable convex optimization layers. Advonces in neural information processing systems, 32, 2019b.
Akshay Agrawal, Shane Barratt, Stephen Boyd, Enzo Busseti, and Walaa M Moursi. Differentiating through cone program. Xiv preprint or
2019c.
14

---

Published in Transactions on Machine Learning Research (07/2024)
Brandon Amos and J Zico Kolter. Optnet: Differentiable optimization layer in neural networks. In International Conference on Machine Learning, pp. 136-145. PML.R, 2017.
Shaojie Bai, Zico Kolter, and Vladlen Koltun. Deep equilibrium models. Advances in Neural Information Processing Systems, 32, 2019.
Jean-Bernard Baillon and Georges Haddad. Quelques propriétés des opérateurs angle-bornés et cycliquement monotones. Israel Journal of Mathematics, 26:137 150, 1977.
Heinz H Bauschke and Patrick L Combettes. arXiv:0906.0 0807, 2009.
The
theorem revisited. arXiv preprint
Yoshua Bengio. Using a financial training criterion rather than a prediction criterion. Internalional journal of neural systems, 8(04):433 443, 1997.
Quentin Berthet, Mathieu Blondel, Olivier Teboul, Marco Cuturi, Jean-Philippe Vert, and Francis Bach. Learning with differentiable perturbed optimizers. Advances in neural information processing systems, 33: 9508-9519, 2020.
David Bertoin, Jérôme Bolte, Sébastien Gerchinovitz, and Edouard Pauwels. Numerical influence of RcLU on backpropagation. Advances in Newral Information Processing Systems, 34:468 114468-77,, 2021.
Dimitri P Bertsekas. Nonlinear programming. Journal of the Operational Research Society, 48(3):334 334, 1997.
Mathieu Blondel, Quentin Berthet, Marco Cuturi, Roy Frostig, Stephan Hoyer, Felipe Llinares-López, Fabian Pedregosa, and Jean-Philippe Vert. Efficient and modular implicit differentiation. Advances in neural information processing systems, 35:5230-5242 2022.
Jérôme Bolte, Tam Le, Edouard Pauwels, and Tony Silveti-Falls. Nonsmooth implicit differentiation for machinc-learning and optimization. Advances in nearal information processing systems, 34:13537-13599, 2021.
Jérôme Bolte, Edouard Pauwels, and Samuel Vaiter. One-step differentiation of iterative algorithms. Ad- vances in Neural Information Processing Systems, 36, 2024.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Mihrn et al. Jax: composable trans- formations of python+ numpy programs. 2018.
Bingqing Chen, Priya L Donti, Kyri Baker, Zico Kolter, and Mario Bergés. Enforcing policy leasibility constraints through differentiable projection for energy optimization. In Proceedings of the Twelfth ACM International Conference on Future Energy Systems, pp. 199 210, 2021.
Tianlong Chen, Xiaohan Chen, Wuyang Chen, Howard Heaton, Jialin Liu, Zhangyang Wang, and Wotao Yin. Learning optimize: primer and benchmark. Journal of Machine Learning Research, 23(189): 1 1-59, 2022.
Xiaohan Chen, Jialin Liu, Zhangyang Wang, and Wotao Yin. Theoretical linear convergence of unfolded ista and its practical weights and thresholds. Advances in Neural Information Processing Systems, 31, 2018.
FH Clarke. Optimization and nonsmooth analysis, wiley-interscience. New York, 1983.
Laurent Condat. Fast projection onto the simplex and the é ball. Mathematical Programming, 158(1): 575-585, 2016.
Guillaume Dalle, Léo Baty, Louis Bouvier, and Axel Parmentier. Learning with combinatorial optimization layers: probabilistic approach. arXív preprint arXto:2207.1351 2022.
Damek Davis and Wotao Yin. A three-operator splitting scheme and its optimization applications. Set-valuek and variational analysis, 25(4):829 858, 2017.
15

---

Published in Transactions on Machine Learning Research (07/2024)
John Duchi, Shai Shalev-Shward Yoram Singer, and Tushar Chandra. Efficient projections onto the f1-ball for learning in high dimensions. In Proceedings of the 25th international conference on Machine learning, pp. 272 279, 2008.
Laurent EI Ghaoui, Fangda Gu, Bertrand Travacca, Armin Askari, and Alicia Tsai. Implicit deep learning. SIAM Journal on Mathematics of Data Science, 333):930-958, 2021.
Adam Elmachtoub and Paul Grigas. Smart "predict, then optimize". Management Science, 68(1):9 26, 2022.
Aaron Ferber, Bryan Wilder, Bistra Dilkina, and Milind Tambe. Mipaal: Mixed integer program as layer. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 34, pp. 1504 1511, 2020.
Samy Wu Fung, Howard Heaton, Qiuwei Li, Daniel McKenzie, Stanley Osher, and Wotao Yin. JFB: Jacobian- free backpropagation for implicit networks. In Proceedings of the AAAI Conference on Artificial Intelli- gence, 2022.
Zhengyang Geng, Xin- Yin-Yu Zhang, Shaojic Bai, Yisen Wang, and Zhouchen Lin. On training implicit models. Advances in Neural Information Processing Systems, 0224447--4260 2021.
Davis Gilton, Gregory Ongie, and Rebecca Willett. Deep equilibrium architectures for inverse problems in imaging. JEEE Transactions on Computational Imaging, 7:1123-111 1133, 2021.
Jean Guyomarch. Warcraft ii open-source map editor, 2017. hilp://gilhub. com/war2/ 2/war2edi..
Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference computer vision and pattern recognition, pp. 770-778, 2016.
Howard Heaton and Samy Wu Fung. Explainable ai via learning optimize. Scientific Reports, 13(1::00103, 2023.
Imed Kacem, Hans Kellerer, and A Ridha Mahjoub. Preface: New trends on combinatorial optimization for network and logistical applications. Annals of Operations Research, 298(1):1-55 1,5, 2021.
Richard M Karp. Reducibility among combinatorial problems. In Complexity of computer compulations, pp- 85 103. Springer, 1972.
James Kotary, Ferdinando Fioretto, Pascal van Hentenryck, and Bryan Wilder. End-to-end constrained optimization learning; A survey, In 30th International Joint Conference on Artificial Intelligence, IJCAI 2021, pp. 4475-4482, International Joint Conferences on Artificial Intelligence, 2021.
Ke Li and Jitendra Malik. Learning to optimize. In International Conference on Learning Representations. 2017.
Qiuwei Li, Daniel McKenzie, and Wotao Yin. From the simplex to the sphere: faster constrained optimization using the hadamard parametrization. Information and Inference: A Journal of the IMA, 12(3):1898 1937,, 2023.
Renjie Liao, Yuwen Xiong, Ethan Fetaya, Lisa Zhang, KiJung Yoon, Xaq Pitkow, Raquel Urtasun, and Richard Zemel. Reviving and improving recurrent back-propagation. In International Conference. Machine Learning, pp. 3082-3091. PMLR, 2018.
Jialin Liu, Xiaohan Chen, Zhangyang Wang, Wotao Yin, and HanQin Cai. Towards constituting mathemati- cal structures for learning to optimize. In International Conference on Machine Learning, pp. 21426 21449.
PMLR, 2023.
Jiaming Liu, Xiaojian Xu, Weijie Gan, Ulugbek Kamilov, et al. Online deep equilibrium learning for regu-
larization by denoising. Advances in Neurol Information Processing Systems, 35:25363-283 -56676, 2022.
Jayanta Mandi and Tias Guns. Interior point solving for lp- based prediction+ optimisation. Advances in Neural Information Processing Systems, 33:7272-7265 2020.
16

---

Published in Transactions on Machine Learning Research (07/2024)
Jayanta Mandi, James Kotary, Senne Berden, Maxime Mulamba, Victor Bucarey, Tias Guns, and Ferdinando Fioretto. Decision-focused learning: Foundations, state of the art, benchmark and future opportunities. arXiv preprint arXic:2307.13563 2023.
Daniel McKenzie, Howard Heaton, Qiuwei Li, Samy Wu Fung, Stanley Osher, and Wotao Yin. Operator splitting for learning predict equilibria in convex games. arXiv e-prints, pp. arXiv 2106, 2021.
Yurii Nesterov and Vladimir Spokoiny. Random gradient-free minimization of convex functions. Foundations of Computational Mathematics, 17:527-566, 2017.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al. Pytorch: An imperative style, high-performance deep learning library. Advances in Neural Information Processing Systems, 32, 2019.
Fabian Pedregosa and Gauthier Gidel. Adaptive three operator splitting. In International Conference on Machine Learning, pp. 4085-4094. PMLR, 2018.
Michal Piorkowski, Natasa Sarafjanovic-Djukii,, and Matthias Grossglauser. Crawdad data set pofl/mooility (v. 2009-02-24), 2009.
Marin Vlastelica Pogančté, Anselm Paulus, Vit Musil, Georg Martius, and Michal Rolinek. Differentiation of blackbox combinatorial solvers. In International Conference Learning Representations, 2019.
Lars Ruthotto, Julianne Chung, and Matthias Chung, Optimal experimental design for inverse problems with state constraints. STAM Journal on Scientific Computing, 40(4):B1000-B1100, 2018.
Ernest Ryu and Wotao Yin. Large-Scale Convex Optimization: Algorithm Designs via Monotone Operators. Cambridge University Press, Cambridge, England, 2022.
Subham Sckhar Sahoo, Anselm Paulus, Marin Vlastelica, Vit Musil, Volodymyr Kuleshov, and Georg Mar- tius. Backpropagation through combinatorial algorithms: Identity with projection works. In The Eleventh International Conference on Learning Representations, 2022.
Abdelkader Sbihi and Richard W Eglese. Combinatorial optimization and green logistics. Annals of Opera- tions Research, 175(1):159 175, 2010.
Bo Tang and Elias Boutros Khalil. Pyepo: A pytorch-based end-to-end predict-then--ptimnie library with linear objective function. In OPT 2022: Optimization for Machine Learning (NeurIPS 2022 Workshop), 2022.
Vladimir Vapnik. The nature of statistical learning theory. Springer science & business media, 1999.
Marin Vlastelica, Anselm Paulus, Vit Musil, Georg Martius, and Michal Rolinek. Differentiation of blackbox combinatorial solvers. In International Conference on Learning Representations, 2019.
Qi Wang and Chunlei Tang. Deep reinforcement learning for transportation network combinatorial opti-
mization: A survey. Knowledge-Bout Based Systems, 233:107526, 2021.
Bryan Wilder, Bistra Dilkina, and Milind Tambe. Melding the data-decisions pipeline: fecision-focused lcarning for combinatorial optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pp. 1658-1665, 2019.
Alp Yurtsever, Varun Mangalick, and Suvrit Sra. Three operator splitting with a nonconvex loss function. In International Conference on Machine Learning, pp. 12267-12277 PMLR, 2021.
Liwei Zhong and Guochun Tang. Preface: Combinatorial optimization drives the future of health care. Journal of Combinatorial Optimization, 42(4):675-677 676, 2021.
Günter M Ziegler. Lectures on polytopes, volume 152. Springer Science & Business Media, 2012.
17

---

Published in Transactions on Machine Learning Research (07 (224))
A Appendix A: Proofs
For the reader's convenience restate cach result given in the main text before proving it.
Theorem Let C1,Co C2 be as in equation 14 and suppose fa d) = we(d) x + 211 for any neural network we d). For any a E (0, 2/7) define the sequence { by:
z+1 = =To(z*) for Kll
(28)
where
=- + Pc, (12 (e-
(29)
If th A (2) then the sequence { converges to De in equation 8 and || x =
First note that (a.ee 7, d) = we(d) + ya. Because fe(ro; v, d) is strongly convex, Te unique and is characterized by the first order optimality condition:
Proof.
VaJe(re;d) (x - xe) 20 for all x ECC.
(30)
Note that equation 30 can equally be viewed
variational inequality with operator Fe(æ; Vzfe(zo;d). We deduce that 2e with ze = (2e( is fixed point of
(= - (z) Pc, (2 (c, (=) (we + (c)()))) Pc) Pc. ((2 an) (2) - awwe(d)
(31a)
(31b)
from McKenzie et 20021
Thcorem 4.2), which itself a standard application of idcas from operator
splitting Davis & Yin
2017
Ryu & Yin 2022).
As v.fe fel (zei7 d) is y- Lipschitz continuous it is also V Vafe is c/o-coceeriive by the Baillon Haddad theorem Bauschke & Combettes 2009). McKenzic et al. (2021 Theorem 3.3) then implics
[Bailln & Haddad 1977
and x= + zee However, this theorem does not give a convergence rate. In fact, the convergence rate can be deduced from known results; because Te is averaged for a < 2/7 the rate |=-+1-==d 0(1 Theorem 1). Thus,
that Rx -
follows from Ryu & Yin
2022
|| k+1 == Pc, (2 Pc, (2*)) s 12 =0=00(1))
(32)
desired.
Next, we state two auxiliary lemmas relating the Jacobian matrices to projections onto linear subspaces. Lemma 8. C A fa: : Ax = b), for full-rank A E Rmxn, b E RR with m < n, and Hi Null(A) then
aPc. d=
Pei for all 2 E ER"
(33)
Proof. Let A T denote the reduced SVD of A, and note that A E mm with m < have U. mm Σ E m and V E Raxm Differentiating the formula for Pc, given in Lemma Cammalyyilda
aPc. 1 - AtA
(34)
Oz
where Al A VT- Note
A'A= (M=(VE-IUT) (UE- VVT.
(35)
which is the orthogonal projection onto Range(V) Range(A) (AT). It follows that I A is the orthogonal projection to Range(AT)e Null(A).
18

---

Published in Transactions on Machine Learning Research (07 (07/2024)
Lemma 9. Define the multi-valued function
if >
da) d maxx(0, a)
0
if 0
(36)
[0
if a a==
and, for 2 E RR,, define H2, _ Span(e; : 24 > 0).
Then
OPc. ReLU(=)) diag(c(=)),
(37)
d=
==
and adopling the convention for choosing an element of (z) stated in the main texl:
d'ca
diag(c(3) PHa.
(38)
dz
Proof. First, suppose 2 E Ik" satisfies žq + 0, for all i E [n], i.e. v is a smooth point of Pc.. Note
d[ReLU(Z.) if and 2; > and d[ReLU(Z)) if or 2; < O.
(39)
dz
dz
Thus, the Jacobian matrix is diagonal with a 1 in the (t, (t,t)-th position whenever 2, > and 0 otherwise, diag(c(Z)). Now suppose 2; = for one For all & E R" with 2; < 0, the Jacobian díe, ds === is well-defined and has O in the (i,t)-th position, while for 2 E R with > 0, the Jacobian d/ca is
L.e. dPe,
de
ds well-defined and has 1 in the (i, i)-th position. Taking the convex hull yields the interval [0, 1] in the (i, i)-h position, as claimed. The case where Σ, = 0 for multiple i is similar.
Consequently, the product of dPc and any vector v R" equals and only if Span(e; : 2, > 0).
dz -==3
This shows the linear operator idempotent with fixed point set H2, dee it is the projection operator PH2.
Theorem 4 Suppose A is full-rank and Hi A Null(A), H2. Span (er : 0). Then for all & E R
iTe
PaL Pug. + (1 (1-w7) Prt. PHo.
(40)
dz
a=a
Proof. Differentiating the expression for Te in equation
16 with respect to yields
8T
dPca
dPc
(2 (2- -7)
dPca dz
1
=I
+
(41a)
0z
1=8 dz
dz (==y(()
==
s=
= PHæ.: + PHt, ((2 P 1 I), for all å IR"
(41b)
where, for notational brevity, we set v1 - (2 - a7) Pc, (3) - - awe (d) in the first line and the second
line follows from Lemmas I and Repeatedly using the fact, for any subspace H. C RR" PH1 = I PR, the derivative 0Te/Oz can be further rewritten:
aT dz
=I PHo.a + (2 (2-0)) PHa,s PH:
(42a)
1==8
= -Put. + +2 - a7): PHt, (1 (1- -PHt.)) - PHt
(42b)
= -PHls + P + (1 - a7): PHs, -Pai PHtle - (1 -ay) PHi PHI - PHi
(42c)
= -(I- (1-Pa PHi)) P + (1-a7)- PHe -P
(42d)
Pe+t PHG. + (1 - an) PHi Pre.es for all E R."
(42e)
completing the proof.
7
19

---

Published in Transactions on Machine Learning Research (07/2024)
We use the following lemma to prove Theorem 6
Lemma 10. If the LICQ condition holds at Te, then Hi n H2 {0}.
Proof. We first rewrite HF and Zizzzz The subspace 22 is spanned by all non-positive coordinates of ze.
By equation 17 [zek max (0, [ze]t) and so i A(re) and only if [ze] S 0. It follows that
H2 Span(e:: [zo] =0) = Span (er : A(re))
(43)
where we enumerate A(ze) via A(re) = {ir....,ddy iet. On the other hand, Hi Range (AT) Span(o)....,.e.l am) where of denotes the i-th row of A.
Let c E H4 n Ho be given. Since w E HT, there are scalars 01,-..,, 0z such that == o1es, +... + Cells Similarly, since Hat there are scalars ,I,...,Bm Bm such that Bial + .... Bm Hence
O=v-v= (a1ei, ... + are,)) ($1a1 + Pmam).
(44)
By the LICQ condition, (eip,...,ein U {ai,..., amd is a linearly independent set of vectors; hence on = ag = 31 = Pm = and, thus, v = O. Since v was arbitrarily chosen in H n H2 the result follows.
Theorem If the LICQ condition holds ne, A is full-rank and a E (0, (0,2/7), then | 1.
Proof. By Lemma 10 the LICQ condition implies Ht n n2 = {0}. This implies that either (i) the first principal angle between these two subspaces is nonzero, and so the cosine of this angle is less than unity, i.e.
>> cos(T)
max HEH, l [|||==1
(u, v),
(45)
or (ii) (at least) one of H Hit is the trivial vector space {0). In either case, let w E R" be given. By Theorem in case (ii)
aTe 0z
Put PutPuut. w + (1-a7) PHù Pu = a1- PMs Pu
(46)
==20
implying that.
0Tg 0z
w
= (1-ar) Prr,Pu w s (1 - ary))//ll]
(47)
where the inequality follows as projection operators are firmly nonexpansive. In case (i), write w = w1 + W22 where wi E H2.2u and W22 E Hz.e0 Appealing to Theorem again
TTe
W
H2..0 + +(1111111 (1- PHù Pu w PHt W2 + (1 ary) PH, W1.
(48)
0z
==
Pythagoras' theorem may be applied to deduce, together with the fact Puc uz and Pr, wi are orthogonal,
aTe
||
Prg wz
+ +1 (1--7))2
. w1 122
(49)
0=
1220
Since Wy E Ho the angle condition
45
implics
P w2 1 (Pxx W2, PHr, W2) = (w2, PH.PPr w2) (w2, Pec w2) S cos(T) ||w,||3,
(50)
where the third equality holds since orthogonal linear projections are symmetric and idempotent. Because projections are non-cxpansive and PH is lincar,
PH, will= 2= [Pr wi PHù 01 || 0112 || 1
(51)
20

---

Published in Transactions on Machine Learning Research (07 (274))
Combining 19
50
and 51
reveals
iTe
cOs(T) /|w2|12 21 + (1 (87))11w11 [/////" 1|2
(52a)
0=
==29
S max (1 a (1|to) 11 + 1|
(52b)
mao((c (1 } |
(52c)
noting the final equality holds since wh and w2 are orthogonal. Because 52
holds for arbitrarily chosen
w E
R",
@Te
{
TTo
sup
0z w
||/||| = 1 =1 V max((cos(T)))) (1-a)) < 1,
(53)
0z
==
==20
where the final incquality holds by 45) and the fact a E (0, (/7) implics 1- ar E (-1,1), as desirod.
Corollary 7 Te continuuusly differentiable with respect to θ at zg, the assumptions in Theorem hold and (iTe/00) (OTo/00) (80) has condition number sufficiently small, then po(zo) a descent direction for l with respect θ.
Proof. From the proof of Theorem we see that Ta is contractive with constant r V (cos a and the main theorem of (Fung al. 2022 guaranteeing pe is descent direction, long the condition number of (0Te/00) (0Te/d0) is less than 1/T.
Remark 11. Similar guarantees, albeit with less restrictive assumptions TTe/0O, can be deduced from the results of the recent work [Bolle et al. 2024).
B Derivation for Canonical Form of Knapsack Prediction Problem
For completeness, we explain how to transform the k-knapsack prediction problem into the canonical form equation g. and how to derive the standardized representation of the constraint polytope C. Rocall that the k-knapsack prediction problem, as originally stated, is
x²
arg max w T x where X {ar E 1), Sa S ef and S == [s] xEX
st] E Rexe
(54)
We introduce slack variables "1,.... Mr so that the inequality constraint Sa S becomes
-Sx +c e 2 - Sxx + c = and 2 0
x
- |s I
We relax the binary constraint 2; {0,1) to 0 I x I 1. account for the upper bound:
We add additional slack variables ,,,..., to
x, 2 and 2; 0 = [lexe Ozxk Lexe]
= 1
(55)
Putting this together, define
S
-Ixk
Oxxt ((++0)x(((+))) and b Ri+e
[-c]
A
(56)
Texe
Oxxk
Texe
1
Finally, redefine = [ an argmin) and obtain:
and w = [-w Ox 0g] (where we're using ww to switch the argmax to
arg min w. rxConv(X)
TT (d)x + 2||0|13 where Conv(X) {o Az b and 2 0}
(57)
21

---

Published in Transactions on Machine Learning Research (07 (07/2024)
C Experimental Details
C.1 Additional Data Details for PyEPO problems
We use PyEPO (Tang & Khalil
2022 specifically the functions shotesshott genData and
data. kwwpsack.packeaac genData .ennaaaaa to generate the training data. Both functions sample E ignxs with Bij = +1
with probability 0.5 and -1 with probability 0.5 and then construct the cost vector w(d) as
(Ba),, +3)
deæ +1 Eij
[w[d)]; 3.5deg
where deg = 4 and eij is sampled uniformly from the interval [0.5, 1.5]. Note that w(d) is, by construction, non-negative.
C.2 Additional Model Details for PyEPO problems
For PertOpt-net and BBOpt-net we use the default hyperparameter values provided by PyEPO, namely X = 5 for BB0pt-net, number of samples equal 3, e = 1 and Gunbel noise for PertOpt-net, See Pogančic el al. (2019) and Berthet et al. (020)) respectively for descriptions of these hyperparameters. We do so as these are selected via hyperparameter tuning Tang & Khalil (2022). As DYS-net and CVX-net solve a regularized form of the underlying LP equation 7 the regularization parameter y needs be set. After experimenting with different values, we select 7 5 x 10-1 for DYS-net and ? x 10-1 for CVX-net.
C.3 Additional Training Details for Knapsack Problem
For all models we use an initial learning rate of 10-3 and a scheduler that reduces the learning rale whenever the validation loss plateaus. We also used weight decay with a parameter of 5 10-4
C.4 Additional Model Details for Large-Scars Scale Shortest Path Problem
@ of the original TensorFlow code associated to the paper Berthet et al. 2020). We experimented with various hyperparameter settings for grids and found setting the number of samples equal 3, the temperature (i.e. E) to 1 and using Gumbel noise lo work best, so we used these values for all other shortest path experiments. Our implementation of BB-net uses the blackbox-backpron package associated to the paper Pogančić et al. (2019). We found setting À 100 to work best for 5-by-5 grids,
Our implementation of PertOpt-net used a PyTTrcch implementation
use this value for all other shortest path experiments. Again, we select 5 x 10-4 for DYS-net and 7 = 5 x 10 for CVX-net
C.5 Additional Training Details for Large-Scale Shortest Path Problem
We use the MSE loss to train DYS-net, PertOpt-net, and CVX-net. We tried using the MSE loss with BB-net but this did not work well, we used the Hamming (also known as 0 1) loss, done in Poganči et al. 2019).
To train DYS-net and CVX-net, we use an initial learning rate of 10-2 and use a scheduler that reduces the learning rate whenever the test loss plateaus we found this perform the best for these two models. For PertOpt-net we found that using a fixed learning rate of 10-2 performed the best. For BB-net, we performed logarithmic grid-search on the learning rate between 10-1 to 10-4 and found that 10-3 performed best; also attempted adaptive learning rate schemes such reducing learning rates on platean but did not obtain improved performance.
6Scc code at github. giterrocotttrreerrrrrrooooosssee
7See code htw.://wwwccc..eeee.eeeeeeeeeeee
88ee code at http::///ghhhtllccoo/aaaaaaaaaaa
22

---

Published in Transactions on Machine Learning Research (07 (07/2024)
10- 10-by-10
20-by-20
30-by-30
Test MSE Loss
10a
DYS
DVS
10a
DYS
10
CVX
CVX
CVX
10a
PertOpt
PertOpt
10-
PertOp
BB
BB
10a
10a
10
10
10
10
0
25
50 epochs
75
100
0
25
50
75
100
25
50 epochs
75
100
epochs
Training Time (in minutes)
80
DYS
DYS
DYS
60
CVX
300
CVX
800
CVX
PertOpt
PerOOpt
600
PertOpt
40
200
BB
BB
400
20
100
0
25
50
75
100
25
50
100
0
50
100
epochs
epochs
epochs
Figure Comparison of of DYS-Net, cyxpylayers (Agrawal 2020), and Blackbox Backpropaggaton-nnt (BB-Net) Pogančić al.
PertOptNet (Berthet et al, 2019 for three different grid sizes: 10 x 10 (first column), 20 20 (second columm), and 30 30 (third column). The first row shows the MSE loss vs. epochs of the lesting dataset. The second row shows the training time epochs.
2019ah.
C.6 Additional Model Details for Warcraft Experiment
For PertOpt-net and Bet we again use the default hyperparametor values provided by PyEPO, namoly = for B80pt-net, number of samples equal 3, = 1 and Gumbel noise for PertOpt-net, We set = 0.5 for both CVX-net and Dys-net.
C.7 Hardware
All networks were trained using a AMD Threadripper Pro 3955WX: 16 cores, 3.90 GHz, 64 MB cache, PCLe 4.0 CPU and an NVIDIA RTX A6000 GPU.
Additional Experimental Results
In Figure we show the test loss and training time per epoch for three architectures: DYS-net, CVX-net, and PertOpt-net for 10-by-10, 20-by-20, and 30-by-30 grids. In terms of MSE loss, CVX-net and DYS-net lead to comparable performance. In the second row of Figure 0 we observe the benefits of combining the threc-operator splitting with JFB (Fung al.) 2022); in particular, DYS-net trains much faster. Figure shows some randomly selected outputs for the three architectures once fully trained.
23

---

Published in Transactions on Machine Learning Research (07 (07/202) /024)
True Path
DYS-net
CVX-net
PertOpt-net
10-by-10
20-by-20
30-by-30
Figure 7: True paths (column 1), paths predicted by DYS-net (column 2), CVX-net (column 3), and PertOpt-net (column 4). Samples are taken from different grid sizes: 10-by-10 (row 1), 20-by-20 (row 2), and 30-by-30 (row 3).
24


Reviews:
[{"review_id": "BCHq8UFPUZ", "reviewer": "Reviewer_UhYG", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is well-written with good organization, related work and technical descriptions that make the paper easy to read. Overall a high-quality presentation. \n- The paper convinces that its approach is sound and well-justified. It  makes use of several very recent methodologies as part of its overall approach. I enjoyed reading it and did not have any major objections on a technical level. \n- The technical development and design choices make sense and are sufficiently explained, and the theory is thoroughly described, though I may have overlooked some details. The proofs, based on the gradients' validity as a descent direction, are a nice inclusion.\n- The experimental evaluations are clear, well-illustrated and decently thorough, on the test problems that were chosen. \n\n\nWeaknesses:\n- The title is misleading, implying that the paper will be about “Learning to Solve ILP problems”. Ultimately, the proposed method is evaluated on “Predict-Then-Optimize (PtO)” tasks against methods in that domain. The paper makes no contribution in learning to solve ILP problems, but rather in learning the uncertain coefficients which may appear in a ILP problem. In this paper's framework, the ILP problem must be solved (at test time) by conventional methods, and not a learned model, unless I am mistaken. \n\n- Further, the Related Work is framed in a misleading way, which contributes to this paper's confusing framing of its own contributions. First, I think it's generally accepted that \"Learning to Optimize (L2O)\" is about speeding up the solution of optimization problems. This paper expands the notion of L2O to include learning the uncertain coefficients of an optimization problem, but this setting already has a name (Predict-Then-Opimize, or Decision-Focused Learning). Moreover, some works cited in the subsection \"Differentiable Combinatorial Optimizers\" are miscategorized. For example, Berthet 2020 is explicitly only developed for linear programs with continuous variables, and not ILPs. Elmachtoub 2022 does not give a differentiable optimizer, but rather only a surrogate subgradient for its resulting regret loss. \n\n- Also in terms of its framing of technical contributions: not only does the paper not show how to learn to solve ILPs (as implied by the title, see the first bullet), it also does not really even show how to learn the coefficients of ILPs. It only shows how to learn the coefficients of convex continuous Quadratic Programming problems, which is a much lesser feat. It's true that a couple other works have claimed to learn ILP coefficients in Predict-Then-Optimize, when what they really do is form a QP approximation to the ILP and learn its coefficients instead. So there is precedent for what the authors do here. However, I think it's an unjustified practice, since it assumes that QP relaxations of ILPs are valid in general, and this has never been sufficiently demonstrated in the literature (unless I'm mistaken). This may seem like a nitpick, but I think it's important, because we are now seeing it become \"folklore\" that ILPs can simply be replaced with QPs in the general case. Unless this can be justified, it seems like a misleading trend in the literature motivated by \"up-selling\" of the papers' contributions, where any QP method can be sold as an ILP method. I'm open to push-back from other reviewers on this view.\n\n- Because the paper claims to handle ILP problems, the experimental evaluation is arguably weak. We are only evaluating the method on two problem types, knapsack and shortest path. Of those, shortest path is not properly a ILP problem, it is a continuous LP with totally unimodular constraints which guarantee integral solutions. But it is solved as an LP, or even n log n methods like Dijstra's. So the shortest path experiment does not demonstrate the premise of the paper, which promises to learn to solve ILPs. The knapsack problem has been shown in empirically in other works to admit a good QP relaxation for Predict-Then-Optimize, so it's not surprising to see as the other experiment. For the paper to justify its claims, it needs to evaluate on another NP-hard, proper ILP problem (in my opinion).\n\n- About the choice of forward pass. The QP approximation is solved by Davis-Yin splitting. I think the paper could do a better job of motivating this. Given that the paper is about solving large-scale problems, this cannot be the fastest way to solve the QP in general. In fact, we seem to have the option to use any black-box optimization solver to replace this forward pass, and proceed directly to the backward pass in Step 4. Is this correct? If so, it should be made clear. This is important to question, since the forward pass is part of an overall method which claims efficiency advantages over other methods. The results of Figure 2 do not indicate whether the efficiency gains stem from the forward pass or the backward pass.\n\n- This DYS-net method requires pre-computation of optimal solutions (x*) to each instance of the NP-Hard MILP problems that make up the training set. They are needed for its loss function. Importantly, none of the evaluated baseline methods require this step. If this is correct, I think it presents a major challenge to the claimed efficiency gains over other methods, which use the ILP's empirical objective value as their loss function.\n\n- Beyond the efficiency drawback of requiring precomputed optimal solutions, I expect that this requirement will also lead to issues in training. Since the MILP problems are NP-Hard, we can not expect to find their true optimal solutions in general. If we settle for sub-optimal precomputed solutions, you cannot expect the training set to preserve the true relationship between unknown problem parameters and the optimal solutions. For this reason and the one above, the loss function seems flawed.\n\n- Building on the above three points, I would conjecture that the reason for using the regression loss (rather than minimizing the objective as a loss) is because the Davis-Yin forward pass will tend not to converge in finite iterations in the latter case. However, replacing the Davis-Yin method with a black-box forward pass optimization would seem to resolve this issue, guaranteeing feasibility so that we can use f as a loss function and avoid the issue of precomputed optimal solutions.  Is this correct, and if so, why not take this type of approach?\n\n-It's not clear what is the significance of the theoretical results. Are they special cases of previously established results, or is this theory novel? This is a question, not necessarily a weakness.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Summary of the Paper: \nThe paper proposes a method for ‘learning to solve’ mixed-integer linear programs, with uncertain coefficients. It is based on a  convex relaxation of the MILP problem’s feasible set, combined with quadratic regularization to allow differentiation through the approximated problem. The paper claims competitive accuracy with other methods on its test problems, with a decided advantage in training efficiency. \n\nSummary Review:\nThe paper, to the best of my evaluation, is technically solid and sound. Its claims seem well-justified, though I did not check the proofs. The proposed method uses an interesting mix of very modern techniques to try and minimize computational effort in solving a difficult problem. On the other hand, it suffers from a lack of clarity in its key insights and contributions, so that I am ultimately unconvinced of its overall significance, despite being interesting. The “Strengths” and “Weaknesses” below elaborate on these points.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Many of my criticisms (the first four Weaknesses bullets) would be addressed if the authors were willing to reframe their paper as one which learns the coefficients of QP problems rather than ILPs. This would require a change of title and significant reframing of the intro sections. Depending on the opinion of other reviewers, I may see this as a major issue with the paper, because it means the paper is not supporting its own claims. I also welcome the authors to give their contrasting view on the issue.\n\nThe authors should also better motivate their use of Davis-Yin splitting in the forward pass (bullet 5 in Weaknesses), and clarify where their computational advantage actually stems from (forward or backward pass).", "broader_impact_concerns": "I don't foresee any significant ethical concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "ZnC3Tamo4p", "reviewer": "Reviewer_4Zt6", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "- Strengths:  \n\n  \n\n1. The paper presents its claims with proof.  \n\n  \n\n2. The paper proposes a novel descent direction for differentiating combinatorial optimization problems.   \n\n  \n\n- Weakness:  \n\n  \n\nThe experimental evaluation is relatively weak. Especially, in the warcraft dataset, it performs significantly worse than BB-net.   \n\n  \n\n- Questions: \n\n1. In Algorithm 1 Line 3, should it not be $P_{C_2} (Z^K)\\approx x_{\\Theta}$? Because in Equation 17 , $x_{\\Theta} = P_{C_2} (Z_\\Theta)$. \n\nSimilarly, it's not clear to me how  $P_{C_1}$ appear in Eq. (19), to me it seems it would be   $P_{C_2}$. \n\n2. Can you comment on how to set the value of hyperparameter $\\alpha$ and $K$. Obviously they could be selected by hyperparameter tuning; but my question is what would be a suitable range to look for?  \n\n3. Both Figure 2 (a) (b) and Figure 4 are for the shortest path problem. However, there are some discrepancies while comparing them. For example, the regret vs grid size plot for BB in Fig 2 (a) does not match with Fig 4 (c) till Grid size 30. In Figure 2 (a)regret is increasing with grid size; which is not the case in Figure 4 (c).  \n\n4. Training time of PertOpt is sometime higher than BB, sometime lower. Why would this happen? How many Monte Carlo samples have been used for PertOpt? \n\n5. Do you have any intuition why BB performs better than DYS  by such a significant margin in the Warcraft problem? I am not sure what do you mean by \"discrete\" nature of cost vector?\n\n6. Finally, in my opinion, DYS should be compared with CVX, as both of them are implementation of differentiable optimization layers, which do not require another solver (Gurobi for example). BB and PertOpt, on the other hand, improve finite difference approximation and require an optimization solver. Can you replace the Warcraft experiment with another problem, where you can run CVX?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The problem considered in this paper is contextual optimization, where a parameter of an optimization problem is unknown but depends on the context. The goal is to learn a neural network to map from the context to the parameter, minimizing the final discrepancy of the solution to the optimization problem. As this requires embedding the optimization problem into the neural network training, an implementation of a differentiable optimization layer is necessary. Moreover, it is assumed that in the training data, the true parameter vector is not recorded; only the solution and the context are present in the training data. While this problem has been covered in many recent works, as noted in Section3, the focus of this paper is on efficiency.   \n\nImplementing a projection gradient descent for this problem is difficult, as the projection must not only ensure $Ax=b$ but also $x \\geq 0$. Put simply, finding even a *feasible* solution for a combinatorial optimization problem is inherently challenging. The projections $P_{C_1}$ and $P_{C_2}$ are associated with satisfying the constraints $Ax=b$ and $x \\geq 0$, respectively. While combining these projections, $P_{C_1} \\cap P_{C_2}$, presents difficulties, implementing each projection individually, either $P_{C_1}$ or $P_{C_2}$, is relatively straightforward. Theorem 2 addresses this by posing the projection as a fixed-point iteration, where $P_{C_1}$ and $P_{C_2}$ appear separately. The proposed technique is supported by experimental evaluation.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- How we can come from Equation 19 from Equation 18? [critical]\n- Elaborate the explanation of the result (5.1.3 contains only one sentence) [strengthen]\n- Consistency between Figure 2 and Figure 4 [crucial]", "broader_impact_concerns": "I do not see any immediate requirement of mentioning broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "W7lADKGa9s", "reviewer": "Reviewer_DjLc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n**Would some individuals in TMLR's audience be interested in the findings of this paper? YES**\n\nThe paper is really well-written. It is easy to follow despite the intrication of several different ideas, and explains in detail the motivation for `DYS-Net`. I have enjoyed it very much, and will use it in my research once it is published.\n\nThe proposed method `DYS-Net` is straightforward to implement and fully GPU-compatible. It is designed to scale to large instances thanks to its use of projected gradient, instead of costly primal-dual methods that differentiate through KKT conditions.\n\nThe theoretical analysis is sound, and shows that the counter-intuitive Jacobian-free approach still works in this non-smooth, non-contractive case. \n\n## Weaknesses\n\n**Are the claims made in the submission supported by accurate, convincing and clear evidence? NOT YET**\n\nMy main criticism is that `DYS-Net` claims to allow high-dimensional settings... but still requires computing a full SVD for the constraint matrix $A$. The SVD is a very expensive (cubic) operation, so even if you just perform it once at the beginning, I'm unsure how much you gain on truly large instances.\nFurthermore, the constraint matrix of real-life problems has a lot of structure (sparsity, blocks, bands) which the SVD does not exploit. Typically, in graph optimization, the constraint matrix is some variant of the adjacency matrix. And unfortunately, I think (but I'm not 100% sure) that the SVD of a sparse matrix is not sparse in general.\n\nLuckily, I think there is a solution. Indeed, the SVD $A = U \\Sigma V^\\top$ is only used to derive the pseudoinverse $A^\\dagger = V \\Sigma^{-1} U^\\top$, and the pseudoinverse is only used for solving linear systems in $P_{\\mathcal{C}_1}(z) = z - A^\\dagger (Az - b)$.\nGiven this observation, why not solve the linear system directly with the matrix $A$, using an iterative Krylov solver like GMRES which only requires matrix-vector products?\nThe same remark applies to the solution of a linear system involving the Jacobian $\\mathcal{J}$: no need to compute a (pseudo)inverse if a Krylov solver is available. \n\nOf course, the previous criticism is only valid if the SVD is actually used in the `DYS-Net` algorithm. If I understand correctly, the Jacobian-free approach dispenses us from ever computing $P_{\\mathcal{C}_1}$, which means we never need the SVD in the code? I am not completely certain because the notations $\\mathcal{C}_1$ and $\\mathcal{C}_2$ are sometimes switched (probably typos).\nIn any case, if the SVD is only an ingredient of the theoretical proofs, then it must be clarified, and there is no need to insist so much on the performance of the split projection operator.\n\nMy other criticism concerns the experimental part. First, I am not sure that the encoder $w(d)$ is guaranteed to yield positive weights, which would be a problem both for shortest paths (with Dijkstra's algorithm in particular) and for knapsack.\nSecond, I do not understand how `DYS-Net` can be more precise than `CVX`. Faster, for sure, but as far as precision is concerned, `DYS-Net` can only be worse than the true convex program, right?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a method called `DYS-Net` for automatic differentiation of Integer Linear Programs (ILPs) with respect to a cost vector. The goal is to facilitate the use of ILPs in deep learning pipelines, especially when the dimension of the problem gets large.\n\nTo approximate the Jacobian of the ILP\n$$\nx(w) = \\mathrm{argmin} \\\\{w^\\top x : x \\in \\mathcal{C} \\cap \\mathbb{Z}^n \\\\}\n\\quad \\text{where} \\quad\n\\mathcal{C} = \\\\{x \\in \\mathbb{R}^n : Ax = b, x \\geq 0 \\\\},\n$$\n\n`DYS-Net` combines known ideas from convex optimization and decision-focused learning in a novel way:\n\n1. Relax the ILP and add a quadratic regularization to the objective, yielding a Linearly Constrained Quadratic Program (LCQP) [Wilder et al., 2019]\n$$\nx(w) = \\mathrm{argmin} \\\\{w^\\top x + \\frac{\\gamma}{2} \\lVert x \\rVert^2 : x \\in \\mathcal{C} \\\\}\n$$\n2. Solve the LCQP with a projected gradient algorithm, using Davis-Yin splitting to efficiently compute the projection onto the polynomial $\\mathcal{C}$ [McKenzie et al., 2023]. The key ingredients are the decomposition\n$$\n\\mathcal{C} = \\mathcal{C}_1 \\cup \\mathcal{C}_2\n\\quad \\text{where} \\quad\n\\mathcal{C}_1 = \\\\{ x \\in \\mathbb{R}^n : Ax = b \\\\}\n\\quad \\text{and} \\quad\n\\mathcal{C}_2 = \\\\{ x \\in \\mathbb{R}^n : x \\geq 0 \\\\}\n$$\nand the use of a precomputed SVD for the fixed constraint matrix $A$.\n3. Avoid backpropagating through the iterative procedure, and instead leverage implicit differentiation on the fixed point condition of the projected gradient [Blondel et al., 2022]\n4. Approximate the implicit function theorem by replacing the Jacobian of the fixed point condition with the identity [Fung et al, 2022]\n\nA theorem justifies that the Jacobian-free heuristic finds a descent direction, which is non-trivial to obtain from previous results in the literature.\nNumerical experiments showcase the efficiency and precision of `DYS-Net` compared with competing approaches on two standard benchmarks: the shortest path problem and the knapsack problem.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "## Important requests (necessary for acceptance)\n\n- Clarify the performance impact of the SVD step, and whether the SVD is used in practice or only for proof purposes\n- Discuss the possible use of iterative linear solvers in lieu of (pseudo)inversion\n- Comment on what happens in realistic scenarios where the constraint matrix $A$ is hard to write down explicitly, but instead given implicitly through a domain-specific language (like [cvxpy](https://www.cvxpy.org/) or [JuMP](https://jump.dev/))\n- Fix the crucial notation switches between $\\mathcal{C}_1$ and $\\mathcal{C}_2$ that occurred e.g. in Equation (19) and the one before.\n- The paper mentions a \"provided code\" but I couldn't find it, did I just miss the line? Or was it removed due to anonymity concerns? If it is the latter, tools such as https://anonymous.4open.science/ allow for anonymous copies of existing GitHub repos\n- Make sure that the cost vector $w(d)$ for both categories of experiments remains nonnegative. For instance in the shortest paths, $w(d) = Wd$ with $d_k \\sim \\mathcal{U}[0, 1]$ but the sign of $W$ is not specified, and the issue is similar for the knapsack (appendix C.1)\n\n## Additional suggestions and questions\n\n- Shorten the introductive part of the abstract, and put more focus on the actual contributions of the paper\n- Could your framework be used when there are no \"ground truth\" optimal solutions $x^\\star$, that is, when learning by experience instead of imitation?\n- In your Theorem 2, the result $\\lVert x^{k+1} - x^k \\rVert_2^2 = O(1/k)$ does not seem to imply that the sequence converges. Is that truly what you wanted to write?\n- In Corollary 7, how do you ensure that the condition number of $(\\partial_{\\Theta} T)^\\top \\partial_{\\Theta} T$ is small enough? Is this something the user can influence? How do we check it?\n- Ideally, keep the notations for $\\mathcal{C}_1$ and $\\mathcal{C}_2$ coherent with [McKenzie et al., 2023]: when I compare your equation (16) with their equation (4.5), they are reversed. This might be the cause of the switches mentioned above\n- Fix the references to the theorems in [McKenzie et al., 2023]: I think what you call 3.2 and 3.3 are actually 4.2 and 4.3, and one of them is a Lemma in the published version\n- In the proof of Theorem 6, can you give a reference for the result on the first principal angle between subspaces with zero intersection? It was not obvious to me intuitively, although I don't doubt it\n- You claim that you cannot use PyEPO for grid sizes larger than 30 because it automatically calls Gurobi. However, you can actually plug any algorithm into PyEPO with a [user-defined model](https://khalil-research.github.io/PyEPO/build/html/content/examples/model.html#user-defined-models-from-scratch), so a Dijkstra implementation is also possible instead of Gurobi.\n- I'm very surprised that CVX is unable to load when the grid size exceeds 30. Is the constraint matrix $A$ encoded in a sparse manner?\n- Even if the cost vector given by $w(d)$ is nonnegative to begin with, the added random vector of `PertOpt` or the black box interpolation of `BB` might change $w$ so much that one of its components becomes negative. This is the reason behind the multiplicative perturbation introduced in [[Dalle et al., 2022](https://arxiv.org/abs/2207.13513)], which preserves the sign. Of course such a behavior is rarely observed if the perturbation is small enough, and the original code from Berthet et al implicitly handles it with (what I think is) an incorrect Dijkstra implementation. I'm not asking you to change your experiments because of that, but I would be curious to see what you think, and if it is worthy of discussion. \n\n## Typos and style\n\n- Page 2: \"L20\" (L twenty) $\\to$ \"L2O\"\n- Page 6: \"Lipshitz\" $\\to$ \"Lipschitz\"\n- Page 10, Table 1: \"network size\" $\\to$ \"neural network size\" (not to be confused with the graph) \n- Page 16, proof of Theorem 2: $\\nabla_x f(x_{\\Theta}; \\gamma, d)$ $\\to$ $\\nabla_x f(x; \\gamma, d)$", "broader_impact_concerns": "There are no ethical implications that would require a broader impact statement.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Reviewers were unanimous in recommending acceptance of the paper, and generally feel excited about the new proposed approach for differentiating through LPs (with quadratic regularization).

There was also an internal discussion on whether we should recommend this submission for any additional certifications and/or ICLR presentations. While we are excited about the paper, the consensus is that additional recognition would require further revisions on the paper to strengthen some points. We detail our concrete comments below, and the authors can decide whether they want to make the requested revisions. If the requested changes are made satisfactorily, I will ask the Editors in Chief to add relevant certifications/recommendations, who have the ultimate authority to grant these recognitions.

Minor revisions (required for the camera ready version for acceptance to TMLR):

- Clarify with a clear acknowledgement that the focus on (features,solution) training data instead of (features,true cost parameters) is due to not being able to outperform SPO+. Please also clarify what "outperforming" means, in terms of training time and/or predictive accuracy for regret.

- Add discussion on forward vs backward pass efficiency gains. My understanding is that by using fast Jacobian-Free backpropagation for the backward pass, the forward pass no longer needs primal-dual solving (via not differentiating the KKT conditions at LP optimality), which further allows for the fast Davis-Yin solving in the forward pass. This is an interesting point, and the discussion would serve to strengthen the paper.

- Update the experimental result tables according to the rebuttal discussion.

More major revisions (required for my recommendation as AE for additional certification and/or ICLR journal-to-conference track):

- Reframe the introduction to emphasize that the main contribution is a novel generic differentiable LP layer approach, before saying that the focus on (features,solution) is due to comparison with SPO+ etc. (see above), and before mentioning that the approach is then also applied to LP relaxations of MILPs.

- Add at least some smaller experiments demonstrating the approach on (feature,true cost parameter) data, comparing with SPO+ and perhaps also Mandi and Guns (Interior Point Solving for LP-based prediction+optimisation, NeurIPS 2020), demonstrating that the proposed approach might not be fully competitive with other methods in the context of (feature,true cost) training data. Reviewers and I believe that the transparency will make the work more convincing.

- For the runtimes reported in experiments, separate out the times for the forward pass and backward pass for both cvx and DYS-net, in the context of the discussion outlined above.

- At least one more benchmark on a (mixed) *integer* linear program problem, ideally on a problem that has some decent integrality gap, to further demonstrate the viability of integrality relaxation+quadratic regularization for MILP problems. I think some non-decision-focused "Baseline" method is still needed here to make this point.
Decision: Accept (poster)