Guillaume Dalle

Léo Baty

Louis Bouvier

Axel Parmentier
email: axel.parmentier@enpc.fr

Learning with Combinatorial Optimization Layers: a Probabilistic Approach

Keywords: combinatorial optimization, machine learning, automatic differentiation, graphs, Julia programming language

Combinatorial optimization (CO) layers in machine learning (ML) pipelines are a powerful tool to tackle data-driven decision tasks, but they come with two main challenges. First, the solution of a CO problem often behaves as a piecewise constant function of its objective parameters. Given that ML pipelines are typically trained using stochastic gradient descent, the absence of slope information is very detrimental. Second, standard ML losses do not work well in combinatorial settings. A growing body of research addresses these challenges through diverse methods. Unfortunately, the lack of well-maintained implementations slows down the adoption of CO layers.

In this paper, building upon previous works, we introduce a probabilistic perspective on CO layers, which lends itself naturally to approximate differentiation and the construction of structured losses. We recover many approaches from the literature as special cases, and we also derive new ones. Based on this unifying perspective, we present InferOpt.jl, an open-source Julia package that 1) allows turning any CO oracle with a linear objective into a differentiable layer, and 2) defines adequate losses to train pipelines containing such layers. Our library works with arbitrary optimization algorithms, and it is fully compatible with Julia's ML ecosystem. We demonstrate its abilities using a pathfinding problem on video game maps as guiding example, as well as three other applications from operations research.

Introduction

Machine learning (ML) and combinatorial optimization (CO) are two essential ingredients of modern industrial processes. While ML extracts meaningful information from noisy data, CO enables decision-making in high-dimensional constrained environments. But in many situations, combining both of these tools is necessary: for instance, we might want to generate predictions from data, and then use those predictions to make optimized decisions. To do that, we need pipelines that contain two types of layers: ML layers and CO layers.

Due to their many possible applications, hybrid ML-CO pipelines currently attract a lot of research interest. The recent reviews by [START_REF] Bengio | Machine Learning for Combinatorial Optimization: A Methodological Tour d'horizon[END_REF] and [START_REF] Kotary | End-to-End Constrained Optimization Learning: A Survey[END_REF] are excellent resources on this topic. Unfortunately, relevant software implementations are scattered across paper-specific repositories, with few tests, minimal documentation and sporadic code maintenance. Not only does this make comparison and evaluation difficult for academic purposes, it also hurts practitioners wishing to experiment with such techniques on real use cases.

Let us discuss a generic hybrid ML-CO pipeline, which includes a CO oracle amid several ML layers:

Input x -----→ § ¦ ¤ ¥ ML layers Objective θ -------→ CO oracle Solution y ------→ § ¦ ¤ ¥ More ML layers Output ----→ (1)
The inference problem consists in predicting an output from a given input. It is solved online, and requires the knowledge of the parameters (weights) for each ML layer. On the other hand, the learning problem aims at finding parameters that lead to "good" outputs during inference. It is solved offline based on a training set that contains several inputs, possibly complemented by target outputs.

In Equation (1), we use the term CO oracle to emphasize that any algorithm may be used to solve the optimization problem, whether it relies on an existing solver or a handcrafted implementation. Conversely, when we talk about a layer, it is implied that we can compute meaningful derivatives using automatic differentiation (AD). Since it may call black box subroutines, an arbitrary CO oracle is seldom compatible with AD. And even when it is, its derivatives are zero almost everywhere, which gives us no exploitable slope information. Therefore, according to our terminology, a CO oracle is not a layer (yet), and the whole point of this paper is to turn it into one.

Modern ML libraries provide a wealth of basic building blocks that allow users to assemble and train complex pipelines. We want to leverage these libraries to create hybrid ML-CO pipelines, but we face two main challenges. First, while ML layers are easy to construct, it is not obvious how to transform a CO oracle into a usable layer. Second, standard ML losses are ill-suited to our setting, because they often ignore the underlying optimization problem.

Our goal is to remove these difficulties. We introduce InferOpt.jl1 , a Julia package which 1) can turn any CO oracle into a layer with meaningful derivatives, and 2) provides structured loss functions that work well with the resulting layers. It contains several state-of-the-art methods that are fully compatible with Julia's AD and ML ecosystem, making CO layers as easy to use as any ML layer. To describe the available methods in a coherent manner, we leverage the unifying concept of probabilistic CO layer, hence the name of our package.

Motivating example

Let us start by giving an example of hybrid ML-CO pipeline. Suppose we want to find shortest paths on a map, but we do not have access to an exact description of the underlying terrain. Instead, all we have are images of the area, which give us a rough idea of the topography and obstacles. To solve our problem, we need a pipeline comprising two layers of very different natures. First, an image processing layer, which is typically implemented as a convolutional neural network (CNN). The CNN is tasked with translating the images into a weighted graph. Second, a CO layer performing shortest path computations on said weighted graph (e.g. using Dijkstra's algorithm).

Map image

Cell costs Shortest path

Figure 1: Pipeline for computing shortest paths on Warcraft maps -data from [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF] This pipeline is exactly the one considered by [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF] and [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF] for pathfinding on video game maps. We illustrate it on Figure 1, and we describe it in more detail in Section 2.3. The goal is to learn appropriate weights for the CNN, so that it feeds accurate cell costs to Dijkstra's algorithm. This is done by minimizing a loss function, such as the distance between the true optimal path and the one we predict.

Our setting

In our hybrid ML-CO pipelines, we consider CO oracles f that solve the following kind of problem:

f : θ -→ argmax v∈V θ ⊤ v (2)
Here, the input θ ∈ R d is the objective direction. Meanwhile, V ⊂ R d (for vertices) denotes a finite set of feasible solutions -which may be exponentially large in d -among which the optimal solution f (θ) shall be selected. For simplicity, we assume that f is single-valued, i.e. that the optimal solution is unique.

The feasible set V and its dimension may depend on the instance. For instance, if Equation ( 2) is a shortest path problem, the underlying graph may change from one input to another. If we wanted to remain generic, we should therefore write V(x) ⊂ R d(x) . To keep notations simple, we omit the dependency in x whenever it is clear from the context. Note that we could also study more general CO oracles given by

argmax v∈V θ ⊤ g(v) + h(v)
where g is any function from an arbitrary finite set V to R d , and h is any function from V to R. As long as the objective is linear in θ, the theory we present generalizes seamlessly. However, for ease of exposition, we keep g(v) = v and h(v) = 0. In this case, Equation ( 2) is equivalent to

argmax v∈conv(V) θ ⊤ v.
Indeed, when the objective is linear in v, it makes no difference to optimize over the convex hull conv(V) instead of optimizing over V.

From an optimization problem to an oracle

It is important to note that the formulation argmax v∈V θ ⊤ v is very generic. Any linear program (LP) or integer linear program (ILP) can be written this way, as long as its feasible set is bounded. Indeed, the optimum of an LP is always reached at a vertex of the polytope, of which there are finitely many. Meanwhile, the optimum of an ILP is always reached at an integral point, or more precisely, at a vertex of the convex hull of the integral points. As a result, Equation (2) encompasses a variety of well-known CO problems related to graphs (paths, flows, spanning trees, coloring), resource management (knapsack, bin packing), scheduling, etc. See [START_REF] Korte | Combinatorial Optimization: Theory and Algorithms[END_REF] for an overview of CO and its applications.

For every one of these problems, dedicated algorithms have been developed over the years, which sometimes exploit the domain structure better than a generic ILP solver (such as Gurobi or SCIP). Thus, we have no interest in restricting the procedure used to compute an optimal solution: we want to pick the best algorithm for each application. That is why the methods discussed in this paper only need to access the CO oracle f as a black box function, without making assumptions about its implementation.

From an oracle to a probability distribution

When using CO oracles within ML pipelines, the first challenge we face is the lack of useful derivatives.

Training often relies on stochastic gradient descent (SGD), so we need to be able to backpropagate loss gradients onto the weights of the ML layers. Unfortunately, since the feasible set V of Equation ( 2) is finite, the CO oracle is a piecewise constant mapping and its derivatives are zero almost everywhere. To recover useful slope information, we seek approximate derivatives, which is where the probabilistic approach comes into play.

To describe it, we no longer think about a CO oracle as a function returning a single element f (θ) from V. Instead, we use it to define a probability distribution p(•|θ) on V. The naive choice would be the Dirac mass p(v|θ) = δ f (θ) (v), but it shares the lack of differentiability of the oracle itself. Thus, our goal is to spread out the distribution p into an approximation p, such that the probability mapping θ -→ p(•|θ) becomes smooth with respect to θ. If we can do that, then the expectation mapping

f : θ -→ E p(•|θ) [V ] = v∈V v p(v|θ), (3)
where it is understood that V ∼ p(•|θ), will be just as smooth. This expectation mapping f is what we take to be our probabilistic CO layer : see Section 3 for detailed examples. In what follows, plain letters (p, f ) always refer to the initial CO oracle, while letters with a hat ( p, f ) refer to the probabilistic CO layer that we wrap around it.

From a probability distribution to a loss function

The presence of CO oracles in ML pipelines gives rise to a second challenge: the choice of an appropriate loss function to learn the parameters. As highlighted by [START_REF] Bengio | Machine Learning for Combinatorial Optimization: A Methodological Tour d'horizon[END_REF], this choice heavily depends on the data at our disposal. They distinguish two main paradigms, which we illustrate using the pipeline of Figure 1.

If our dataset only contains the map images, then we are in a weakly supervised setting, which they call learning by experience (see Section 4). In that case, the loss function will evaluate the solutions computed by our pipeline using the true cell costs. On the other hand, if our dataset happens to contain precomputed targets such as the true shortest paths, then we are in a fully supervised setting, which they call learning by imitation (see Section 5). In that case, the loss function will compare the paths computed by our pipeline with the optimal ones, hoping to minimize the discrepancy. For both of these cases, the probabilistic perspective plays an important role in ensuring smoothness of the loss.

Complete pipeline

The typical pipeline we will focus on is a special case of Equation ( 1), which we now describe in more detail:

Input x -----→ § ¦ ¤ ¥ ML layer φ w Objective θ = φw(x) -------------→ CO oracle f Solution y = f (θ) -----------→ £ ¢ ¡ Loss function L (4)
Our pipeline starts with an ML layer φ w , where w stands for the vector of weights. Its role is to encode the input x into an objective direction θ = φ w (x), which is why we often refer to it as the encoder. Then, the CO oracle defined in Equation ( 2) returns an optimal solution y = f (θ). Finally, the loss function L is used during training to evaluate the quality of the solution. If our dataset contains N input samples x (1) , ..., x (N ) , training occurs by applying SGD to the following loss minimization problem:

min w 1 N N i=1 L y (i) f φ w (x (i) ) θ (i) , . . . . (5)
The dots . . . correspond to additional arguments that may be used by the loss. For instance, the loss may depend on the input x (i) itself, or require targets t(i) for comparison (see Section 5).

Remark 1.1. Our work focuses on individual layers and loss functions. Although we present various concrete examples, we do not give generic advice on how to build the whole pipeline for a specific application. If the use case corresponds to a "predict, then optimize" setting such as the one from Figure 1, then [START_REF] Elmachtoub | Smart "Predict, Then Optimize[END_REF] give a few useful pointers. If the goal is to approximate hard optimization problems with easier ones, the reader can refer to Parmentier (2021a) for a general methodology.

Contributions

Our foremost contribution is the open-source package InferOpt.jl, which is written in the Julia programming language [START_REF] Bezanson | Julia: A Fresh Approach to Numerical Computing[END_REF]. Given a CO oracle provided as a callable object, our package wraps it into a probabilistic CO layer that is compatible with Julia's AD and ML ecosystem. This is achieved thanks to the ChainRules.jl2 interface [START_REF] White | [END_REF]. Moreover, InferOpt.jl defines several structured loss functions, both for learning by experience and for learning by imitation.

On top of that, we present theoretical insights that fill some gaps in previous works. In addition to the framework of probabilistic CO layers, we propose:

• A new perturbation technique designed for CO oracles that only accept objective vectors with a certain sign (such as Dijkstra's algorithm, which fails on graphs with negative edge costs): see Section 3.3.2.

• A way to differentiate through a large subclass of probabilistic CO layers (those that rely on convex regularization) by combining the Frank-Wolfe algorithm with implicit differentiation: see Section 3.3.3.

• A probabilistic regularization of the regret for learning by experience: see Section 4.

• A generic decomposition framework for imitation losses, which subsumes most of the literature so far and suggests ways to build new loss functions: see Section 5.1.

Finally, we describe numerical experiments on our motivating example of Warcraft shortest paths, as well as three combinatorial optimization problems from operations research: the stochastic vehicle scheduling problem, the single-machine scheduling problem, and the two-stage stochastic minimum weight spanning tree problem. Here are a few highlights:

• We benchmark and show the strengths of the different learning methods proposed in the literature on hybrid ML-CO pipelines.

• We use the pipeline of Figure 1 for learning by experience on the Warcraft shortest paths problem, even though the CNN encoder has tens of thousands of parameters. To the best of our knowledge, previous attempts to learn such pipelines by experience were restricted to ML layers with fewer than 100 parameters.

• We tackle the stochastic vehicle scheduling problem efficiently by approximating it with its deterministic counterpart.

• We obtain a fast heuristic with state-of-the-art performance for the single machine scheduling problem, which has been a focus of the by the scheduling community for several decades.

• We show that our library enables the use of graph neural networks (GNNs) instead of generalized linear models (GLMs) for instance encoding, bringing additional performance on the two-stage spanning tree problem.

Notations

We write 1 for the vector with all components equal to 1, and e i for the basis vector corresponding to dimension i. The notation 1{E} corresponds to the indicator function of the set (or event) E. The operator ⊙ denotes the Hadamard (componentwise) product between vectors of the same size. We use ∆ d to refer to the unit simplex of dimension d, and E p to denote an expectation with respect to the distribution p. If S is a set, we write

conv(S) = i p i s i : s i ∈ S, p i ≥ 0, i p i = 1 = {E p [S] : p ∈ ∆ S } (6)
for its convex hull and proj S for the orthogonal projection onto S. If h is a real-valued function, we denote by ∇ a h(x) the gradient of h with respect to parameter a at point x and by ∂ a h(x) its convex subdifferential (set of subgradients). The notation dom(h) stands for the domain of h, i.e. the set on which it takes finite values. If h is a vector-valued function, we denote by J a h(x) its Jacobian matrix.

Outline

In Section 2, we review the literature on differentiable optimization layers, before focusing on the Warcraft example. Section 3 introduces the family of probabilistic CO layers by splitting it into perturbed and regularized approaches. Then, Section 4 gives tools for learning by experience, while Section 5 discusses loss functions for learning by imitation. Practical applications of our package are presented in Section 6, before we conclude in Section 7. Proofs for our main theoretical results can be found in Appendix A.

2 Related work

Optimization layers in ML

A significant part of modern ML relies on AD: see [START_REF] Baydin | Automatic Differentiation in Machine Learning: A Survey[END_REF] for an overview and [START_REF] Griewank | Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation[END_REF] for an in-depth treatment. In particular, AD forms the basis of the backpropagation algorithm used to train neural networks.

The notion of implicit layer

Standard neural architectures draw from a small collection of explicit layers [START_REF] Goodfellow | Deep Learning[END_REF]. Whatever their connection structure (dense, convolutional, recurrent, etc.) and regardless of their activation function, these layers all correspond to input-output mappings that can be expressed using an analytic formula. This same formula is then used by AD to compute gradients. On the other hand, the layers defined by InferOpt.jl are of the implicit kind, which means they can contain arbitrarily complex iterative procedures. While we focus here on optimization algorithms, those are not the only kind of implicit layers: fixed point iterations and differential equation solvers are also widely used, depending on the application at hand. See the tutorial by [START_REF] Kolter | Deep Implicit Layers -Neural ODEs[END_REF] for more thorough explanations.

Due to the high computational cost of unrolling iterative procedures, efficient AD of implicit layers often relies on the implicit function theorem. As long as we can specify a set of conditions satisfied by the inputoutput pair, this theorem equates differentiation with solving a linear system of equations. See the Python package jaxopt3 for an example implementation, and its companion paper for theoretical details [START_REF] Blondel | Efficient and Modular Implicit Differentiation[END_REF]. The recent Python package theseus4 showcases an application of this technique to robotics and vision [START_REF] Pineda | Theseus: A Library for Differentiable Nonlinear Optimization[END_REF].

Convex optimization layers

Among the early works on optimization layers for deep learning, the seminal OptNet paper by [START_REF] Amos | OptNet: Differentiable Optimization as a Layer in Neural Networks[END_REF] stands out. It describes a way to differentiate through quadratic programs (QPs) by using the Karush-Kuhn-Tucker (KKT) optimality conditions and plugging them into the implicit function theorem.

More sophisticated tools exist for disciplined conic programs, such as the Python package cvxpylayers5 [START_REF] Agrawal | Differentiable Convex Optimization Layers[END_REF]. The recent Julia package DiffOpt.jl6 [START_REF] Sharma | Flexible Differentiable Optimization via Model Transformations[END_REF]) extends these ideas beyond the conic case to general convex programs. Note that both libraries only accept optimization problems formulated in a domain-specific modeling language, as opposed to arbitrary oracles.

Strong convexity makes differentiation easier because the solutions evolve smoothly as a function of the constraints and objective parameters. In particular, this means the methods listed above return exact derivatives and do not rely on approximations. Regrettably, this nice behavior falls apart as soon as we enter the combinatorial world.

Linear optimization layers

Let us consider an LP whose feasible set is a bounded polyhedron, also called polytope. It is well-known that for most objective directions, the optimal solution will be unique and located at a vertex of the polytope. Even though LPs look like continuous optimization problems, this property shows that they are fundamentally combinatorial. Indeed, a small change in the objective direction can cause the optimal solution to suddenly jump to another vertex, which results in a discontinuous mapping from objectives to solutions. In fact, this mapping is piecewise constant, which means no useful differential information can come from it: its Jacobian is undefined at the jump points and zero everywhere else.

Therefore, when differentiating LPs with respect to their objective parameters, we need to resort to approximations. [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF] use interpolation to turn a piecewise constant mapping into a piecewise linear and continuous one. However, the dominant approximation paradigm in the literature is regularization, as formalized by [START_REF] Blondel | Learning with Fenchel-Young Losses[END_REF].

For instance, [START_REF] Wilder | Melding the Data-Decisions Pipeline: Decision-Focused Learning for Combinatorial Optimization[END_REF] add a quadratic penalty to the linear objective, which allows them to reuse the QP computations of [START_REF] Amos | OptNet: Differentiable Optimization as a Layer in Neural Networks[END_REF]. Mandi and Guns (2020) propose a log-barrier penalty, which lets them draw a connection with interior-point methods. [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF] suggest perturbing the optimization problem by adding stochastic noise to the objective direction, which is a form of implicit regularization.

When LP layers are located at the end of a pipeline, a clever choice of loss function can also simplify differentiation. This is illustrated by the structured support vector machine (S-SVM) loss [START_REF] Nowozin | Structured Learning and Prediction in Computer Vision[END_REF], smart "predict, then optimize" (SPO+) loss [START_REF] Elmachtoub | Smart "Predict, Then Optimize[END_REF] and Fenchel-Young (FY) loss [START_REF] Blondel | Learning with Fenchel-Young Losses[END_REF].

Integer and combinatorial optimization layers

In theory, the methods from the previous section still work in the presence of integer variables, that is, for ILPs. To apply them, we only need to consider the polytope defined by the convex hull of integral solutions.

Alas, in the general case, there is no concise way to describe this convex hull. This is why many authors decide to differentiate through the continuous relaxation of the ILP instead [START_REF] Mandi | Smart Predict-and-Optimize for Hard Combinatorial Optimization Problems[END_REF]: it is an outer approximation of the integral polytope, but it can be sufficient for learning purposes. Some suggest taking advantage of techniques specific to integer programming, such as integrality cuts [START_REF] Ferber | MIPaaL: Mixed Integer Program as a Layer[END_REF] or a generalization of the notion of active constraint [START_REF] Paulus | CombOptNet: Fit the Right NP-Hard Problem by Learning Integer Programming Constraints[END_REF].

There has also been significant progress on finding gradient approximations for combinatorial problems such as ranking [START_REF] Blondel | Fast Differentiable Sorting and Ranking[END_REF]) and shortest paths (Parmentier 2021b). Yet these techniques are problem-specific, and therefore hard to generalize, which is why we leave them aside.

