11affiliationtext: TU Wien, Vienna, Austria. 22affiliationtext: Institut
Camille Jordan, École Centrale Lyon, CNRS UMR 5208, France.33affiliationtext:
Institut Universitaire de France (IUF)44affiliationtext: CERMICS, École des
Ponts, France.55affiliationtext: SIERRA, INRIA Paris, France.

# Generalization Bounds of Surrogate Policies
for Combinatorial Optimization Problems

Pierre Cyril Aubin-Frankowski  Yohann De Castro  Axel Parmentier  Alessandro
Rudi

###### Abstract

A recent stream of structured learning approaches has improved the practical
state of the art for a range of combinatorial optimization problems with
complex objectives encountered in operations research. Such approaches train
policies that chain a statistical model with a surrogate combinatorial
optimization oracle to map any instance of the problem to a feasible solution.
The key idea is to exploit the statistical distribution over instances instead
of dealing with instances separately. However learning such policies by risk
minimization is challenging because the empirical risk is piecewise constant
in the parameters, and few theoretical guarantees have been provided so far.
In this article, we investigate methods that smooth the risk by perturbing the
policy, which eases optimization and improves generalization. Our main
contribution is a generalization bound that controls the perturbation bias,
the statistical learning error, and the optimization error. Our analysis
relies on the introduction of a uniform weak property, which captures and
quantifies the interplay of the statistical model and the surrogate
combinatorial optimization oracle. This property holds under mild assumptions
on the statistical model, the surrogate optimization, and the instance data
distribution. We illustrate the result on a range of applications such as
stochastic vehicle scheduling. In particular, such policies are relevant for
contextual stochastic optimization and our results cover this case.

##  1 Introduction

###  1.1 Learning policies instead of minimizing

Operations Research (OR) aims to optimize industrial processes by allocating
resources efficiently. Often, these resources are indivisible, such as a truck
that either delivers a request or not. Indivisible ressources lead to
combinatorial optimization problems, where larger instances offer more
flexibility and economy of scale. Consider for instance a vehicle routing
problem: if there are more requests to deliver in the same neighborhood, the
distance travelled between two requests is smaller, and the marginal costs per
request decreases. However, achieving economy of scale requires algorithms
that can handle huge instances. Exploring large combinatorial spaces (with up
to $10^{100}$ solutions) is extremely challenging and most of the research
focus has been on scalable algorithms. This scalability has been achieved at
the price of simplifying the objective function, a fact reflected in the OR
practitioners’ tools, where mixed integer linear programming (MILP) dominates.
When it comes to making algorithms scale, focusing on MILP has been a success
story, with solvers in 2015 being $540$ billion times faster than in 1991 as
recalled by Bertsimas et al. (2016), and instances with millions of variables
and constraints being solved routinely in the industry.

Consider now the following combinatorial optimization problem where we have
made the dependency on the instance ${\bm{x}}$ in $\mathcal{X}$ explicit

$\min_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}f^{0}({\bm{y}},{\bm{x}})\,$ (1.1)

where $\mathcal{Y}({\bm{x}})$ is the finite set of feasible solutions for
instance ${\bm{x}}\in\mathcal{X}$. The objective $f^{0}$ may model any complex
phenomena. We assume that $f^{0}$ comes with an oracle, that can be mildly
expansive to compute, i.e., can take a few seconds to evaluate and for
instance be the result of a simulation. We focus on settings where the problem
defined by $f^{0}$ is challenging in two ways: i) its modeling is too complex
to be cast as a mathematical programming formulation that can be solved with
off the shelf algorithms; ii) its computational cost prohibits the use of a
local descent heuristic. Since combinatorial algorithms must scale to be
useful in OR, practitioners generally simplify (1.1) into

$\min_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}\tilde{f}^{0}({\bm{y}},{\bm{x}})\,$
(1.2)

where $\tilde{f}^{0}$ a possibly-over-simplified deterministic objective
amenable to classical algorithms. Such algorithms are implicitly supposed to
work on any instance of the problem, and distributional aspects over the
instances, i.e., that instances that have to be solved in practice have a
specific structure, have been overlooked in the literature. This results in
worst case analysis to be the main setting considered by complexity theory.

On the other hand, a recent line of works detailed below has proposed with
some successes a learning methodology that directly address the objective of
(1.1) in its full complexity. Instead of solving each instance independently,
we consider a family $\mathcal{H}$ of policies
$h:{\bm{x}}\in\mathcal{X}\mapsto{\bm{y}}\in\mathcal{Y}({\bm{x}})$. Our goal is
to learn a mapping $h^{\star}$ that minimizes the risk

$\min_{h\in\mathcal{H}}\mathcal{R}(h):=\mathds{E}_{X}\big{[}f^{0}(h(X),X)\big{]},$
(1.3)

where $X$ is a random instance distributed according to some distribution
$\mathds{P}_{X}$. If the set $\mathcal{H}$ corresponds to all the possible
measurable mappings
${\bm{x}}\in\mathcal{X}\to(\mathcal{Y}({\bm{x}}))_{{\bm{x}}\in\mathcal{X}}$,
then we get back to the classical approach where each instance is solved
independently.

We instead restrict ourselves to a family of tractable policies parametrized
by a vector $\bm{w}$ belonging to a set $\mathcal{W}$, and we learn the
parameter $\bm{w}$ that minimizes the risk. For such an approach to work, the
policy mappings of $\mathcal{H}$ need to be able to properly explore the large
combinatorial set of solutions $\mathcal{Y}({\bm{x}})$. One approach has
proved successful on a variety of applications such as portfolio optimization
(Elmachtoub and Grigas, 2021), path finding in images (Pogancic et al., 2020),
vehicle scheduling (Parmentier, 2021a), machine scheduling (Parmentier and
T’Kindt, 2021), dynamic vehicle routing (Baty et al., 2024) in the last few
years. It consists in delegating the exploration of $\mathcal{Y}({\bm{x}})$ to
the linear optimization problem

$\hat{y}(\bm{\theta})\in\operatorname*{arg\,max}_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}\langle{\bm{y}},\bm{\theta}\rangle\,,$
(1.4)

which we assume tractable, and where
$\hat{\bm{y}}\,:\,\theta\in\mathcal{Y}({\bm{x}})\mapsto\hat{\bm{y}}(\bm{\theta})\in\mathcal{Y}({\bm{x}})$
is an oracle returning an optimal solution of (1.4). We omit the dependency of
$\hat{\bm{y}}$ in ${\bm{x}}$ throughout the paper, for sake of readability. In
Equation (1.4), vector $\bm{\theta}$ is a direction in
$\mathds{R}^{d({\bm{x}})}$ and $\langle\cdot,\cdot\rangle$ denotes the
standard inner product. A statistical model ${\psi}_{\bm{w}}$ is used to find
a direction vector $\bm{\theta}$ that leads to promising solutions. In other
words, we consider the following family of parameterized policy mappings

$\mathcal{H}:=\big{\\{}h_{\bm{w}}\,:\,{\bm{x}}\in\mathcal{X}\mapsto\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))\in\mathcal{Y}({\bm{x}})\,,\
{\bm{w}}\in\mathcal{W}\big{\\}}\,,$ (1.5)

Note that the objective of (1.4) nevertheless relies on an embedding of
$\mathcal{Y}({\bm{x}})$ in $\mathds{R}^{d({\bm{x}})}$, whose dimension
$d({\bm{x}})$ may depend on the instance considered.

In a nutshell, instead of approximating the objective function $f^{0}$ in
(1.1) by a simpler $\tilde{f}^{0}$ as in (1.2), we study an average error
$\mathds{E}_{X}\big{[}f^{0}(h(X),X)\big{]}$ as in (1.3). The latter involves a
statistical model $\mathcal{H}$ defined in (1.5) and derived from the linear
combinatorial optimization problem (1.4). The role of the model is to explore
efficiently the solution set $\mathcal{Y}({\bm{x}})$ and we do not seek to
approximate $f^{0}$ in any way. In practice, as usual in machine learning, it
is the distribution $\mathds{P}_{X}$ in (1.3) that is approximated through
samples $(X_{1},\ldots,X_{n})$, in an empirical risk minimization fashion. Our
Theorem 3 provides error guarantees when $n$ goes to $+\infty$ to show the
convergence to (1.3).

The idea of transforming an optimization problem into a learning one is
motivated by industrial practice. In many real-world scenarios, the problem
instances that practitioners in the industry need to solve come from a
specific but unknown distribution. For example, an airline optimizes the
routes of its airplanes based on a schedule and some operational constraints
that do not vary much from one day to another. The traditional approach does
not take advantage of this distributional information. Learning is a sensible
strategy as it allows the algorithm to adapt to the problem structure.
Moreover, it enables to focus on optimizing the true objective $f^{0}$ instead
of pushing the optimization of (1.2) beyond the statistical accuracy of
$\tilde{f}$. Furthermore, the choice of $\mathcal{H}$ allows the combinatorial
aspects to be handled only through the linear optimization oracle. This
enables the use of very efficient algorithms, developed by the combinatorial
optimization community for problems with linear objectives. Finally, most of
the computational burden is pushed into the learning problem, which is solved
offline: only the linear surrogate (1.4) needs to be solved online. This is
helpful in industrial settings, where the full information on the instance is
often known at the last minute, which strongly constrains the computing time
and forces to use only objective functions simple enough to allow for fast
optimization algorithms.

###  1.2 Learning problem

Finding a suitable policy in $\mathcal{H}$ requires designing a learning
algorithm that selects a parameter ${\bm{w}}$ that leads to a good performance
in practice. Two typologies of learning strategies have been proposed in the
literature: supervised learning and regret minimization, a.k.a. imitation
learning and learning by experience (Bengio et al., 2021). Supervised learning
approaches assume that the training set contains instances
$X_{1},\ldots,X_{n}$ of (1.1) and their solutions
$\bar{\bm{y}}_{1},\ldots,\bar{\bm{y}}_{n}$. In this setting, the objective
$f^{0}$ is dismissed, and surrogate losses evaluate how far the model output
$\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))\in\mathcal{Y}({\bm{x}})$ of the model
is from the training set target $\bar{\bm{y}}_{i}$. Regret minimization
approaches instead only require instances $X_{1},\ldots,X_{n}$ and being able
to evaluate $f^{0}$, by focusing on minimizing the empirical regret

$\min_{{\bm{w}}\in\mathcal{W}}\mathcal{R}_{n}(h_{\bm{w}}),\quad\text{where}\quad\mathcal{R}_{n}(h_{\bm{w}})=\frac{1}{n}\sum_{i=1}^{n}f^{0}\big{(}h_{\bm{w}}(X_{i}),X_{i}\big{)}.$
(1.6)

Unfortunately, due to the combinatorial nature of (1.4), this objective is
piecewise constant, which has two negative impacts: it makes the optimization
problem intractable, and leads to poor generalization properties. To that
purpose, we rely on a perturbation strategy developed in Berthet et al. (2020)
to turn the linear optimization problem (1.4) into a family of tractable
distributions $p_{\lambda}({\bm{y}}|{\bm{\theta}})$, with $\lambda>0$ a
parameter controlling the smoothness of ${\bm{\theta}}\mapsto
p_{\lambda}({\bm{y}}|{\bm{\theta}})$. And we formulate the learning problem as

$\min_{{\bm{w}}\in\mathcal{W}}\Bigg{\\{}\frac{1}{n}\sum_{i=1}^{n}\sum_{{\bm{y}}\in\mathcal{Y}(X_{i})}p_{\lambda}({\bm{y}}|\psi_{{\bm{w}}}(X_{i}))\,f^{0}\big{(}{\bm{y}},X_{i}\big{)}\Bigg{\\}}\,.$
(1.7)

###  1.3 Related works

#### Relevant applications

A first natural application of our approach is _black-box optimization_ for
combinatorial problems, see for instance Vu et al. (2017). When modeling
complex phenomena, the cost function oracle $f^{0}({\bm{y}},{\bm{x}})$
typically involves running a simulation, which can be computationally
expensive. When the set of feasible solutions $\mathcal{Y}({\bm{x}})$ is
large, it is intractable to directly optimize the objective function $f^{0}$
in (1.1). It is then natural to replace the black-box $f^{0}$ by a tractable
surrogate in a parametric family. As a variant, $f^{0}$ may be sufficiently
simple to be optimized directly, but the algorithm solving (1.1) is too slow
or does not scale to large instances. In that setting, using a surrogate $f$
can be a way to _obtain a fast heuristic for large instances_. Parmentier and
T’Kindt (2021) obtained the best fast heuristic for large instances of a
classic single machine scheduling problem using such an approach. The need for
such algorithms is notably important in the context of _stochastic
optimization_ (Shapiro et al., 2021), where the objective
$f^{0}({\bm{y}},{\bm{x}})=\mathds{E}_{\xi}g({\bm{y}},\xi,{\bm{x}})$ involves
an expectation with respect to an uncertain parameter. There exists plenty of
algorithms to solve (1.1), but these are typically able to solve only
instances of moderate size, while algorithms for deterministic versions can
address large industrial instances. Such techniques have been applied to a
stochastic vehicle scheduling problem (Parmentier, 2021a). They can also be
used to approximate of two stage stochastic problems by single stage
deterministic surrogates (Dalle et al., 2022; Parmentier, 2021b). In that
context,
$f^{0}({\bm{y}},{\bm{x}})=\mathds{E}_{\xi,{\bm{z}}}g({\bm{y}},{\bm{z}},\xi,{\bm{x}})$,
where $\xi$ is again an uncertainty, ${\bm{y}}$ is a decision taken before
knowing $\xi$, and ${\bm{z}}$ is a decision taken after observing $\xi$. The
surrogate ${\psi}_{\bm{w}}({\bm{x}})^{\top}{\bm{y}}$ is then a deterministic
approximation of the function
$\min_{\bm{z}}\mathds{E}_{\xi,{\bm{z}}}g({\bm{y}},{\bm{z}},\xi,{\bm{x}})$.
Moreover, these methods are particularly well suited for _contextual
stochastic optimization problems_ (Sadana et al., 2024), see Example 1 below.
Finally, one context where such surrogates are very useful but which is beyond
the scope of this paper is _multistage stochastic optimization_ (Shapiro et
al., 2021) for combinatorial problems, as illustrated by the winning approach
for the 2022 EURO-NeurIPS challenge on dynamic vehicle routing (Baty et al.,
2024; EUR, ).

To illustrate our theoretical guarantees, we present three examples, one
generic and two more industrial.

######  Example 1.

(Contextual Stochastic Optimization) Consider a stochastic optimization
problem where the cost over a solution is influenced by some noise $\xi$,
(Sadana et al., 2024). The decision maker does not observe $\xi$, but has
access to some contextual data $\tilde{X}$ that is correlated to $\xi$. A
solution is a policy $h$ that maps any context $\tilde{\bm{x}}$ to a feasible
decision $\bm{y}\in\mathcal{Y}(\tilde{\bm{x}})$ given $\tilde{\bm{x}}$. Given
a joint distribution over $(\tilde{X},\xi)$ and a hypothesis class
$\mathcal{H}$ for $h$, the goal of contextual stochastic optimization is to
find the policy that minimizes the expected risk

$\min_{h\in\mathcal{H}}\mathds{E}\Big{[}f^{\mathrm{c}}\big{(}h(\tilde{X}),\tilde{X},\xi\big{)}\Big{]}$
(1.8)

An optimal policy is provided by the conditional distribution

$h^{\star}:\tilde{\bm{x}}\longmapsto\operatorname*{arg\,min}_{{\bm{y}}\in\mathcal{Y}(\tilde{\bm{x}})}\mathds{E}\Big{[}c\big{(}\delta(\tilde{X}),\tilde{X},\xi\big{)}\Big{|}\tilde{X}=\tilde{\bm{x}}\Big{]}$
(1.9)

