Title: DIFFERENTIATION OF BLACKBOX COMBINATORIAL SOLVERS

Abstract: Achieving fusion of deep learning with combinatorial algorithms promises transformative changes to artificial intelligence. One possible approach is to introduce combinatorial building blocks into neural networks. Such end-to-end architectures have the potential to tackle combinatorial problems on raw input data such as ensuring global consistency in multi-object tracking or route planning on maps in robotics. In this work, we present a method that implements an efficient backward pass through blackbox implementations of combinatorial solvers with linear objective functions. We provide both theoretical and experimental backing. In particular, we incorporate the Gurobi MIP solver, Blossom V algorithm, and Dijkstra's algorithm into architectures that extract suitable features from raw inputs for the traveling salesman problem, the min-cost perfect matching problem and the shortest path problem. The code is available at https://github.com/martius-lab/blackbox-backprop. * These authors contributed equally.
Accordingly, we will not rely on standard techniques for gradient estimation (see (Mohamed et al., 2019) for a comprehensive survey).

1. INTRODUCTION:

The toolbox of popular methods in computer science currently sees a split into two major components. On the one hand, there are classical algorithmic techniques from discrete optimization -graph algorithms, SAT-solvers, integer programming solvers -often with heavily optimized implementations and theoretical guarantees on runtime and performance. On the other hand, there is the realm of deep learning allowing data-driven feature extraction as well as the flexible design of end-to-end architectures. The fusion of deep learning with combinatorial optimization is desirable both for foundational reasons -extending the reach of deep learning to data with large combinatorial complexity -and in practical applications. These often occur for example in computer vision problems that require solving a combinatorial sub-task on top of features extracted from raw input such as establishing global consistency in multi-object tracking from a sequence of frames.

The fundamental problem with constructing hybrid architectures is differentiability of the combinatorial components. State-of-the-art approaches pursue the following paradigm: introduce suitable approximations or modifications of the objective function or of a baseline algorithm that eventually yield a differentiable computation. The resulting algorithms are often sub-optimal in terms of runtime, performance and optimality guarantees when compared to their unmodified counterparts. While the sources of sub-optimality vary from example to example, there is a common theme: any differentiable algorithm in particular outputs continuous values and as such it solves a relaxation of the original problem. It is well-known in combinatorial optimization theory that even strong and practical convex relaxations induce lower bounds on the approximation ratio for large classes of problems (Raghavendra, 2008; Thapper & Živný, 2017) which makes them inherently sub-optimal. This inability to incorporate the best implementations of the best algorithms is unsatisfactory.

In this paper, we propose a method that, at the cost of one hyperparameter, implements a backward pass for a blackbox implementation of a combinatorial algorithm or a solver that optimizes a linear objective function. This effectively turns the algorithm or solver into a composable building block of neural network architectures, as illustrated in Fig. 1 . Suitable problems with linear objective include classical problems such as SHORTEST-PATH, TRAVELING-SALESMAN (TSP), MIN-COST-PERFECT-MATCHING, various cut problems as well as entire frameworks such as integer programs (IP), Markov random fields (MRF) and conditional random fields (CRF).

The main technical challenge boils down to providing an informative gradient of a piecewise constant function. To that end, we are able to heavily leverage the minimization structure of the underlying combinatorial problem and efficiently compute a gradient of a continuous interpolation.

While the roots of the method lie in loss-augmented inference, the employed mathematical technique for continuous interpolation is novel. The computational cost of the introduced backward pass matches the cost of the forward pass. In particular, it also amounts to one call to the solver.

In experiments, we train architectures that contain unmodified implementations of the following efficient combinatorial algorithms: general-purpose mixed-integer programming solver Gurobi (Gurobi Optimization, 2019), state-of-the-art C implementation of MIN-COST-PERFECT-MATCHING algorithm -Blossom V (Kolmogorov, 2009 ) and Dijkstra's algorithm (Dijkstra, 1959) for SHORTEST-PATH. We demonstrate that the resulting architectures train without sophisticated tweaks and are able to solve tasks that are beyond the capabilities of conventional neural networks.

2. RELATED WORK:

Multiple lines of work lie at the intersection of combinatorial algorithms and deep learning. We primarily distinguish them by their motivation.

Motivated by applied problems. Even though computer vision has seen a substantial shift from combinatorial methods to deep learning, some problems still have a strong combinatorial aspect and require hybrid approaches. Examples include multi-object tracking (Schulter et al., 2017) , semantic segmentation (Chen et al., 2018) , multi-person pose estimation (Pishchulin et al., 2016; Song et al., 2018) , stereo matching (Knöbelreiter et al., 2017) and person re-identification (Ye et al., 2017) . The combinatorial algorithms in question are typically Markov random fields (MRF) (Chen et al., 2015) , conditional random fields (CRF) (Marin et al., 2019 ), graph matching (Ye et al., 2017) or integer programming (Schulter et al., 2017) . In recent years, a plethora of hybrid end-to-end architectures have been proposed. The techniques used for constructing the backward pass range from employing various relaxations and approximations of the combinatorial problem (Chen et al., 2015; Zheng et al., 2015) over differentiating a fixed number of iterations of an iterative solver (Paschalidou et al., 2018; Tompson et al., 2014; Liu et al., 2015) all the way to relying on the structured SVM framework (Tsochantaridis et al., 2005; Chen et al., 2015) .