Instead, we want to allow implicit manipulation of the integral polytope itself, without making assumptions on its structure. To achieve that, we can only afford to invoke CO oracles as black boxes [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF][START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF]. Compatibility with arbitrary algorithms is one of the fundamental tenets of InferOpt.jl. The recent Python package PyEPO7 (Tang and Khalil 2022) shares our generic perspective, but it only implements a subset of InferOpt.jl and does not address the notion of probabilistic CO layer.

This probabilistic aspect is central to our proposal, and it is mostly inspired by the works of [START_REF] Blondel | Learning with Fenchel-Young Losses[END_REF] and [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF]. Another related approach is the Implicit Maximum Likelihood Estimation of [START_REF] Niepert | Implicit MLE: Backpropagating Through Discrete Exponential Family Distributions[END_REF], whose probabilistic formulation is slightly different and based on exponential families.

Finally, note that CO layers can be computationally heavy due to the frequent need to recompute optimal solutions at each epoch. Mulamba et al. (2021) address this issue by proposing a noise contrastive approach with solution caching, while [START_REF] Mandi | Decision-Focused Learning: Through the Lens of Learning to Rank[END_REF] build upon this method with a focus on learning to rank.

Similarities and differences with reinforcement learning

As we will see in Section 4, learning by experience can be reminiscent of reinforcement learning (RL), which also relies on a reward or cost signal given by the environment [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]. Furthermore, the encoder layer of Equation ( 4) is similar to a parametric approximation of the value function, which forms the basis of deep RL approaches. This prompts us to discuss a few differences between RL and the framework we study.

The standard mathematical formulation of RL is based on Markov decision processes (MDPs), where a reward and state transition are associated with each action. Usually, the available actions are elementary decisions: pull one lever of a multi-armed bandit, cross one edge on a graph, select one move in a board game, etc. The resulting value or policy update is local: it is specific to both the current state and the action taken. The more reward information we gather, the more efficient learning becomes.

In our framework, the basic step is a call to the optimizer. But combinatorial algorithms can go beyond simple actions: they often output a structured and high-dimensional solution, which aggregates many elementary decisions. This in turn triggers a global update in our knowledge of the system, whereby the final reward is redistributed between all elementary decisions. In a way, backpropagation through the optimizer enables efficient credit assignment, even for sparse reward signals.

Another difference is related to the Bellman fixed point equation. In an RL setting, the Bellman equation is used explicitly to derive parameter updates. In our setting, the Bellman equation is used implicitly within optimizers such as Dijkstra's algorithm.

To conclude, while standard RL decomposes a policy into elementary decisions, the pipelines we study here are able to look directly for complex multistep solutions. Note that a similar concept of option exists in hierarchical RL [START_REF] Barto | Recent Advances in Hierarchical Reinforcement Learning[END_REF]: comparing both perspectives in detail would no doubt be fruitful, and we leave it for future work.

Our guiding example: shortest paths on Warcraft maps

As a way to clarify the concepts introduced in this paper, we illustrate them on the problem of Warcraft shortest paths [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF][START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF], which was already introduced on Figure 1. The associated dataset 8 , assembled by [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF], contains randomly-generated maps similar to those from the Warcraft II video game. Each of these maps is a red-green-blue (RGB) image of size ks × ks containing k × k square cells of side length s. Every cell has its own terrain type (grass, forest, water, earth, etc.) which incurs a specific cost when game characters cross it.

The goal is to find the shortest path from the top left corner of the map to the bottom right corner. At prediction time, cell costs are unknown, which means they must be approximated from the image alone. Solving the problem of Warcraft shortest paths thus requires adapting the pipeline of Equation (4) as follows:

Map image --------→ x∈R ks×ks×3 © Conv. Neural Network φ w Negative cell costs -----------→ θ∈R k×k Dijkstra's algorithm argmax v∈P k θ ⊤ v Shortest path --------→ y∈P k © Loss function L (7)
where we have defined the set of feasible paths

P k = {v ∈ {0, 1} k×k : v represents a path from (1, 1) to (k, k)}.
Training proceeds based on the map images and (possibly) the true shortest paths or cell costs provided in the dataset.

As suggested by [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF], we design the CNN based on the first few layers of a ResNet18 [START_REF] He | Deep Residual Learning for Image Recognition[END_REF]. We also append a negative softplus activation, in order to make sure that all outputs are negative. The sign constraint is there to ensure that Dijkstra's algorithm will terminate, but it is not obvious why we require negative costs instead of positive ones. The reason behind this sign switch is that Dijkstra's algorithm is a minimization oracle, whereas by convention InferOpt.jl works with maximization oracles.

We now show how to implement this pipeline in Julia. Throughout the paper, in addition to InferOpt.jl, the following packages are used: Flux.jl 9 [START_REF] Innes | Fashionable Modelling with Flux[END_REF][START_REF] Innes | Flux: Elegant Machine Learning with Julia[END_REF] Finally, Code sample 3 demonstrates the full prediction and optimization pipeline. It assumes that we have already parsed the Warcraft dataset into three vectors:

• images, whose elements are three-dimensional arrays representing map images;

• cells, whose elements are floating-point matrices representing true cell costs;

• paths, whose elements are binary matrices representing true shortest paths.

These functions do not rely on InferOpt.jl, but they will be used throughout the paper inside the differentiable wrappers provided by the package. Note that some snippets shown here have been shortened In this section, we focus on the CO oracle f defined in Equation (2), which is piecewise constant. By adopting a probabilistic point of view, we construct several smooth approximations f , which can be computed and differentiated based solely on calls to f .

The expectation of a differentiable probability distribution

As announced in Section 1.2.2, a probabilistic CO layer works in two steps. First, it constructs a probability distribution p( 2 illustrates this behavior on a two-dimensional polytope, with a maximization problem defined by the vector θ (black arrow). While the CO oracle outputs a single optimal vertex (red square), the probabilistic CO layer defines a distribution on all the vertices (light blue circles). Its output (dark blue hexagon) is a convex combination of the vertices with nonzero weights, which belongs to the convex hull of V (gray surface).

•|θ) ∈ ∆ V . Second, it returns the expectation f (θ) = E p(•|θ) [V ] ∈ conv(V). Figure
For such a layer to be useful in our setting, we impose several conditions. First, all computations must only require calls to the CO oracle f . Second, the expectation E p(•|θ) [V ] must be tractable (whether it is with an explicit formula, Monte-Carlo sampling, variational inference, etc.). Third, the mapping θ -→ p(•|θ) must be differentiable. If the last condition is satisfied, then the Jacobian of f is easily deduced from Equation (3):