However, contextual stochastic optimization considers the setting where the
joint distribution over $(\tilde{X},\xi)$ is unknown, but we have access to a
training set
$\tilde{\bm{x}}_{1},\bm{\xi}_{1},\ldots,{\bm{x}}_{n},\bm{\xi}_{n}$. This makes
the optimal policy (1.9) impractical. Sadana et al. (2024) survey this
emergent field. Contextual stochastic optimization can be considered as a
special case of our setting. It suffices to define $X$ as $(\tilde{X},\xi)$,
and $f^{0}(\bm{y},\bm{x})$ as $f^{\mathrm{c}}(\bm{y},\tilde{\bm{x}},\bm{\xi})$
for $\bm{x}=(\tilde{\bm{x}},\bm{\xi})$. We then restrict ourselves to
statistical models ${\psi}_{\bm{w}}(\bm{x})$ of the form
$\tilde{\psi}_{\bm{w}}(\tilde{\bm{x}})$ that exploit only the contextual
information and not the noise.

######  Example 2.

(Stochastic Vehicle Scheduling (Parmentier, 2021a)) Consider a set of tasks
with fixed begin and end times operated by a fleet of vehicles. When a task is
late, the assigned vehicle may propagate the delay to the next tasks. The goal
is to build the sequence of tasks operated by each vehicle in such a way that
each task is operated by exactly one vehicle and the expected cost of delay
$f^{0}({\bm{y}},{\bm{x}})$ is minimized. The instance ${\bm{x}}$ can typically
be encoded as a digraph $D=(V,A)$ where $V$ is the set of tasks, and there is
an arc $a$ in $A$ between $v$ and $v^{\prime}$ if these two tasks can be
operated by the same vehicle in a sequence, and a vector
$\bm{\phi}(a,{\bm{x}})$ of additional information for each arc $a$. Since a
solution corresponds to a partition of the vertices of an acyclic digraph into
paths, the set of feasible solutions $\mathcal{Y}({\bm{x}})$ can be identified
with the vertices of the flow polytope over $D$. As statistical model, one can
take $\psi_{\bm{w}}:\bm{x}\mapsto\bm{\theta}=(\theta_{a})_{a\in A}$ with
$\bm{\theta}_{a}=\bm{w}^{\top}\bm{\phi}(a,{\bm{x}})$ for each $a$. Aircraft
routing problems solved routinely by airlines are variants of this problem.
The stochastic version of it aims for its part at controlling delay
propagation. Solvers used in production typically solve the linear oracle.
Switching from the deterministic and delay-agnostic approaches in production
to our policy that avoids delay propagation does not impact the algorithm used
online in production. It amounts to changing the parametrization of the online
solver. Although, computing this new parametrization requires solving a
computationally intensive learning problem, this can be done offline.

######  Example 3.

(Single Machine Scheduling (Parmentier and T’Kindt, 2021)) In single machine
scheduling problems, $n$ jobs have to be performed on a given machine. A
solution is a permutation of $\\{1,\ldots,n\\}$ giving the order in which the
jobs are processed. Several objectives can be used. In the reference above,
jobs have release times, before which they cannot be started, and processing
times, which indicate the time needed by the machine to operate the job. The
vector $\bm{y}$ encodes a permutation, and the goal here is to minimize
$f^{0}(\bm{y},{\bm{x}})$ which is equal to the sum (the average) of the job
completion times, taking into account the release time. Permutations $\bm{y}$
can be encoded by vectors $\bm{y}_{i}$ giving the position of job $i$ in the
permutation. Optimizing a linear objective $\bm{\theta}^{\top}\bm{y}$ amounts
to sorting $\bm{\theta}$, and we can therefore use a sorting algorithm as
linear minimization oracle.

#### Learning algorithms and generalization bounds.

In order to train policies that embed a combinatorial optimization oracle
using a supervised learning approach, recent advances have focused on Fenchel-
Young losses (Blondel et al., 2020; Berthet et al., 2020; Dalle et al., 2022).
In contrast, the regret minimization paradigm has been overlooked (Parmentier,
2021b). As far as we know, the only well-established for regret is the “smart
predict then optimize” (Elmachtoub and Grigas, 2021), which relies on a
surrogate Structured Hinge Loss (Nowozin and Lampert, 2011, Chapter 6). It
applies only to the case of an $f^{0}$ that is linear but noisy. In this
paper, we focus on the theoretical analysis of the regret minimization
approach (1.6). One natural estimator consists in solving the empirical
learning problem. To the best of our knowledge, generalization bounds for
structured regret minimization have been studied in Parmentier (2021b), which
bounds the empirical process error on the regularized risk when parameters
${\bm{w}}$ are restricted to some $\ell_{\infty}$ ball. A similar result has
been shown in Qi et al. (2021) which also used Rademacher complexity to derive
the statistical excess risk in contextual stochastic optimization.

Our problem could also be seen as a special case of data-driven algorithm
design, in a setting introduced by Gupta and Roughgarden (2016). We refer to
Balcan (2021) for a summary of the literature on this topic. Given a training
set of instances of a problem coming from an unknown distribution, a regret
minimization approach is used to select the hyperparameters. Balcan et al.
(2021) consider generalization guarantees for the optimum of piecewise
decomposable risk such as our piecewise constant non-regularized risk (1.6).
This framework could apply in our context when $\psi_{w}$ is linear, but we do
not have an algorithm to optimize such a piecewise constant risk. Since the
regularized risk is smooth and therefore no more piecewise decomposable, their
result cannot be applied in our context. The counterpart of their result is
our Theorem 3 and relies on tools more appropriate to our smoothed context.

Such questions of generalization are also studied in the structured learning
literature (Ciliberto et al., 2016, 2020), where the problem consist in
finding the best function $h_{w}$ minimizing
$\mathds{E}_{X,Y}f^{0}(h_{w}(X),Y)$ for some probability distribution over
$(X,Y)$, given examples sampled from the same probability. The results are
very general, holding under mild assumptions on $f^{0}$, and have been
specialized to discrete spaces, in Nowak-Vila et al. (2019). In particular,
the results in the latter paper allow deriving automatically functions $h_{w}$
with a method that avoids the curse of dimensionality in the dimension of the
discrete space. However, it is unclear how to generalize this technique to an
arbitrary $f^{0}$ since the result depends on some internal properties of such
function that are not easy to check in a general case.

###  1.4 Contributions

In this paper, we build upon the normal fan geometry (see Section 2) involved
by the linear optimization program (1.4) and upon a perturbation approach to
deliver a generalization bound that controls the perturbation bias, the
statistical learning error, and the optimization error (Theorem 1).

Perturbation makes the learning problem smooth and amenable to gradient
descent methods. We choose to study the so-called kernel Sum-of-Squares method
(Rudi et al., 2024) in this paper. Furthermore, we prove non-asymptotic
guarantees on the performance of the learned policy, as well as convergence
rates toward the best policy encoded by (1.5).

To obtain our performance guarantees, we introduce the uniform weak (UW)
moment property. It corresponds to controlling uniformly over ${\bm{w}}$ the
moments of the distance of ${\psi}_{\bm{w}}(X)$ and the boundaries of the
cones of the normal fan of $\mathcal{Y}({\bm{x}})$. Hence, it quantifies the
interplay of the statistical model and (1.4). From this property, we prove an
explicit control of the variation of the risk when varying the scale
$\lambda>0$ of the perturbation. Furthermore, we show that the UW property
holds under mild assumptions on the statistical model $\psi_{w}$ on and on
$\mathds{P}_{X}$. We illustrate this on several examples in Section 3.2.

###  1.5 Outline

The next section (Section 2) presents our surrogate policy model (Section
2.1), the different kind of risks (Section 2.2) and the main theorem on
generalization guarantees (Section 2.3). Section 3 proves the main theorems.
It starts with a presentation the conditions required by our analysis (Section
3.1) and a discussion on the UW property (Section 3.2). In Section 3.3 we show
a control on the regularized risk under UW property. Section 3.4 focuses on
the discrepancy between the risk and its empirical version. Finally, we
provide in Section 3.5 the optimization error bound when using kernel Sum-of-
Squares.

###  1.6 Notation

We denote by $\mathcal{L}(U)$ the law of some random variable $U$ and by
$\mathds{P}_{U}$ (resp. $\mathds{E}_{U}$) the probability (resp. expectation)
with respect to the law $\mathcal{L}(U)$. We denote the stochastic boundedness
by $\mathcal{O}_{\mathds{P}}$: for a sequence of random variables
$(U_{n})_{n}$, notation $U_{n}=\mathcal{O}_{\mathds{P}}(1)$ means that for any
$\varepsilon>0$, there exists a finite $M>0$ and a finite $m\geq 0$ such that
for every $n\geq m$, $\mathds{P}(|U_{n}|>M)\leq\varepsilon$. For a sequence of
non-zero reals $(a_{n})_{n}$, $U_{n}=\mathcal{O}_{\mathds{P}}(a_{n})$ denotes
that $(U_{n}/a_{n})=\mathcal{O}_{\mathds{P}}(1)$.

The set of instances is denoted by $\mathcal{X}$, and the set of feasible
solutions for an instance ${\bm{x}}\in\mathcal{X}$ by
$\mathcal{Y}({\bm{x}})\subseteq\mathds{R}^{d({\bm{x}})}$. Given an instance
${\bm{x}}\in\mathcal{X}$, the objective function is $f^{0}(\cdot,{\bm{x}})$
whose optimal solution is denoted by ${\bm{y}}^{0}({\bm{x}})$ or simply
${\bm{y}}^{0}$. The oscillation of $f^{0}$ is defined by

$\operatorname{osc}(f^{0}):=\sup_{{\bm{y}},{\bm{x}}}\big{\\{}f^{0}({\bm{y}},{\bm{x}})\big{\\}}-\inf_{{\bm{y}},{\bm{x}}}\big{\\{}f^{0}({\bm{y}},{\bm{x}})\big{\\}}\,.$

The set of parameters is denoted by $\mathcal{W}$ and the statistical model by
${\psi}_{\bm{w}}\,:\,{\bm{x}}\in\mathcal{X}\mapsto\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$.
The solution given by parameter ${\bm{w}}\in\mathcal{W}$ is given by
$\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))$.

##  2 The surrogate policy model and its guarantees

###  2.1 Surrogate policy model

Statistical model ${\psi}_{\bm{w}}({\bm{x}})$CO oracle $\hat{\bm{y}}$
$\displaystyle\min_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}\langle{\bm{y}}|{\bm{\theta}}\rangle$Instance${\bm{x}}\in\mathcal{X}$Direction
vector$\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$Surrogate
Policy${\bm{y}}\in\mathcal{Y}({\bm{x}})$ Figure 1: Surrogate policy encoded by
the statistical model
${\psi}_{\bm{w}}\,:\,{\bm{x}}\in\mathcal{X}\mapsto\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$
with combinatorial optimization (CO) layer given by a linear program over
solutions ${\bm{y}}\in\mathcal{Y}$.

Figure 1 illustrates the policies we consider. We start with a statistical
model ${\bm{x}}\mapsto{\psi}_{\bm{w}}({\bm{x}})\in\mathds{R}^{d({\bm{x}})}$
that embeds an instance ${\bm{x}}\in\mathcal{X}$ into
$\mathds{R}^{d({\bm{x}})}$ and we set
$\bm{\theta}:={\psi}_{\bm{w}}({\bm{x}})\in\mathds{R}^{d({\bm{x}})}$. We build
our surrogate policy
$\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))\in\mathcal{Y}({\bm{x}})$ by
considering a solution to the linear program

$\max_{{\bm{y}}\in\mathcal{C}({\bm{x}})}\langle{\bm{y}},\bm{\theta}\rangle\text{
where }\bm{\theta}:={\psi}_{\bm{w}}({\bm{x}}).$ (2.1)

We minimize over the convex hull

$\mathcal{C}({\bm{x}}):=\operatorname*{conv}(\mathcal{Y}({\bm{x}})).$

which is a polytope since $\mathcal{Y}({\bm{x}})$ is finite. This surrogate
optimization problem is tractable in practice. It is given by the family of
surrogate maps
$f({\bm{y}},{\bm{x}};{\bm{w}})=\langle{\bm{y}},{\psi}_{\bm{w}}({\bm{x}})\rangle$
linear in $\bm{\theta}={\psi}_{\bm{w}}({\bm{x}})$.

By construction, the surrogate policy $\hat{\bm{y}}(\bm{\theta})$ is an
extreme point of $\mathcal{C}({\bm{x}})$. For each
${\bm{y}}\in\mathcal{Y}({\bm{x}})$ one can define the normal cone
$\mathcal{F}_{{\bm{y}}}$ as the set of directions $\bm{\theta}$ such that
${\bm{y}}$ is a solution to (2.1), equivalently

$\mathcal{F}_{{\bm{y}}}:=\Big{\\{}\bm{\theta}\in\mathds{R}^{d({\bm{x}})}\
\Big{|}\
\langle{\bm{y}},\bm{\theta}\rangle\geq\langle{\bm{y}}^{\prime},\bm{\theta}\rangle\text{
for all }{\bm{y}}^{\prime}\in\mathcal{C}({\bm{x}})\Big{\\}}$

The normal fan $\mathcal{N}(\mathcal{C}({\bm{x}}))$ is the collection of all
such cones, illustrated on Figure 2.

${\bm{y}}_{5}$${\bm{y}}_{6}$${\bm{y}}_{1}$${\bm{y}}_{2}$${\bm{y}}_{3}$${\bm{y}}_{4}$$\bm{\theta}$ |  | $\mathcal{F}_{{\bm{y}}_{1}}$$\mathcal{F}_{{\bm{y}}_{2}}$$\mathcal{F}_{{\bm{y}}_{3}}$$\mathcal{F}_{{\bm{y}}_{4}}$$\mathcal{F}_{{\bm{y}}_{5}}$$\mathcal{F}_{{\bm{y}}_{6}}$$\bm{\theta}$$\rho$
---|---|---
Figure 2: Normal cone at point ${\bm{y}}_{1}$ to the polytope (left) and
normal fan with internal radius $\rho$ at point $\bm{\theta}$ (right).

######  Example 2 (continued).

In the context of stochastic vehicle scheduling, $\mathcal{Y}(\bm{x})$
corresponds to the partitions of the vertices $V$ into paths. Solutions are
encoded by indicator vectors $\bm{y}=(y_{a})_{a\in A}$ where $y_{a}$ is a
binary variable equal to $1$ if arc $a$ is in the solution, and to $0$
otherwise. Its convex hull, $\mathcal{C}(\bm{x})$, is known to be a flow
polytope and has a compact formulation.

######  Example 3 (continued).

The feasible solutions $\mathcal{Y}(\bm{x})$ of single machine scheduling are
the vectors $\bm{y}$ encoding a permutation of $\\{1,\ldots,n\\}$. Their
convex hull $\mathcal{C}(\bm{x})$ is a polytope known as the permutahedron.

For a given ${\bm{x}}\in\mathcal{X}$, it can occur that several
${\bm{y}}\in\cal Y$ are optimal, i.e., there exists several
${\bm{y}}\in\mathcal{Y}({\bm{x}})$ such that
${\psi}_{\bm{w}}({\bm{x}})\in\mathcal{F}_{\bm{y}}$. In this case,
$\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))$ is defined as a probability
distribution on such ${\bm{y}}$’s. A formal definition of the surrogate policy
$\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))$ is given as follows.

######  Definition 1 (Surrogate policy).

Let $\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$, and define the internal cone
radius as

$\rho(\bm{\theta}):=\sup\Big{\\{}r>0\ \Big{|}\ \forall{\bm{u}}\in
B(0,1)\,,\exists{\bm{y}}\in\mathcal{Y}\text{ s.t.
}\bm{\theta}\in\mathcal{F}_{y}\text{ and
}\bm{\theta}+r{\bm{u}}\in\mathcal{F}_{y}\Big{\\}}\,,$