Motivated by "bridging the gap". Building links between combinatorics and deep learning can also be viewed as a foundational problem; for example, (Battaglia et al., 2018) advocate that "combinatorial generalization must be a top priority for AI". One such line of work focuses on designing architectures with algorithmic structural prior -for example by mimicking the layout of a Turing machine (Sukhbaatar et al., 2015; Vinyals et al., 2015; Graves et al., 2014; 2016) or by promoting behaviour that resembles message-passing algorithms as it is the case in Graph Neural Networks and related architectures (Scarselli et al., 2009; Li et al., 2016; Battaglia et al., 2018) . Another approach is to provide neural network building blocks that are specialized to solve some types of combinatorial problems such as satisfiability (SAT) instances (Wang et al., 2019) , mixed integer programs (Ferber et al., 2019) , sparse inference (Niculae et al., 2018) , or submodular maximization (Tschiatschek et al., 2018) . A related mindset of learning inputs to an optimization problem gave rise to the "predict-and-optimize" framework and its variants (Elmachtoub & Grigas, 2017; Demirovic et al., 2019; Mandi et al., 2019) . Some works have directly addressed the question of learning combinatorial optimization algorithms such as the TRAVELING-SALESMAN-PROBLEM in (Bello et al., 2017) or its vehicle routing variants (Nazari et al., 2018) . A recent approach also learns combinatorial algorithms via a clustering proxy (Wilder et al., 2019) .

There are also efforts to bridge the gap in the opposite direction; to use deep learning methods to improve state-of-the-art combinatorial solvers, typically by learning (otherwise hand-crafted) heuristics. Some works have again targeted the TRAVELING-SALESMAN-PROBLEM (Kool et al., 2019; Deudon et al., 2018; Bello et al., 2017) as well as other NP-Hard problems (Li et al., 2018) . Also, more general solvers received some attention; this includes SAT-solvers (Selsam & Bjørner, 2019; Selsam et al., 2019) , integer programming solvers (often with learning branch-and-bound rules) (Khalil et al., 2016; Balcan et al., 2018; Gasse et al., 2019) and SMT-solvers (satisfiability modulo theories) (Balunovic et al., 2018) .

3. METHOD:

Let us first formalize the notion of a combinatorial solver. We expect the solver to receive continuous input w ∈ W ⊆ R N (e.g. edge weights of a fixed graph) and return discrete output y from some finite set Y (e.g. all traveling salesman tours on a fixed graph) that minimizes some cost c(w, y) (e.g. length of the tour). More precisely, the solver maps w → y(w) such that y(w) = arg min y∈Y c(w, y).

(1)

We will restrict ourselves to objective functions c(w, y) that are linear , namely c(w, y) may be represented as c(w, y) = w • φ(y) for w ∈ W and y ∈ Y (2)

in which φ : Y → R N is an injective representation of y ∈ Y in R N . For brevity, we omit the mapping φ and instead treat elements of Y as discrete points in R N .

Note that such definition of a solver is still very general as there are no assumptions on the set of constraints or on the structure of the output space Y .

Example 1 (Encoding shortest-path problem). If G = (V, E) is a given graph with vertices s, t ∈ V , the combinatorial solver for the (s, t)-SHORTEST-PATH would take edge weights w ∈ W = R |E| as input and produce the shortest path y(w) represented as φ(y) ⊆ {0, 1} |E| an indicator vector of the selected edges. The cost function is then indeed the inner product c(w, y) = w • φ(y).

The task to solve during back-propagation is the following. We receive the gradient dL/dy of the global loss L with respect to solver output y at a given point ŷ = y( ŵ). We are expected to return dL/dw, the gradient of the loss with respect to solver input w at a point ŵ.

Since Y is finite, there are only finitely many values of y(w). In other words, this function of w is piecewise constant and the gradient is identically zero or does not exist (at points of jumps). This should not come as a surprise; if one does a small perturbation to edge weights of a graph, one usually does not change the optimal TSP tour and on rare occasions alters it drastically. This has an important consequence:

The fundamental problem with differentiating through combinatorial solvers is not the lack of differentiability; the gradient exists almost everywhere. However, this gradient is a constant zero and as such is unhelpful for optimization.

f (y(w))

f λ (w) w W λ eq W λ dif g 1 g 2 (a) δ f (y(w)) f λ (w) w f (y 1 )

f (y 2 )

g 3 P 1 G (b)

Figure 2 : Continuous interpolation of a piecewise constant function. (a) f λ for a small value of λ; the set W λ eq is still substantial and only two interpolators g 1 and g 2 are incomplete. Also, all interpolators are 0-interpolators. (b) f λ for a high value of λ; most interpolators are incomplete and we also encounter a δ-interpolator g 3 (between y 1 and y 2 ) which attains the value f (y 1 ) δ-away from the set P 1 . Despite losing some local structure for high λ, the gradient of f λ is still informative.

First, we simplify the situation by considering the linearization f of L at the point ŷ. Then for

f (y) = L(ŷ) + dL dy (ŷ) • (y -ŷ) we have df y(w) dw = dL dw

and therefore it suffices to focus on differentiating the piecewise constant function f y(w) .

If the piecewise constant function at hand was arbitrary, we would be forced to use zero-order gradient estimation techniques such as computing finite differences. These require prohibitively many function evaluations particularly for high-dimensional problems.

However, the function f y(w) is a result of a minimization process and it is known that for smooth spaces Y there are techniques for such "differentiation through argmin" (Schmidt & Roth, 2014; Samuel & Tappen, 2009; Foo et al., 2008; Domke, 2012; Amos et al., 2017; Amos & Kolter, 2017) . It turns out to be possible to build -with different mathematical tools -a viable discrete analogy. In particular, we can efficiently construct a function f λ (w), a continuous interpolation of f y(w) , whose gradient we return (see Fig. 2 ). The hyper-parameter λ > 0 controls the trade-off between "informativeness of the gradient" and "faithfulness to the original function".

Before diving into the formalization, we present the final algorithm as listed in Algo. 1. It is simple to implement and the backward pass indeed only runs the solver once on modified input. Providing the justification, however, is not straightforward, and it is the subject of the rest of the section. 

:= ŵ + λ • dL dy (ŷ) // Calculate perturbed weights y λ := Solver(w ) return ∇ w f λ ( ŵ) := -1 λ ŷ -y λ // Gradient of continuous interpolation

3.1. CONSTRUCTION AND PROPERTIES OF f λ:

Before we give the exact definition of the function f λ , we formulate several requirements on it. This will help us understand why f λ (w) is a reasonable replacement for f y(w) and, most importantly, why its gradient captures changes in the values of f . Property A1. For each λ > 0, f λ is continuous and piecewise affine.