J θ f (θ) = J θ E p(•|θ) [V ] = v∈V v∇ θ p(v|θ) ⊤ (8)
Let us give an example where analytic formulas exist. If V = {e 1 , ..., e d } is the set of basis vectors, then its convex hull conv(V) = ∆ d is the unit simplex of dimension d. Given an objective direction θ, solving argmax v∈V θ ⊤ v yields the basis vector f (θ) = e i where i is the index maximizing θ i .

Figure 2: Effect of a probabilistic CO layer

We want a probability distribution that evolves smoothly with θ, so we need to spread out the naive Dirac mass δ f (θ) (•) by putting weight on several vertices instead of just one. Let us assign to each vertex a probability that depends on its level of optimality in the optimization problem argmax v∈V θ ⊤ v, that is, on its inner product with θ. The Boltzmann distribution is a natural candidate, leading to p(e i |θ) ∝ e θ ⊤ ei = e θi . Computing the expectation reveals a well-known operation:

f (θ) = E p(•|θ) [V ] = d i=1 e θi d j=1 e θj e i = softmax(θ)
Unlike the "hardmax" function f , the softmax function f is differentiable, which justifies its frequent use as an activation function in classification tasks.

While the Boltzmann distribution can be used for specific sets V, in the general case, it has an intractable normalizing constant. This means we would need Markov chain Monte-Carlo (MCMC) methods to compute derivatives, which undermines the simplicity we are looking for. Fortunately, Sections 3.3.1 and 3.3.2 present other probability distributions which can be easily approximated through sampling.

Regularization as another way to define a distribution

Although this probabilistic point of view was recently put forward by [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF], the most popular paradigm in the literature remains regularization [START_REF] Blondel | Learning with Fenchel-Young Losses[END_REF]. Instead of using the CO oracle (2), regularization solves a different problem:

f Ω : θ -→ argmax µ∈dom(Ω) θ ⊤ µ -Ω(µ) (9)
where Ω : R d → R is a smooth and convex function that penalizes the output µ. Usually, Ω is chosen to enforce Ω(µ) = +∞ whenever µ / ∈ conv(V), which means dom(Ω) ⊆ conv(V). Remember that since V is finite, conv(V) is a polytope whose vertices form a subset of V.

The change of notation from v (vertex) to µ (moment) stresses the fact that v is an element of V, while µ belongs to dom(Ω) ⊆ conv(V). By Equation ( 6), any feasible µ is the expectation of some distribution over V, hence our choice of letter. This means we can write f Ω (θ) as a convex combination of the elements of V, whose weights are then interpreted as probabilities p Ω (•|θ). In other words, the two perspectives are not opposed: regularization is just another way to define a probability distribution.

For instance, if we go back to the concrete example from Section 3.1 and select Ω(µ) = i µ i log µ i (the negative Shannon entropy), we find once again that f Ω (θ) = softmax(θ). But the case of the unit simplex is very peculiar, because for any µ ∈ conv(V), the convex decomposition of µ onto the vertices V is unique. In other words, the correspondence between regularizations Ω and probability mappings p Ω is one-to-one.

This does not hold for arbitrary polytopes. As a result, we need to be more specific in how we choose the convex decomposition. In particular, we need the weights to be differentiable, in order to compute the Jacobian with Equation ( 8). Section 3.3.3 describes one possible approach, which relies on the Frank-Wolfe algorithm and implicit differentiation.

Conversely, the probability distributions given in Sections 3.3.1 and 3.3.2 also give rise to an implicit regularization, which can be expressed using Fenchel conjugates. This point of view is especially useful when we want to combine these layers with Fenchel-Young losses (Section 5.2.3). In essence, we claim that probabilistic CO layers and regularization are two sides of the same coin.

Collection of probabilistic CO layers

Our package implements various flavors of probabilistic CO layers, which are summed up in Table 1. Code sample 4 displays the operations they all support.

Remark 3.1. We also implement an Interpolation layer, corresponding to the piecewise linear interpolation of [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF]. However, to the best of our knowledge it cannot be cast as a probabilistic CO layer, so it does not support as many operations, and we only mention it for benchmarking purposes.

We now present each row of Table 1 in more detail. The main goal of Sections 3.3.1, 3.3.2 and 3.3.3 is to explain how p(•|θ) and f (θ) are computed, as well as their derivatives. We also draw connections with the regularization paradigm, which will ease the introduction of Fenchel-Young losses in Section 5.2.3. A hasty reader can safely skip to Section 4.

Layer Notations Probability p(•|θ) Regularization

PerturbedAdditive

p + ε , f + ε Explicit: f (θ + εZ) Implicit: Fenchel conjugate PerturbedMultiplicative p ⊙ ε , f ⊙ ε Explicit: f (θ ⊙ e εZ-ε 2 1/2 ) Implicit: Fenchel conjugate RegularizedGeneric p FW Ω , f FW Ω Implicit: Frank-Wolfe weights Explicit: function Ω
Table 1: Probabilistic CO layers and their defining features

Additive perturbation

A natural way to define a distribution on V is to solve (2) with a stochastic perturbation of the objective direction θ. [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF] suggest the following additive perturbation mechanism: ¦ ¥ Code sample 4: Supported operations for a probabilistic CO layer smoothness (large ε) and accuracy of the approximation (small ε). The associated probability distribution on V can be described explicitly:

f + ε (θ) = E argmax v∈V (θ + εZ) ⊤ v = E [f (θ + εZ)] (10
f + ε (θ) = v∈V v p + ε (v|θ) with p + ε (v|θ) = P (f (θ + εZ) = v) . (11)
Meanwhile, Proposition 3.1 allows us to compute differentials. Although the expectations cannot be expressed in closed form, they can be estimated using M Monte-Carlo samples Z 1 , ..., Z M ∼ N (0, I). Increasing M yields smoother approximations but makes the complexity grow linearly.

Proposition 3.1 (Differentiating through an additive perturbation [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF]). We have:

∇ θ p + ε (v|θ) = 1 ε E [1{f (θ + εZ) = v}Z] J θ f + ε (θ) = 1 ε E f (θ + εZ)Z ⊤ Proof. See Appendix A.1.1.
In order to recover the regularization associated with p + ε , we leverage convex conjugation. Let F + ε be the function defined by

F + ε (θ) = E max v∈V (θ + εZ) ⊤ v and let Ω + ε = (F + ε ) * denote its Fenchel conjugate.
Proposition 3.2 (Regularization associated with an additive perturbation [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF]). The function

Ω + ε is convex, it satisfies dom(Ω + ε ) ⊂ conv(V) and f + ε (θ) = argmax µ∈conv(V) θ ⊤ µ -Ω + ε (µ) = f Ω + ε (θ). Proof. See Appendix A.1.2.
Code sample 5 shows how this translates into InferOpt.jl syntax.

Multiplicative perturbation

Since the Gaussian distribution puts mass on all of R d , it can happen that some components of θ + εZ switch their sign with respect to θ. This may cause problems whenever the CO oracle for f has sign-dependent behavior. For instance, Dijkstra's algorithm for shortest paths requires all the edges of a graph to have a positive cost. In those cases, we need a sign-preserving kind of perturbation. Changing the distribution of Z to make it positive almost surely is not the right answer because it would bias the pipeline, leading to E[θ + εZ] > θ for the componentwise order. So instead of being additive, the perturbation becomes multiplicative:

f ⊙ ε (θ) = E argmax v∈V θ ⊙ e εZ-ε 2 1/2 ⊤ v = E f θ ⊙ e εZ-ε 2 1/2 (12)
Here, ⊙ denotes the Hadamard product, and the exponential is taken componentwise. Since E[e εZ ] = e ε 2 1/2 ̸ = 1, we add a correction term in the exponent to remove any bias: E[θ ⊙ e εZ-ε 2 1/2 ] = θ. As before, the associated probability distribution is easy to describe:

f ⊙ ε (θ) = v∈V v p ⊙ ε (v|θ) with p ⊙ ε (v|θ) = P f θ ⊙ e εZ-ε 2 1/2 = v . (13)
And Proposition 3.3 provides differentiation formulas that are very similar to the additive case.

Proposition 3.3 (Differentiating through a multiplicative perturbation). We have:

∇ θ p ⊙ ε (v|θ) = 1 εθ ⊙ E 1 f θ ⊙ e εZ-ε 2 1/2 = v Z J θ f ⊙ ε (θ) = 1 εθ ⊙ E f θ ⊙ e εZ-ε 2 1/2 Z ⊤ Proof. See Appendix A.2.1.
As far as regularization is concerned, we need a slight tweak compared to the additive case. Let F ⊙ ε be the function defined by

F ⊙ ε (θ) = E max v∈V θ ⊙ e εZ-ε 2 1/2 ⊤ v and let Ω ⊙ ε = (F ⊙ ε ) * denote its Fenchel conjugate. We define f ⊙scaled ε (θ) = E e εZ-ε 2 1/2 ⊙ f θ ⊙ e εZ-ε 2 1/2
Proposition 3.4 (Regularization associated with a multiplicative perturbation). The function

Ω ⊙ ε is convex and satisfies f ⊙scaled ε (θ) = argmax µ∈dom(Ω ⊙ ε ) θ ⊤ µ -Ω ⊙ ε (µ) = f Ω ⊙ ε (θ).
Unlike in the additive case, it is not f ⊙ ε itself that can be viewed as the product of regularization with

Ω ⊙ ε , but f ⊙scaled ε . Furthermore, this time we have dom(Ω ⊙ ε ) ̸ ⊆ conv(V).
Proof. See Appendix A.2.2.

Code sample 5 shows how this translates into InferOpt.jl syntax.

Generic regularization

We now switch our focus to the case of an explicit regularization Ω. Provided the regularization is convex and smooth, approximate computation of f Ω (θ) is made possible by the Frank-Wolfe algorithm [START_REF] Frank | An Algorithm for Quadratic Programming[END_REF]). This algorithm is interesting for two reasons. First, it only requires access to the CO oracle f and the gradient of Ω. Second, its output is expressed as a convex combination of only a few polytope vertices [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-Free Sparse Convex Optimization[END_REF]. In other words, the Frank-Wolfe algorithm does not just return a single point f Ω (θ) ∈ conv(V): it also defines a sparse probability distribution p FW Ω (•|θ) over the vertices V such that This distribution is called sparse because most of the weights are actually zero. Note that p FW Ω (•|θ) is not uniquely specified by the regularization Ω, but instead depends on the precise implementation of the Frank-Wolfe algorithm (initialization, step size, convergence criterion, etc.). In particular, the number of atoms in the distribution is upper-bounded by the number of Frank-Wolfe iterations.

f Ω (θ) = v∈V v p FW Ω (
As pointed out by Blondel, Berthet, et al. (2022, Appendix C), there exists a function g(p, θ) defined on ∆ V ×R d such that p FW Ω (•|θ) is a fixed point of its projected gradient operator p -→ proj ∆ V (p-∇ p g(p, θ)). Since the orthogonal projection onto the simplex ∆ V is itself differentiable [START_REF] Martins | From Softmax to Sparsemax: A Sparse Model of Attention and Multi-Label Classification[END_REF], we can apply the implicit function theorem to this fixed point equation. Doing so yields gradients ∇ θ p Ω (v|θ) that we use to compute a Jacobian for f Ω (θ). Again, by sparsity, this sum only has a few non-zero terms, which makes it tractable:

J θ f Ω (θ) = v∈V v∇ θ p FW Ω (v|θ) ⊤ . (14)
Among all the possible functions Ω, the quadratic penalty Ω(µ) = 1 2 ∥µ∥ 2 is particularly interesting. It gives rise to the SparseMAP method [START_REF] Niculae | SparseMAP: Differentiable Sparse Structured Inference[END_REF], whose name comes from the sparsity of the Euclidean projection onto a polytope:

f Ω (θ) = argmax µ∈conv(V) θ ⊤ µ - 1 2 ∥µ∥ 2 = argmin µ∈conv(V) ∥µ -θ∥ 2 .
This is the one we used for the example of Code sample 6. Our implementation relies on the recent package FrankWolfe.jl14 [START_REF] Besançon | FrankWolfe.jl: A High-Performance and Flexible Toolbox for Frank-Wolfe Algorithms and Conditional Gradients[END_REF].

Remark 3.2. Blondel, Martins, and Niculae (2020) also suggest distribution regularization, whereby Ω(µ) is defined through a generalized entropy H(p) on ∆ V :

Ω(µ) = -max p∈∆ V H(p) s.t. E p [V ] = µ.
Distribution regularization can only be computed explicitly for certain entropies H (Shannon entropy, Gini index) and certain polytopes conv(V) (unit simplex, permutahedron, spanning trees, etc.). In each case, a custom combinatorial algorithm is required. Since we aim for a generic approach, we only consider mean regularization, which is defined directly on the expectation µ.

The case of inexact CO oracles

In our discussion so far, an implicit assumption was that the CO oracle f returns an exact solution to Equation (2). For most polynomial problems (as well as some NP-hard problems which are tractable in practice), this is perfectly reasonable. But in some cases, exact solutions are too expensive to compute. Then, our CO oracle may only be able to return an inexact solution, for instance because branch & bound has to be interrupted before the whole tree can be explored. What kind of impact does this have on the precision of the computed Jacobian? Let us denote by f a hypothetical exact oracle, and by g an inexact oracle.

Proposition 3.5 (Jacobian precision for inexact oracles -perturbed case). Suppose we use g instead of f with additive (resp. multiplicative) perturbation. Then the error on the Jacobian of the probabilistic CO layer satisfies:

J θ g + ε (θ) -J θ f + ε (θ) 2 ≤ √ d ε ∥g -f ∥ ∞ J θ g ⊙ ε (θ) -J θ f ⊙ ε (θ) 2 ≤ √ d ε min i |θ i | ∥g -f ∥ ∞ .
While requiring the inexact oracle g to be uniformly close to f is quite restrictive, this result does provide heuristic justification for the use of inexact oracles in practice.

Proof. See Appendix A.3.

Learning by experience

Now that we have seen several ways to construct probabilistic CO layers, we turn to the definition of an appropriate loss function. Let us start with learning by experience, which takes place when we only have access to input samples without target outputs. In that case, Equation (5) simplifies as

min w 1 N N i=1 L f φ w (x (i) ) . ( 15
)
As we will see below, the regret, which is the natural choice of loss, does not yield interesting gradients.

That is why we propose a family of smooth regret surrogates derived from our probabilistic CO layers, and explain how to differentiate them. While similar losses have been hinted at in previous works, to the best of our knowledge, our general point of view is new.

To make notations lighter, we restrict ourselves to a single input x. Furthermore, we write losses as functions of θ instead of w. Indeed, our losses do not just rely on y = f (θ): they use f as an ingredient internally. In practice, we leave it to AD to exploit the relation θ = φ w (x) in order to compute gradients with respect to w.

Minimizing a smooth regret surrogate

When we learn by experience, the problem statement usually includes a cost function c : V → R, and we want our pipeline to generate solutions that are as cheap as possible. Internally, this cost function may use parameters that are unknown to us at prediction time: typically, it may assess the quality of our solution using the true objective direction θ. It may be useful to think about c as the feedback provided by an outside evaluator, rather than a function we implement ourselves.

The natural loss to minimize is the cost incurred by our prediction pipeline, also called regret: This function relies on the CO oracle f , which is piecewise constant. Our spontaneous impulse would be to replace the CO oracle f with a probabilistic CO layer f , thus minimizing c( f (θ)). Unfortunately, the cost function c is not necessarily smooth either. To make matters worse, c may only be defined on vertices v ∈ V, and not on general convex combinations µ ∈ conv(V).

R(θ) = c(f (θ)). (
The solution we propose relies on the pushforward measure (also called image measure) of p(•|θ) with respect to the function c. Recall that a probabilistic CO layer is defined by

f (θ) = E p(•|θ) [V ].
To compose it with an arbitrary cost, instead of applying c outside the expectation, we apply it inside the expectation. In other words, we first push the measure p(•|θ) forward through the function c, before taking the expectation. This gives rise to the notion of expected regret:

R p (θ) = E p(•|θ) [c(V )] (17)
By integration, this loss is just as smooth as the probability mapping θ -→ p(•|θ), which means we can compute its gradient easily. We therefore suggest using the expected regret R p that stems from the probabilistic CO layers p + ε , p ⊙ ε , and p FW Ω defined in Section 3. Note that if c is linear and defined on all of conv(V), then Since the learning problem is non-convex, we may also try to minimize the (non-smooth) regret R using global optimization algorithms such as DIRECT [START_REF] Jones | Lipschitzian Optimization without the Lipschitz Constant[END_REF]. Perhaps surprisingly, this has been shown to yield good results when φ w is a generalized linear model and the dimension of the weights w is not too large, i.e., no greater that 100 (Parmentier 2021a). When the ML layer φ w is a large neural network, we cannot use this approach anymore.

E p(•|θ) [c(V )] = c E p(•|θ) [V ] and the two quantities coincide. Furthermore, if c is convex, then E p(•|θ) [c(V )] ≥ c E p(•|θ) [V ]

Derivatives of the regret for learning by experience

When p comes from a random perturbation, we can formulate the gradient of the expected regret as an expectation too, and approximate it with Monte-Carlo samples.

Proposition 4.1 (Gradient of the expected regret, perturbation setting). We have:

∇ θ R p + ε (θ) = 1 ε E [(c • f )(θ + εZ)Z] ∇ θ R p ⊙ ε (θ) = 1 εθ ⊙ E (c • f ) θ ⊙ e εZ-ε 2 1/2 Z .
Proof. For the additive perturbation, it is a consequence of Proposition 3.1. For the multiplicative perturbation, it is a consequence of Proposition 3.3.

These gradients obey a simple logic: the more the perturbation Z increases the cost of a solution, the more positive weight it gets, and vice versa. To see it, we remember that E[Z] = 0 for a standard Gaussian, and rewrite the regret gradients as follows:

∇ θ R p + ε (θ) = 1 ε E [(c • f )(θ + εZ)Z -(c • f )(θ)Z] ∇ θ R p ⊙ ε (θ) = 1 εθ ⊙ E (c • f ) θ ⊙ e εZ-ε 2 1/2 Z -(c • f )(θ)Z .
On the other hand, when p is derived from an explicit regularization Ω, the expected regret is amenable to implicit differentiation of the Frank-Wolfe algorithm. Once more, the sparsity property makes exact computation tractable by reducing the number of terms in the sum:

∇ θ R p FW Ω (θ) = v∈V c(v)∇ θ p FW Ω (v|θ)
Remark 4.2. Although the previous discussion focuses on a scalar-valued cost, it actually applies to any pushforward function c, even with vector values. The formulas for the generic Jacobian are given below:

J θ E p + ε (•|θ) [c(V )] = 1 ε E (c • f )(θ + εZ)Z ⊤ J θ E p ⊙ ε (•|θ) [c(V )] = 1 εθ ⊙ E (c • f ) θ ⊙ e εZ-ε 2 1/2 Z ⊤ J θ E p FW Ω (•|θ) [c(V )] = v∈V c(v)∇ θ p FW Ω (v|θ) ⊤
At the moment, InferOpt.jl only handles the case where c is a fully-defined function without free parameters. In the near future, we will add support for the case where c is itself an ML layer with learnable weights.

Learning by imitation

We now move on to learning by imitation, where additional information is used to guide the training procedure. For each input sample x (i) , we assume we have access to a target t(i) . In that case, Equation ( 5) simplifies as

min w 1 N N i=1 L f φ w (x (i) ) , t(i) , (18)
and we can see that the loss takes the target as an additional argument. In this section, we introduce imitation losses that are well-suited to hybrid ML-CO pipelines, and explain how to compute their gradients. As in Section 4, we only consider a single input x, and we write losses as L(θ, t).

A loss that takes the optimization layer into account

There are two main kinds of target. The first one is a good quality solution t = ȳ. The second one is the true objective direction θ, from which we can also deduce ȳ = f ( θ), so that t = ( θ, ȳ). When learning by imitation, it is tempting to focus only on reproducing the targets, but this would be misguided. To explain why, we revisit the pipeline of Equation ( 4).

Remember that we may have access to the true objective direction θ during training, but at prediction time, the CO oracle f is applied to the encoder output θ = φ w (x) instead. Minimizing a naive square loss like ∥φ w (x) -θ∥ 2 completely neglects the asymmetric impacts of the prediction errors on θ: for example, overestimating or underestimating θ may have very different consequences on the quality of the downstream solution. That is why, according to [START_REF] Elmachtoub | Smart "Predict, Then Optimize[END_REF], we need a loss function that takes the optimization step into account. The same holds true when we have access to a precomputed solution ȳ. [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF] present experiments showing that the naive square loss ∥ f (φ w (x))-ȳ∥ 2 performs poorly compared with more refined approaches. Our own numerical findings (Section 6) support their conclusion.

To sum up, we want a loss that does not neglect the optimization step. Let y temporarily denote the output of our pipeline. When surveying the literature, we realized that most flavors of imitation learning use losses that combine the same components:

L aux ℓ,Ω (θ, t, y) = ℓ(y, t) base loss + θ ⊤ (y -ȳ)
gap between y and ȳ for the CO problem ( 2)

-(Ω(y) -Ω(ȳ)) regularization term (19)
Here is another way to write it:

L aux ℓ,Ω (θ, t, y) = ℓ(y, t) base loss + θ ⊤ y -Ω(y) -θ ⊤ ȳ -Ω(ȳ)
gap between y and ȳ for the regularized CO problem (9)

The base loss ℓ(y, t) is similar in spirit to the cost function c(y) from Section 4. But it is the gap term that truly makes it possible for the optimization problem to play a role in the loss. Indeed, minimizing the gap encourages the (regularized) CO problem to output a solution y that is close to the target ȳ.

Putting these components together yields a linear function of θ, and we can remove the dependency in y by maximizing over y:

L gen ℓ,Ω (θ, t) = max y∈dom(Ω) L aux ℓ,Ω (θ, t, y) = max y∈dom(Ω) ℓ(y, t) + θ ⊤ (y -ȳ) -(Ω(y) -Ω(ȳ)) . ( 20
)
The following result justifies why this is an interesting loss.

Proposition 5.1 (Properties of the generic loss for learning by imitation). The function L gen ℓ,Ω (θ, t) is convex with respect to θ, and a subgradient is given by

argmax y∈dom(Ω) L aux ℓ,Ω (θ, t, y) -ȳ ∈ ∂ θ L gen ℓ,Ω (θ, t). (21)
Proof. As a pointwise maximum of affine functions, θ -→ L gen ℓ,Ω (θ, t) is convex. Its subgradient is obtained using Danskin's theorem [START_REF] Danskin | The Theory of Max-Min and Its Application to Weapons Allocation Problems[END_REF].

Method Notation Target Base loss Regul. The idea is that solving argmax y∈dom(Ω) L aux ℓ,Ω (θ, t, y) should not be much harder than the regularized CO problem (9). Therefore, using such a loss function dispenses us from differentiating through the probabilistic CO layer: most of the time, we only need to compute the layer output in order to obtain a loss subgradient for free.

Loss formula

S-SVM

L S-SVM ℓ ȳ ℓ(y, ȳ) No max y ℓ(y, ȳ) + θ ⊤ (y -ȳ) SPO+ L SPO+ ( θ, ȳ) θ⊤ (ȳ -y) No max y θ⊤ (ȳ -y) + 2θ ⊤ (y -ȳ) FY L FY Ω ȳ 0 Yes max y θ ⊤ (y -ȳ) -(Ω(y) -Ω(ȳ)) Generic L gen ℓ,Ω t ℓ(y, t) Yes max y ℓ(y, t) + θ ⊤ (y -ȳ) -(Ω(y) -Ω(ȳ))

Collection of losses for learning by imitation

Several prominent loss functions from the literature are special cases of our decomposition (20): we gather them in Table 2. Code sample 10 clarifies their construction, while Code sample 11 displays supported operations. Finally, the entire program necessary for learning by imitation is shown on Code sample 12.

In Sections 5.2.1, 5.2.2, 5.2.3 and 5.2.4, we go over these special cases to explain how to compute each loss and its subgradient using Equation (21). They can be skipped without danger.

Remark 5.1. While the S-SVM and SPO+ losses do not fall within the framework of probabilistic CO layers (due to the absence of regularization), we still include them for benchmarking purposes.

Structured support vector machines

The structured support vector machine (S-SVM) was among the first methods introduced for learning in structured spaces (Nowozin and Lampert 2010, Chapter 6). Given a target solution ȳ and an underlying distance function ℓ(y, ȳ) on V, the S-SVM loss is computed as follows:

L S-SVM ℓ (θ, ȳ) = max y∈V {ℓ(y, ȳ) + θ ⊤ (y -ȳ)}. ( 22
)
The subgradient formula (21) becomes

argmax y∈V {ℓ(y, ȳ) + θ ⊤ (y -ȳ)} -ȳ ∈ ∂ θ L S-SVM ℓ .
Note that due to the presence of ℓ, computing a subgradient requires an auxiliary solver that is different from the linear oracle f . This is why we do not illustrate the S-SVM with a code sample. In InferOpt.jl, we only implement this auxiliary solver for the unit simplex, in the case where ℓ is the Hamming distance. However, we also provide a generic layer where the user can plug in the relevant auxiliary solver.

Smart "predict, then optimize"

The smart "predict, then optimize" (SPO) paradigm is applicable when the true objective direction θ is known (remember that in this case, we have ȳ = f ( θ)). Elmachtoub and Grigas (2022) define the SPO+ loss function as follows:

L SPO+ (θ, θ) = (2θ -θ) ⊤ f (2θ -θ) + ( θ -2θ) ⊤ ȳ (23) = max y∈V θ⊤ (ȳ -y) + 2θ ⊤ (y -ȳ) .
It can be seen as a special case of S-SVM. But this time, computing the loss and its subgradient only requires calling f twice: 2f (2θ -θ) -2ȳ ∈ ∂ θ L SPO+ (θ, θ).

Fenchel-Young losses

The framework of Fenchel-Young losses is built on the theory of convex conjugates, in particular the Fenchel-Young inequality [START_REF] Blondel | Learning with Fenchel-Young Losses[END_REF]. Starting from a target solution ȳ and a regularization Ω, a loss is constructed as follows:

L FY Ω (θ, ȳ) = Ω * (θ) + Ω(ȳ) -θ ⊤ ȳ (24) = max y∈conv(V) θ ⊤ y -Ω(y) -θ ⊤ ȳ -Ω(ȳ)
This time, the loss and subgradient require access to f Ω :

f Ω (θ) -ȳ ∈ ∂ θ L FY Ω (θ, ȳ).
As can be inferred from the expression above, there are deep connections between Fenchel-Young losses and the regularization paradigm of Section 3.2. In particular, it is also possible to use implicit regularization by perturbation [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF]. The fact that we cannot compute Ω + ε (y) or Ω ⊙ ε (y) is not a real obstacle: since those terms do not depend on θ, we can just drop them from the loss during training. We end up with the following estimators for the loss and its subgradient:

L FY Ω + ε (θ, ȳ) = F + ε (θ) -θ ⊤ ȳ f + ε (θ) -ȳ ∈ ∂ θ L FY Ω + ε (θ, ȳ) L FY Ω ⊙ ε (θ, ȳ) = F ⊙ ε (θ) -θ ⊤ ȳ f ⊙scaled ε (θ) -ȳ ∈ ∂ θ L FY Ω ⊙ ε (θ, ȳ)

Generic imitation loss

Of course, it is tempting to fill in the blanks of Table 2 by combining every single term of the loss decomposition (20). To the best of our knowledge, this has not been done before in the literature, but there is no theoretical obstacle. If we use this generic loss together with regularization, then it is interesting to remark that L gen ℓ,Ω (θ, t) acts as a convex upper bound on the base loss ℓ( f Ω (θ), t). Indeed, since ȳ is a worse solution than f Ω (θ) for ( 9), we have

ℓ( f Ω (θ), t) ≤ ℓ( f Ω (θ), t) + θ ⊤ f Ω (θ) -Ω( f Ω (θ)) -θ ⊤ ȳ -Ω(ȳ) ≤ max y∈conv(V) ℓ(y, t) + θ ⊤ (y -ȳ) -(Ω(y) -Ω(ȳ)) = L gen ℓ,Ω (θ, t).
Therefore, our generic loss can be seen as a crossover between the Fenchel-Young loss and a problem-specific base loss. It is not yet implemented in InferOpt.jl, and we leave its thorough testing for future work.

Applications

In this section, we define specific ML-CO pipelines, in both learning settings (by experience and by imitation), to address four problems: our shortest path problem on Warcraft maps, a stochastic vehicle scheduling problem, a single-machine scheduling problem, and a two-stage stochastic minimum weight spanning tree problem. The first one stems from the ML community. The other three are classics from the field of operations research. They illustrate the idea of approximating hard optimization problems with simpler ones using ML-CO pipelines.

Shortest paths on Warcraft maps

First we come back to our guiding example of Section 2.3. Our aim is to illustrate the various learning settings introduced in this paper, and to evaluate their relative performance. We do so with two kinds of shortest path (SP) oracles. The first one uses Dijkstra's algorithm. The second one uses the Ford-Bellman algorithm with a bounded number of iterations.

Experimental setting

In every experiment presented here, we only consider a sub-dataset, made up of 1% of the original Warcraft dataset from [START_REF] Vlastelica | Differentiation of Blackbox Combinatorial Solvers[END_REF]. It contains 200 samples, or maps, which we split into 80 training samples, 100 validation samples (for hyperparameter tuning) and 20 test samples (for performance evaluation). For each learning setting, the train, test and validation sets remain the same, and we individually tune a subset of the hyperparameters stated in Table 3. Our motivation for reducing the dataset is to show that we can still obtain convincing results with a limited amount of computation. We use the Metalhead.jl package to build a truncated ResNet18 CNN, Flux.jl to train our pipelines with the Adam optimizer [START_REF] Kingma | Adam: A Method for Stochastic Optimization[END_REF], and GridGraphs.jl to compute shortest paths. Our code is available in the WarcraftShortestPaths.jl15 repository. The experiments are conducted on a MacBook Pro with 2,3 GHz Intel Core i9, 8 cores and 16 Go 2667 MHz DDR4 RAM.

To obtain a precise learning setting, we need to define:

1. The combinatorial problem we need to solve, along with an appropriate oracle.

2. The probabilistic CO layer used to wrap said oracle.

3. The data we have at our disposal to train our pipeline. 4. The loss function we want to minimize.

Hyperparameter

No matter the setting, the data always contains a list of RGB map images. When we learn by experience, we also have access to a black box cost function, which evaluates paths based on the true cell costs (see the beginning of Section 4). On the other hand, when we learn by imitation, we add targets to the maps (as defined in Section 5). The target in our case always includes the optimal path, with or without the true cell costs.

All those ingredients are detailed in Table 4 for each learning setting we consider. The names given in the first column are reused in the legends of Figure 3. Most of the probabilistic CO layers considered in this paper do not prevent the objective vector θ from changing its sign, and the same goes for the losses As a result, we need oracles able to accommodate negative cell costs. That is why we use the Ford-Bellman algorithm, while limiting the number of iterations to the number of nodes in the grid graph (to ensure termination even with negative cycles). Our multiplicative perturbation is the only approach that preserves non-negative costs. It enables us to apply Dijkstra's algorithm, which has a smaller time and space complexity.

Results

In Figure 3, we show the average train (Figure 3a) and test (Figure 3b) optimality gaps, computed using the true cell costs. We compare all the settings detailed in Table 4, with the exception of MSE base loss + additive noise (we could not get satisfactory results using only 80 training samples). To quantify training effort, instead of counting epochs (i.e. passes through the dataset), we use the number of optimizer calls, because these calls are the truly time-consuming part. This aims at comparing learning settings which involve different amounts of computation per gradient step. For instance, using SPO+, we need 2 optimizer calls to compute the loss gradient for one sample. On the other hand, we need M optimizer calls if we choose Fenchel-Young perturbed additive noise.

When learning by imitation, SPO+ reaches almost zero average gap both on train and test sets after very few optimizer calls, even though we only kept 1% of the initial dataset. This impressive result can be understood since, with SPO+, we have access to the true cell costs during training, and we leverage the problem structure within the loss.

Assuming we only have access to target paths, we obtain better results with Fenchel-Young losses than with MSE losses. The train and test average gaps are lower than 5% with the former, and we observe good generalization performance. This may be explained by the use of the optimization problem in the Fenchel-Young loss definition. On the contrary, in the MSE setting, although we have access to target paths, we only seek to imitate them without truly accounting for solution cost.

Perhaps surprisingly, we also manage to learn by experience with our small sub-dataset. Indeed, using the techniques introduced in Section 4, we reach 7% average test gaps in the cost perturbed multiplicative noise setting, which is better than learning by imitation with an MSE loss. To the best of our knowledge, it is the first time that learning by experience (as defined in Section 4) is combined with CNNs.

Approximating hard optimization problems

Solving hard optimization problems is one of the main applications of hybrid CO-ML pipelines. The principle is to build a pipeline which approximates a hard combinatorial problem, for which we do not have an efficient algorithm. The approximation relies on a similar but easier problem, for which an efficient algorithm exists. This "easy" algorithm is then used as a CO layer. In order to obtain the best possible approximation (in terms of closeness to the original problem), this CO layer is complemented with encoding and decoding layers, whose weights can be learned.

In the remaining of this section, we use InferOpt.jl in this framework, and apply it to three different problems. In Section 6.3, we focus on the stochastic vehicle scheduling problem with similar experiments as in Parmentier (2021a) for the learning by imitation, and Parmentier (2021b) for the learning by experience. Then, in Section 6.4, we study the single machine scheduling problem and compare our results to those in [START_REF] Parmentier | Learning to Solve the Single Machine Scheduling Problem with Release Times and Sum of Completion Times[END_REF]. Finally, in Section 6.5, we look at the two-stage stochastic minimum weight spanning tree problem, and demonstrate the use of a GNN in combination with InferOpt.jl. Source code for these three applications can be found in their respective satellite packages.

Stochastic vehicle scheduling problem

We use InferOpt.jl to solve the stochastic vehicle scheduling problem (StoVSP), by learning a transformation that approximates its instances as instances of the easier to solve vehicle scheduling problem (VSP). Source code can be found in the satellite package StochasticVehicleScheduling.jl16 .

Problem formulation

Vehicle scheduling problem The (deterministic) VSP consists in assigning vehicles to cover a set of scheduled tasks, while minimizing the total cost. Let V be the set of tasks. Each task v ∈ V has a scheduled beginning time t b v , and a scheduled end time t e v , such that t e v > t b v . We denote t tr (u,v) the travel time from task u to task v. A task v can be scheduled after another task u only if we can reach it in time, before it starts:

t b v ≥ t e u + t tr (u,v) (25)
An instance of VSP can be modeled with a directed graph D = (V, A), with V = V ∪ {o, d}, and o, d origin and destination dummy nodes. For all task v ∈ V , (o, v) and (v, d) are arcs in A. Additionally, there is an arc between two tasks u and v only if ( 25) is satisfied. The resulting graph D is acyclic.

A solution of the VSP problem is a set of od paths partitioning D, such that all tasks are covered exactly once. The objective is to minimize the sum of path edge costs. This can be formulated as an LP (see Appendix B.1.1), and can be solved either using a flow algorithm, or using a general purpose LP solver.

Stochastic Vehicle Scheduling

In the stochastic VSP, we consider the same setting as the deterministic version, to which we add the following. Once the scheduling decision is set, we observe random delays, which propagate along vehicle paths. The objective is to minimize the sum of vehicle costs and expected total delay costs.

We consider a finite set of scenarios s ∈ S. For each task v ∈ V , we denote γ s v the intrinsic delay of v in scenario s, and d s v its total delay. We also denote δ s u,v the slack between tasks u and v. These quantities follow the delay propagation equation when u and v are consecutively operated by the same vehicle:

d s v = γ s v + max(d s u -δ s u,v , 0) (26)
This leads to a much more difficult problem to solve. In Appendix B.1.2, we provide a compact MILP formulation, which enables to easily solve optimally instances with up to 25 tasks using commercial MILP solvers.

Datasets

To generate our instance datasets, we use a generator similar to the one used by Parmentier (2021b). More details are given in Appendix B.1.3. We consider 3 training/validation datasets used to train our models. Each dataset contains 100 instances, and is divided into 50 training instances and 50 validation instances. The first dataset contains only instances with 25 tasks and 10 scenarios, the second one only instances with 50 tasks and 50 scenarios, and the last one only instances with 100 tasks and 50 scenarios. The 25 tasks dataset contains small instances for which we can easily compute optimal solutions to use as target solutions during the learning. For the two other datasets, we instead used a local search heuristic to compute "good" (non necessarily optimal) solutions of each instance. For testing purposes, we consider several additional test datasets we only use for evaluating the final performances and generalization abilities of our learned models. These test datasets contain larger instances with up to 1000 tasks. Datasets content is summarized in Table 5. Table 5: Summary of the number of instances in each used dataset

InferOpt pipeline

We use the deterministic VSP as a CO layer in the following pipeline to solve pipeline to solve the stochastic VSP. The goal is to train a model able to generalize well enough to instances it has not seen before, especially bigger instances, for which we cannot compute a good solution in reasonable time. For this, we consider the pipeline presented in Figure 4. a∈δ -(v)

© GLM φ w VSP flow Linear Program © Loss function StoVSP instance Edge weights θ a , ∀a ∈ A Vehicle routes o v 1 v 2 v 3 v 4 v 5 v 6 v 7 v 8 d time o v 1 v 2 v 3 v 4 v 5 v 6 v 7 v 8
y a = a∈δ + (v) y a , ∀v ∈ V (flow polytope) a∈δ -(v) y a = 1, ∀v ∈ V (task covering) y a ≥ 0, ∀a ∈ A
Training this pipeline means finding parameter vector w that minimizes the chosen loss.

Experimental setting

For each of the three datasets considered, we train two models: we learn the first one by imitation thanks to a Fenchel-Young loss, and the second one by experience thanks to the expected regret. In both cases, we use a PerturbedAdditive regularization with parameters ε and nb_samples. We use the Flux.jl Julia library and the Adam optimizer for training. Hyperparameters used for each model/dataset pair can be found in Table 6.

Hyperparameters

Learn by imitation models Learn by experience models 25 tasks 50 tasks 100 tasks 25 tasks 50 tasks 100 tasks Datasets features are normalized during training by dividing each feature by its standard deviation in the train dataset. At the end of training, the obtained w is renormalized accordingly. We optimize hyperparameters using metric values on validation datasets. The experiments are conducted on a 2.6 GHz Intel Core i9-11950H, 16 cores, with 64 Go RAM. We faced one practical difficulty during the training: the calibration of hyperparameter ε when learning by experience. Indeed, training performance is quite sensible to its value, which lead to a lot of hyperparameter tuning on the validation set. Detailed training plots can be found in Appendix B.1.4.

ε 0.1 0.1 0.

Results

Finally, we can look at the performance of the learned pipelines on all the different datasets, and compare them with each others. In Figure 5, we observe the resulting solution of three different algorithms on the same input instance. Interestingly, the model learned by imitation outputs a different solution than the model learned by experience. In Table 7, we compare the average execution time between the local search heuristic and one of our learned pipelines. We see that, once the model is learned offline, online prediction of a solution through the pipeline is very fast.

Test dataset

tasks

50 tasks 100 tasks heuristic learned pipeline heuristic learned pipeline heuristic learned pipeline 0.12s 0.004s 1.47s 0.01s 5.28s 0.043s

Table 7: Learning by imitation: comparison of average execution time on an instance between the local search heuristic (10000 iterations) a learned pipeline

Learning by imitation We train three models in the learning by imitation setting, one for each training dataset. Final results on test data are summarized in Tables 8 and9. First, we observe the average and maximum cost gaps on test datasets with instances of the same size as in training datasets, between predicted solutions and labeled solutions. We recall that labeled solutions are optimal for 25 tasks, and heuristic (obtained from local search) for 50 and 100 tasks. That is why gaps should stay positive when predicting on the 25 tasks test dataset, but can be negative when predicting on the two other datasets.

We can see that all three models obtain very good gap values, with less than 1% on 25 tasks instances on average, and average negative gaps on 50 and 100 tasks instances, which means better solutions on average than the local search heuristic used as labels for the learning.

Train dataset

Test dataset 25 tasks 50 tasks 100 tasks mean max mean max mean max 25 tasks 0.68% 9.46% -0.41% 4.26% -1.02% 2.4% 50 tasks 0.49% 3.01% -0.46% 2.34% -1.6% 0.62% 100 tasks 0.62% 3.36% -0.14% 9.9% -1.2% 0.11% Then, we measure the performance of our models on test datasets with larger instances than training ones. We do not necessarily have label solutions available for these big instances, but can evaluate the average cost per task. This metric is very useful to assess how well a model generalizes. Indeed, since all instances share the same map size regardless of the number of tasks, the average cost per task should decrease when the number of tasks increases, because optimal paths contain more tasks. It is the case for the model trained on 25 tasks instances, but not for the other two, which deteriorate on instances larger than 300 tasks. This is probably due to the fact that 50 and 100 tasks datasets are not labeled with optimal solutions, but only with solutions from the local search heuristic. Indeed, the 100-task model performs even worse than the 50-task model, because the local search heuristic performance decreases when the number of tasks increases.

Train dataset

Test dataset (number of tasks in each instance) 25 50 100 200 300 500 750 1000 25 tasks 274.72 225.29 207.14 194.46 186.68 182.56 178.57 177.3 50 tasks 274.27 225.23 205.97 195.78 193.12 194.48 196.99 199.38 100 tasks 274.61 225.87 206.8 197.97 195.53 207.02 219.34 227.14 Table 9: Learning by imitation: average cost per task

Learning by experience In the learning by experience setting, we also train three models on the same three training datasets, and summarize results in Tables 10 and11. First, we observe that average cost gaps and maximum cost gaps are lower on all models, respect to the corresponding models learned by imitation.

Train dataset

Test dataset 25 tasks 50 tasks 100 tasks mean max mean max mean max 25 tasks 0.45% 4.2% -0.77% 0.63% -2.11% -0.14% 50 tasks 0.43% 3.04% -0.78% 0.74% -2.06% -0.22% 100 tasks 0.43% 3.28% -0.83% 0.97% -2.06% -0.29% Finally, unlike the learning by imitation setting, we observe decreasing average cost per task on all three models for instances up to 1000 tasks. This means models learned by experience generalize much better, because they are not biased trying to imitate non-optimal solutions. In the stochastic VSP, learning by experience seems more performant than learning by imitation, the only downside being the practical difficulty of the training. Indeed, results are more sensible to hyperparameters and the model needs to be trained for more epochs because the loss convergence is slower.

Single-machine scheduling

We solve a single-machine scheduling problem with release times and sum of completion times. Given an instance of this NP-hard problem, we use an ML predictor to define the input of a ranking problem, which is an LP on the permutahedron and can be solved instantaneously. Its optimal solution (permutation) can then be decoded into an approximate solution of the hard problem. This follows the recent work of Parmentier and T'Kindt (2021), re-coded in Julia in the satellite package SingleMachineScheduling.jl17 .

Problem formulation

We need to schedule n jobs on a single machine. Each job j ∈ [n] has an associated processing time p j , and release date r j . Before r j , the job j cannot start. The machine can process only one job at a time. Preemption is not allowed, which means whenever a job has started, it must be completed, and cannot be paused to be finished later. Our aim is to find a schedule (permutation) s = (j 1 , ..., j n ) of [n], such that the sum of completion times j C j (s) is minimal. In this sum, C j (s) is the completion time of job j in the schedule s. More precisely, given the schedule s, we have C j1 (s) = r j1 + p j1 , and for k > 1, we have C j k (s) = max r j k , C j k-1 (s) + p j k .

Existing algorithms for the hard problem

Two kinds of algorithms are considered, both to benchmark our learning pipelines, and to provide solutions for the learning dataset. The first kind includes several heuristics: the RDI/APRTF heuristic [START_REF] Chand | An Iterative Heuristic for the Single Machine Dynamic Total Completion Time Scheduling Problem[END_REF], referred to as RDIA, the RBS heuristic with beam width w = 2 (Della Croce, Salassa, and T'kindt 2014), referred to as RBS, and the matheuristic (Della Croce, Salassa, and T'kindt 2014), referred to as MATH. The second kind is an exact algorithm, used to solve instances with up to n = 110 jobs, the branch-and-memorize algorithm (Shang, T'Kindt, and Della Croce 2021).

Datasets

We consider two types of datasets generated randomly: Table 12: Hyperparameters for scheduling experiments.

InferOpt pipeline

In the pipeline ( 27), an instance of the hard problem x is passed through an ML layer involving the GLM φ w , leading to the parameter θ = φ w (x). We then sort the components of θ using a ranking algorithm. It leads to a permutation s. This permutation is then passed through a decoder. The output of the decoder y is the candidate solution to the hard problem.

Hard problem instance x -----------→ § ¦ ¤ ¥ ML layer φ w Parameter θ = φ w (x) ---------→ Ranking f Permutation s = f (θ) ----------→ § ¦ ¤ ¥ Decoder ψ Solution y = ψ(s) -------→ (27)
Parmentier and T'Kindt (2021) consider two possibilities for the decoder ψ. First, we can use the permutation s = f (θ) as the candidate solution of the hard problem. We call this algorithm pure machine learning heuristic, referred to as PMLH. It corresponds to ψ = id. Second, in the improved machine learning heuristic, referred to as IMLH, we combine a local search and the RDI procedure in the decoder ψ. We consider PMLH in the following.

Experimental setting

A fast ML-CO pipeline. We focus on a learning by imitation setting, with ψ = id. Since our train dataset is made of small instances (with up to 110 jobs), and our test dataset has large instances (with up to 2550 jobs), we want to see if our model generalizes well. Since ranking is quasi-instantaneous, and the algorithms defined in Section 6.4.2 require much more computational effort, small average gaps on the test set would lead to an ML-CO algorithm that saves a lot of computation time on large-scale instances.

Learning process. In the learning phase, we use the permutation s = f (θ) as candidate solution to be compared with the associated target solution in the dataset. We use three types of probabilistic CO layers and associated Fenchel-Young loss functions: additive perturbation, multiplicative perturbation, and half square norm regularization. Our hyperparameters are detailed in Table 12. The value of ε is tuned with a grid search on a validation dataset for each probabilistic CO layer, selecting the one leading to the smallest average gap. We use a batch size of 1 sample and learn through 300 epochs in each case. Experiments are done on a computing cluster with processor 32 x Intel Xeon E5-2667, 3.30GHz (hyperthreading), and 192 Go of RAM.

Results

In Table 13, we compare the average train and test gaps given by each of the three probabilistic CO layers detailed in Table 12, as well as the GLM-based pipeline from [START_REF] Parmentier | Learning to Solve the Single Machine Scheduling Problem with Release Times and Sum of Completion Times[END_REF]. In each case, we derive solutions in the following way: given an instance x, we follow the pipeline (27), where the weights of the GLM w are fixed thanks to the learning phase based on the specific probabilistic CO layer (or by [START_REF] Parmentier | Learning to Solve the Single Machine Scheduling Problem with Release Times and Sum of Completion Times[END_REF]), and where the decoder is ψ = id. We recall that on the train set, gaps are computed with respect to the optimal solutions. On the test set, we consider the best solutions found by the algorithms detailed in Section 6.4.2. Overall, we manage to obtain low test gaps. We even get better solutions than the heuristics with our ML-CO pipelines on the test set. In Figure 6, we show the histograms of train and test gaps given by the pipeline trained with the additive perturbation, which gives the best average results in Table 13. It emphasizes a good generalization performance, even trained on small instances. Since ranking is very fast, it highlights the practical interest of this kind of ML-CO algorithm. For further studies in a learning by experience setting, as well as additional details on the decoder part, we refer to Parmentier (2021a) and [START_REF] Parmentier | Learning to Solve the Single Machine Scheduling Problem with Release Times and Sum of Completion Times[END_REF].

Figure 6: Histograms of train and test gaps in % given by the pipeline trained with the additive perturbation CO layer and decoder ψ = id.

Two-stage stochastic minimum weight spanning tree

Let us now focus on the two-stage minimum weight spanning tree problem. This last subsection has two main purposes: first, having an application on a two stage stochastic optimization problem where we use the single stage problem as a CO layer, and second, we show how to use InferOpt.jl with a pipeline containing a graph neural network (GNN), and compare it to a basic generalized linear model (GLM). Source code for this application can be found in the corresponding satellite package TwoStageSpanningTree.jl18 .

Problem formulation

Let G = (V, E) be an undirected graph, and S be a (finite) set of scenarios. The goal is to build a spanning tree on G over two stages at minimum cost, the second-stage building costs being unknown when first-stage decisions are taken.

y e + z es ≤ |Y | -1, ∀Y, ∅ ⊊ Y ⊊ V, ∀s ∈ S, y e ∈ {0, 1}, ∀e ∈ E, z es ∈ {0, 1}, ∀e ∈ E, ∀s ∈ S, (28a) (28b)
(28c) (28d) (28e)
where y e is a binary variable indicating if e is selected in the first-stage solution, and z es is a binary variable indicating if e is selected in the second-stage solution for scenario s. The two constraints (28b) and (28c) jointly ensures that y e + z es belongs to the spanning tree polytope. Formulation ( 28) contains an exponential number of constraints. Several mathematical programming approaches enable to solve it. Instances with up to 50 vertices can be solved to optimality using a cutgeneration approach or a Benders decomposition. A Lagrangian relaxation and a Lagrangian heuristic can provide solutions within a 2% optimality gap for instances with up to 10000 vertices in approximately one hour. A detailed description of this heuristic can be found in Appendix B.2.1.

Datasets

We build three datasets, one for training, one for validation, and one for testing. Each dataset contains 600 grid graphs instances whose widths range from 10 to 60 (that is, with 100 to 3600 vertices), and from 5 to 20 scenarios. For all instances, first stage costs are drawn uniformly between 0 and 20. For the second stage costs, depending on the instances, they are drawn uniformly between 0 and n, with n varying between 10 and 30. Each instance is labeled with a solution given by the Lagrangian heuristic.

InferOpt pipeline

The learning pipeline is presented in Figure 8.

ML predictor

The feature matrix is given as input to an ML layer φ w with learnable parameters w, which predicts edge weights θ e . We tested and compared two different predictors: first a GLM such that θ e = w ⊤ ϕ(x, e), and then a GNN, in particular the GraphConv architecture [START_REF] Morris | Weisfeiler and Leman Go Neural: Higher-Order Graph Neural Networks[END_REF]. 1

Combinatorial layer We use the predicted edge weights θ as the objective of a maximum weight forest problem layer. A maximum weight forest can be efficiently found using Kruskal's algorithm.

Decoder As a post-processing, the output forest of the training pipeline is interpreted as a first stage solution of instance x. For each scenario, we can complete it into a minimum cost spanning tree containing this forest, using Kruskal's algorithm.

Experimental setting

We train each model by imitation using a FenchelYoungLoss with PerturbedAdditive regularization.

Hypeparameters used for each model can be found in table 14. We use the Flux.jl Julia library and the Adam optimizer for training. We construct our GNN with the GraphNeuralNetworks.jl19 Julia library [START_REF] Lucibello | GraphNeuralNetworks.jl: A Geometric Deep Learning Library for the Julia Programming Language[END_REF]. The GNN architecture contains 3 GraphConv layers with hidden size 50.

ε
The GLM and GNN models are both trained on each of the two feature sets. Features are normalized by dividing each feature by its train dataset standard deviation. We optimize hyperparameters and perform model selection by using the average optimality gap values on the validation dataset. The experiments are conducted on a 2.6 GHz Intel Core i9-11950H, 16 cores, with 64 Go RAM. Epochs with best gaps on the validation dataset are selected for each experiment.

Results

The CO layer takes the same CPU time as one iteration of the Lagrangian relaxation, and the feature computation takes about the same time as 2 iterations. Therefore, resulting pipelines are roughly 15000 times faster than the Lagrangian heuristic.

We evaluate the performance of the trained pipelines on the test dataset. For this we compute the average, minimum, and maximum gaps respect to the lower bound and upper bounds respectively given by the Lagrangian relaxation and the Lagrangian heuristic. Results are gathered in Table 15 On both feature sets we observe low gaps for both GLM and GNN, with the GNN performing better. We also see that the GLM trained on advanced features performs better that the GNN on basic features, which means that feature engineering is more effective than replacing the linear model with a more complex neural network. Advanced features combined with the GNN performs best. As can be seen in section 6.3 on the stochastic VSP, we could maybe further improve performance if we learn by experience instead of learning by imitation. Indeed, the labels computed by the Lagrangian heuristic are not necessarily optimal and probably introduce bias in the learning.

Conclusion

In this paper, we introduce InferOpt.jl, a Julia package that seamlessly turns CO algorithms into layers of ML pipelines. The package is compatible with most approximate differentiation techniques and structured losses from the literature, but also includes novel ones. Its theoretical foundation is a probabilistic point of view on CO layers, which unifies existing and new methods. Our probabilistic framework notably provides natural "learning by experience" counterparts to the more common "learning by imitation" approaches.

InferOpt.jl also unlocks the use of large neural networks in hybrid pipelines, even when learned by experience. We demonstrate it with a hybrid ML-CO pipeline containing a CNN with tens of thousands of parameters, which to the best of our knowledge, has never been achieved. Using InferOpt.jl, we benchmark pipeline architectures and learning algorithms. In general, we observe that hybrid pipelines can yield fast algorithms with state-of-the-art performance on well-studied applications such as the single machine scheduling problem. Depending on the learning setting, we may require more or less information in the training set. Unsurprisingly, the more information we have in the training set, the better our predicted solutions become. This is notably clear on the Warcraft shortest path problem. However, on very hard combinatorial tasks where no algorithm scales on larges instances, learning by experience can outperform learning by imitation, as can be seen on the stochastic vehicle scheduling problem.

When we learn by imitating precomputed solutions, Fenchel-Young losses seem to outperform alternatives. If we use a CO algorithm that is sensitive to the sign of the input, the multiplicative perturbation seems the best approach. Concerning the pipeline architecture, the most important aspect seems to be the choice of the combinatorial optimization layer and of the features. Using GNNs instead of GLMs for encoding can be beneficial, as we remarked on the two-stage spanning tree problem.

We also highlight some perspectives. In this paper, we learn the objective function of CO problems based on the notion of probabilistic CO layers. We do not consider the case when the desired flexibility concerns the constraints, as studied by [START_REF] Paulus | CombOptNet: Fit the Right NP-Hard Problem by Learning Integer Programming Constraints[END_REF] for instance, even though this perspective would be equally natural. Furthermore, we make the assumption that the CO oracle derives optimal solutions within our layer. We briefly mention the case of inexact CO oracles, but additional research would be necessary in this direction. Finally, a large part of the ML literature is dedicated to reinforcement learning. The connection between this field and structured learning deserves further exploration.

A Proofs

In this Section, the dominated convergence theorem is used implicitly to justify any differentiation under the integral sign. When dealing with polyhedral functions such as θ -→ max v∈V θ ⊤ v, we often write ∇ θ for simplicity even though they are only subdifferentiable, because the set of non-differentiability has measure zero. We denote the standard Gaussian density by ν.

A.1 Additive perturbation

These proofs are already given by Berthet et al. ( 2020), but we include them for comparison purposes.

A.1.1 Derivatives

Proof of Proposition 3.1.

Proof. The following change of variable is a diffeomorphism:

u = θ + εz ⇐⇒ u -θ ε = z.
We apply it to the definition of p + ε (v|θ).

p + ε (v|θ) = R d 1{f(θ + εz) = v} ν(z) dz = R d 1{f(u) = v} ν u -θ ε du ε d .
We now differentiate with respect to θ before applying the reverse change of variable:

∇ θ p + ε (v|θ) = R d 1{f(u) = v} -1 ε ∇ν u -θ ε du ε d = -1 ε R d 1{f(θ + εz) = v} ∇ν(z) dz = 1 ε R d 1{f(θ + εz) = v} zν(z) dz.
The last equality holds because the standard Gaussian density satisfies ∇ν(z) = -zν(z). From there, we deduce the Jacobian of f + ε (θ):

J θ f + ε (θ) = v∈V v∇ θ p + ε (θ, y) ⊤ = 1 ε R d v∈V v 1{f(θ + εz) = v} f (θ+εz) z ⊤ ν(z) dz 40
We arrive at the following simple expression, which was already given by [START_REF] Berthet | Learning with Differentiable Perturbed Optimizers[END_REF]:

J θ f + ε (θ) = 1 ε E f (θ + εZ)Z ⊤ A.1.2 Regularization Proof of Proposition 3.2. Proof. Because Ω + ε = (F + ε ) * is a Fenchel conjugate, it is automatically convex. Furthermore, Ω + ε (µ) = sup θ∈R d θ ⊤ µ -F + ε (θ) = sup θ∈R d θ ⊤ µ -E max v∈V (θ + εZ) ⊤ v .
We consider µ / ∈ conv(V). By convex separation, there exists θ ∈ R d and α > 0 such that θ⊤ µ ≥ α + θ⊤ v for all v ∈ V. This implies that, for all t > 0,

Ω + ε (µ) ≥ t θ⊤ µ -E max v∈V t θ + εZ ⊤ v ≥ t θ⊤ µ -tE max v∈V θ⊤ v -εE max v∈V Z ⊤ v ≥ tα -εE max v∈V Z ⊤ v ----→ t→+∞ +∞
We have shown that µ / ∈ dom(Ω + ε ), and therefore dom(Ω + ε ) ⊂ conv(V). We define

f Z (θ, v) = (θ + εZ) ⊤ v so that F + ε (θ) = E max v∈V f Z (θ, v) .
Danskin's theorem helps us compute the gradient of F + ε :

∇ θ F + ε (θ) = E ∇ θ max v∈V f Z (θ, v) = E ∇ 1 f Z θ, argmax v∈V f Z (θ, v) = E argmax v∈V f Z (θ, v) = E argmax v∈V (θ + εZ) ⊤ v = f + ε (θ).
As shown by Berthet et al. (2020, Proposition 2.2), the function Ω + ε is a Legendre type function, which means that

∇ θ F + ε = ∇ θ (Ω + ε ) * = (∇ θ Ω + ε ) -1 . From this, we deduce ∇ θ F + ε (θ) = argmax µ∈R d θ ⊤ µ -Ω + ε (µ) = argmax µ∈dom(Ω + ε ) θ ⊤ µ -Ω + ε (µ) = f Ω + ε (θ).
Hence, we can conclude:

f + ε (θ) = ∇ θ F + ε (θ) = argmax µ∈dom(Ω + ε ) {θ ⊤ µ -Ω + ε (µ)} = f Ω + ε (θ).
To recover the formula given in Proposition 3.2, we simply remember that dom(Ω + ε ) ⊆ conv(V).

A.2 Multiplicative perturbation

A.2.1 Derivatives

Proof of Proposition 3.3.

Proof. Suppose θ ∈ R d only has positive components. Then the following change of variable is a diffeomorphism:

u = θ ⊙ e εz-ε 2 1/2 ⇐⇒ log(u) -log(θ) ε + ε1 2 = z
We apply it to the definition of p ⊙ ε (v|θ).

p ⊙ ε (θ, y) = R d 1 f θ ⊙ e εz-ε 2 1/2 = v ν(z) dz = (0,+∞) d 1{f(u) = v} ν log(u) -log(θ) ε + ε1 2 du ε d i u i
.

We now differentiate with respect to θ before applying the reverse change of variable:

∇ θ p ⊙ ε (θ, y) = (0,+∞) d 1{f(u) = v} -1 εθ ⊙ ∇ν log(u) -log(θ) ε + ε1 2 du ε d i u i = -1 εθ ⊙ R d 1 f θ ⊙ e εz-ε 2 1/2 = v ∇ν(z) dz = 1 εθ ⊙ R d 1 f θ ⊙ e εz-ε 2 1/2 = v zν(z) dz.
From there, we deduce the Jacobian of f ⊙ ε (θ):

J θ f ⊙ ε (θ) = v∈V v∇ θ p ⊙ ε (θ, y) ⊤ = 1 εθ ⊙ R d v∈V v 1 f θ ⊙ e εz-ε 2 1/2 = v f (θ⊙e εz-ε 2 1/2 ) z ⊤ ν(z) dz
We arrive at a simple variant of the previous expression:

J θ f ⊙ ε (θ) = 1 εθ ⊙ E f θ ⊙ e εz-ε 2 1/2 Z ⊤ A.2.2 Regularization Proof of Proposition 3.4. Proof. Because Ω ⊙ ε = (F ⊙ ε ) * is a Fenchel conjugate, it is automatically convex. Furthermore, Ω ⊙ ε (µ) = sup θ∈R d θ ⊤ µ -F ⊙ ε (θ) = sup θ∈R d θ ⊤ µ -E max v∈V θ ⊙ e εZ-ε 2 1/2 ⊤ v = sup θ∈R d θ ⊤ µ -E max v∈V θ ⊤ v ⊙ e εZ-ε 2 1/2 .
This last expression shows why we do not have dom(Ω ⊙ ε ) ⊂ conv(V) (unlike in the additive case). Indeed, even when µ / ∈ conv(V), the multiplicative scaling of v might allow it to compensate the inner product θ ⊤ µ and stop Ω ⊙ ε (µ) from going to +∞. We define

f Z (θ, v) = θ ⊙ e εZ-ε 2 1/2 ⊤ v so that F ⊙ ε (θ) = E max v∈V f Z (θ, v) .
Danskin's theorem helps us compute the gradient of F ⊙ ε :

∇ θ F ⊙ ε (θ) = E ∇ θ max v∈V f Z (θ, v) = E ∇ 1 f Z θ, argmax v∈V f Z (θ, v) = E e εZ-ε 2 1/2 ⊙ argmax v∈V f Z (θ, v) = E e εZ-ε 2 1/2 ⊙ argmax v∈V θ ⊙ e εZ-ε 2 1/2 ⊤ v = f ⊙scaled ε (θ) ̸ = f ⊙ ε (θ)
We could prove in a way similar to Berthet et al. (2020, Proposition 2.2) that Ω ⊙ ε is a Legendre type function, which means that

∇ θ F ⊙ ε = ∇ θ (Ω ⊙ ε ) * = (∇ θ Ω ⊙ ε ) -1 . From this, we deduce ∇ θ F ⊙ ε (θ) = argmax µ∈R d θ ⊤ µ -Ω ⊙ ε (µ) = argmax µ∈dom(Ω ⊙ ε ) θ ⊤ µ -Ω ⊙ ε (µ) = f Ω + ε (θ).
This time we cannot replace dom(Ω ⊙ ε ) by conv(V), but we still obtain a similar conclusion:

f ⊙scaled ε (θ) = ∇ θ F ⊙ ε (θ) = argmax µ∈dom(Ω ⊙ ε ) {θ ⊤ µ -Ω ⊙ ε (µ)} = f Ω ⊙ ε (θ).

A.3 Inexact oracles

Proof of Proposition 3.5.

Proof. We start with additive perturbation. By Proposition 3.1, we have:

J θ g + ε (θ) -J θ f + ε (θ) = 1 ε E g(θ + εZ)Z ⊤ - 1 ε E f (θ + εZ)Z ⊤ = 1 ε E (g(θ + εZ) -f (θ + εZ)) Z ⊤ .
We bound the spectral norm of the error using Jensen's inequality:

J θ g + ε (θ) -J θ f + ε (θ) 2 ≤ 1 ε 2 E (g(θ + εZ) -f (θ + εZ)) Z ⊤ 2 = 1 ε 2 E ∥g(θ + εZ) -f (θ + εZ)∥ 2 ∥Z∥ 2 ≤ 1 ε 2 ∥g -f ∥ 2 ∞ E ∥Z∥ 2 = d ε 2 ∥g -f ∥ 2 ∞ .
We now move on to multiplicative perturbation. For multiplicative perturbation, following the same proof starting from Proposition 3.3 yields a similar inequality:

J θ g ⊙ ε (θ) -J θ f ⊙ ε (θ) 2 ≤ d ε 2 min i |θ i | 2 ∥g -f ∥ 2 ∞ .
B More details on applications Since the constraint matrix is totally unimodular, the constraint polytope gives a perfect formulation, therefore binary constraints can be relaxed, and this problem can be solved using either a flow algorithm or a linear programming solver. Since there is an exponential number of variables, this formulation can be solved using a column generation algorithm. We denote λ v the dual variables of (31). The associated sub-problem is a constrained shortest path problem : min

P ∈P (c P - v∈P λ v ) (32)
The sub-problem (32) can be solved using a stochastic constrained shortest path algorithm (see Parmentier 2019 for in-depth details about these algorithms, and ConstrainedShortestPaths.jl20 for a Julia implementation). Solving (31) using column generation gives good bounds on instances with up to 50 tasks. A Branchand-price could be used to try to find the optimal solution.

B.1.3 Data generation

Instance generator To generate our dataset, we use a similar approach as by Parmentier 2021b. We define a squared map of 50 minutes width, divided in 25 squared (10 × 10) districts. Each task v has a (uniformly drawn) start point, a start time t b v , end point, and end time t e v . Let v be a task corresponding to a trip between district d 1 and d 2 . We compute the perturbed start and end times are:

ξ b v = t b v + β v ξ e v = ξ b v + t e v -
with ξ 1 = ξ e u , ξ 2 = ξ 1 + ζ dis d1,h(ξ1) , and ξ 3 = ξ 2 + t tr a + ξ inter h(ξ2) .

Encoding features For an instance, we compute a matrix of 20 features for every arc of the corresponding graph. Let a = (u, v) ∈ A. The first feature is the length of a in minutes, i.e. the deterministic travel time between u and v. The second feature is equal to c vehicle if a is connected to the source (i;.e if u = o), else 0. Then, we add the 9 deciles of ξ b v -(ξ e u + ξ tr a ) which represents the slack between u and v. Finally, we add the values of the cumulative probability distribution of the slack, evaluated in -100, -50, -20, -10, 0, 50, 200, and 500.

Label solutions We label each instance with a local search heuristic initialized with the solution of the deterministic problem. Its implementation can be found on the GitHub repository21 .

B.2.1 Lagrangian heuristic

We present the heuristic based on Lagrangian relaxation used to label our datasets.

Lagrangian relaxation An option to compute a lower bound on the optimal value of the two-stage minimum weight spanning tree is to apply a Lagrangian relaxation. For this, we duplicate complicating variables y for each scenario and introduce associated θ es dual variable for each new constraint:

where M is a large constant. In theory M is not needed, we can take M = +∞, but taking a finite M leads to more informative gradients. G can be maximized by gradient ascent:

(∇ θ G(θ)) es = 1 |S| (y ey es ) (39)

Lagrangian heuristic Once we find a good θ with the Lagrangian relaxation, we can retrieve the corresponding y es solution, and use it as input of the Lagrangian heuristic 1. The general idea is to rebuild a forest from the edges that are selected in the first stage for most scenarios. In the numerical experiments, we use the Lagrangian heuristic to build the training set. Practically, we stop the Lagrangian relaxation after 50.000 iterations or if the gap between the lower bound provided by the G(θ) and the upper bound obtained from the heuristic is non-grater than 0.1%.

  , Graphs.jl 10 (Fairbanks et al. 2021), GridGraphs.jl 11 , Metalhead.jl 12 and Zygote.jl 13 (Innes 2019), along with the base libraries LinearAlgebra and Statistics. Code sample 1 creates the CNN encoder layer. Meanwhile, Code sample 2 shows how to define the Dijkstra oracle and the true cost function.

  for clarity. Please refer to the documentation of InferOpt.jl and satellite packages for actual runnable examples.

  )where ε > 0 controls the amplitude of the perturbation, Z ∼ N (0, I) is a standard Gaussian vector and the expectation is taken with respect to Z unless otherwise specified. Choosing ε is a trade-off between § ¤

  v|θ).

  by Jensen's inequality, which means the expected regret is an upper bound.Code sample 7 demonstrates how to define an expected regret from probabilistic CO layers, while Code sample 8 shows that we can compute and differentiate it automatically. Finally, Code sample 9 displays a complete program for learning by experience. The rest of this section explains how to compute derivatives of the expected regret R p and can be skipped without danger. Code sample 9: Learning with an expected regret Remark 4.1.

  Code sample 12: Learning with an imitation loss

Figure 4 :

 4 Figure 4: Pipeline for the stochastic VSP

Figure 5 :

 5 Figure 5: Comparison of predictions of three different models on the same 50-task instance. Each circle represents a task, with color value and size corresponding to its total delay (intrinsic + propagated) in the predicted solution. Tasks are ordered respect to start time along the x axis. Each arrow path represents a vehicle path in the predicted solution.

•

  A train dataset made of 420 instances with n ∈ {50, 70, 90, 110} jobs, and corresponding solutions derived by means of the branch-and-memorize algorithm (Shang, T'Kindt, and Della Croce 2021). Each value of n is equally represented. • A test dataset made of 1820 instances with n ∈ {50 + k × 100, k ∈ [25]} jobs, and the associated best solutions derived by the algorithms noted above. Each value of n is equally represented.

Figure 7 :

 7 Figure 7: Example of a two stage spanning tree on a 3 × 4 grid graph. Blue edges correspond to first-stage forest, and red dashed edges correspond to second-stage forest.

Figure 8 :

 8 Figure 8: Two-stage minimum spanning tree pipeline

B. 1

 1 Stochastic vehicle scheduling problemB.1.1 Linear program formulation of (VSP):The deterministic Vehicle Scheduling Problem (VSP), can be formulated as the following linear program:

Figure 9 :

 9 Figure 9: Example of a map with 10 tasks. Task 0 represents the depot, in the middle of the map. Square with task indices are starting points, and tasks without are end points.

Figure 11 :

 11 Figure 11: Learning by experience metric evolution during learning

d

  z es = |V | -1 ∀s ∈ S e∈E(Y ) y e + z es ≤ |Y | -1, ∀Y, ∅ ⊊ Y ⊊ E, ∀s ∈ S y es = y e , ∀e ∈ E, ∀s ∈ S (θ es dual variables) y, z ∈ {0, 1}(37)Lagrangian dual problem can be formulated as follows: es z esθ es y es s.t. y s + z s ∈ spanning tree polytope ,

Algorithm 1 :

 1 Lagrangian heuristic Data: y es ∈ {0, 1}, ∀e ∈ E, ∀s ∈ S forall e ∈ E do if 1 |S| s∈S y es > 0.5 then w e ← 1; else w e ← -1; end end Find maximum weight forest for weights w; Result: y e ∈ {0, 1}, ∀e ∈ E

Figure 3 :

 3 Figure 3: Train and test optimality gaps along training in the Warcraft application

Table 2 :

 2 A common decomposition for loss functions in imitation learning § ¤

			§	¤
	using InferOpt		using Zygote
	fyl_pert = FenchelYoungLoss(perturbed_add)		L = loss(theta, y_ref)
	fyl_reg = FenchelYoungLoss(regularized) spol = SPOPlusLoss(warcraft_maximizer) ¦	¥	Zygote.gradient(loss, theta, y_ref) ¦ Code sample 11: Supported operations for an ¥
	Code sample 10: Example imitation losses		imitation loss

Table 3 :

 3 Hyperparameters for learning Warcraft shortest paths.

		Description
	epsilon	Scale of the noise for perturbation.
	nb_samples	Number of noise samples M for perturbation.
	batch_size	Size of the batches to compute gradients.
	lr_start	Starting learning rate.

Table 4 :

 4 Learning settings for Warcraft shortest paths.

	Seeting name	CO problem (CO oracle)	Probabilistic CO layer	Exp./Imit. Target	Loss
	Cost perturbed multiplicative noise	SP with non-negative costs (Dijkstra)	Multiplicative perturbation	Experience No target	Perturbed cost
	Cost perturbed additive noise	SP on an extended acyclic graph (Ford-Bellman)	Additive perturbation	Experience No target	Perturbed cost
	Cost regularized half square norm	SP on an extended acyclic graph (Ford-Bellman)	Half square norm	Experience No target	Regularized cost
	SPO+	SP on an extended acyclic graph (Ford-Bellman)	No regularization	Imitation Cost and path	SPO+ loss
	MSE perturbed multiplicative noise	SP with non-negative costs (Dijkstra)	Multiplicative perturbation	Imitation Path	Mean squared error
	MSE regularized half square norm	SP on an extended acyclic graph (Ford-Bellman)	Half square norm	Imitation Path	Mean squared error
	Fenchel-Young perturbed multiplicative noise	SP with non-negative costs (Dijkstra)	Multiplicative perturbation	Imitation Path	Fenchel-Young
	Fenchel-Young perturbed additive noise	SP on an extended acyclic graph (Ford-Bellman)	Additive perturbation	Imitation Path	Fenchel-Young
	Fenchel-Young regularized half square norm	SP on an extended acyclic graph (Ford-Bellman)	Half square norm	Imitation Path	Fenchel-Young

Table 6 :

 6 Hyperparameter for Stochastic Vehicle Scheduling experiments

	1	50	100	300

Table 8 :

 8 Learning by imitation: cost gap

Table 10 :

 10 Learning by experience: cost gap

Table 11 :

 11 Learning by experience: average cost per task

	Train dataset	25	Test dataset (number of tasks in each instance) 50 100 200 300 500 750	1000
	25 tasks	274.19 224.55 204.9 191.86 184.71 181.29 178.0 177.02
	50 tasks	274.12 224.51 205.0 191.85 184.3 180.48 176.96 176.0
	100 tasks	274.13 224.41 205.0 191.85 184.63 181.08 177.81 176.74

Table 13 :

 13 Train and test average gaps in % for pipelines with decoder ψ = id.

	Train set gap (%) Test set gap (%)

  For each edge e ∈ E and scenario and s ∈ S, we denote by c e in R the scenario-independent first stage cost of building e, and by d es the scenario-dependent second stage cost. The two-stage spanning tree problem can be formulated as the following ILP:

	min y,z	e∈E	c e y e +	1 |S| e∈E s∈S	d es z es ,
	s.t.	e∈E	y e + z es = |V | -1,	∀s ∈ S,
		e∈E(Y )

Table 14 :

 14 Hypeparameters for two-stage spanning tree experiments

		nb_samples learning rate batch size training epochs
	GLM 1	20	0.01	1	20
	GNN 1	20	0.001	1	20

Table 15 :

 15 below. Test gaps in % for models trained with the basic features

		Optimality gap (%)	Lagrangian heuristic gap (%)
		Average Min Max Average Min	Max
	GLM on basic features	4.77	1.08 12.33	2.91	0.71	8.81
	GNN on basic features	3.06	0.09 12.09	1.23	-0.24	7.07
	GLM on advanced features	2.21	0.09	7.66	0.41	-0.73	3.12
	GNN on advanced features	2.11	0.09	7.91	0.31	-0.73	2.86

  B.1.2 Two MILP formulations for the StoVSPMILP formulation We denote respectively c vehicle and c delay the cost of one vehicle, and the cost of one minute of delay. The StoVSP problem can be formulated as the following MIP with quadratic constraints:These quadratic delay constraints (30b) can be linearized using McCormick inequalities, and we can then solve the resulting MILP with industrial solvers. This approach does not scale well on instances with large number of tasks and scenarios.Column generation formulation Another possible formulation for StoVSP is the column generation approach. Let P the set of feasible vehicle routes. Let us define the cost of a vehicle route P ∈ P by c s P = c vehicle + c delay × v∈P d s v . We introduce decision variables y P ∈ {0, 1} that equals 1 if the vehicle route P is chosen. We can then formulate StoVSP as follows:

	min d,y	c delay	1 |S| s∈S v∈V \{o,d}	d s v + c vehicle	a∈δ + (o)	y a	(30a)
	s.t.	a∈δ -(v)	y a =	a∈δ + (v)	y a	∀v ∈ V \{o, d}	(30b)
		a∈δ -(v)	y a = 1		∀v ∈ V \{o, d}	(30c)
	d s v ≥ γ s v +	a∈δ -(v)	(d s u -s s u,v )y a	∀v ∈ V \{o, d}, ∀s ∈ S	(30d)
					a=(u,v)
	d s v ≥ γ s v y a ∈ {0, 1}			∀v ∈ V \{o, d}, ∀s ∈ S ∀a ∈ A	(30e) (30f)
					min	1 |S| s∈S P ∈P	c s P y P
						s.t.		(31)
						y P ∈ {0, 1}	∀p ∈ P

p∋v y P = 1 ∀v ∈ V \{o, d} (λ v ∈ R)

  = ξ b v , ξ 2 = ξ 1 + ζ dis d1,h(ξ1), andξ 3 = ξ 2 + t e vt b v + ξ inter h(ξ2). And the perturbed travel times along arcs a = (u, v):

	t b v + ζ dis d1,h(ξ1) + ζ inter h(ξ2) + ζ dis d2,h(ξ3)	(35)
	with ξ 1 ξ tr a = ξ b v + t tr a + ζ dis d1,h(ξ1) + ζ inter h(ξ2) + ζ dis d2,h(ξ3)

https://github.com/axelparmentier/InferOpt.jl

https://github.com/JuliaDiff/ChainRules.jl

https://github.com/google/jaxopt

https://github.com/facebookresearch/theseus

https://github.com/cvxgrp/cvxpylayers

https://github.com/jump-dev/DiffOpt.jl

https://github.com/khalil-research/PyEPO

https://github.com/ZIB-IOL/FrankWolfe.jl

https://github.com/LouisBouvier/WarcraftShortestPaths.jl

https://github.com/BatyLeo/StochasticVehicleScheduling.jl

https://github.com/axelparmentier/SingleMachineScheduling.jl

https://github.com/axelparmentier/MinimumWeightTwoStageSpanningTree.jl

https://github.com/CarloLucibello/GraphNeuralNetworks.jl

https://github.com/BatyLeo/ConstrainedShortestPaths.jl

https://github.com/BatyLeo/StochasticVehicleScheduling.jl

Acknowledgements

The authors want to thank Mathieu Besançon for his support with various Julia packages, as well as his expert proofreading. We also thank Éloïse Berthier, Alexandre Forel, Jérôme Malick, Clément Mantoux, and Pierre Marion for their time and advice.

B.1.4 Training plots

B.2.2 Instance features

We detail here all the features used for training.

Basic feature set Basic features only contain twelve features. Let e be an edge. The first one is the first stage cost c e , and the other 11 are the quantiles of the second stage costs (d es ) s∈S for values {0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9, 1}. All quantiles in the advanced feature set use the same values.

Advanced feature set Advanced features contain 79 features in total. Let e be an edge. The first 12 features are the same as the basic features. To these, we add quantiles of best stage costs (i.e for each scenario s the value min(c e , d es )), quantiles of neighbors first stage costs, and quantiles of neighbors second stage costs. Then, we compute several single stage minimum spanning trees with Kruskal's algorithm: one with first stage costs, one for each scenario with corresponding second stage costs, and one for each scenario with corresponding best stage costs min(c e , d es ). We use these spanning trees to compute new useful features: a binary features that tells if e is in the first stage spanning tree, and quantile features on its presence in second stage trees. Finally, we use best stage spanning trees to compute quantiles for the cost of e in the best stage costs spanning tree: one set of quantiles for first edge costs, and another set for second stage costs.