where $B(0,1)\subset\mathds{R}^{d({\bm{x}})}$ is the standard Euclidean ball.
We define the surrogate policy measure

$\hat{\bm{y}}(\bm{\theta}):=\sum_{{\bm{y}}\in\mathcal{Y}}p_{0}({\bm{y}}|\bm{\theta})\delta_{\bm{y}}$

where $\delta_{\bm{y}}$ is the Dirac mass at point ${\bm{y}}$ and

$p_{0}({\bm{y}}|\bm{\theta}):=\mathds{P}_{U}\big{[}\bm{\theta}+\rho(\bm{\theta})U\in\mathcal{F}_{\bm{y}}\big{]}$

where $U$ is uniformly distributed on $B(0,1)$.

Note that $p_{0}({\bm{y}}|\bm{\theta})$ measures the proportion of normal cone
$\mathcal{F}_{\bm{y}}$ locally around $\bm{\theta}$.

######  Remark 2 (Generic case and an abuse of notation).

Generically, a direction $\bm{\theta}$ is in the interior of a normal cone
$\mathcal{F}_{\bm{y}}$ and ${\bm{y}}^{\prime}\mapsto
p_{0}({\bm{y}}^{\prime}|\bm{\theta})$ is zero except at point ${\bm{y}}$. In
this case, the measure $\hat{\bm{y}}(\bm{\theta})$ is a Dirac mass at point
${\bm{y}}$ and, by an abuse of notation, we denote by
$\hat{\bm{y}}(\bm{\theta})$ the vector ${\bm{y}}\in\mathds{R}^{d({\bm{y}})}$.

To ease notation, we will consider that, by convention,

$f(\hat{\bm{y}}(\bm{\theta})):=\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{0}({\bm{y}}|\bm{\theta})f({\bm{y}})$

for all functions $f$ and all $\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$.

###  2.2 The regularized and empirical risks

In our algorithms, we optimize an empirical and regularized version of the
risk $\mathcal{R}$, which we now introduce. We observe $n$ instances
$X_{1},\ldots,X_{n}$ drawn independently according to a (possibly unknown)
distribution $\mathds{P}_{X}$. For all ${\bm{w}}\in\mathcal{W}$, we denote the
policy by
$h_{{\bm{w}}}\,:\,{\bm{x}}\mapsto\hat{\bm{y}}(\psi_{{\bm{w}}}({\bm{x}}))$ and
introduce the regularized population risk $\mathcal{R}_{t}$ and its empirical
version $\mathcal{R}_{n,t}$ by

$\displaystyle\mathcal{R}_{t}(h_{\bm{w}})$
$\displaystyle:=\mathds{E}_{X,Z}\big{[}f^{0}(\hat{\bm{y}}({\psi}_{\bm{w}}(X)+tZ(X)),X)\big{]}\,,$
(2.2) $\displaystyle\mathcal{R}_{n,t}(h_{\bm{w}})$
$\displaystyle:=\frac{1}{n}\sum_{i=1}^{n}\mathds{E}_{Z}\Big{\\{}\big{[}f^{0}(\hat{\bm{y}}({\psi}_{\bm{w}}(X_{i})+tZ(X_{i})),X_{i})\big{]}\Big{\\}}\,,$
(2.3)

where $t\geq 0$ is some regularization parameter and $Z(\cdot)$ is a
perturbation defined as follows.

######  Definition 3 (Law of the perturbation).

Given ${\bm{x}}\in\mathcal{X}$, $Z({\bm{x}})$ has the same law as $R\times U$
where the random vector $U$ is uniformly distributed on the Euclidean sphere
of $\mathds{R}^{d({\bm{x}})}$ and the random positive scalar $R\in(0,\infty)$
is independent of $U$.

######  Remark 4.

The perturbation $Z({\bm{x}})$ is a random vector in
$\mathds{R}^{d({\bm{x}})}$ that is uniformly distributed on the Euclidean
sphere of $\mathds{R}^{d({\bm{x}})}$ and is scaled by a random positive scalar
$R$. As such, the perturbed vector ${\psi}_{\bm{w}}({\bm{x}})+tZ({\bm{x}})$
belongs to the normal cone of ${\psi}_{\bm{w}}({\bm{x}})$ whenever
$tR\leq\rho({\psi}_{\bm{w}}({\bm{x}}))$, where
$\rho({\psi}_{\bm{w}}({\bm{x}}))$ is the internal radius of the normal cone at
${\psi}_{\bm{w}}({\bm{x}})$. Note that this latter condition does not depend
on the direction $U(G)$ of the perturbation $Z({\bm{x}})$.

For later use in our analysis we introduce the following quantity.

######  Definition 5 (Perturbed surrogate policy probabilities).

Let $\lambda>0$ and ${\bm{x}}\in\mathcal{X}$. Let
$\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$ and
${\bm{y}}\in\mathcal{Y}({\bm{x}})$. We define the perturbed surrogate policy
probability at solution point ${\bm{y}}$ for the direction $\bm{\theta}$ as

$p_{\lambda}({\bm{y}}|\bm{\theta}):=\mathds{E}_{Z}[p_{0}({\bm{y}}|\bm{\theta}+\lambda
Z({\bm{x}}))]\,.$

The (non-perturbed) _surrogate policy_ of ${\bm{x}}\in\mathcal{X}$ is defined
as a solution
$\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}}))\in\mathds{R}^{d({\bm{x}})}$ to (1.4),
and the surrogate policy model is defined by (1.5).

Now, let ${\varepsilon_{0}}\geq 0$ and consider the risk

$\displaystyle\mathcal{R}_{\varepsilon_{0}}(\mathcal{H})$
$\displaystyle:=\min_{{\bm{w}}\in\mathcal{W}}\big{\\{}\mathds{E}_{X,Z}\big{[}f^{0}(\hat{\bm{y}}({\psi}_{\bm{w}}(X)+{\varepsilon_{0}}Z(X)),X)\big{]}\big{\\}}\,,$
(2.4)
$\displaystyle=\min_{{\bm{w}}\in\mathcal{W}}\big{\\{}\mathds{E}_{X,Z}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{0}({\bm{y}}|{\psi}_{\bm{w}}(X)+{\varepsilon_{0}}Z(X))\,f^{0}({\bm{y}},X)\big{]}\big{\\}}\,,$
$\displaystyle=\min_{{\bm{w}}\in\mathcal{W}}\big{\\{}\mathds{E}_{X}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{\varepsilon_{0}}({\bm{y}}|{\psi}_{\bm{w}}(X))\,f^{0}({\bm{y}},X)\big{]}\big{\\}}\,,$
$\displaystyle=\inf_{{\bm{w}}\in\mathcal{W}}\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})\,,$
$\displaystyle=\mathcal{R}_{\varepsilon_{0}}(h_{{\bm{w}}^{\star}})\,.$ (2.5)

We aim at finding the optimal surrogate policy
$h_{{\bm{w}}^{\star}}\,:\,{\bm{x}}\mapsto\hat{\bm{y}}(\psi_{{\bm{w}}^{\star}}({\bm{x}}))$,
where ${\bm{w}}^{\star}$ is a solution to (2.4), assumed to exist.

######  Remark 6.

(cases $\varepsilon_{0}=0$ and $\varepsilon_{0}>0$) In the case where
$\varepsilon_{0}=0$, note that
$\mathcal{R}_{0}(h_{{\bm{w}}})=\mathds{E}_{X}\big{[}f^{0}(\hat{\bm{y}}({\psi}_{\bm{w}}(X)),X)\big{]}$
and it may fail to be continuous in ${\bm{w}}$. Hence it may not be clear that
a solution to (2.4) exists.

In the case $\varepsilon_{0}>0$, we will prove that the risk is continuous and
hence a solution always exists, under mild assumptions (Proposition 1).
Furthermore, we will also prove that the Property
$($$\mathrm{UW}_{\varepsilon_{0}}$$)$, introduced later, holds, under mild
assumptions (Proposition 2). So rather than studying a perturbation scaling
$\lambda$ such that $\lambda\in(0,+\infty)$, we suggest instead to fix
$\varepsilon_{0}>0$ and focus on $\lambda\in[\varepsilon_{0},+\infty)$. The
small constant $\varepsilon_{0}$ should be then understood as a minimal amount
of regularization to introduce, playing the role of machine precision.

###  2.3 Guarantees

#### Policy model misspecification error.

We define the policy model misspecification error as

$\displaystyle\mathcal{E}_{\varepsilon_{0}}(\mathcal{H})$
$\displaystyle:=\mathcal{R}_{\varepsilon_{0}}(\mathcal{H})-\mathds{E}_{X}\big{[}f^{0}({\bm{y}}^{0}(X),X)\big{]}$
$\displaystyle=\mathds{E}_{X,Z}\big{[}f^{0}(\hat{\bm{y}}(\psi_{{\bm{w}}^{\star}}(X)+{\varepsilon_{0}}Z(X)),X)\big{]}-\mathds{E}_{X}\big{[}f^{0}({\bm{y}}^{0}(X),X)\big{]}\geq
0$

where ${\bm{y}}^{0}({\bm{x}})$ minimizes (1.1).

In $\mathds{P}_{X}$-average, it measures, for the loss given by $f^{0}$, how
the policy $\hat{\bm{y}}(\psi_{{\bm{w}}^{\star}}(X)+\varepsilon_{0}Z(X))$
fails to approximate ${\bm{y}}^{0}(X)$ or equivalently how
$\psi_{{\bm{w}}^{\star}}(X)+\varepsilon_{0}Z(X)$ fails to point towards normal
cones $\mathcal{F}_{{\bm{y}}}$ at point ${\bm{y}}$ such that
$f^{0}({\bm{y}},X)-f^{0}({\bm{y}}^{0}(X),X)$ is small. To apprehend the term
$\mathcal{E}_{\varepsilon_{0}}(\mathcal{H})$, one needs some regularity-
related notion for ${\bm{y}}^{0}$ and $f^{0}$. This study is out of the scope
of this paper, and we focus on the optimal risk defined in (2.4).

#### Error bounds

We will introduce an estimator $h_{{\bm{w}}_{M,n,\lambda}}$ which depends on
some complexity integer parameter $M\geq 1$, the sample size $n$, and a real
valued tuning parameter $\lambda>0$ such that $\lambda\geq\varepsilon_{0}$.
Assume that ${\bm{w}}^{\star}$ exists, then one has the following
decomposition of the error, with each term handled separately,

$\displaystyle 0\leq\;$
$\displaystyle\mathcal{R}_{\varepsilon_{0}}(h_{{\bm{w}}_{M,n,\lambda}})-\mathds{E}_{X}\big{[}f^{0}({\bm{y}}^{0}(X),X)\big{]}$
$\displaystyle=\underbrace{\mathcal{R}_{\varepsilon_{0}}(h_{{\bm{w}}_{M,n,\lambda}})-\mathcal{R}_{\lambda}(h_{{\bm{w}}_{M,n,\lambda}})}_{\text{Theorem\leavevmode\nobreak\
\ref{thm:reg_excess}}}+\underbrace{\mathcal{R}_{\lambda}(h_{{\bm{w}}_{M,n,\lambda}})-\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{M,n,\lambda}})}_{\text{Theorem\leavevmode\nobreak\
\ref{thm:A3_statistical}}}$
$\displaystyle+\underbrace{\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{M,n,\lambda}})-\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{n,\lambda}})}_{\text{Theorem\leavevmode\nobreak\
\ref{thm:kSoS_certificate}
}}+\underbrace{\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{n,\lambda}})-\mathcal{R}_{n,\lambda}(h_{{\bm{w}}^{\star}})}_{\leq
0}$
$\displaystyle+\underbrace{\mathcal{R}_{n,\lambda}(h_{{\bm{w}}^{\star}})-\mathcal{R}_{\lambda}(h_{{\bm{w}}^{\star}})}_{\text{Theorem\leavevmode\nobreak\
\ref{thm:A3_statistical}}}+\underbrace{\mathcal{R}_{\lambda}(h_{{\bm{w}}^{\star}})-\mathcal{R}_{\varepsilon_{0}}(h_{{\bm{w}}^{\star}})}_{\text{Theorem\leavevmode\nobreak\
\ref{thm:reg_excess}}}$
$\displaystyle+\mathcal{E}_{\varepsilon_{0}}(\mathcal{H})$

where, for any ${\bm{w}}\in\mathcal{W}$, we define $h_{{\bm{w}}_{n,\lambda}}$
such that

$\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{n,\lambda}})=\min_{{\bm{w}}\in\mathcal{W}}\mathcal{R}_{n,\lambda}(h_{{\bm{w}}})\,.$

In practice, we will make some working assumptions, to obtain the following,
possibly probabilistic, bounds on each of the terms.

######  Theorem 1.

Under the conditions given in Section 3.1, the following holds true. Let
$\varepsilon_{0}\geq 0$ and $\lambda>0$ be such that
$\lambda\geq\varepsilon_{0}$. Let $\tau\in(0,1)$. There exists a constant
$C>0$ that depends only on $\varepsilon$, $\tau$ and $f^{0}$ such that for any
${\bm{w}}\in\mathcal{W}$ and $n\geq 1$, one has

$\displaystyle|\mathcal{R}_{\varepsilon_{0}}(h_{{\bm{w}}})-\mathcal{R}_{\lambda}(h_{{\bm{w}}})|$
$\displaystyle=C\lambda^{\tau}\mathrm{polylog}(\lambda)$ (Theorem 2,
Perturbation bias)
$\displaystyle|\mathcal{R}_{\lambda}(h_{{\bm{w}}})-\mathcal{R}_{n,\lambda}(h_{{\bm{w}}})|$
$\displaystyle=\mathcal{O}_{\mathds{P}}\Big{(}\frac{1}{\lambda\sqrt{n}}\Big{)}$
(Theorem 3, Empirical process) and, for $h_{{\bm{w}}_{n,\lambda}}$ given by
the kernel Sum-of-Squares estimate solution to
$\displaystyle\min_{{\bm{w}}\in\mathcal{W}}\mathcal{R}_{n,\lambda}$,
$\displaystyle|\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{M,n,\lambda}})-\mathcal{R}_{n,\lambda}(h_{{\bm{w}}_{n,\lambda}})|$
$\displaystyle=\mathcal{O}_{\mathds{P}}\Bigg{[}\bigg{(}\frac{1}{\lambda(M/\log
M)^{\frac{1}{d}}}\bigg{)}^{s-\frac{d}{2}}\Bigg{]}$ by Theorem 4 (K-SoS)

where $\mathrm{polylog}(\lambda)$ is a polynomial logarithm term and $s>d/2$
is some tuning parameter on the order of regularity of the admissible
functions. It yields the following bound on the error between the risk and the
regularized risk

$0\leq\mathcal{R}_{\varepsilon_{0}}(h_{{\bm{w}}_{M,n,\lambda}})-\mathds{E}_{X}\big{[}f^{0}({\bm{y}}^{0}(X),X)\big{]}\leq
C\lambda^{\tau}\mathrm{polylog}(\lambda)+\mathcal{O}_{\mathds{P}}\Bigg{[}\Big{(}\frac{1}{\lambda\sqrt{n}}\Big{)}+\bigg{(}\frac{1}{\lambda(M/\log
M)^{\frac{1}{d}}}\bigg{)}^{s-\frac{d}{2}}\Bigg{]}+\mathcal{E}_{\varepsilon_{0}}(\mathcal{H})\,,$

where the K-SoS bound is given up to a small a posteriori error, see Remark
11.