The second property describes the trade-off induced by changing the value of λ. For λ > 0, we define sets W λ eq and W λ dif as the sets where f y(w) and f λ (w) coincide and where they differ, i.e.

W λ eq = w ∈ W : f λ (w) = f y(w) and W λ dif = W \ W λ eq .

Property A2. The sets W λ dif are monotone in λ and they vanish as λ → 0 + , i.e.

W λ1 dif ⊆ W λ2 dif for 0 < λ 1 ≤ λ 2 and W λ dif → ∅ as λ → 0 + .

In other words, Property A2 tells us that λ controls the size of the set where f λ deviates from f and where f λ has meaningful gradient. This behaviour of f λ can be seen in Fig. 2 .

In the third and final property, we want to capture the interpolation behavior of f λ . For that purpose, we define a δ-interpolator of f . We say that g, defined on a set G ⊂ W , is a δ-interpolator of f between y 1 and y 2 ∈ Y , if

• g is non-constant affine function;

• the image g(G) is an interval with endpoints f (y 1 ) and f (y 2 );

• g attains the boundary values f (y 1 ) and f (y 2 ) at most δ-far away from where f (y(w)) does.

In particular, there is a point

w k ∈ G for which g(w k ) = f (y k ) and dist(w k , P k ) ≤ δ, where P k = {w ∈ W : y(w) = y k }, for k = 1, 2.

In the special case of a 0-interpolator g, the graph of g connects (in a topological sense) two components of the graph of f y(w) . In the general case, δ measures displacement of the interpolator (see also Fig. 2 for some examples). This displacement on the one hand loosens the connection to f y(w) but on the other hand allows for less local interpolation which might be desirable.

Property A3. The function f λ consists of finitely many (possibly incomplete) δ-interpolators of f on W λ dif where δ ≤ Cλ for some fixed C. Equivalently, the displacement is linearly controlled by λ.

Intuitively, the consequence of Property A3 is that f λ has reasonable gradients everywhere since it consists of elementary affine interpolators.

For defining the function f λ , we need a solution of a perturbed optimization problem y λ (w) = arg min y∈Y {c(w, y) + λf (y)}.

(3)

Theorem 1. Let λ > 0. The function f λ defined by

EQUATION

satisfies Properties A1, A2, A3.

Let us remark that already the continuity of f λ is not apparent from its definition as the first term f y λ (w) is still a piecewise constant function. Proof of this result, along with geometrical description of f λ , can be found in section A.2. Fig. 3 visualizes f λ for different values if λ. Now, since f λ is ensured to be differentiable, we have

EQUATION

The second equality then holds due to (2). We then return ∇f λ as a loss gradient.

Remark 1. The roots of the method we propose lie in loss-augmented inference. In fact, the update rule from (5) (but not the function f λ or any of its properties) was already proposed in a different context in (Hazan et al., 2010; Song et al., 2016) and was later used in (Lorberbom et al., 2018; Mohapatra et al., 2018) . The main difference to our work is that only the case of λ → 0 + is recommended and studied, which in our situation computes the correct but uninformative zero gradient. Our analysis implies that larger values of λ are not only sound but even preferable. This will be seen in experiments where we use values λ ≈ 10 -20. In other words, y λ is the output of calling the solver on input w .

4. EXPERIMENTS:

In this section, we experimentally validate a proof of concept: that architectures containing exact blackbox solvers (with backward pass provided by Algo. 1) can be trained by standard methods. To that end, we solve three synthetic tasks as listed in Tab. 1. These tasks are designed to mimic practical examples from Section 2 and solving them anticipates a two-stage process: 1) extract suitable features from raw input, 2) solve a combinatorial problem over the features. The dimensionalities of input and of intermediate representations also aim to mirror practical problems and are chosen to be prohibitively large for zero-order gradient estimation methods. Guidelines of setting the hyperparameter λ are given in section A.1.

We include the performance of ResNet18 (He et al., 2016) as a sanity check to demonstrate that the constructed datasets are too complex for standard architectures. Remark 2. The included solvers have very efficient implementations and do not severely impact runtime. All models train in under two hours on a single machine with 1 GPU and no more than 24 utilized CPU cores. Only for the large TSP problems the solver's runtime dominates.

4.1. WARCRAFT SHORTEST PATH:

Problem input and output. The training dataset for problem SP(k) consists of 10000 examples of randomly generated images of terrain maps from the Warcraft II tileset (Guyomarch, 2017) . The maps have an underlying grid of dimension k × k where each vertex represents a terrain with a fixed cost that is unknown to the network. The shortest (minimum cost) path between top left and bottom right vertices is encoded as an indicator matrix and serves as a label (see also Fig. 4 ). We consider datasets SP(k) for k ∈ {12, 18, 24, 30}. More experimental details are provided in section A.3. Architecture. An image of the terrain map is presented to a convolutional neural network which outputs a k × k grid of vertex costs. These costs are then the input to the Dijkstra algorithm to compute the predicted shortest path for the respective map. The loss used for computing the gradient update is the Hamming distance between the true shortest path and the predicted shortest path. Test % 12 99.7 ± 0.0 96.0 ± 0.3 100.0 ± 0.0 23.0 ± 0.3 18 98.9 ± 0.2 94.4 ± 0.2 99.9 ± 0.0 0.7 ± 0.3 24 97.8 ± 0.2 94.4 ± 0.6 100.0 ± 0.0 0.0 ± 0.0 30 97.4 ± 0.1 94.0 ± 0.3 95.6 ± 0.5 0.0 ± 0.0

Input → Label      1 0 • • • 0 1 0 • • • 0 . . . . . . . . . . . . 0 0 • • • 1      k × k indicator matrix of shortest path

Results. Our method learns to predict the shortest paths with high accuracy and generalization capability, whereas the ResNet18 baseline unsurprisingly fails to generalize already for small grid sizes of k = 12. Since the shortest paths in the maps are often nonunique (i.e. there are multiple shortest paths with the same cost), we report the percentage of shortest path predictions that have optimal cost. The results are summarized in Tab. 2.