Our excess risk bound is of the type statistical to computational showing the
effect of regularisation $\lambda$, sampling $n$ and the complexity of the
algorithm $M$. It shows that, when $\lambda$ goes to zero, the excess risk
goes to zero if $\sqrt{n}$ dominates $1/\lambda$ and $M$ dominates
$1/\lambda^{d}$. The first shows quantitatively how many sample points the
practitioner needs when diminishing the regularization $\lambda$. Our
statistical bound matches the parametric rate $1/\sqrt{n}$ and $1/\lambda$ can
be interpreted as a standard error. The latter is a so-called curse of
dimension, the complexity of the algorithm is exponential in the dimension $d$
and $1/\lambda$ can be interpreted as a smoothness parameter.

######  Example 1 (continued).

In our reduction of contextual stochastic optimization to the setting of this
article, the definition of $X$ as $(\tilde{X},\xi)$ and of
$f^{0}(\bm{y},\bm{x})$ as $f^{\mathrm{c}}(\bm{y},\tilde{\bm{x}},\bm{\xi})$ may
have seemed surprising. The common practice in stochastic optimization would
have suggested defining $X$ as $\tilde{X}$ and $f^{0}(\bm{y},\bm{x})$ as
$\mathds{E}\big{[}f^{\mathrm{c}}(\bm{y},X,\xi)\big{|}X=\bm{x}\big{]}$. Our
reduction has however two advantages.

First, we do not know the distribution over $(\tilde{X},\xi)$, but only have
access to a training set
$\tilde{\bm{x}}_{1},\bm{\xi}_{1},\ldots,\tilde{\bm{x}}_{n},\bm{\xi}_{n}$.
Hence, we have all the data needed to formulate the empirical risk
minimization problem with our reduction, while the alternative reduction would
require the conditional distribution on $\xi$ given $\tilde{X}={\bm{x}}_{i}$,
which we do not know.

Second, with our reduction, Theorem 1 applies out of the box to contextual
stochastic optimization. In stochastic optimization, we select a decision by
coordinating the different scenarios for a given instance/context. In our
case, we select instead a policy by coordinating scenarios across different
instances/context. Using a machine learning terminology, our reduction is a
discriminative approach, while the other one, which would be closer to the
common practice in stochastic optimization, is a generative approach. Theorem
1 provides as theoretical basis for this discriminative approach.

##  3 Main results

We present in this section the results we use to obtain Theorem 1. As
discussed above, our analysis involves a parameter $\varepsilon_{0}\geq 0$
that accounts for a small perturbation defining the risk
$\mathcal{R}_{\varepsilon_{0}}$. Our main result is an upper bound on the
excess risk when Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$, defined
later, is met. We show that this property is always satisfied when
$\varepsilon_{0}>0$. In the case $\varepsilon_{0}=0$ there is no
regularization on the risk, but one has to check Property
$($$\mathrm{UW}_{\varepsilon_{0}}$$)$ to apply our result. We will give four
cases in which Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$ is met when
$\varepsilon_{0}=0$.

###  3.1 The conditions

#### Partition of the instances

On many combinatorial optimization problems, instances ${\bm{x}}$ are pairs
$(G,\xi)$, where $G$ is a combinatorial object, and $\xi_{G}$ is a vector of
parameters whose dimensions depend only on $G$, e.g., $G$ is typically a graph
and $\xi_{G}$ contains a vector of parameters (a label) for each vertex or
edge of $G$. Furthermore, we know an embedding of $\mathcal{Y}({\bm{x}})$ in
the Euclidean space $\mathds{R}^{d({\bm{x}})}$ with $d({\bm{x}})$ of
polynomial size in ${\bm{x}}$. For instance, solutions of the shortest path
problems are the vertices of the path polytope, solutions of the minimum
weight spanning tree problem are the vertices of the spanning tree polytope,
and so on. More generally, most combinatorial optimization problems admit a
mixed integer linear programming formulation. As illustrated by these
examples, the set of feasible solutions $\mathcal{Y}({\bm{x}})$ frequently
depends only on $G$. To avoid technicalities, we omit $G$ and $\xi_{G}$ in the
analysis but make the following assumption about a partition of the space of
instances.

######  Assumption $($Part$)$.

The following holds.

  * •

There is a finite partition $\mathcal{G}$ of $\mathcal{X}$ into
$(\mathcal{X}_{G})_{G\in\mathcal{G}}$ such that $\mathcal{Y}({\bm{x}})$ are
constant and finite for all ${\bm{x}}\in\mathcal{X}_{G}$;

  * •

The absolute value of the target function $|f^{0}|$ is uniformly bounded on
$\mathcal{X}_{G}$ for each $G\in\mathcal{G}$ and, hence, its oscillation is
finite, $\operatorname{osc}(f^{0})<\infty$, with

$\operatorname{osc}(f^{0}):=\sup f^{0}-\inf f^{0}\,.$

######  Example 2 (continued).

In the context of the stochastic vehicle scheduling, $G$ corresponds to the
digraph $D$ we introduced, and $\mathcal{G}$ to all possible digraphs"?
Assuming that $\mathcal{G}$ is finite is therefore natural: since no firm is
going to serve an infinite number of requests, $\mathcal{G}$ is a subset of
the set of acyclic digraphs with at most $n^{\mathrm{v}}$ vertices for some
$n^{\mathrm{v}}$. Since we are working with the flow polytope, $d(G)$
corresponds to the number of arcs $|A|$ in digraph $D$. The second part of the
assumption holds if costs are positive and if there is a maximum cost that can
be encountered per task, which is natural in any real life setting.

######  Example 3 (continued).

In the context of single machine scheduling, the set of feasible solutions
depends only on the number $n$ of tasks. As a consequence, $G$ can be
identified with $n$, and $\mathcal{G}$ is the subset of $\mathds{Z}_{+}$
corresponding to instance sizes. Supposing that it is finite only amounts to
having an upper bound on the number of tasks. Finally, $d(G)=n$.

######  Remark 7.

Under Assumption $($Part$)$, note that, for each $G\in\mathcal{G}$, and all
${\bm{x}}\in\mathcal{X}_{G}$, the embedding of $\mathcal{Y}({\bm{x}})$ into
$\mathds{R}^{d({\bm{x}})}$ and the dimension ${\bm{x}}\to d({\bm{x}})$ are
constant. We denote by $d(G)$ the value of $d({\bm{x}})$ for
${\bm{x}}\in\mathcal{X}_{G}$. Without loss of generality, we assume that
$\mathcal{Y}({\bm{x}})\subseteq\mathds{R}^{d({\bm{x}})}$ and
$\mathrm{Vect}(\mathcal{Y}({\bm{x}}))=\mathds{R}^{d({\bm{x}})}$.

Under Assumption $($Part$)$, we have the following decomposition of the law
$\mathds{P}_{X}$ of $X$.

$\mathds{E}_{X}[t(X)]:=\int_{\cal
X}t({\bm{x}})\mathrm{d}\mathds{P}_{X}({\bm{x}})=\sum_{G\in\mathcal{G}}\mathds{P}[X\in\mathcal{X}_{G}]\int_{\mathcal{X}_{G}}t({\bm{x}})\mathrm{d}\mathds{P}_{\xi_{G}}({\bm{x}})\,,$
(3.1)

where $\mathds{P}_{\xi_{G}}$ is the law of $\xi_{G}$, namely the conditional
law of $\mathcal{L}(X|{X\in\mathcal{X}_{G}})$, and $t$ is any bounded
continuous function. Denote by $\mathds{P}_{X}^{(n)}$ its empirical version,
given by

$\mathds{E}_{X}^{(n)}[t(X)]:=\frac{1}{n}\sum_{i\in[n]}t(X_{i})=\frac{1}{n}\sum_{i\in[n]}\sum_{G\in\mathcal{G}}\mathds{1}_{\\{X_{i}\in\mathcal{X}_{G}\\}}t(X_{i}).$
(3.2)

#### The Gaussian perturbation assumption:

Throughout this paper, we consider that the perturbation random field
$Z(\cdot)$ is a Gaussian random field satisfying the following assumption.

######  Assumption $($Gauss$)$.

We assume that $\sqrt{d(G)}R(G)$ has chi-distribution with $d(G)$ degrees of
freedom which is equivalent to

$\forall{\bm{x}}\in\mathcal{X}_{G}\,,\quad\sqrt{d(G)}Z({\bm{x}})\sim\mathcal{N}(0,\mathrm{Id}_{d(G)})\,.$
(3.3)

Under Assumptions $($Part$)$ and $($Gauss$)$, we have the following
decomposition of the law $\mathcal{L}(X,Z)$,

$\mathds{E}_{X,Z}\big{[}t(X,Z)\big{]}:=\sum_{G\in\mathcal{G}}\mathds{P}[X\in\mathcal{X}_{G}]\int_{\mathcal{X}_{G}}\Big{(}\int_{\mathds{R}^{d(G)}}t({\bm{x}},{\bm{z}})\varphi_{d(G)}({\bm{z}})\mathrm{d}{\bm{z}}\Big{)}\mathrm{d}\mathds{P}_{\xi_{G}}({\bm{x}})\,,$
(3.4)

where $\mathds{P}_{\xi_{G}}$ is the law of $\xi_{G}$, namely the conditional
law of $\mathcal{L}(X|{X\in\mathcal{X}_{G}})$, $\varphi_{d(G)}$ is the
Gaussian density of $\mathcal{N}(0,(1/d(G))\mathrm{Id}_{d(G)})$, and
$t(\cdot)$ is any bounded continuous function.

#### Lipschitz feature map:

The following assumption is about the Lipschitz continuity of the feature map.

######  Assumption $($Lip$)$.

For all ${\bm{x}}\in\mathcal{X}$, the function
${\bm{w}}\in\mathcal{W}\mapsto{\psi}_{\bm{w}}({\bm{x}})\in\mathds{R}^{d({\bm{x}})}$
is $L_{\mathcal{W}}$-Lipschitz continuous with a constant $L_{\mathcal{W}}$
which does not depend on ${\bm{x}}$.

######  Example 2 (continued).

In the stochastic vehicle scheduling problem example, we use a generalized
linear model with $\theta_{a}=\bm{w}^{\top}\bm{\phi}_{a}$, which
$\bm{\phi}_{a}$ being some features of the connection $a$, related to the
slack between the tasks that are the tails and the head of $a$. In
applications, it is natural to assume that this slack is bounded, which gives
our Lipschitz assumption.

Under these assumptions

######  Remark 8.

When ${\bm{w}}$ encodes the parameters of a neural network with smooth
activation functions one can compute the differential with respect to
${\bm{w}}$ at input point ${\bm{x}}$ of neural network, namely the gradient of
${\bm{w}}\mapsto{\psi}_{\bm{w}}({\bm{x}})$. By the chain rule, one obtain that
partial gradients $\partial{\psi}_{\bm{w}}({\bm{x}})/\partial w_{i}$ are
product of derivatives of the activation functions evaluated at layers values,
matrices of some of the parameters ${\bm{w}}$, and input vector ${\bm{x}}$.
Since the layers values are themselves composition of activation functions and
product of matrices of some parameters ${\bm{w}}$ and input vector ${\bm{x}}$,
they are bounded by compactness of $\mathcal{W}$ and $\mathcal{X}$. Hence,
partial gradients $\partial{\psi}_{\bm{w}}({\bm{x}})/\partial w_{i}$ are also
bounded by compactness. We deduce that $($Lip$)$ is satisfied.

#### Existence of solutions:

We assume that the optimization problem (2.4) has a solution and that the
minimum of the risk has a minimizer.

######  Assumption $($Sol$)$.

There exists ${\bm{w}}^{\star}\in\mathcal{W}$ solution to (2.4) and
$\displaystyle\min_{{\bm{w}}\in\mathcal{W}}\mathcal{R}_{n,\varepsilon_{0}}(h_{{\bm{w}}})$
has a minimizer.

######  Proposition 1.

Let $\varepsilon_{0}$ be positive and let Assumptions $($Gauss$)$ and
$($Lip$)$ hold. Then, Assumption $($Sol$)$ is met.

###### Proof.

When $\varepsilon_{0}>0$, the risk
${\bm{w}}\mapsto\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})$ is continuous with
respect to ${\bm{w}}$ (see Lemma 7) and the set $\mathcal{W}$ is compact.
Therefore, the infimum in (2.4) is attained. ∎

#### The uniform weak (UW) moment property:

We introduce a property that is useful for the proofs, and use it to discuss
the impact of the perturbation on the policy.

Given a regularization parameter $\lambda>\varepsilon_{0}\geq 0$, we will see
later that the uniform weak moment property, given below, implies that
${\bm{w}}\mapsto|\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})|$
is uniformly bounded by $\mathcal{O}(\lambda^{\tau}\mathrm{polylog}(\lambda))$
where $\mathrm{polylog}(\lambda)$ is a polynomial logarithm term and
$\tau\in(0,1)$, see Theorem 2.

######  Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$.

For all $\tau\in(0,1)$, it exists a positive constant
$C_{\varepsilon_{0},\tau}>0$ such that

$\forall{\bm{w}}\in\mathcal{W}\,,\
\mathds{E}_{X,Z}\bigg{[}\Big{(}\frac{\rho(\psi_{{\bm{w}}}(X)+\varepsilon_{0}Z(X))}{\sqrt{d(X)}}\Big{)}^{-\tau}\bigg{]}\leq
C_{\varepsilon_{0},\tau}\,.$ (3.5)

###  3.2 Discussion on the UW property

#### UW is always met for positive regularization:

We have the following proposition showing that, when $\varepsilon_{0}>0$,
Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$ holds with a constant
$C_{\varepsilon_{0},\tau}$ that depends only on $\varepsilon_{0}$ and $\tau$.

######  Proposition 2.

If $\varepsilon_{0}>0$ then, under Assumptions $($Part$)$ and $($Gauss$)$,
Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$ holds and

$C_{\varepsilon_{0},\tau}\leq\big{(}\int_{\mathds{R}}|t|^{-\tau}\varphi(t)\mathrm{d}t\big{)}\times\varepsilon_{0}^{-\tau}\times\sum_{G\in\mathcal{G}}\mathds{P}[X\in\mathcal{X}_{G}]\,|\mathcal{Y}(G)|^{2}\,d(G)^{\tau}\,\,,$

where $\varphi$ is the standard Gaussian density on $\mathds{R}$.

###### Proof.

In view of (3.4) and without loss of generality, we consider that
$X\in\mathcal{X}_{G}$ almost surely, for $G$ fixed. Let $H$ be an hyperplane
and $d(\theta,H)$ be the Euclidean distance between $\theta$ and $H$. By
(3.4), note that

$\displaystyle\mathds{E}_{X,Z}\bigg{[}\Big{(}\frac{d(\psi_{{\bm{w}}}(X)+\varepsilon_{0}Z(X))}{\sqrt{d(X)}}\Big{)}^{-\tau}\bigg{]}$
$\displaystyle=\mathds{E}_{\xi_{G},Z}\bigg{[}\Big{(}\frac{d(G)^{\frac{\tau}{2}}}{d(\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G}),H)^{\tau}}\Big{)}\bigg{]}$
$\displaystyle=d(G)^{\frac{\tau}{2}}\,\mathds{E}_{\xi_{G}}\bigg{[}\mathds{E}_{Z}\Big{(}\frac{1}{d(\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G}),H)^{\tau}}\Big{|}\xi_{G}\Big{)}\bigg{]}\,.$

Conditionally to $\xi_{G}$, observe that
${d(\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G}),H)^{\tau}}$ is the
distance between the hyperplane $H$ and a Gaussian with mean
$\psi_{{\bm{w}}}(\xi_{G})$ and variance
$(\varepsilon_{0}^{2}/d(G))\mathrm{Id}_{d(G)}$. This distance is

${d(\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G}),H)^{\tau}}=|\langle
n_{H}\,,\,\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G})\rangle|^{\tau}$

where $n_{H}$ is a unitary orthogonal vector to $H$. Conditionally to
$\xi_{G}$, $\langle
n_{H}\,,\,\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G})\rangle$ is a
Gaussian random variable with mean $\langle
n_{H}\,,\,\psi_{{\bm{w}}}(\xi_{G})\rangle$ and variance
$\varepsilon_{0}^{2}/d(G)$.

Observe that the level sets of $t\mapsto|t|^{-\tau}$ are symmetric convex
bodies, namely intervals of the form $[-a,a]$ for $a>0$. By Anderson’s lemma
(Giné and Nickl, 2021, Theorem 2.4.4), it holds that

$\mathds{E}_{Z}\Big{(}\frac{1}{d(\psi_{{\bm{w}}}(\xi_{G})+\varepsilon_{0}Z(\xi_{G}),H)^{\tau}}\Big{|}\xi_{G}\Big{)}\leq\mathds{E}_{Z}\Big{(}\frac{1}{d(\varepsilon_{0}Z(\xi_{G}),H)^{\tau}}\Big{|}\xi_{G}\Big{)}=\frac{d(G)^{\frac{\tau}{2}}}{\varepsilon_{0}^{\tau}}\int_{\mathds{R}}|t|^{-\tau}\varphi(t)\mathrm{d}t\,.$

We deduce that

$\mathds{E}_{X,Z}\bigg{[}\Big{(}\frac{d(\psi_{{\bm{w}}}(X)+\varepsilon_{0}Z(X))}{\sqrt{d(X)}}\Big{)}^{-\tau}\bigg{]}\leq\frac{d(G)^{\tau}}{\varepsilon_{0}^{\tau}}\int_{\mathds{R}}|t|^{-\tau}\varphi(t)\mathrm{d}t<\infty$

Then observe that for $Q$ a polyhedral cone with non-empty interior we have
that $Q$ is the intersection of halfspaces delimited by $H_{1},\ldots,H_{D}$.
Note also that
$\frac{1}{d(\theta,Q)^{\tau}}\leq\frac{1}{d(\theta,H_{1})^{\tau}}+\ldots+\frac{1}{d(\theta,H_{D})^{\tau}}$,
and thus $\theta\mapsto\frac{1}{d(\theta,Q)^{\tau}}$ is integrable with
respect to the law of $\psi_{{\bm{w}}}(X)+\varepsilon_{0}Z(X)$ conditional to
$X\in\mathcal{X}_{G}$. Finally, observe that $\theta\mapsto\rho(\theta)$ is
the distance to the boundaries of polyhedral cones. ∎

#### The no perturbation case:

We would like to show that Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$ is
met when $\varepsilon_{0}=0$ in some cases. Under Assumption $($Part$)$ and
without loss of generality, we consider that $d(X)$ is constant almost surely,
say $d(X)=d$. We consider four scenarios on the family of laws

$\mathds{L}:=\big{\\{}\mathcal{L}(\psi_{{\bm{w}}}(X))\ :\
{\bm{w}}\in\mathcal{W}\big{\\}}$

where $\mathcal{L}(\psi_{{\bm{w}}}(X))$ denotes the law of
$\psi_{{\bm{w}}}(X)$, given by

  * (i)

There exists a compact $K\subset\mathds{R}^{d}$ such that the supports of the
laws of $\mathds{L}$ are included in $K$ and these laws are absolutely
continuous, with respect to the Lebesgue measure, with uniformly bounded
density;

  * (ii)

There exists a finite set $K\subset\mathds{R}^{d}$ such that the laws of
$\mathds{L}$ are finite with support included in $K$;

  * (iii)

There exists a compact Riemannian submanifold $K\subset\mathds{R}^{d}$ (with
metric induced by the Euclidean metric) such that the supports of the laws of
$\mathds{L}$ are included in $K$ and these laws are absolutely continuous,
with respect to the Riemannian measure, with uniformly bounded density;

  * (iv)

The laws of $\mathds{L}$ are a mixture of an absolutely continuous law
satisfying (i), a discrete law satisfying (ii), and a singular continuous law
satisfying (iii).

Now, using that

  * •

$\theta\mapsto\rho(\theta)$ is the distance to the boundaries of polyhedral
cones;

  * •

$\frac{1}{d(\theta,Q)^{\tau}}\leq\frac{1}{d(\theta,H_{1})^{\tau}}+\ldots+\frac{1}{d(\theta,H_{D})^{\tau}}$
where $Q$ a polyhedral cone wich is the intersection of halfspaces delimited
by $H_{1},\ldots,H_{D}$;

  * •

$\theta\mapsto\frac{1}{d(\theta,Q)^{\tau}}$ is integrable with respect to the
uniform measure on compacts, hence with respect to any bounded density with
compact support;

we deduce that Property $($$\mathrm{UW}_{\varepsilon_{0}}$$)$ holds for the
four scenarios above.

######  Remark 9.

Let us first consider a _generalized linear model_ , such as in our stochastic
vehicle scheduling of Example 2, where $\bm{\theta}=\bm{\Phi}(x)^{\top}\bm{w}$
for a given matrix of features. Then Case (i) is satisfied when every feature
is continuous and upper bounded and the number of features is rich enough. If
the number of features is small, then the image of the instances through the
model may be included in a hyperplane, which leads to Case (iii). Case (ii)
corresponds to features that take integer values and are finite. Using a
submanifold instead of a hyperplane, Case (iii) enables us to deal with neural
networks in $\psi_{\bm{w}}$.

######  Example 2 (continued).

In the stochastic vehicle scheduling problem example, continuous features
would typically correspond to features reflecting properties of the delay or
the distance, which are by nature continuous and upper bounded. Discrete
features can for instance indicate combinatorial properties such as the number
of arcs departing from the tail of an arc $a$.

###  3.3 Regularization by perturbation and perturbation bias

We start by introducing

$V_{\bm{w}}(\lambda):=\mathds{E}_{X}\mathds{P}_{Z}\Big{[}\|Z(X)\|_{2}>\frac{\rho(\psi_{{\bm{w}}}(X))}{\lambda}\,\big{|}\,X\Big{]}=\mathds{E}_{X}\mathds{P}_{R}\Big{[}{R(X)}>\frac{\rho(\psi_{{\bm{w}}}(X))}{\lambda}\,\big{|}\,X\Big{]}\,,$

for $\lambda>0$ and ${\bm{w}}\in\mathcal{W}$. By dominated convergence, note
that $V_{\bm{w}}(\lambda)$ tends to zero as $\lambda$ goes to zero and that
$V_{\bm{w}}(\lambda)$ tends to one as $\lambda$ goes to infinity. Note also
that the function $\lambda\mapsto V_{\bm{w}}(\lambda)\in[0,1]$ is non-
decreasing.

Our first contribution is a bound on the risk perturbation.

######  Proposition 3.

Let $\varepsilon_{0}\geq 0$ and $\lambda>0$ be such that
$\lambda\geq\varepsilon_{0}$. It holds that, for all ${\bm{w}}\in\mathcal{W}$,

$\big{|}\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}(h_{\bm{w}})\big{|}\leq
2\operatorname{osc}(f^{0})V_{\bm{w}}(\lambda)\,\quad\text{and}\quad\big{|}\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})\big{|}\leq
4\operatorname{osc}(f^{0})V_{\bm{w}}(\lambda)\,.$

###### Proof.

Recall that
$p_{0}({\bm{y}}|\psi_{{\bm{w}}}({\bm{x}}))=\mathds{P}[\psi_{{\bm{w}}}({\bm{x}})+\rho(\psi_{{\bm{w}}}({\bm{x}}))U\in\mathcal{F}_{\bm{y}}]$
where $U$ is uniformly on $B(0,1)$ and

$p_{\lambda}({\bm{y}}|\psi_{{\bm{w}}}({\bm{x}}))=\mathds{P}_{Z}[\psi_{{\bm{w}}}({\bm{x}})/\lambda+Z({\bm{x}})\in\mathcal{F}_{\bm{y}}]=\mathds{P}_{Z}[\psi_{{\bm{w}}}({\bm{x}})/\lambda+R({\bm{x}})U({\bm{x}})\in\mathcal{F}_{\bm{y}}]\,.$

Observe that
$p_{0}({\bm{y}}|\psi_{{\bm{w}}}({\bm{x}}))=\mathds{P}[\psi_{{\bm{w}}}({\bm{x}})+\rho(\psi_{{\bm{w}}}({\bm{x}}))U\in\mathcal{F}_{\bm{y}}]=\mathds{P}_{U}[\psi_{{\bm{w}}}({\bm{x}})+\rho(\psi_{{\bm{w}}}({\bm{x}}))U({\bm{x}})\in\mathcal{F}_{\bm{y}}]$
where $U({\bm{x}})$ is uniformly distributed on the Euclidean sphere
$\mathds{R}^{d({\bm{x}})}$.

Note that

$\displaystyle\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}(h_{\bm{w}})$
$\displaystyle=\mathds{E}_{X}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}p_{\lambda}({\bm{y}}|\psi_{{\bm{w}}}(X))f^{0}({\bm{y}},X)-\sum_{{\bm{y}}\in\mathcal{Y}(X)}p_{0}({\bm{y}}|\psi_{{\bm{w}}}(X))f^{0}({\bm{y}},X)\big{]}$
$\displaystyle=\mathds{E}_{X,R,U}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}\underbrace{\Big{(}\mathbbm{1}\big{(}\psi_{{\bm{w}}}({\bm{x}})/\lambda+R({\bm{x}})U({\bm{x}})\in\mathcal{F}_{\bm{y}}\big{)}-\mathbbm{1}\big{(}\psi_{{\bm{w}}}({\bm{x}})+\rho(\psi_{{\bm{w}}}({\bm{x}}))U({\bm{x}})\in\mathcal{F}_{\bm{y}}\big{)}\Big{)}}_{A(X,{\bm{y}},R,U)}f^{0}({\bm{y}},X)\big{]}$
$\displaystyle=\mathds{E}_{X,R,U}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}\Big{(}\mathbbm{1}\big{(}R(x)\leq\rho(\psi_{{\bm{w}}}(x))/\lambda\big{)}+\mathbbm{1}\big{(}R(x)>\rho(\psi_{{\bm{w}}}(x))/\lambda\big{)}\Big{)}A(X,{\bm{y}},R,U)f^{0}({\bm{y}},X)\big{]}$

where we have introduced the random variable $A(X,{\bm{y}},R,U)$. By a zero-
sum argument, one can replace $f^{0}({\bm{y}},X)$ by $f^{0}({\bm{y}},X)-\inf
f^{0}$ in the calculation above.

By definition of $\rho$, if $R(x)\leq\rho(\psi_{{\bm{w}}}(x))/\lambda$ then
$\mathbbm{1}\big{(}\psi_{{\bm{w}}}(x)/\lambda+R(x)U(x)\in\mathcal{F}_{\bm{y}}\big{)}$
is equal to
$\mathbbm{1}\big{(}\psi_{{\bm{w}}}({\bm{x}})+\rho(\psi_{{\bm{w}}}({\bm{x}}))U\in\mathcal{F}_{\bm{y}}\big{)}$.
Hence
$A(X,{\bm{y}},R,U)\mathbbm{1}\big{(}R(x)\leq\rho(\psi_{{\bm{w}}}(x))/\lambda\big{)}=0$.

We therefore have

$\displaystyle\big{|}\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}(h_{\bm{w}})\big{|}$
$\displaystyle=\big{|}\mathds{E}_{X,R,U}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}\mathbbm{1}\big{(}R(X)>\rho(\psi_{{\bm{w}}}(X))/\lambda\big{)}A(X,{\bm{y}},R,U)(f^{0}({\bm{y}},X)-\inf
f^{0})\big{]}\big{|}$
$\displaystyle\leq\mathds{E}_{X,R,U}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}\mathbbm{1}\big{(}R(x)>\rho(\psi_{{\bm{w}}}(x))/\lambda\big{)}|A(X,{\bm{y}},R,U)|(f^{0}({\bm{y}},X)-\inf
f^{0})\big{]}$
$\displaystyle\leq\mathds{E}_{X,R,U}\big{[}\operatorname{osc}(f^{0})\mathbbm{1}\big{(}R(X)>\rho(\psi_{{\bm{w}}}(X))/\lambda\big{)}\underbrace{\sum_{{\bm{y}}\in\mathcal{Y}(X)}|A(X,{\bm{y}},R,U)|}_{\leq
2}\big{]}$ $\displaystyle\leq
2\operatorname{osc}(f^{0})\mathds{E}_{X}\mathds{P}_{R}\big{(}R(X)>\rho(\psi_{{\bm{w}}}(X))/\lambda|X\big{)}=2\operatorname{osc}(f^{0})V_{\bm{w}}(\lambda).$

Using that the function $\lambda\mapsto V_{\bm{w}}(\lambda)\in[0,1]$ is non-
decreasing, we have

$\displaystyle\big{|}\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})\big{|}$
$\displaystyle\leq\big{|}\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}(h_{\bm{w}})\big{|}+\big{|}\mathcal{R}(h_{\bm{w}})-\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})\big{|}$
$\displaystyle\leq
2\operatorname{osc}(f^{0})V_{\bm{w}}(\lambda)+2\operatorname{osc}(f^{0})V_{\bm{w}}(\varepsilon_{0})$
$\displaystyle\leq 4\operatorname{osc}(f^{0})V_{\bm{w}}(\lambda)\,,$

which concludes the proof. ∎

######  Corollary 4.

Under Condition $($Sol$)$, it holds that

$\mathcal{R}_{\lambda}(\mathcal{H})-\mathcal{R}(\mathcal{H})\leq
2\operatorname{osc}(f^{0})V_{{\bm{w}}^{\star}}(\lambda)\,.$

###### Proof.

Observe that
$\mathcal{R}(\mathcal{H})=\mathds{E}_{X}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}p_{0}({\bm{y}}|\psi_{{\bm{w}}^{\star}}(X))f^{0}({\bm{y}},X)\big{]}$
and remark that it holds

$\mathcal{R}_{\lambda}(\mathcal{H})\leq\mathds{E}_{X}\big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X)}p_{\lambda}({\bm{y}}|\psi_{{\bm{w}}^{\star}}(X))f^{0}({\bm{y}},X)\big{]}\,,$

which gives the result by substracting the first equality. ∎

Combining Proposition 3 and Proposition 10 we get the following theorem.

######  Theorem 2.

Under $($$\mathrm{UW}_{\varepsilon_{0}}$$)$ and $($Gauss$)$, the following
holds true. Let $\varepsilon_{0}\geq 0$ and $\lambda>0$ be such that
$\lambda\geq\varepsilon_{0}$. Let $\tau\in(0,1)$. There exists a constant
$C>0$ $($that may depend on the constant $C_{\varepsilon_{0},\tau}$, appearing
in $($$\mathrm{UW}_{\varepsilon_{0}}$$)$, and $\operatorname{osc}(f^{0}))$
such that, for all ${\bm{w}}\in\mathcal{W}$,

$\big{|}\mathcal{R}_{\lambda}(h_{\bm{w}})-\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})\big{|}\leq
C\,\lambda^{\tau}\mathrm{polylog}(\lambda)\,,$

where $\mathrm{polylog}(\lambda)$ is a polynomial logarithm term.

###### Proof.

Assumption (Mom.) is implied by Assumption
$($$\mathrm{UW}_{\varepsilon_{0}}$$)$, see Proposition 3 and Proposition 10 of
Appendix B. ∎

###  3.4 Empirical process control

We consider the empirical regularized risk (2.3) and

$\mathcal{R}_{n,\lambda}(\mathcal{H}):=\min_{{\bm{w}}\in\mathcal{W}}\Big{\\{}\frac{1}{n}\sum_{i=1}^{n}\mathds{E}_{Z}\big{[}f^{0}\big{(}\hat{\bm{y}}({\psi}_{\bm{w}}(X_{i})+\lambda
Z(X_{i})),X_{i}\big{)}\big{]}\Big{\\}}\,.$ (3.6)

We would like to control the random variable