4.2. GLOBE TRAVELING SALESMAN PROBLEM:

Problem input and output. The training dataset for problem TSP(k) consists of 10000 examples where the input for each example is a k-element subset of fixed 100 country flags and the label is the shortest traveling salesman tour through the capitals of the corresponding countries. The optimal tour is represented by its adjacency matrix (see also Fig. 5 ). We consider datasets TSP(k) for k ∈ {5, 10, 20, 40}. Architecture. Each of the k flags is presented to a convolutional network that produces k threedimensional vectors. These vectors are projected onto the unit sphere in R 3 ; a representation of the globe. The TSP solver receives a matrix of pairwise distances of the k computed locations. The loss of the network is the Hamming distance between the true and the predicted TSP adjacency matrix. The architecture is expected to learn the correct representations of the flags (i.e. locations of the respective countries' capitals on Earth, up to rotations of the sphere). The employed Gurobi solver optimizes a mixed-integer programming formulation of TSP using the cutting plane method (Marchand et al., 2002) for lazy sub-tour elimination. Train % Test % Train % Test % 5 99.8 ± 0.0 99.2 ± 0.1 100.0 ± 0.0 1.9 ± 0.6 10 99.8 ± 0.1 98.7 ± 0.2 99.0 ± 0.1 0.0 ± 0.0 20 99.1 ± 0.1 98.4 ± 0.4 98.8 ± 0.3 0.0 ± 0.0 40 97.4 ± 0.2 96.7 ± 0.4 96.9 ± 0.3 0.0 ± 0.0 Results. This architecture not only learns to extract the correct TSP tours but also learns the correct representations. Quantitative evidence is presented in Tab. 3, where we see that the learned locations generalize well and lead to correct TSP tours also on the test set and also on somewhat large instances (note that there are 39! ≈ 10 46 admissible TSP tours for k = 40). The baseline architecture only memorizes the training set. Additionally, we can extract the suggested locations of world capitals and compare them with reality. To that end, we present Fig. 5b , where the learned locations of 10 capitals in Southeast Asia are displayed.

EQUATION

4.3. MNIST MIN-COST PERFECT MATCHING:

Problem input and output. The training dataset for problem PM(k) consists of 10000 examples where the input to each example is a set of k 2 digits drawn from the MNIST dataset arranged in a k × k grid. For computing the label, we consider the underlying k × k grid graph (without diagonal edges) and solve a MIN-COST-PERFECT-MATCHING problem, where edge weights are given simply by reading the two vertex digits as a two-digit number (we read downwards for vertical edges and from left to right for horizontal edges). The optimal perfect matching (i.e. the label) is encoded by an indicator vector for the subset of the selected edges, see example in Fig. 6 .

Architecture. The grid image is the input of a convolutional neural network which outputs a grid of vertex weights. These weights are transformed into edge weights as described above and given to the solver. The loss function is Hamming distance between solver output and the true label. Table 4 : Results for MNIST Min-cost perfect matching. Reported is the accuracy of predicting an optimal matching. Standard deviations are over five restarts.

Embedding Blossom V ResNet18 k:

Train % Test % Train % Test % 4 99.97 ± 0.01 98.32 ± 0.24 100.0 ± 0.0 92.5 ± 0.3 8 99.95 ± 0.04 99.92 ± 0.01 100.0 ± 0.0 8.3 ± 0.8 16 99.02 ± 0.84 99.06 ± 0.57 100.0 ± 0.0 0.0 ± 0.0 24 95.63 ± 5.49 92.06 ± 7.97 96.1 ± 0.5 0.0 ± 0.0 Results. The architecture containing the solver is capable of good generalizations suggesting that the correct representation is learned.

The performance is good even on larger instances and despite the presence of noise in supervision -often there are many optimal matchings. In contrast, the ResNet18 baseline only achieves reasonable performance for the simplest case PM(4). The results are summarized in Tab. 4.

5. DISCUSSION:

We provide a unified mathematically sound algorithm to embed combinatorial algorithms into neural networks. Its practical implementation is straightforward and training succeeds with standard deep learning techniques. The two main branches of future work are: 1) exploring the potential of newly enabled architectures, 2) addressing standing real-world problems. The latter case requires embedding approximate solvers (that are common in practice). This breaks some of our theoretical guarantees but given their strong empirical performance, the fusion might still work well in practice. 

Input →        0 1 . . . 1 0        (a) (b)

A APPENDIX:

A.1 GUIDELINES FOR SETTING THE VALUES OF λ.

In practice, λ has to be chosen appropriately, but we found its exact choice uncritical (no precise tuning was required). Nevertheless, note that λ should cause a noticeable disruption in the optimization problem from equation ( 3), otherwise it is too likely that y(w) = y λ (w) resulting in a zero gradient.

In other words, λ should roughly be of the magnitude that brings the two terms in the definition of w in Prop. 1 to the same order:

λ ≈ w dL dy

where • stands for the average. This again justifies that λ is a true hyperparameter and that there is no reason to expect values around λ → 0 + .

A.2 PROOFS:

Proof of Proposition 1. Let us write L = L(ŷ) and ∇L = dL dy (ŷ), for brevity. Thanks to the linearity of c and the definition of f , we have

c( ŵ, y) + λf (y) = ŵy + λ L + ∇L(y -ŷ) = ( ŵ + λ∇L)y + λL -λ∇Lŷ = c(w , y) + c 0 ,

where c 0 = λL -λ∇Lŷ and w = ŵ + λ∇L as desired. The conclusion about the points of minima then follows.

Before we prove Theorem 1, we make some preliminary observations. To start with, due to the definition of the solver, we have the fundamental inequality c(w, y) ≥ c w, y(w) for every w ∈ W and y ∈ Y .

(6)

Observation 1. The function w → c w, y(w) is continuous and piecewise linear.

Proof. Since c's are linear and distinct, c w, y(w) , as their pointwise minimum, has the desired properties.