$\Delta_{n}=\sup_{{\bm{w}}\in\mathcal{W}}\big{|}\mathcal{R}_{n,\lambda}({\bm{w}})-\mathcal{R}_{\lambda}({\bm{w}})\big{|}\,.$

We are going to use a Bernstein inequality for bounded sub-Gaussian random
variables. In this analysis, we use the oscillation of $f^{0}$, recalled in
Assumption $($Part$)$, and the Dudley’s entropy integral, defined as

$\displaystyle\mathcal{I}_{\mathcal{W}}$
$\displaystyle:=\int_{0}^{\infty}\\!\\!\sqrt{\log\mathcal{N}(\mathcal{W},\|\cdot\|,\varepsilon)}\,\,\mathrm{d}\varepsilon\,,$
(3.7)

where $\mathcal{N}(\mathcal{W},\|\cdot\|,\varepsilon)$ is
$\varepsilon$-covering number of $\mathcal{W}$ with respect to $\|\cdot\|$,
for further details the reader may refer to (Wainwright, 2019, Chapter 5).

######  Remark 10 (Example).

Let $\mathds{B}$ be the unit ball of $\|\cdot\|$ and assume that there exists
$R>0$ such that $\mathcal{W}\subset R\mathds{B}$. The volume ratios lemma
(Wainwright, 2019, Lemma 5.7) states that

$\mathcal{I}_{\mathcal{W}}\leq Cd\log(R+1/R)\,,$

where $C>0$ is a universal constant.

######  Theorem 3.

Let Assumptions $($Lip$)$ and $($Gauss$)$ hold. Assume that $\displaystyle
d(\mathcal{X}):=\max_{{\bm{x}}\in\mathcal{X}}d({\bm{x}})<\infty$ $($which is
implied by $($Part$)$$)$. Then, for all $\delta\in(0,1)$, it holds that

$\Delta_{n}=\sup_{{\bm{w}}\in\mathcal{W}}\big{|}\mathcal{R}_{n,\lambda}({\bm{w}})-\mathcal{R}_{\lambda}({\bm{w}})\big{|}\leq\frac{\operatorname{osc}(f^{0})}{\lambda\sqrt{n}}\Bigg{(}{(\ln{2})^{\nicefrac{{-3}}{{4}}}}\,L_{\mathcal{W}}\,\mathcal{I}_{\mathcal{W}}\,\sqrt{d(\mathcal{X})}+4\sqrt{\ln{\dfrac{8}{\delta}}}\Bigg{)}\,,$
(3.8)

with probability higher than $1-\delta$.

###### Proof.

Bounding $\Delta_{n}$ is related to studying the zero-mean random variables

$V_{i,{\bm{w}}}:=\frac{1}{n}\Big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X_{i})}p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}(X_{i}))f^{0}({\bm{y}},X_{i})-\mathds{E}_{(X,Y)}\big{[}p_{\lambda}(Y|{\psi}_{\bm{w}}(X))f^{0}(Y,X)\big{]}\Big{]}\,,$
(3.9)

which depend on the random samples $X_{i}\sim\mathds{P}_{X}$. Indeed
$\Delta_{n}=\sup_{{\bm{w}}\in\mathcal{W}}\big{|}\sum_{i\in[n]}V_{i,w}\big{|}$
for a given draw. Rather than studying the absolute value, we focus on
$Z:=\sup_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}V_{i,{\bm{w}}}$ and use a
symmetrization argument. Since
$\Delta_{n}=\max(\sup_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}V_{i,{\bm{w}}},-\inf_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}V_{i,{\bm{w}}})$,
focusing on the first term, the second one being dealt similarly, we will show
that $Z$ is close to its mean $\mathds{E}Z$ with high probability (through a
Bernstein inequality) and then bound the mean $\mathds{E}Z$ by
$\mathcal{O}(\nicefrac{{1}}{{\sqrt{n}}})$ (through Hoeffding and Dudley’s
inequalities).

Let $(X^{\prime}_{i})\sim\mathds{P}_{X}^{\otimes n}$ be drawn independently
from $(X_{i})_{i\in[n]}$ and define $V^{\prime}_{i,{\bm{w}}}$ as in Equation
(3.9). Then $(V_{i,{\bm{w}}})_{i\in[n]}$ and
$(V^{\prime}_{i,{\bm{w}}})_{i\in[n]}$ are independent with the same
distribution.

######  Lemma 5.

We have that almost surely

$\sup_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}(V_{i,{\bm{w}}}-V^{\prime}_{i,{\bm{w}}})^{2}\leq\frac{(\operatorname{osc}(f^{0}))^{2}}{n}\,.$
(3.10)

###### Proof.

We use the fact that $p_{\lambda}\in[0,1]$ to obtain that

$\displaystyle V_{i,{\bm{w}}}-V^{\prime}_{i,{\bm{w}}}$
$\displaystyle=\frac{1}{n}\Big{[}\sum_{{\bm{y}}\in\mathcal{Y}(X_{i})}p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}(X_{i}))f^{0}({\bm{y}},X_{i})-\sum_{{\bm{y}}\in\mathcal{Y}(X^{\prime}_{i})}p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}(X^{\prime}_{i}))f^{0}(y,x^{\prime}_{i})\big{]}\Big{]}$
$\displaystyle\leq\frac{1}{n}\Big{[}\sup_{{\bm{x}},{\bm{x}}^{\prime}\in\mathcal{X},\,{\bm{y}},{\bm{y}}^{\prime}\in\mathcal{Y}}f^{0}({\bm{y}},{\bm{x}})-f^{0}({\bm{y}}^{\prime},{\bm{x}}^{\prime})\big{]}\Big{]}$
$\displaystyle=\frac{\operatorname{osc}(f^{0})}{n}\,,$

Since $(V^{\prime}_{i,{\bm{w}}}-V_{i,{\bm{w}}})$ is equal in law to
$(V_{i,{\bm{w}}}-V^{\prime}_{i,{\bm{w}}})$, this yields (3.10). ∎

Owing to Lemma 5, we can thus apply Boucheron et al. (2013, Theorem 12.3, Page
333) since

$\mathds{E}_{X^{\prime}_{i}\sim\mathds{P}_{X}}[\sup_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}(V_{i,{\bm{w}}}-V^{\prime}_{i,{\bm{w}}})^{2}]\leq\frac{\operatorname{osc}(f^{0})}{n}\,.$

Consequently, for all $t\geq 0$, we have

$\mathds{P}[Z\geq\mathds{E}Z+t]\leq
4e^{-\frac{nt^{2}}{16(\operatorname{osc}(f^{0}))^{2}}}.$ (3.11)

We are now going to show that
$\mathds{E}Z=\mathcal{O}(\nicefrac{{1}}{{\sqrt{n}}})$ through a chaining
argument. Fix ${\bm{w}},{\bm{w}}^{\prime}\in\mathcal{W}$. By Lemma 9, note
that
${\bm{w}}\mapsto\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}({\bm{x}}))$
is Lipschitz with constant $L_{\mathcal{W}}\sqrt{d({\bm{x}})}/\lambda$. We
obtain that for all $i\in[n]$, as
$\frac{1}{n}\sum_{{\bm{y}}\in\mathcal{Y}(X_{i})}p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}(X_{i}))=1$,

$\displaystyle V_{i,{\bm{w}}}-V_{i,{\bm{w}}^{\prime}}$
$\displaystyle=\frac{1}{n}\sum_{{\bm{y}}\in\mathcal{Y}(X_{i})}(p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}(X_{i}))-p_{\lambda}({\bm{y}}|\psi_{{\bm{w}}^{\prime}}(X_{i}))f^{0}({\bm{y}},X_{i})$
$\displaystyle=\frac{1}{n}\sum_{{\bm{y}}\in\mathcal{Y}(X_{i})}(p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}(X_{i}))-p_{\lambda}({\bm{y}}|\psi_{{\bm{w}}^{\prime}}(X_{i}))(f^{0}({\bm{y}},X_{i})-\inf_{{\bm{x}}^{\prime},{\bm{y}}^{\prime}}f^{0}({\bm{y}}^{\prime},{\bm{x}}^{\prime}))$
$\displaystyle\leq\frac{1}{n}L_{\mathcal{W}}\frac{\sqrt{d(\mathcal{X})}}{\lambda}\operatorname{osc}(f^{0})\|{\bm{w}}-{\bm{w}}^{\prime}\|=:a_{n}\,,$

and similarly $V_{i,{\bm{w}}^{\prime}}-V_{i,{\bm{w}}}\leq a_{n}$. For
$V_{\bm{w}}:=\sum_{i\in[n]}V_{i,{\bm{w}}}$, we deduce that the random variable
$(V_{\bm{w}}-V_{{\bm{w}}^{\prime}})$ is sub-Gaussian as a sum of zero-mean
independent sub-Gaussians (Vershynin, 2018, Proposition 2.6.1, Page 29), and
it holds that for the associated $\psi_{2}$-norm (Vershynin, 2018, Example
2.5.8, Page 28),

$\|V_{\bm{w}}-V_{{\bm{w}}^{\prime}}\|_{\psi_{2}}\leq\Big{[}\frac{1}{\sqrt{\ln{2}}}\sum_{i\in[n]}\|V_{i,{\bm{w}}}-V_{i,{\bm{w}}^{\prime}}\|^{2}_{\psi_{2}}\Big{]}^{\nicefrac{{1}}{{2}}}\\!\\!\leq\frac{a_{n}\sqrt{n}}{(\ln{2})^{\nicefrac{{1}}{{4}}}}=\frac{L_{\mathcal{W}}\sqrt{d(\mathcal{X})}}{(\ln{2})^{\nicefrac{{1}}{{4}}}\lambda\sqrt{n}}\operatorname{osc}(f^{0})\|{\bm{w}}-{\bm{w}}^{\prime}\|\,,$

with, for a sub-Gaussian $X$, the $\psi_{2}$-norm is defined as
$\|X\|_{\psi_{2}}:=\inf\\{t\geq
0\,|\,\mathds{E}e^{\nicefrac{{X^{2}}}{{t^{2}}}}\leq 2\\}$. In other words, we
have shown that the zero-mean random variable $V_{\bm{w}}$ has sub-Gaussian
increments. We can thus apply Dudley’s inequality (Vershynin, 2018, Theorem
8.1.3, Page 188),

$\displaystyle\mathds{E}V=\mathds{E}[\sup_{{\bm{w}}\in\mathcal{W}}V_{\bm{w}}]$
$\displaystyle\leq\frac{L_{\mathcal{W}}\sqrt{d(\mathcal{X})}}{(\ln{2})^{\nicefrac{{3}}{{4}}}\lambda\sqrt{n}}\operatorname{osc}(f^{0})\,\mathcal{I}_{\mathcal{W}}\,.$

We proceed similarly with
$\inf_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}V_{i,{\bm{w}}}$. We obtain that,
setting
$t=\sqrt{\frac{16(\operatorname{osc}(f^{0}))^{2}}{n}\ln{\frac{4}{\delta}}}$ in
(3.11), with probability $1-2\delta$,

$\Delta_{n}=\max\big{\\{}\sup_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}V_{i,{\bm{w}}},-\inf_{{\bm{w}}\in\mathcal{W}}\sum_{i\in[n]}V_{i,{\bm{w}}}\big{\\}}\leq\mathds{E}V+t$

Using the bound on $\mathds{E}V$, this yields (3.8). ∎

###  3.5 Optimization error bound on the Kernel-SoS case

The last part of the proof consists in bounding the optimization error when
minimizing $\mathcal{R}_{n,\lambda}$. Up to our knowledge, many global
optimization algorithms converge because they provably do an exhaustive
search, as is the case for DIRECT for instance (Jones et al., 1993, Section
5). Consequently, their convergence bounds are most often affected by the
curse of dimensionality. Here we have an advantage since
$\mathcal{R}_{n,\lambda}$ is a smooth function thanks to the parameter
$\lambda>0$. This incites us to use a specific kernel-based approach recently
suggested in Rudi et al. (2024) to leverage the smoothness of the function
$\mathcal{R}_{n,\lambda}$. To alleviate notation, we set

$R_{n,\lambda}({\bm{w}}):=\mathcal{R}_{n,\lambda}(h_{{\bm{w}}}).$

The key idea will be to find an (approximate) Sum-of-Squares representation
for $R_{n,\lambda}$ while minimizing $R_{n,\lambda}$. Fix $s>1+d/2$ and let
$\mathcal{H}_{\phi}$ be the Sobolev space over $\mathcal{W}$ of smoothness
$s$. Denote by $k$ its reproducing kernel, i.e.
$k(\cdot,{\bm{w}})\in\mathcal{H}_{\phi}$ such that for all
$g\in\mathcal{H}_{\phi},{\bm{w}}\in\mathcal{W}$, we have
$g({\bm{w}})=\left\langle{g},{k(\cdot,{\bm{w}})}\right\rangle_{\mathcal{H}_{\phi}}$.
Let $\phi:\mathcal{W}\to\mathcal{H}_{\phi}$ be the kernel embedding, i.e.
$\phi({\bm{w}})=k(\cdot,{\bm{w}})$. Set $\tilde{s}=s-d/2\geq 1$. Let
$S^{+}(\mathcal{H}_{\phi})$ be the set of self-adjoint positive semidefinite
operators from $\mathcal{H}_{\phi}$ to $\mathcal{H}_{\phi}$.

######  Assumption $($kSoS$)$.

We assume that $\mathcal{W}$ is the closure of a union of closed balls
$\cup_{{\bm{w}}\in\mathcal{S}}\mathds{B}({\bm{w}},r)$ for some $r>0$ and a
bounded set $\mathcal{S}\subset\mathds{R}^{d}$. Assume furthermore that
$R_{n,\lambda}\in C^{s+3}(\mathcal{W},\mathds{R})$ and that $R_{n,\lambda}$
has a kernel Sum-of-Square representation, i.e. there exists $A_{n,\lambda}\in
S^{+}(\mathcal{H}_{\phi})$ such that

$R_{n,\lambda}({\bm{w}})-\min_{{\bm{w}}^{\prime}\in\mathcal{W}}R_{n,\lambda}({\bm{w}}^{\prime})=\left\langle{\phi({\bm{w}})},{A_{n,\lambda}\phi({\bm{w}})}\right\rangle_{\mathcal{H}_{\phi}},\,\forall{\bm{w}}\in\mathcal{W}.$
(3.12)

For the following, we fix the $A_{n,\lambda}$ with minimal trace satisfying
(3.12).

The requirement on $\mathcal{W}$ is to use scattering bounds and is quite mild
for compact sets. We mostly ensure $R_{n,\lambda}\in
C^{s+3}(\mathcal{W},\mathds{R})$ by regularizing $R_{n,0}$ through $\lambda>0$
and using a smooth feature map ${\bm{w}}\mapsto{\psi}_{\bm{w}}$, e.g.
${\psi}_{\bm{w}}({\bm{x}})=\Phi_{x}{\bm{w}}$ with
$\Phi_{x}\in\mathds{R}^{d(x)\times d}$. By Rudi et al. (2024, Theorem 3) (see
also Marteau-Ferey et al. (2024)), the existence of $A_{n,\lambda}$ is for
instance guaranteed if $R_{n,\lambda}$ and its second order partial
derivatives belong to $\mathcal{H}_{\phi}$, and $R_{n,\lambda}$ has a finite
number of minimizers, all with positive Hessian. We will see an explicit
example below.

Recall that a global optimization problem can always be cast as a convex one
under the form

$\displaystyle\underset{\displaystyle\begin{subarray}{c}c\in\mathds{R},\end{subarray}}{\mathrm{max}}\quad
c$ $\displaystyle c$ $\displaystyle\leq
R_{n,\lambda}({\bm{w}}),\,\forall\,{\bm{w}}\in\mathcal{W}.$

We consider the kernel Sum-Of-Squares approximation of it, namely

$\displaystyle\underset{\displaystyle\begin{subarray}{c}c\in\mathds{R},\\\
A\in
S^{+}(\mathcal{H}_{\phi})\end{subarray}}{\mathrm{min}}\quad-c+\lambda_{\phi}\operatorname{Tr}(A)$$\displaystyle{R_{n,\lambda}(\tilde{{\bm{w}}}_{m})-c}{=\left\langle{\phi(\tilde{{\bm{w}}}_{m})},{A\phi(\tilde{{\bm{w}}}_{m})}\right\rangle_{\mathcal{H}_{\phi}},\,\forall\,m\in[M]}$
(3.13a)

where $\lambda_{\phi}>0$ is some parameter, the
$(\tilde{{\bm{w}}}_{m})_{m\in[M]}$ are sampled uniformly at random in
$\mathcal{W}$.

######  Remark 11.

A few words are in order. We want a $c\in\mathds{R}$ such that $w\mapsto
R_{n,\lambda}({\bm{w}})-c\geq 0$. On the other hand, for any $A\in
S^{+}(\mathcal{H}_{\phi})$, the function
$w\mapsto\left\langle{\phi({\bm{w}})},{A\phi({\bm{w}})}\right\rangle_{\mathcal{H}_{\phi}}$
is nonnegative, and under some assumptions on $\phi$ and $\mathcal{H}_{\phi}$
(Marteau-Ferey et al., 2020), these functions are pointwise dense in the
nonnegative continuous functions. So we look for one such $A$ for which the
two functions are equal at least at the sampled points and we penalize the
extra variable $A$ by its trace. For (3.13), the Newton-based algorithm
detailed in Rudi et al. (2024, Section 6) outputs the infimum value estimate
$\hat{R}=\hat{c}$, the corresponding $\hat{A}$ and a candidate
$\hat{\bm{w}}_{n,\lambda}=\sum_{m\in[M]}\alpha_{m}\tilde{{\bm{w}}}_{m}$ where
each $\alpha_{m}$ is the Lagrange multiplier of the equality constraint
(3.13a) at $\tilde{{\bm{w}}}_{m}$. Contrary to Rudi et al. (2024), an extra
term $|R_{n,\lambda}(\hat{\bm{w}}_{n,\lambda})-\hat{R}|$ appears in our
analysis. One expects that this error goes to zero as $M$ grows large, and we
referred to it as a posteriori error.

######  Theorem 4.

Fix $\delta\in(0,1)$. With the above notation, under Assumption $($kSoS$)$,
there exists explicit constants $\bar{M}\in\mathds{N}$ and
$\bar{C}\in\mathds{R}_{+}$ depending on
$s,d,\delta,r,\operatorname{diam}(\mathcal{W})$ such that for $M\geq\bar{M}$
and
$\lambda_{\phi}\geq\bar{C}M^{-\tilde{s}/d}\left(\log\frac{M}{\delta}\right)^{\tilde{s}/d}$
we have, with probability at least $1-\delta$,

$|R_{n,\lambda}(\hat{\bm{w}}_{n,\lambda})-R_{n,\lambda}({\bm{w}}_{n,\lambda})|\leq
R_{n,\lambda}(\hat{\bm{w}}_{n,\lambda})-\hat{R}+\lambda_{\phi}\left(\operatorname{Tr}\left(A_{n,\lambda}\right)+|R_{n,\lambda}|_{\mathcal{W},\lceil\tilde{s}\rceil}\right),$
(3.14)

depending on the Sobolev norm
$|R_{n,\lambda}|_{\mathcal{W},\lceil\tilde{s}\rceil}=\sup_{|\alpha|=\tilde{s}}\sup_{{\bm{w}}^{\prime}\in\mathcal{W}}|D^{\alpha}R_{n,\lambda}({\bm{w}}^{\prime})|$.

The estimate (3.14) has the advantage that if $s\gg d$ then the smoothness
mitigates the curse of dimensionality. Indeed, to achieve an a priori error
$\epsilon$ on the last term of (3.14), which we denote
$\lambda_{\phi}C_{n,\lambda}$, we can take $M\propto\epsilon^{-d/\tilde{s}}$
rather than the more usual, and much larger, $M\propto\epsilon^{-d}$. Note
however that the constants $\bar{M}$ and $\bar{C}$ can still depend
exponentially in the dimension, their closed-form expressions can be found in
Rudi et al. (2024, Theorem 5). Furthermore the Sobolev norm
$|R_{n,\lambda}|_{\mathcal{W},\lceil s-d/2\rceil}$ and the cost of the Sum-of-
Square representation $\operatorname{Tr}\left(A_{n,\lambda}\right)$ depend
crucially on the noise parameter $\lambda$ since the latter plays a role on
the regularity of $R_{n,\lambda}$, and we expect both terms to diverge when
$\lambda$ goes to $0$. On the other hand, as just said, $\lambda_{\phi}$ can
be taken very small when the number of sampling points $M$ is large, so that
it compensates the other terms. In other words the computational complexity,
which is polynomial in $M$, of solving larger SDP problems will offset the
cases when $\lambda$ is taken small and the function $R_{n,\lambda}$ is less
smooth. We face thus a classical tradeoff between having a better estimate of
$R_{n,0}$ or doing fewer computations.

### The case of generalized linear embeddings

Assume $\bm{\theta}={\psi}_{\bm{w}}({\bm{x}})=\Phi_{x}{\bm{w}}$ with
$\Phi_{x}\in\mathds{R}^{d(x)\times d}$ and the added noise $Z(x)$ is Gaussian
with zero mean, then we can provide more explicit estimates for
$\operatorname{Tr}\left(A_{n,\lambda}\right)$ and
$|R_{n,\lambda}|_{\mathcal{W},\lceil s-d/2\rceil}$.

Decompose $Z(x)=\Phi_{x}Z^{\prime}(x)+Z^{\prime\prime}(x)$ with $Z^{\prime}$
Gaussian over $\mathds{R}^{d}$ with covariance $\Sigma_{x}$ and
$Z^{\prime\prime}$ Gaussian with values in the complement
$\operatorname{Im}(\Phi_{x})^{\perp}$ to the range of $\Phi_{x}$. Then

$\mathds{E}_{Z(x)}\big{[}f^{0}(\hat{\bm{y}}({\psi}_{\bm{w}}({\bm{x}})+\lambda
Z({\bm{x}})),{\bm{x}})\big{]}=\mathds{E}_{Z^{\prime}(x)}\mathds{E}_{Z^{\prime\prime}(x)}\big{[}f^{0}(\hat{\bm{y}}(\Phi_{x}({\bm{w}}+Z^{\prime}(x))+\lambda
Z^{\prime\prime}(x)),{\bm{x}})\big{]}$ (3.15)

Let
$\tilde{R}_{n,\lambda,x}=\mathds{E}_{Z^{\prime\prime}(x)}\big{[}f^{0}(\hat{\bm{y}}(\Phi_{x}{\bm{w}}+\lambda
Z^{\prime\prime}(x)),{\bm{x}})\big{]}$. We then observe that the objective
function is a sum of Gaussian convolutions

$R_{n,\lambda}({\bm{w}})=\frac{1}{n}\sum_{i=1}^{n}\tilde{R}_{n,\lambda,x_{i}}\star\mathcal{N}_{\lambda\Sigma_{x_{i}}}({\bm{w}})$
(3.16)

so $R_{n,\lambda}$ is smooth. Moreover the convolution commutes with the
derivation, and $|g\star g^{\prime}|_{\infty}\leq|g|_{\infty}|g^{\prime}|_{1}$
for any functions $g,g^{\prime}$. Moreover
$|R_{n,\lambda,x_{i}}|_{\infty}\leq|f^{0}|_{\infty}$, whence for
$\tilde{s}=s-d/2$,

$\displaystyle|R_{n,\lambda}|_{\mathcal{W},\lceil\tilde{s}\rceil}\leq$
$\displaystyle\frac{1}{n}\sum_{i=1}^{n}\sup_{|\alpha|=\tilde{s}}\sup_{\bm{\theta}\in\Omega}|\underbrace{D^{\alpha}R_{n,\lambda,x_{i}}\star\mathcal{N}_{\lambda\Sigma_{x_{i}}}({\bm{w}})}_{=R_{n,\lambda,x_{i}}\star
D^{\alpha}\mathcal{N}_{\lambda\Sigma_{x_{i}}}}|$
$\displaystyle\leq|f^{0}|_{\infty}\sup_{i\in[n],|\alpha|=\tilde{s}}|D^{\alpha}\mathcal{N}_{\lambda\Sigma_{x_{i}}}|_{1}$
$\displaystyle=\mathcal{O}\big{(}\lambda^{-\tilde{s}}\big{)}\,.$

Moreover, defining
$\mathcal{M}_{\lambda,i}=\sqrt{\mathcal{N}_{\lambda\Sigma_{x_{i}}}}$, clearly
$\mathcal{M}_{\lambda,i}\in\mathcal{H}_{\phi}$ and we can obtain a closed-form
sum-of-squares expression of $R_{n,\lambda}$ for the chosen kSoS feature map
$\phi$, so that we justify the assumption (3.12). Indeed by the reproducing
property, we have

$\displaystyle R_{n,\lambda}({\bm{w}})$
$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\left\langle{\phi({\bm{w}})},{\left[\int_{\mathds{R}^{d}}\mathcal{M}_{\lambda,i}({\bm{w}}^{\prime},\cdot)\otimes\mathcal{M}_{\lambda,i}({\bm{w}}^{\prime},\cdot)R_{n,0}({\bm{w}}^{\prime})d{\bm{w}}^{\prime}\right]\phi({\bm{w}})}\right\rangle_{H_{\phi}}$
$\displaystyle=:\left\langle{\phi({\bm{w}})},{\tilde{A}_{n,\lambda}\phi({\bm{w}})}\right\rangle_{H_{\phi}}\,.$

and since $A_{n,\lambda}$ was chosen as the operator with minimal trace

$\operatorname{Tr}(A_{n,\lambda})\leq\operatorname{Tr}(\tilde{A}_{n,\lambda})=\int_{\mathds{R}^{d}}\left\langle{\mathcal{M}_{\lambda,i}({\bm{w}}^{\prime},\cdot)},{\mathcal{M}_{\lambda,i}({\bm{w}}^{\prime},\cdot)}\right\rangle_{H_{\phi}}R_{n,0}({\bm{w}}^{\prime})d{\bm{w}}^{\prime}\,.$

Since the RKHS associated with $\phi$ is the Sobolev space of order $s$, then
we have a characterization of the norm through the Fourier transform

$\left\langle{\mathcal{M}_{\lambda,i}({\bm{w}}^{\prime},\cdot)},{\mathcal{M}_{\lambda,i}({\bm{w}}^{\prime},\cdot)}\right\rangle_{H_{\phi}}=\int_{\mathds{R}^{d}}(1+|\xi|^{2})^{s}|\hat{\mathcal{M}}_{\lambda,i}({\bm{w}}^{\prime},\xi)|^{2}d\xi=\mathcal{O}(\lambda^{-\tilde{s}}).$

######  Example 2 (continued).

The stochastic vehicle scheduling problem (StoVSP) is an instance of
generalized linear embedding. As such, the previous calculation holds.

## References

  * [1] EURO Meets NeurIPS 2022 Vehicle Routing Competition.  https://euro-neurips-vrp-2022.challenges.ortec.com/.
  * Balcan [2021] Maria-Florina Balcan.  Data-driven algorithm design.  In Tim Roughgarden, editor, _Beyond the Worst-Case Analysis of Algorithms_ , chapter 28, pages 626–645. Cambridge University Press, 2021.
  * Balcan et al. [2021] Maria-Florina Balcan, Dan DeBlasio, Travis Dick, Carl Kingsford, Tuomas Sandholm, and Ellen Vitercik.  How much data is sufficient to learn high-performing algorithms? generalization guarantees for data-driven algorithm design.  In _Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing_ , pages 919–932. ACM, 2021.  ISBN 978-1-4503-8053-9.
  * Baty et al. [2024] Léo Baty, Kai Jungel, Patrick S. Klein, Axel Parmentier, and Maximilian Schiffer.  Combinatorial Optimization-Enriched Machine Learning to Solve the Dynamic Vehicle Routing Problem with Time Windows.  _Transportation Science_ , 2024.
  * Bengio et al. [2021] Yoshua Bengio, Andrea Lodi, and Antoine Prouvost.  Machine learning for combinatorial optimization: A methodological tour d’horizon.  _European Journal of Operational Research_ , 290(2):405–421, 2021.  ISSN 0377-2217.
  * Berthet et al. [2020] Quentin Berthet, Mathieu Blondel, Olivier Teboul, Marco Cuturi, Jean-Philippe Vert, and Francis Bach.  Learning with differentiable perturbed optimizers.  In _Advances in Neural Information Processing Systems 33 (NeurIPS)_ , pages 9508–9519, 2020.
  * Bertsimas et al. [2016] Dimitris Bertsimas, Angela King, and Rahul Mazumder.  Best subset selection via a modern optimization lens.  _The Annals of Statistics_ , 44(2):813–852, 2016\.  ISSN 00905364.
  * Blondel et al. [2020] Mathieu Blondel, André F. T. Martins, and Vlad Niculae.  Learning with fenchel-young losses.  _J. Mach. Learn. Res._ , 21:35:1–35:69, 2020.
  * Boucheron et al. [2013] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  _Concentration inequalities: a nonasymptotic theory of independence_.  Oxford University Press, 1st ed edition, 2013.  ISBN 9780199535255.
  * Chafaï et al. [2012] Djalil Chafaï, Olivier Guédon, Guillaume Lecué, and Alain Pajor.  _Interactions between compressed sensing random matrices and high dimensional geometry_ , volume 37.  Société Mathématique de France Paris, 2012.
  * Ciliberto et al. [2016] Carlo Ciliberto, Lorenzo Rosasco, and Alessandro Rudi.  A consistent regularization approach for structured prediction.  In Daniel D. Lee, Masashi Sugiyama, Ulrike von Luxburg, Isabelle Guyon, and Roman Garnett, editors, _Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain_ , pages 4412–4420, 2016.
  * Ciliberto et al. [2020] Carlo Ciliberto, Lorenzo Rosasco, and Alessandro Rudi.  A general framework for consistent structured prediction with implicit loss embeddings.  _J. Mach. Learn. Res._ , 21:98:1–98:67, 2020.
  * Dalle et al. [2022] Guillaume Dalle, Léo Baty, Louis Bouvier, and Axel Parmentier.  Learning with combinatorial optimization layers: a probabilistic approach.  _ArXiv preprint_ , abs/2207.13513, 2022.
  * Elmachtoub and Grigas [2021] Adam N. Elmachtoub and Paul Grigas.  Smart “Predict, then Optimize”.  _Management Science_ , 2021.  ISSN 0025-1909.
  * Giné and Nickl [2021] Evarist Giné and Richard Nickl.  _Mathematical foundations of infinite-dimensional statistical models_.  Cambridge university press, 2021.
  * Gupta and Roughgarden [2016] Rishi Gupta and Tim Roughgarden.  A pac approach to application-specific algorithm selection.  In _Proceedings of the 2016 ACM Conference on Innovations in Theoretical Computer Science_ , pages 123–134, 2016.
  * Jones et al. [1993] D. R. Jones, C. D. Perttunen, and B. E. Stuckman.  Lipschitzian optimization without the lipschitz constant.  _Journal of Optimization Theory and Applications_ , 79(1):157–181, 1993.  ISSN 1573-2878.
  * Laurent and Massart [2000] Beatrice Laurent and Pascal Massart.  Adaptive estimation of a quadratic functional by model selection.  _Annals of statistics_ , pages 1302–1338, 2000.
  * Marteau-Ferey et al. [2020] Ulysse Marteau-Ferey, Francis R. Bach, and Alessandro Rudi.  Non-parametric models for non-negative functions.  In Hugo Larochelle, Marc’Aurelio Ranzato, Raia Hadsell, Maria-Florina Balcan, and Hsuan-Tien Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.
  * Marteau-Ferey et al. [2024] Ulysse Marteau-Ferey, Francis Bach, and Alessandro Rudi.  Second order conditions to decompose smooth functions as sums of squares.  _SIAM Journal on Optimization_ , 34(1):616–641, 2024.
  * Nowak-Vila et al. [2019] Alex Nowak-Vila, Francis Bach, and Alessandro Rudi.  A general theory for structured prediction with smooth convex surrogates.  _ArXiv preprint_ , abs/1902.01958, 2019.
  * Nowozin and Lampert [2011] Sebastian Nowozin and Christoph H. Lampert.  Structured Learning and Prediction in Computer Vision.  _Foundations and Trends® in Computer Graphics and Vision_ , 6(3–4):185–365, 2011.  ISSN 1572-2740, 1572-2759.
  * Parmentier [2021a] Axel Parmentier.  Learning to Approximate Industrial Problems by Operations Research Classic Problems.  _Operations Research_ , 2021a.  ISSN 0030-364X.
  * Parmentier [2021b] Axel Parmentier.  Learning structured approximations of operations research problems.  _ArXiv preprint_ , abs/2107.04323, 2021b.
  * Parmentier and T’Kindt [2021] Axel Parmentier and Vincent T’Kindt.  Learning to solve the single machine scheduling problem with release times and sum of completion times.  _ArXiv preprint_ , abs/2101.01082, 2021.
  * Pogancic et al. [2020] Marin Vlastelica Pogancic, Anselm Paulus, Vít Musil, Georg Martius, and Michal Rolinek.  Differentiation of blackbox combinatorial solvers.  In _8th International Conference on Learning Representations, ICLR 2020, Addis Ababa, Ethiopia, April 26-30, 2020_. OpenReview.net, 2020.
  * Qi et al. [2021] Meng Qi, Paul Grigas, and Zuo-Jun Max Shen.  Integrated conditional estimation-optimization.  _arXiv preprint arXiv:2110.12351_ , 2021.
  * Rudi et al. [2024] Alessandro Rudi, Ulysse Marteau-Ferey, and Francis Bach.  Finding global minima via kernel approximations.  _Mathematical Programming_ , 2024.  ISSN 1436-4646.
  * Sadana et al. [2024] Utsav Sadana, Abhilash Chenreddy, Erick Delage, Alexandre Forel, Emma Frejinger, and Thibaut Vidal.  A survey of contextual optimization methods for decision-making under uncertainty.  _European Journal of Operational Research_ , 2024.
  * Shapiro et al. [2021] Alexander Shapiro, Darinka Dentcheva, and Andrzej Ruszczynski.  _Lectures on stochastic programming: modeling and theory_.  SIAM, 2021.
  * Vershynin [2018] Roman Vershynin.  _High-dimensional probability: an introduction with applications in data science_.  Cambridge series in statistical and probabilistic mathematics. Cambridge University Press, 2018.  ISBN 9781108231596.
  * Vu et al. [2017] Ky Khac Vu, Claudia d’Ambrosio, Youssef Hamadi, and Leo Liberti.  Surrogate-based methods for black-box optimization.  _International Transactions in Operational Research_ , 24(3):393–424, 2017.
  * Wainwright [2019] Martin J Wainwright.  _High-dimensional statistics: A non-asymptotic viewpoint_ , volume 48.  Cambridge university press, 2019.