Analogous fundamental inequality c(w, y) + λf (y) ≥ c w, y λ (w) + λf y λ (w) for every w ∈ W and y ∈ Y (7)

follows from the definition of the solution to the optimization problem (3).

A counterpart of Observation 1 reads as follows.

Observation 2. The function w → c w, y λ (w) + λf y λ (w) is continuous and piecewise affine.

Proof. The function under inspection is a pointwise minimum of distinct affine functions w → c(w, y) + λf (y) as y ranges Y .

As a consequence of above-mentioned fundamental inequalities, we obtain the following two-sided estimates on f λ . Observation 3. The following inequalities hold for w ∈ W f y λ (w) ≤ f λ (w) ≤ f y(w) .

Proof. Inequality (6) implies that c w, y(w) -c w, y λ (w) ≤ 0 and the first inequality then follows simply from the definition of f λ . As for the second one, it suffices to apply (7) to y = y(w).

Now, let us introduce few notions that will be useful later in the proofs. For a fixed λ, W partitions into maximal connected sets P on which y λ (w) is constant (see Fig. 7 ). We denote this collection of sets by W λ and set W = W 0 .

For λ ∈ R and

y 1 = y 2 ∈ Y , we denote F λ (y 1 , y 2 ) = w ∈ W : c(w, y 1 ) + λf (y 1 ) = c(w, y 2 ) + λf (y 2 ) .

We write F (y 1 , y 2 ) = F 0 (y 1 , y 2 ), for brevity. For technical reasons, we also allow negative values of λ here.

(a) The situation for λ = 0. We can see the polytope P on which y(w) attains y1 ∈ Y . The boundary of P is composed of segments of lines F (y1, y k ) for k = 2, . . . , 5.

(b) The same situation is captured for some relatively small λ > 0. Each line F λ (y1, y k ) is parallel to its corresponding F (y1, y k ) and encompasses a convex polytope in W λ .

Figure 7 : The family W λ of all maximal connected sets P on which y λ is constant.

Note, that if W = R N , then F λ is a hyperplane since c's are linear. In general, W may just be a proper subset of R N and, in that case, F λ is just the restriction of a hyperplane onto W . Consequently, it may happen that F λ (y 1 , y 2 ) will be empty for some pair of y 1 , y 2 and some λ ∈ R.

To emphasize this fact, we say "hyperplane in W ". Analogous considerations should be taken into account for all other linear objects. The note "in W " stands for the intersection of these linear object with the set W .

Observation 4. Let P ∈ W λ and let y λ (w) = y for w ∈ P . Then P is a convex polytope in W , where the facets consist of parts of finitely many hyperplanes F λ (y, y k ) in W for some {y k } ⊂ Y .

Proof. Assume that W = R N . The values of y λ may only change on hyperplanes of the form F λ (y, y ) for some y ∈ Y . Then P is an intersection of corresponding half-spaces and therefore P is a convex polytope. If W is a proper subset of R N the claim follows by intersecting all the objects with W .

Observation 5. Let y 1 , y 2 ∈ Y be distinct. If nonempty, the hyperplanes F (y 1 , y 2 ) and F λ (y 1 , y 2 ) are parallel and their distance is equal to |λ|K(y 1 , y 2 ), where

K(y 1 , y 2 ) = |f (y 1 ) -f (y 2 )| y 1 -y 2 .

Proof. If we define a function c(w) = c(w, y 1 ) -c(w, y 2 ) = w(y 1 -y 2 ) and a constant C = f (y 2 ) -f (y 1 ), then our objects rewrite to

F (y 1 , y 2 ) = {w ∈ W : c(w) = 0} and F λ (y 1 , y 2 ) = {w ∈ W : c(w) = λC}.

Since c is linear, these sets are parallel and F (y 1 , y 2 ) intersects the origin. Thus, the required distance is the distance of the hyperplane F λ (y 1 , y 2 ) from the origin, which equals to |λC|/ y 1y 2 .

As the set Y is finite, there is a uniform upper bound K on all values of K(y 1 , y 2 ). Namely

EQUATION

A.2.1 PROOF OF THEOREM 1

Proof of Property A1. Now, Property A1 follows, since

f λ (w) = 1 λ c w, y λ (w) + λf y λ (w) - 1 λ c w, y(w)

and f λ is a difference of continuous and piecewise affine functions.

Proof of Property A2. Let 0 < λ 1 ≤ λ 2 be given. We show that W λ2 eq ⊆ W λ1 eq which is the same as showing W λ1 dif ⊆ W λ2 dif . Assume that w ∈ W λ2 eq , that is, by the definition of W λ2 eq and f λ ,

EQUATION

in which we denoted y 2 = y λ2 (w). Our goal is to show that c w, y(w)

EQUATION

where y 1 = y λ1 (w) as this equality then guarantees that w ∈ W λ1 eq . Observe that (7) applied to λ = λ 1 and y = y(w), yields the inequality "≥" in (10).

Let us show the reversed inequality. By Observation 3 applied to λ = λ 1 , we have f y(w) ≥ f (y 1 ).

(11)

We now use (7) with λ = λ 2 and y = y 1 , followed by equality (9) to obtain c(w, y 1 ) + λ 1 f (y 1 ) = c(w, y 1 ) + λ 2 f (y

1 ) + (λ 1 -λ 2 )f (y 1 ) ≥ c(w, y 2 ) + λ 2 f (y 2 ) + (λ 1 -λ 2 )f (y 1 ) = c w, y(w) + λ 2 f y(w) + (λ 1 -λ 2 )f (y 1 ) = c w, y(w) + λ 1 f y(w) + (λ 2 -λ 1 ) f y(w) -f (y 1 ) ≥ c w, y(w) + λ 1 f y(w)

where the last inequality holds due to (11).

Next, we have to show that W λ dif → ∅ as λ → 0 + , i.e. that for almost every w ∈ W , there is a λ > 0 such that w / ∈ W λ dif . To this end, let w ∈ W be given. We can assume that y(w) is a unique solution of solver (1), since two solutions, say y 1 and y 2 , coincide only on the hyperplane F (y 1 , y 2 ) in W , which is of measure zero. Thus, since Y is finite, the constant c = min 

For the remaining y's, (13) holds trivially for every λ > 0. Therefore, y(w) is a solution of the minimization problem (3), whence y λ (w) = y(w). This shows that w ∈ W λ eq as we wished. If d = 0, then f y(w) ≤ f (y) for every y ∈ Y and (13) follows again. The condition on c tells us that h is a non-constant affine function. It follows by the definition of F (y 1 , y 2 ) and F λ (y 1 , y 2 ) that h(w) = 0 if and only if w ∈ F (y 1 , y 2 ) (15) and h(w) = λ f (y 2 ) -f (y 1 ) if and only if w ∈ F λ (y 1 , y 2 ).

(16) By Observation 5, the sets F and F λ are parallel hyperplanes. Denote by G the nonempty intersection of their corresponding half-spaces in W . We show that g is a δ-interpolator of f on G between y 1 and y 2 , with δ being linearly controlled by λ.

We have already observed that g is the affine function ranging from f (y 1 ) -on the set F λ (y 1 , y 2 )to f (y 2 ) -on the set F (y 1 , y 2 ). It remains to show that g attains both the values f (y 1 ) and f (y 2 ) at most δ-far from the sets P 1 and P 2 , respectively, where

P k ∈ W denotes a component of the set {w ∈ W : y(w) = y k }, k = 1, 2.

Consider y 1 first. By Observation 4, there are z 1 , . . . , z ∈ Y , such that facets of P 1 are parts of hyperplanes F (y 1 , z 1 ), . . . , F (y 1 , z ) in W . Each of them separates W into two half-spaces, say W + k and W - k , where W - k is the half-space which contains P 1 and W + k is the other one. Let us denote c k (w) = c(w, y 1 ) -c(w, z k ) for w ∈ W and k = 1, . . . , .

Every c k is a non-zero linear function which is negative on W - k and positive on W + k . By the definition of y 1 , we have c(w, y 1 ) + λf (y 1 ) ≤ c(w, z k ) + λf (z k ) for w ∈ P 1 and for k = 1, . . . , , that is c k (w) ≤ λ f (z k ) -f (y 1 ) for w ∈ P 1 and for k = 1, . . . , . 

Now, denote:

W λ k = w ∈ W : c k (w) ≤ λ f (z k ) -f (y 1 ) for k = 1, . . . , .

Each W λ k is a half-space in W containing W - k and hence P 1 . Let us set P λ 1 = k=1 W λ k . Clearly, P 1 ⊆ P λ 1 (see Fig. 8 ). By Observation 5, the distance of the hyperplane w ∈ W : c k (w) = λ f (z k ) -f (y 1 ) from P 1 is at most λK, where K is given by (8). Therefore, since all the facets of P λ 1 are at most λK far from P 1 , there is a constant C such that each point of P λ 1 is at most Cλ far from P 1 .

Finally, choose any w 1 ∈ P λ 1 ∩ F λ (y 1 , y 2 ). By ( 16), we have g(w 1 ) = f (y 1 ), and by the definition of P λ 1 , w 1 is no farther than Cλ away from P 1 . Now, let us treat y 2 and define the set P λ 2 analogous to P λ 1 , where each occurrence of y 1 is replaced by y 2 . Any w 2 ∈ P λ 2 ∩ F (y 1 , y 2 ) has desired properties. Indeed, (15) ensures that g(w 2 ) = f (y 2 ) and w 2 is at most Cλ far away from P 2 .

A.3 DETAILS OF EXPERIMENTS:



A.3.1 WARCRAFT SHORTEST PATH:

The maps for the dataset have been generated with a custom random generation process by using 142 tiles from the Warcraft II tileset (Guyomarch, 2017) . The costs for the different terrain types range from 0.8-9.2. Some example maps of size 18 × 18 are presented in Fig. 9a together with a histogram of the shortest path lengths. We used the first five layers of ResNet18 followed by a max-pooling operation to extract the latent costs for the vertices. The dataset consists of randomly generated grids of MNIST digits that are sampled from a subset of 1000 digits of the full MNIST dataset. We trained a fully convolutional neural network with two convolutional layers followed by a max-pooling operation that outputs a k × k grid of vertex costs for each example. The vertex costs are transformed into the edge costs via the known cost function and the edge costs are then the inputs to the Blossom V solver (Edmonds, 1965) as implemented in (Kolmogorov, 2009) .

Regarding the optimization procedure, we employed the Adam optimizer along with scheduled learning rate drops dividing the learning rate by 10 at epochs 10 and 20, respectively. Other training details are in Tab. 6. Lower batch sizes were used to reduce GPU memory requirements. For the Globe Traveling Salesman Problem we used a convolutional neural network architecture of three convolutional layers and two fully connected layers. The last layer outputs a vector of dimension 3k containing the k 3-dimensional representations of the respective countries' capital cities. These representations are projected onto the unit sphere and the matrix of pairwise distances is fed to the TSP solver.

The high combinatorial complexity of TSP has negative effects on the loss landscape and results in many local minima and high sensitivity to random restarts. For reducing sensitivity to restarts, we set Adam parameters to β 1 = 0.5 (as it is done for example in GAN training (Radford et al., 2015) ) and = 10 -3 .

The local minima correspond to solving planar TSP as opposed to spherical TSP. For example, if all cities are positioned to almost identical locations, the network can still make progress but it will never have the incentive to spread the cities apart in order to reach the global minimum. To mitigate that, we introduce a repellent force between epochs 15 and 30. In particular, we set

L rep = E i =j e -xi-xj

where x i ∈ R 3 for i = 1, . . . , k are the positions of the k cities on the unit sphere. The regularization constants C k were chosen as 2.0, 3.0, 6.0, and 20.0 for k ∈ {5, 10, 20, 40}.

For fine-tuning we also introduce scheduled learning rate drops where we divide the learning rate by 10 at epochs 80 and 90. [20, 4, 2] , [50, 4, 2] , 500] 100 50 20