Supplemental Materials

##  Appendix A Technical lemmas

######  Lemma 6.

Let $u\,:\,\mathds{R}^{d}\to\mathds{R}$ be an integrable bounded function, $Z$
a standard normal random vector of $\mathds{R}^{d}$, $\sigma>0$ a positive
real number, and $G({\bm{w}}):=\mathds{E}_{Z}[u({\bm{w}}+\sigma Z)]$. Then
${\bm{w}}\mapsto G({\bm{w}})$ is $\sqrt{d}\|u\|_{\infty}/\sigma$-Lipschitz.

###### Proof.

Let $h$ be the density of $\sigma Z$. One has
$G({\bm{w}})=\int_{\mathds{R}^{d}}u({\bm{w}}+{\bm{v}})h({\bm{v}})\mathrm{d}{\bm{v}}=\int_{\mathds{R}^{d}}u({\bm{v}})h({\bm{v}}-{\bm{w}})\mathrm{d}{\bm{v}}$.
By dominated convergence, it holds that

$\nabla G({\bm{w}})=-\int_{\mathds{R}^{d}}u({\bm{v}})\nabla
h({\bm{v}}-{\bm{w}})\mathrm{d}{\bm{v}}=-\int_{\mathds{R}^{d}}u({\bm{w}}+{\bm{v}})\nabla
h({\bm{v}})\mathrm{d}{\bm{v}}\,.$

By Jensen inequality and standard calculation with the Gaussian density, one
has

$\|\nabla G({\bm{w}})\|\leq\|u\|_{\infty}\int_{\mathds{R}^{d}}\|\nabla
h({\bm{v}})\|\mathrm{d}{\bm{v}}\leq\frac{\sqrt{d}\|u\|_{\infty}}{\sigma}\,,$

which gives the result. ∎

######  Lemma 7.

Under Assumptions $($Lip$)$ and $($Gauss$)$, the following holds. When
$\varepsilon_{0}>0$, the risk
${\bm{w}}\mapsto\mathcal{R}_{\varepsilon_{0}}(h_{\bm{w}})$ is continuous with
respect to ${\bm{w}}$.

###### Proof.

Let ${\bm{x}}\in\mathcal{X}$ and
${\bm{y}},\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$. Invoke Lemma 6 with
$u(\bm{\theta})=p_{0}({\bm{y}}|\bm{\theta})\in[0,1]$ to get that
$\bm{\theta}\mapsto\mathds{E}_{Z}[\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{0}({\bm{y}}|\bm{\theta}+{\varepsilon_{0}}Z({\bm{x}}))\,f^{0}({\bm{y}},{\bm{x}})]$
is continuous. Now, note that ${\bm{w}}\mapsto{\psi}_{\bm{w}}({\bm{x}})$ is
$L_{\mathcal{W}}$-Lipschitz by Assumption $($Lip$)$, hence
${\bm{w}}\mapsto\mathds{E}_{Z}[\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{0}({\bm{y}}|{\psi}_{\bm{w}}({\bm{x}})+{\varepsilon_{0}}Z({\bm{x}}))\,f^{0}({\bm{y}},{\bm{x}})]$
is continuous. Taking the expectation, by dominated convergence, we deduce the
result. ∎

######  Lemma 8.

Let $\lambda>0$, ${\bm{x}}\in\mathcal{X}$ and
${\bm{y}}\in\mathcal{Y}({\bm{x}})$. Then, the function
$\bm{\theta}\in\mathds{R}^{d({\bm{x}})}\mapsto\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{\lambda}({\bm{y}}|\bm{\theta})$
is $\sqrt{d({\bm{x}})}/\lambda$-Lipschitz, under Assumption $($Gauss$)$.

###### Proof.

Note that
$\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{\lambda}({\bm{y}}|\bm{\theta})=\mathds{E}_{Z}[u(\bm{\theta}+\lambda
Z)]$ where
$u(\cdot)=\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{0}({\bm{y}}|\cdot)\in[0,1]$
and apply Lemma 6. ∎

######  Lemma 9.

The function
${\bm{w}}\in\mathcal{W}\mapsto\sum_{{\bm{y}}\in\mathcal{Y}({\bm{x}})}p_{\lambda}({\bm{y}}|{\psi}_{\bm{w}}({\bm{x}}))\in[0,1]$
is $L_{\mathcal{W}}\sqrt{d({\bm{x}})}/\lambda$-Lipschitz, under Assumptions
$($Lip$)$ and $($Gauss$)$.

###### Proof.

Let ${\bm{x}}\in\mathcal{X}$ and
${\bm{y}},\bm{\theta}\in\mathds{R}^{d({\bm{x}})}$. Invoke Lemma 8 to get that
$\bm{\theta}\mapsto p_{\lambda}({\bm{y}}|\bm{\theta})$ is
$\sqrt{d({\bm{x}})}/\lambda$-Lipschitz. Now, observe that
${\bm{w}}\mapsto{\psi}_{\bm{w}}({\bm{x}})$ is $L_{\mathcal{W}}$-Lipschitz to
conclude. ∎

##  Appendix B Perturbation bias with remarks

######  Proposition 10.

If Assumption $($Gauss$)$ holds, then, for all $\lambda\in(0,1)$ and all
${\bm{w}}\in\mathcal{W}$,

$V_{\bm{w}}(\lambda)\leq\inf_{q\in(0,1)}\Bigg{\\{}\mathds{P}_{X}\Big{(}\frac{\rho(\psi_{{\bm{w}}}(X))}{\sqrt{d(X)}}<\lambda^{q}\Big{)}+\exp\Big{(}-\frac{1}{10\lambda^{2(1-q)}}\Big{)}\Bigg{\\}}\,.$

Furthermore, under the following Assumption (Orl.),

$\exists\tau>0\text{ s.t.
}\mathds{E}_{X}\bigg{[}\exp\Big{(}\big{(}\frac{\rho(\psi_{{\bm{w}}}(X))}{\sqrt{d(X)}}\big{)}^{-\tau}\Big{)}\bigg{]}<\infty\,,$
(Orl.)

it holds

$V_{\bm{w}}(\lambda)=\mathcal{O}_{\lambda\to
0}\Big{(}\exp(-\lambda^{-\frac{2\tau}{2+\tau}})\Big{)}\,;$

And under the following Assumption (Mom.),

$\exists\tau>0\text{ s.t.
}\mathds{E}_{X}\bigg{[}\Big{(}\frac{\rho(\psi_{{\bm{w}}}(X))}{\sqrt{d(X)}}\Big{)}^{-\tau}\bigg{]}<\infty\,,$
(Mom.)

it holds

$V_{\bm{w}}(\lambda)=\mathcal{O}_{\lambda\to
0}\big{(}\lambda^{\tau}\mathrm{polylog}(\lambda)\big{)}\,,$

where $\mathrm{polylog}(\lambda)$ is a polynomial logarithm term.

###### Proof.

Invoke [Laurent and Massart, 2000, Lemma 1] to get that, for
${\bm{x}}\in\mathcal{X}$, for all $t\geq d({\bm{x}})$,

$\mathds{P}[\|Z({\bm{x}})\|_{2}^{2}\geq t]\leq\exp(-0.1\,t)\,,$

for the Gaussian model (see Assumption ($($Gauss$)$)). Hence, for all
$\lambda\in(0,1)$, all $q\in(0,1)$,

$\forall{\bm{x}}\text{ s.t.
}\rho(\psi_{{\bm{w}}}({\bm{x}}))\geq\lambda^{q}\sqrt{d({\bm{x}})}\,,\quad\mathds{P}\Big{[}\|Z({\bm{x}})\|_{2}\geq\frac{\rho(\psi_{{\bm{w}}}({\bm{x}}))}{\lambda}\Big{]}\leq\exp\Big{(}-0.1\,\frac{\rho^{2}(\psi_{{\bm{w}}}({\bm{x}}))}{\lambda^{2}}\Big{)}\,.$

Hence It follows that, for all $\lambda\in(0,1)$ and all $q\in(0,1)$,

$\displaystyle V(\lambda)$
$\displaystyle=\mathds{E}_{X}\mathds{P}_{Z}\Big{[}\|Z\|_{2}\geq\frac{\rho(\psi_{{\bm{w}}}(X))}{\lambda}\,\big{|}\,X\Big{]}\,,$
$\displaystyle\leq\mathds{P}_{X}\Big{[}\rho(\psi_{{\bm{w}}}(X))<\lambda^{q}\sqrt{d(X)}\Big{]}+\mathds{E}_{X}\Big{[}\exp\Big{(}-0.1\,\frac{\rho^{2}(\psi_{{\bm{w}}}(X))}{\lambda^{2}}\Big{)}\mathbf{1}_{\\{\rho(\psi_{{\bm{w}}}(X))\geq\lambda^{q}\sqrt{d(X)}\\}}\Big{]}\,,$
$\displaystyle\leq\mathds{P}_{X}\Big{[}\rho(\psi_{{\bm{w}}}(X))<\lambda^{q}\sqrt{d(X)}\Big{]}+\mathds{E}_{X}\Big{[}\exp\Big{(}-\frac{d(X)}{10\lambda^{2(1-q)}}\Big{)}\Big{]}\,,$
$\displaystyle\leq\mathds{P}_{X}\Big{[}\rho(\psi_{{\bm{w}}}(X))<\lambda^{q}\sqrt{d(X)}\Big{]}+\exp\Big{(}-\frac{1}{10\lambda^{2(1-q)}}\Big{)}\,,$

the last inequality stemming from $d(X)\geq 1$ $\mathds{P}_{X}$-almost surely.

Assume (Orl.), by Markov’s inequality it holds

$\mathds{P}_{X}\Big{[}\rho(\psi_{{\bm{w}}}(X))<\lambda^{q}\sqrt{d(X)}\Big{]}\leq
C\exp\big{(}-\lambda^{-\tau q}\big{)}\,.$

Taking $q=2/(2+\tau)$ gives the result.

Assume (Mom.), by Markov’s inequality it holds

$\mathds{P}_{X}\Big{[}\rho(\psi_{{\bm{w}}}(X))<\lambda^{q}\sqrt{d(X)}\Big{]}\leq
C\,\lambda^{\tau q}\,.$

Taking $q=1-\log(-10\tau\log\lambda)/(-2\log(\lambda))$ gives the result (note
that $q<1$, and $q>0$ for $\lambda$ sufficient small). ∎

######  Remark 12 (On the conditions (Orl.) and (Mom.)).

These conditions describe the tail behavior of the random variable
$R:={\sqrt{d(X)}}/{\rho(\psi_{{\bm{w}}}(X))}$. Condition (Orl.) is of Orlicz-
norm type, it states that $R$ belongs to $\tau$-Orlicz space and, as such as,
it has $\tau$-exponentially light tail, see for instance [Chafaï et al., 2012,
Section 1.1] for a review on Orlicz spaces. Hence, the distribution of the
random variable ${\rho(\psi_{{\bm{w}}}(X))}/{\sqrt{d(X)}}=1/R$ is
$\tau$-exponentially light at edge support point zero. The random variable
${\rho(\psi_{{\bm{w}}}(X))}/{\sqrt{d(X)}}=1/R$ is close to zero when switching
from a normal cone to another. Roughly speaking, Condition (Orl.) means that
the distribution of $\psi_{{\bm{w}}}(X)$ does not put too much mass around
decision boundaries (the border of the normal cones), the learnt direction
$\psi_{{\bm{w}}}(X)$ points towards the interior of the normal cone far from
its boundaries.

Condition (Mom.) is a moment type assumption on $R$. It is weaker than
Condition (Orl.). Indeed, by [Chafaï et al., 2012, Theoreme 1.1.5], note that
(Orl.) is equivalent to

$\forall
p\geq\tau\,,\quad\bigg{[}\mathds{E}_{X}\Big{(}\frac{\rho(\psi_{{\bm{w}}}(X))}{\sqrt{d(X)}}\Big{)}^{-p}\bigg{]}^{\frac{1}{p}}\leq
Cp^{\tau}\,,$

for some positive constant $C>0$. Hence, Condition (Orl.) implies (Mom.).
Condition (Mom.) gives a slower rate of convergence of $V(\lambda)$ than
Condition (Orl.).

Note that both convergence rates are decreasing in $\tau$ and one cannot have
better rate than the sub-Gaussian type bound $\mathcal{O}_{\lambda\to
0}\big{(}\exp(-\lambda^{-2})\big{)}$.