In Fig. 5b , we compare the true city locations with the ones learned by the hybrid architecture. Due to symmetries of the sphere, the architecture can embed the cities in any rotated or flipped fashion. We resolve this by computing "the most favorable" isometric transformation of the suggested locations.

In particular, we solve the orthogonal Procrustes problem (Gower & Dijksterhuis, 2004 )

R * = arg min R:R T R=I RX -Y 2

where X are the suggested locations, Y the true locations, and R * the optimal transformation to apply. We report the resulting offsets in kilometers in Tab. 8. Since approximate solvers often appear in practice where the combinatorial instances are too large to be solved exactly in reasonable time, we test our method also in this setup. In particular, we use the approximate solver (OR-Tools (ort, 2019)) for the Globe TSP. We draw two conclusions from the numbers presented below in Tab. 9.

(i) The choice of the solver matters. Even if OR-Tools is fed with the ground truth representations (i.e. true locations) it does not achieve perfect results on the test set (see the right column). We expect, that also in practical applications, running a suboptimal solver (e.g. a differentiable relaxation) substantially reduces the maximum attainable performance.

(ii) The suboptimality of the solver didn't harm the feature extraction -the point of our method. Indeed, the learned locations yield performance that is close to the upper limit of what the solver allows (compare the middle and the right column). 

:



Figures and Tables:

Figure 1: Architecture design enabled by Theorem 1. Blackbox combinatorial solver embedded into a neural network.

Forward and Backward Pass function FORWARDPASS( ŵ) ŷ := Solver( ŵ) // ŷ = y( ŵ) save ŵ and ŷ for backward pass return ŷ function BACKWARDPASS( dL dy (ŷ), λ) load ŵ and ŷ from forward pass w

Figure 3: Example f λ for w ∈ R 2 and λ = 3, 10, 20 (left to right). As λ changes, the interpolation f λ is less faithful to the piecewise constant f (y(w)) but provides reasonable gradient on a larger set.

Figure 4: The SP(k) dataset. (a) Each input is a k × k grid of tiles corresponding to a Warcraft II terrain map, the respective label is a the matrix indicating the shortest path from top left to bottom right. (b) is a different map with correctly predicted shortest path.

Figure 5: The TSP(k) problem. (a) illustrates the dataset. Each input is a sequence of k flags and the corresponding label is the adjacency matrix of the optimal TSP tour around the corresponding capitals. (b) displays the learned locations of 10 country capitals in southeast Asia and Australia, accurately recovering their true position.

Figure 6: Visualization of the PM dataset. (a) shows the case of PM(4). Each input is a 4 × 4 grid of MNIST digits and the corresponding label is the indicator vector for the edges in the min-cost perfect matching. (b) shows the correct min-cost perfect matching output from the network. The cost of the matching is 348 (46 + 12 horizontally and 27 + 45 + 40 + 67 + 78 + 33 vertically).

0, set λ < c/d. Then, for every y ∈ Y such that f y(w) > f (y), we have λ < c(w, y) -c w, y(w) f y(w) -f (y) which rewrites c w, y(w) + λf y(w) < c(w, y) + λf (y).

Proof of Property A3. Let y 1 = y 2 ∈ Y be given. We show that on the component of the set{w ∈ W : y(w) = y 1 and y λ (w) = y 2 } (14)the function f λ agrees with a δ-interpolator, where δ ≤ Cλ and C > 0 is an absolute constant. The claim follows as there are only finitely many sets and their components of the form (14) in W λ dif . Let us set h(w) = c(w, y 1 ) -c(w, y 2 ) for w ∈ W and g(w) = f (y 2 ) -

(a) The facets of P1 consist of parts of hyperplanes F (y1, z k ) in W . Each facet F (y1, z k ) has its corresponding shifts F λ and F -λ , from which only one intersects P . The polytope P λ 1 is then bounded by those outer shifts. (b) The interpolator g attains the value f (y1) on a part of F λ (y1, y2) -a border of the domain G. The value f (y2) is attained on a part of F (y1, y2) -the second border of the strip G.

Figure 8: The polytopes P 1 and P λ 1 and the interpolator g.

the shortest path distribution in the training set. All possible path lengths (18-35) occur.

Figure 9: Warcraft SP(18) dataset.

Experiments Overview.

Table content: <table><tr><td>Graph Problem</td><td>Solver</td><td>Solver instance size</td><td>Input format</td></tr><tr><td>Shortest path</td><td>Dijkstra</td><td>up to 900 vertices</td><td>(image) up to 240 × 240</td></tr><tr><td>Min Cost PM</td><td>Blossom V</td><td>up to 1104 edges</td><td>(image) up to 528 × 528</td></tr><tr><td>Traveling Salesman</td><td>Gurobi</td><td>up to 780 edges</td><td>up to 40 images (20 × 40)</td></tr></table>

Results for Warcraft shortest path. Reported is the accuracy, i.e. percentage of paths with the optimal costs. Standard deviations are over five restarts.

Table content: <table><tr><td/><td colspan="2">Embedding Dijkstra</td><td>ResNet18</td></tr><tr><td>k</td><td>Train %</td><td>Test %</td><td>Train %</td></tr></table>

Results for Globe TSP. Reported is the full tour accuracy. Standard deviations are over five restarts.

Table content: <table><tr><td>Embedding TSP Solver</td><td>ResNet18</td></tr><tr><td>k</td><td/></tr></table>

Oriol Vinyals, Meire Fortunato, and Navdeep Jaitly. Pointer networks. In Advances in Neural Information Processing Systems 28 (NIPS'15), pp. 2692-2700. Curran Associates, Inc., 2015. Po-Wei Wang, Priya L. Donti, Bryan Wilder, and Zico Kolter. SATNet: Bridging deep learning and logical reasoning using a differentiable satisfiability solver. arXiv, 1905.12149, 2019. URL http://arxiv.org/abs/1905.12149. Advances in Neural Information Processing Systems 32, pp. 4674-4685. Curran Associates, Inc., 2019. URL http://papers.nips.cc/paper/ 8715-end-to-end-learning-and-optimization-on-graphs.pdf. Mang Ye, Andy J. Ma, Liang Zheng, Jiawei Li, and Pong C. Yuen. Dynamic label graph matching for unsupervised video re-identification. In IEEE International Conference on Computer Vision (ICCV'17). IEEE Computer Society, Oct 2017. Shuai Zheng, Sadeep Jayasumana, Bernardino Romera-Paredes, Vibhav Vineet, Zhizhong Su, Dalong Du, Chang Huang, and Philip H. S. Torr. Conditional random fields as recurrent neural networks. In IEEE International Conference on Computer Vision (ICCV'15), pp. 1529-1537. IEEE Computer Society, 2015.

Table content: <table><tr><td colspan="2">Bryan Wilder, Eric Ewing, Bistra Dilkina, and Milind Tambe. End to end learning and</td></tr><tr><td>optimization on graphs.</td><td>In H. Wallach, H. Larochelle, A. Beygelzimer, F. d Alche-</td></tr><tr><td colspan="2">Buc, E. Fox, and R. Garnett (eds.),</td></tr></table>

Optimization was carried out via Adam optimizer(Kingma & Ba, 2014) with scheduled learning rate drops dividing the learning rate by 10 at epochs 30 and 40. Hyperparameters and model details are listed in Tab. 5 Experimental setup for Warcraft Shortest Path.

Table content: <table><tr><td>k</td><td>Optimizer(LR)</td><td>Architecture</td><td colspan="2">Epochs Batch Size λ</td></tr><tr><td colspan="3">12, 18, 24, 30 Adam(5 × 10 -4 ) subset of ResNet18</td><td>50</td><td>70</td><td>20</td></tr></table>

Experimental setup for MNIST Min-cost Perfect Matching.

Table content: <table><tr><td>k</td><td>Optimizer(LR)</td><td>Architecture [channels, kernel size, stride]</td><td colspan="3">Epochs Batch Size λ</td></tr><tr><td>4, 8</td><td>Adam(10 -3 )</td><td>[[20, 5, 1], [20, 5, 1]]</td><td>30</td><td>70</td><td>10</td></tr><tr><td>16</td><td>Adam(10 -3 )</td><td>[[50, 5, 1], [50, 5, 1]]</td><td>30</td><td>40</td><td>10</td></tr><tr><td>24</td><td>Adam(10 -3 )</td><td>[[50, 5, 1], [50, 5, 1]]</td><td>30</td><td>30</td><td>10</td></tr><tr><td colspan="3">A.3.3 GLOBE TRAVELING SALESMAN PROBLEM</td><td/><td/><td/></tr></table>

Experimental setup for the Globe Traveling Salesman Problem.

Table content: <table><tr><td/><td/><td>Architecture</td><td/><td/><td/></tr><tr><td>k</td><td>Optimizer(LR)</td><td>[channels, kernel size, stride],</td><td colspan="3">Epochs Batch Size λ</td></tr><tr><td/><td/><td>linear layer size</td><td/><td/><td/></tr><tr><td>5, 10, 20</td><td>Adam(10 -4 )</td><td>[[20, 4, 2], [50, 4, 2], 500]</td><td>100</td><td>50</td><td>20</td></tr><tr><td>40</td><td>Adam(5 × 10 -5 )</td><td>[</td><td/><td/><td/></tr></table>

Average errors of city placement on the Earth. 69 ± 11 19 ± 5 11 ± 5 58 ± 7 A.4 TRAVELING SALESMAN WITH AN APPROXIMATE SOLVER

Table content: <table><tr><td>k</td><td>5</td><td>10</td><td>20</td><td>40</td></tr><tr><td>Location offset (km)</td><td/><td/><td/><td/></tr></table>

Perfect path accuracy for Globe TSP using the approximate solver OR-Tools (ort, 2019). The maximal achievable performance is in the right column, where the solver uses the ground truth city locations.

Table content: <table><tr><td/><td colspan="2">Embedding OR-tools</td><td>OR-tools on GT locations</td></tr><tr><td>k</td><td>Train %</td><td>Test %</td><td>Test %</td></tr><tr><td colspan="3">5 99.8 ± 0.0 99.3 ± 0.1</td><td>100.0</td></tr><tr><td colspan="3">10 84.3 ± 0.2 84.4 ± 0.2</td><td>88.6</td></tr><tr><td colspan="3">20 49.2 ± 0.2 48.6 ± 0.8</td><td>54.4</td></tr><tr><td colspan="3">40 14.6 ± 0.1 15.1 ± 0.3</td><td>15.2</td></tr></table>

Peer review (iclr2020):
This paper shows how end-to-end learning can be done through combinatorial solvers by using the derivative of continuous surrogate function in the backward pass. One elegant part of the method is that no modification or relaxation is done to the combinatorial solver in the forward pass and that the backward pass just requires another call to the blackbox solver. The idea of constructing continuous surrogate functions and using them for differentiating through solvers with piecewise-constant output spaces is thought-provoking and I can see it inspiring many new directions of work. For example looking at Figure 2 for intuition, one could imagine other ways of making the solution space continuous. The solution space of linear programs over continuous spaces, as considered in [Elmachtoub & Grigas], the Sudoku example in [Amos & Kolter], and related papers, is also piecewise constant and it seems like a similar method could be used to bring more informative derivative information to linear programs --- have you considered this as a future direction? One of my concerns with this work is that the ResNet baseline in the experimental results seems like too much of a straw man for the tasks. I do not see why they should have the capacity to generalize well. This paper shows the ResNet baseline achieve near-zero test accuracy but doesn*t compare to other relevant baselines that are mentioned in the related work section: for example [Bello et al, Deudon et al., Kool et al.] for the TSP. And one smaller comment: If one wanted to squeeze the performance even more, would starting the training process with a large lamdba and annealing it to zero help? ---- Elmachtoub, A. N., & Grigas, P. Smart *predict, then optimize*. arXiv 2017.