# Learning-to-Optimize with PAC-Bayesian Guarantees: Theoretical
Considerations and Practical Implementation

Michael Sucker<EMAIL_ADDRESS>
Department of Mathematics
University of Tübingen
Tübingen, Germany Jalal Fadili<EMAIL_ADDRESS>
ENSICAEN
Normandie Université
CNRS, GREYC, France Peter Ochs<EMAIL_ADDRESS>
Department of Mathematics and Computer Science
Saarland University
Saarbrücken, Germany

###### Abstract

We use the PAC-Bayesian theory for the setting of learning-to-optimize. To the
best of our knowledge, we present the first framework to learn optimization
algorithms with provable generalization guarantees (PAC-Bayesian bounds) and
explicit trade-off between convergence guarantees and convergence speed, which
contrasts with the typical worst-case analysis. Our learned optimization
algorithms provably outperform related ones derived from a (deterministic)
worst-case analysis. The results rely on PAC-Bayesian bounds for general,
possibly unbounded loss-functions based on exponential families. Then, we
reformulate the learning procedure into a one-dimensional minimization problem
and study the possibility to find a global minimum. Furthermore, we provide a
concrete algorithmic realization of the framework and new methodologies for
learning-to-optimize, and we conduct four practically relevant experiments to
support our theory. With this, we showcase that the provided learning
framework yields optimization algorithms that provably outperform the state-
of-the-art by orders of magnitude.

##  1 Introduction

Typically, optimization algorithms are derived by considering a _class of
problems_ , which consists of three ingredients: First, the _model_ , which
summarizes the imposed assumptions and common properties of the problems
within that class. Second, the _oracle_ , which describes the information the
algorithm actually receives during run-time. And third, the _stopping
criterion_ , that specifies when the algorithm terminates. For example,
consider all convex and continuously differentiable functions with Lipschitz
continuous gradient (model), a first-order method that has access to the
function value and its gradient in each iteration (oracle), with the overall
goal to decrease the residual in function values up to a certain threshold
(stopping-criterion).
Based on this class, one can perform a worst-case analysis to make sure that
_in any case_ a sufficient descent property, for example in function values,
does hold along the iterations, which allows to prove convergence and derive
corresponding convergence rates. Thus, such a worst-case analysis allows for
theoretical guarantees by being applicable to any case in question. However,
often this comes at the cost of slower convergence in the _typical case_ , as
considering the worst-case only is _inevitably_ accompanied by a negligence of
information and thus, an overly conservative estimation in most of the cases.
Hence, there is an inherent trade-off between applicability in terms of
imposed assumptions, theoretical guarantees and practical performance. A
worst-case analysis aims for one of the extremes: Strong theoretical
guarantees at the cost of potentially strong assumptions and practical
performance. Doing so, we get an _uniform upper bound_ on the performance.
Yet, if the worst-case is actually equally likely as any other problem in that
class, this indeed is the best one can hope for, as there are simple examples
for which this upper bound then also matches the best possible lower bound,
yielding the notion of an _optimal_ algorithm. Therefore, to improve the
performance of an optimization algorithm, one needs to have access to more
information about the problems at hand. This can be done _analytically_ , that
is, by imposing further restrictive assumptions, or _statistically_ under the
assumption that the single problem instances are not equally likely to be
observed. Here, the first approach boils down to reducing the class of
functions in question, hence a change of the model. That means that one might
get a faster algorithm, yet at the expense of less applicability. Therefore,
it does not alter the underlying trade-off effectively. The second approach,
however, rather aims for a change of the oracle and tries to use structure in
the given class of functions, which is _explicitly_ not analytically
tractable, that is, not already encoded in the model. This means that, if
there is actually more structure in the class of functions than the model
describes analytically (the model is _larger_ than the problem class at hand),
fast algorithms can be obtained without making further restrictive
assumptions. The premise for this is that one actually manages to access and
use this structure in other ways than by pen and paper, for example,
statistically. One particular way to do this is to _learn_ the optimization
algorithm on a given data set. For this, the oracle gets enriched by access to
a data set that can be used for a training phase. By doing so, the algorithm
becomes _data-dependent_ and can in turn enrich its oracle by the statistical
structure of the class of problems. Hence, learning actually alleviates the
bounds of analytical tractability by introducing data-dependency. Furthermore,
it also allows for new algorithms, which cannot be analyzed (or even derived)
analytically and therefore were not tractable before. The whole process can be
automated, such that, in total, one can access more information and train more
complicated algorithms with actually less hand-crafting. Hence, there is a
huge potential up-side to learning optimization algorithms. However, there is
no free-lunch. Indeed, there is also a down-side: Since one insists on using
structure and algorithms that are explicitly not analytically tractable, one
looses some theoretical guarantees that were built upon exactly these analytic
features. Hence, while learning alleviates the bounds of analytical
tractability, it does it at the cost of loosing the _typical_ theoretical
guarantees. Since the practical usefulness of an optimization algorithm
without convergence guarantees is at least questionable, this is a major
problem, and poses the first central question:

What kind of theoretical guarantees can be given for such a learned
optimization algorithm? Are we able to _ensure_ its usefulness?

Hence, in the first part of this paper we provide a theoretical framework for
learning an optimization algorithm based on its performance on the training
set, that is, based on the empirical risk together with a generalization bound
for the (true) risk. A popular framework that provides such generalization
bounds is the PAC-Bayesian approach to learning, which we apply to the setting
of learning-to-optimize. This yields the following informal version of our
main theoretical result (compare Theorem 4.1 together with Corollary 3.5):

######  Theorem 1.1 (Informal)

Under mild boundedness assumptions on the optimization algorithm, the
$\mathbb{Q}$-average population loss $\mathcal{R}_{g}$ of the last iterate can
be bounded by the $\mathbb{Q}$-average empirical loss $\hat{\mathcal{R}}_{g}$
of the last iterate plus some remainder term $r_{N}$ that vanishes with the
size $N$ of the data set, that is, for all $\varepsilon>0$:

$\displaystyle\mathbb{P}_{\mathbf{\mathfrak{D}}}\left\\{\forall\lambda\in\Lambda,\
\forall\mathbb{Q}\in\mathcal{P}:\mathbb{Q}[\mathcal{R}_{g}]\leq\mathbb{Q}[\hat{\mathcal{R}}_{g}]+r_{N}\right\\}\geq
1-\varepsilon\,.$

This result states that it is unlikely to observe a data set for which the
given bound on the risk of the last iterate would not hold uniformly in
$\lambda$ or $\mathbb{Q}$. Especially, the uniformity in $\mathbb{Q}$ allows
for _learning_ such a distribution $\mathbb{Q}$. To proof this result, in
Theorem 3.3 we provide a general PAC-Bayesian theorem for data-dependent
exponential families, which holds under mild assumptions, thus being widely
applicable. Then, in Corollary 3.5, by specifying the natural parameters and
the sufficient statistics, that is, by choosing a specific exponential family,
we directly deduce a PAC-Bayesian generalization bound. Here, since these two
results are widely applicable and constructive in nature, they might be of
independent interest. Finally, to be able to apply Theorem 3.3 and Corollary
3.5 to the setting of learning-to-optimize we identified suitable properties
of optimization algorithms, such that the corresponding assumptions are
satisfied (compare Theorem 4.1 and Theorem 4.9).
Taken together, this provides one possible answer to the central question
posed above about the theoretical guarantees. However, while being a
generalization bound, these guarantees are a statement about _relative values_
, that is, how the true risk compares to the empirical risk. Though, they do
not directly imply that these have to be small _absolute values_. Thus, one
still has to train the optimization algorithm in such a way that the empirical
risk is indeed small enough to be worth the effort. Hence, the second central
question that arises is of a more practical nature and pertains to the actual
training of such an algorithm:

How do we learn an optimization algorithm, so that its performance is
_clearly_ superior to the one achieved by a worst-case analysis?

Figure 1: Some numerical results: Loss over iterations (mean as dashed and
median as dotted line) of the learned algorithm compared to a standard choice.

Therefore, in the second part of this work, we develop a concrete algorithmic
realization, which allows for learning an optimization algorithm and
evaluating the corresponding theoretical guarantee. This involves several
design choices and tricks that have not been used before and which are of
interest in their own right. Furthermore, as empirical evaluation of our
claims, we conduct four practically relevant experiments, all dealing with
very different classes of functions, thereby demonstrating the wide
applicability and strong practical performance of our approach. Examples are
shown in Figure 1. Each subplot compares the performance of the learned
algorithm to that of a standard algorithm on different problems: a smooth and
strongly convex quadratic toy problem (upper left), a high-dimensional, smooth
and convex image processing problem (upper right), a non-smooth and convex
LASSO problem (lower left), and a non-smooth and non-convex problem of
training a neural network (lower right). Since the learned algorithm is
clearly superior in each case, this provides a possible answer to the question
about how to train optimization algorithms. Taken together, we provide a
complete framework to train optimization algorithms with theoretical
guarantees that are (in a certain sense) provably faster than their worst-case
optimal counterparts. With this, we get closer to the ideal case of wide
applicability, strong practical performance and sufficient theoretical
guarantees. In particular, this work is a far reaching extension of our
conference paper (Sucker and Ochs, 2023) by extending and clarifying the
theoretical results in Sections 3, and by the algorithmic realization together
with its evaluation in Sections 5 and 6, which additionally includes a
probabilistic constraining procedure for sampling algorithms in Sections 5.2.

###  1.1 Related Work

The literature on both learning-to-optimize and the PAC-Bayes learning
approach is vast. Hence, the discussion of learning-to-optimize will mainly
focus on learning procedures for optimization algorithms and on approaches
that provide some theoretical guarantees. Especially, this excludes many
model-free approaches, which replace the whole update step with a learnable
mapping such as a neural network. Chen et al. (2021) provide a good overview
about the variety of approaches in learning-to-optimize, and good introductory
references for the PAC-Bayesian approach are given by Guedj (2019) and Alquier
(2021).

##### Broader Context of Learning-to-Optimize.

Solving optimization problems is an integral part of machine learning. Thus,
learning-to-optimize has significant overlap with the areas of meta-learning
(or “learning-to-learn”) and AutoML. The first one is a subset of learning-to-
optimize, as, while learning-to-optimize applies to general optimization
problems, meta-learning is mostly concerned with determining parameters of
machine learning models, that is, minimizing a training loss (Vilalta and
Drissi, 2002; Hospedales et al., 2021). AutoML, however, more broadly refers
to automating all steps necessary to create a machine learning application,
which thus also involves the choice of an optimization algorithm and its
hyperparameters (Yao et al., 2018; Hutter et al., 2019; He et al., 2021).
Hence, most relevant for this work is model selection.

##### Learning-to-Optimize with Guarantees.

Chen et al. (2021) point out that learned optimization methods may lack
theoretical guarantees for the sake of convergence speed. That being said,
there are applications where a convergence guarantee is of highest priority.
To underline this, Moeller et al. (2019) provide an example where a purely
learning-based approach fails to reconstruct the crucial details in a medical
image. Also, they prove convergence of their method by restricting the output
to descent directions, for which mathematical guarantees exist. The basic idea
is to trace the learned object back to, or constrain it to, a mathematical
object with convergence guarantees. Similarly, Sreehari et al. (2016) provide
sufficient conditions under which the learned mapping is a proximal mapping.
Related schemes, under different assumptions and guarantees, are given by Chan
et al. (2016), Teodoro et al. (2017), Tirer and Giryes (2018), Buzzard et al.
(2018), Ryu et al. (2019), Sun et al. (2019), Terris et al. (2021) and Cohen
et al. (2021). A major advantage of these methods is the fact that the number
of iterations is not restricted a priori. However, a major drawback is their
restriction to specific algorithms and problems. Another approach, which
limits the number of iterations, yet in principle can be applied to every
iterative optimization algorithm, is unrolling, pioneered by Gregor and LeCun
(2010). Xin et al. (2016) study the IHT algorithm and show that it is, under
some assumptions, able to achieve a linear convergence rate. Likewise, Chen et
al. (2018) establish a linear convergence rate for the unrolled ISTA. However,
a difficulty in the theoretical analysis of unrolled algorithms is actually
the notion of convergence itself, such that one rather has to consider the
generalization performance. Only few works have addressed this: Either
directly by means of Rademacher complexity (Chen et al., 2020b), or indirectly
in form of a stability analysis (Kobler et al., 2020), as algorithmic
stability is linked to generalization (Bousquet and Elisseeff, 2000, 2002;
Shalev-Shwartz et al., 2010). We consider a mixture of the iterative approach
and the unrolling approach: The theoretical analysis corresponds to the
approach of unrolling, that is, a fixed number of iterations, and the
resulting PAC-Bayesian generalization bound applies to this particular number
osubsubsectionf iterations only. However, in our implementation and
experiments, we stay more closely to the iterative approach of learning an
update step that can be applied for an arbitrary number of iterations.
Further, while the resulting algorithms do include neural networks in their
update step, the construction is closely related to the model-based approach,
as it is inspired by analytical algorithms.

##### Typical Problems and Tricks in Learning-to-Optimize.

The literature contains a vast amount of approaches, tricks and problems about
learning procedures for optimization algorithms. We only touch briefly upon
the ones that are most relevant for the present work: Besides the actual
optimizer, a typical design choice is that of the loss function. Here, either
the final loss or a weighted sum of the losses incurred along the iterations
is typically being used (Chen et al., 2021). This extends to choices about the
computational graph, as one has to decide which gradients are computed during
training. An assumption already made by Andrychowicz et al. (2016) and widely
adopted thereafter, is to neglect higher-order derivatives. Importantly, if
the theoretical analysis is coherent with the computational graph, one is
allowed to do whatever works here.
A major problem of many learned optimization algorithms, especially the ones
based on recurrent neural networks, is their restriction to a certain number
of iterations, that is, they cannot be trained for an arbitrary number of
iterations due to instabilities, such as vanishing or exploding gradients, and
memory bottlenecks. Further, often they do not generalize well to more
iterations (Andrychowicz et al., 2016; Chen et al., 2017; Lv et al., 2017;
Chen et al., 2021). A typical way to mitigate this problem is to split the
whole trajectory into smaller parts, for which the vanishing/exploding
gradient phenomenon does not occur, which effectively also solves the memory
problem. However, this often does not lead to fully satisfactory results
either, such that other approaches have been proposed: Lv et al. (2017)
introduce random scaling of the coordinates and the addition of a convex
function to the objective. Metz et al. (2019) propose not to model the
recurrent nature of optimization algorithms directly, and therefore to replace
the recurrent neural network with a multilayer perceptron, to “smooth” the
loss, and to dynamically use two unbiased gradient estimators instead of one.
Doing so they manage to train algorithms that are faster in wall-clock time
than standard ones like Adam. Wichrowska et al. (2017) introduce a
hierarchical RNN architecture consisting of three RNNs, and additionally draw
the number of unrollings and the unrolling length from a heavy-tailed
exponential distribution. While achieving the needed generalization, this
approach is computationally costly and does not achieve the same wall-clock
time as simple optimization algorithms. Chen et al. (2020a) consider improved
training techniques in general, and they introduce a progressive scheme that
gradually increases the unrolling length, and an imitation learning approach
by learning from analytic optimizers first. We will introduce a new loss
function for training optimization algorithms, motivated by an intuitive
theoretical argument. To our best knowledge, this loss function is new. Here,
we will typically also ignore higher-order derivatives in the computational
graph, as the algorithms performance is analyzed independently of this step in
the learning procedure. Further, we found that a combination of many of the
above mentioned approaches works well: 1) We use a single learned update based
on MLPs and analytic algorithms instead of an RNN. 2) We split the trajectory
into subtrajectories and randomize its total length, however, in a different
and new way, again motivated by theoretical considerations. 3) We start the
learning procedure with imitation learning to find an initialization that
performs similar to analytic algorithms.

##### PAC-Bounds through Change-of-Measure.

The PAC-Bayesian framework allows for giving high probability bounds on the
risk, either as an oracle or as an empirical bound. The key ingredients are
so-called change-of-measure inequalities. The choice of such an inequality
strongly influences the corresponding bound. The one used most often is based
on a variational representation of the Kullback–Leibler divergence due to
Donsker and Varadhan (1975), employed, for example, by Catoni (2004, 2007).
Yet, not all bounds are based on a variational representation, that is,
holding uniformly over all posterior distributions (Rivasplata et al., 2020).
However, most bounds involve the Kullback–Leibler divergence as measure of
proximity, for example, those by McAllester (2003a, b), Seeger (2002),
Langford and Shawe-Taylor (2002), or the general PAC-Bayes bound of Germain et
al. (2009). More recently, other divergences have been used: Honorio and
Jaakkola (2014) prove an inequality for the $\chi^{2}$-divergence, which is
also used by London (2017). Bégin et al. (2016) and Alquier and Guedj (2018)
use the Renyi-divergence ($\alpha$-divergence). Ohnishi and Honorio (2021)
propose several PAC-bounds based on general f-divergences, which include the
Kullback–Leibler-, $\alpha$\- and $\chi^{2}$-divergences. Even more recently,
Amit et al. (2022) propose to replace the Kullback-Leibler divergence by so-
called “integral probability metrics”, which encompass, for example, the
Wasserstein distance that obeys many favorable properties and also captures
the geometry of the underlying space (see Villani et al., 2009). Motivated by
this, Haddouche and Guedj (2023) also investigate PAC-Bayesian generalization
bounds for the Wasserstein distance and their interplay with the output of
optimization algorithms. A major advantage of using the Wasserstein distances
instead of the Kullback-Leibler divergence is the fact that it does not
constrain the support of the distribution a-priori through the choice of the
prior. On the other hand, it demands assumptions on the loss function, which
are not necessarily satisfied in learning-to-optimize. We give a general PAC-
Bayesian theorem based on exponential families. Here, the role of prior,
posterior, divergence and data dependence will be given naturally. Moreover,
this approach allows us to implement a more abstract learning framework that
can be applied to a wide variety of algorithms.

##### Boundedness of the Loss Function.

A major drawback of many of the existing PAC-Bayes bounds is the assumption of
a bounded loss-function. However, this assumption is mainly there to apply
some exponential moment inequality like the Hoeffding- or Bernstein-inequality
(Rivasplata et al., 2020; Alquier, 2021). Several ways have been developed to
solve this problem: Germain et al. (2009) explicitly include the exponential
moment in the bound, Alquier et al. (2016) use so-called Hoeffding- and
Bernstein-assumptions, Catoni (2004) restricts to the sub-Gaussian or sub-
Gamma case. Another possibility to ensure the generalization or exponential
moment bounds is to use properties of the algorithm in question. London (2017)
uses algorithmic stability to provide PAC-Bayes bounds for SGD. We consider
suitable properties of optimization algorithms aside from algorithmic
stability to ensure the exponential moment bounds. To the best of our
knowledge, this is new.

##### Minimization of the PAC-Bound.

The PAC-bound relates the true risk to other terms such as the empirical risk.
Yet, it does not directly say anything about the absolute numbers. Thus, in
learning procedures based on the PAC-Bayesian theory one typically aims to
minimize this bound: Langford and Caruana (2001) compute non-vacuous numerical
generalization bounds through a combination of PAC-bounds with a sensitivity
analysis. Dziugaite and Roy (2017) extend this by minimizing the PAC-bound
directly. Pérez-Ortiz et al. (2021) also consider learning as minimization of
the PAC-Bayes bound and provide tight generalization bounds. Thiemann et al.
(2017) are able to solve the minimization problem resulting from their PAC-
bound by alternating minimization. Further, they provide sufficient conditions
under which the resulting minimization problem is quasi-convex. We also follow
this approach and consider learning as minimization of the PAC-Bayesian upper-
bound, however, applied to the context of learning-to-optimize.

##### Choice of the Prior.

A common difficulty in learning with PAC-Bayesian bounds is the choice of the
prior distribution, as it heavily influences the performance of the learned
models and the generalization bound (Catoni, 2004; Dziugaite et al., 2021;
Pérez-Ortiz et al., 2021). In part, and especially for the Kullback-Leibler
divergence, this is due to the fact that the divergence term can dominate the
bound, keeping the posterior close to the prior. This leads to the idea of
choosing a data- or distribution-dependent prior (Seeger, 2002; Parrado-
Hernández et al., 2012; Lever et al., 2013; Dziugaite and Roy, 2018; Pérez-
Ortiz et al., 2021), which, by using an independent subset of the data set,
gets optimized to yield a good performance. The choice of the prior
distribution is crucial for the performance of our learned algorithms. Thus,
we use a data-dependent prior. Further, we show how the prior is essential in
preserving crucial properties during learning. It is key to control the trade-
off between convergence guarantee and convergence speed. Based on these
insights, a major part of our algorithmic realization is concerned with
constructing such a prior distribution.

##### More Generalization Bounds.

There are many areas of machine learning research that study generalization
bounds and have not been discussed here. Importantly, the vast field of
“stochastic optimization” (SO) provides generalization bounds for specific
algorithms. The main differences to our setting are the learning approach and
the assumptions made:

  * •

In most of the cases, the concrete algorithms studied in SO generate a single
point estimate by either minimizing the (regularized) empirical risk
functional over a possibly large data set, or by repeatedly updating the point
estimate based on a newly drawn (small) batch of samples. Then, one studies
the properties of this point in terms of the stationarity measure of the true
risk functional (Bottou et al., 2018; Davis and Drusvyatskiy, 2022; Bianchi et
al., 2022). In our case, as we consider the PAC-Bayesian learning approach,
the final object to be studied is a distribution over hyperparameters.

  * •

As the setting in SO is more explicit, more assumptions have to be made.
Typical assumptions are (weak) convexity (Shalev-Shwartz et al., 2009; Davis
and Drusvyatskiy, 2019), bounded gradients (Défossez et al., 2022), bounded
noise (Davis and Drusvyatskiy, 2022), or smoothness (Kavis et al., 2022).
Since we consider an abstract optimization algorithm, and the problem of
finding a distribution over its hyperparameters, all of these assumptions
cannot be made without severely limiting the applicability of our results.

_We provide a principled way to learn a distribution over general
hyperparameters of an abstract algorithm under weak assumptions and go
explicitly beyond analytically tractable quantities. Therefore, the
methodology is independent of the chosen implementation._

##  2 Problem Setup & Assumptions

In this section we establish the notation, formalize the setting, and state
the main assumptions that are used throughout the remainder of the text.

##### Notation.

We will endow every topological space $X$ with the corresponding
Borel-$\sigma$-algebra $\mathcal{B}(X)$, and, given a product space $X\times
Y$ of two measurable spaces $(X,\mathcal{A})$ and $(Y,\mathcal{B})$, we endow
it with the product-$\sigma$-algebra $\mathcal{A}\otimes\mathcal{B}$. To
denote the product space of a larger number of spaces ${X}_{1},...,{X}_{n}$,
we will use $\prod_{i=1}^{n}{X}_{i}$ as the shorthand for
${X}_{1}\times...\times{X}_{n}$ and $\bigotimes_{i=1}^{n}\mathcal{B}({X}_{i})$
for the product-$\sigma$-algebra. Further, if $X=X_{1}=...=X_{n}$, we also use
$X^{n}=\prod_{i=1}^{n}X_{i}$. If we are given a multivariate function
$f:{X}\times{Y}\to{Z}$, then $f(x,\cdot):{Y}\to{Z}$ denotes the function
$y\mapsto f(x,y)$ with fixed element $x\in{X}$. Similarly, for a set
$\mathsf{C}\subset{X}\times{Y}$ we denote the section of $\mathsf{C}$ for
fixed $x\in{X}$ by $\mathsf{C}_{x}:=\\{y\in{Y}\ :\ (x,y)\in\mathsf{C}\\}$. In
general, generic sets are denoted in the typewriter font, for example
$\mathsf{A}$, and $\mathds{1}_{\mathsf{A}}$ denotes the function that is equal
to one for $x\in\mathsf{A}$ and zero else, while $\iota_{\mathsf{A}}$ denotes
the function that is equal to zero for $x\in\mathsf{A}$ and $+\infty$
else.111We omit the name here, as both $\mathds{1}_{\mathsf{A}}$ and
$\iota_{\mathsf{A}}$ are called “indicator function”. The former in
probability theory, the latter in optimization. For a measure $\mu$ on a
measurable space ${X}$, and a measurable function $f$ on ${X}$, $\mu[f]$
denotes the integral of $f$ w.r.t. $\mu$, while, if $f\geq 0$, $f\cdot\mu$
denotes the measure given by $(f\cdot\mu)[\mathsf{A}]=\int_{\mathsf{A}}f(x)\
\mu(dx)$, that is, $(f\cdot\mu)[{X}]=\mu[f]$ and
$(f\cdot\mu)[\mathsf{A}]=\mu[f\cdot\mathds{1}_{\mathsf{A}}]$. Hence, by
definition, the measure $f\cdot\mu$ is absolutely continuous w.r.t. $\mu$,
written as $f\cdot\mu\ll\mu$, with $f$ being the corresponding density. Here,
the set of all measures on a measurable space $({X},\mathcal{X})$ will be
denoted by $\mathcal{M}({X}):=\\{\mu:\mathcal{X}\to[0,\infty]\ :\ \text{$\mu$
is a measure}\\}$, and the set of all probability measures that are absolutely
continuous w.r.t. $\mu\in\mathcal{M}({X})$ are denoted by
$\mathcal{P}(\mu):=\\{\nu\in\mathcal{M}({X})\ :\ \text{$\|\nu\|=\nu[{X}]=1$
and $\nu\ll\mu$}\\}$. Here, the Kullback-Leibler divergence between two
measures $\nu$ and $\mu$ is defined as

$D_{\rm{KL}}(\mu\ \|\ \nu)=\begin{cases}\mu[\log(f)]=\int_{X}\log(f(x))\
\mu(dx),&\text{$\mu\ll\nu$ with density $f$}\,,\\\
+\infty,&\text{otherwise}\,.\end{cases}$

If $\mu=\mathbb{P}$ is actually the specific probability measure of the
underlying probability space $\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}$,
the corresponding expectation is denoted by
$\mathbb{E}[f]:=\mathbb{P}[f]=\int_{\Omega}f(\omega)\ \mathbb{P}(d\omega)$.
Random variables, that is, measurable functions on
$\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}$, are written in Fraktur-font,
for example $\mathbf{\mathfrak{X}}$. Given two random variables, say
$\mathbf{\mathfrak{X}}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\to{X}$
and
$\mathbf{\mathfrak{Y}}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\to{Y}$,
integration of a function $f$ on ${X}\times{Y}$ w.r.t. the induced probability
measure $\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}$ is
specified by the subscript $(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})$,
that is:

$\mathbb{E}\left[f(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})\right]=\int_{\Omega}f\left(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}}\right)(\omega)\
\mathbb{P}(d\omega)=\int_{{X}\times{Y}}f(x,y)\
\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}(dx,dy)=\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}[f]=\mathbb{E}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}\left[f\right]\,.$

If we are given a regular version of the conditional distribution of
$\mathbf{\mathfrak{Y}}$, given $\mathbf{\mathfrak{X}}$, denoted by
$\mathbb{P}_{\mathbf{\mathfrak{Y}}|\mathbf{\mathfrak{X}}}$, their joint
distribution $\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}$ can
be disintegrated into the product
$\mathbb{P}_{\mathbf{\mathfrak{X}}}\otimes\mathbb{P}_{\mathbf{\mathfrak{Y}}|\mathbf{\mathfrak{X}}}$
of the marginal $\mathbb{P}_{\mathbf{\mathfrak{X}}}$ and the probability
kernel
$(x,\mathsf{B})\mapsto\mathbb{P}_{\mathbf{\mathfrak{Y}}|\mathbf{\mathfrak{X}}=x}[\mathsf{B}]$,
which allows us to use the notation:

$\displaystyle\mathbb{E}[f(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})]=\int_{X}\int_{Y}f(x,y)\
\mathbb{P}_{\mathbf{\mathfrak{Y}}|\mathbf{\mathfrak{X}}=x}(dy)\
\mathbb{P}_{\mathbf{\mathfrak{X}}}(dx)=\mathbb{E}_{\mathbf{\mathfrak{X}}}\left[\mathbb{E}_{\mathbf{\mathfrak{Y}}|\mathbf{\mathfrak{X}}=x}\left[f(x,\cdot)\right]\right]\,.$

Note that, in this case, naively changing the order of integration is not
allowed. However, if $\mathbf{\mathfrak{X}}$ and $\mathbf{\mathfrak{Y}}$ are
independent, their joint distribution
$\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}$ is given by the
product measure
$\mathbb{P}_{\mathbf{\mathfrak{X}}}\otimes\mathbb{P}_{\mathbf{\mathfrak{Y}}}$
for which Fubini’s theorem is applicable, and the iterated integration is
clarified by the corresponding subscripts
$\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}}$:

$\displaystyle\mathbb{E}\left[f(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})\right]=\left(\mathbb{P}_{\mathbf{\mathfrak{X}}}\otimes\mathbb{P}_{\mathbf{\mathfrak{Y}}}\right)[f]=\int_{X}\int_{Y}f(x,y)\
\mathbb{P}_{\mathbf{\mathfrak{Y}}}(dy)\
\mathbb{P}_{\mathbf{\mathfrak{X}}}(dx)=\mathbb{E}_{\mathbf{\mathfrak{X}}}\left[\mathbb{E}_{\mathbf{\mathfrak{Y}}}\left[f(x,\cdot)\right]|_{x=\mathbf{\mathfrak{X}}}\right]\,.$

Furthermore, we will denote the extended real numbers by
$\bar{\mathbb{R}}:=\mathbb{R}\cup\\{-\infty,+\infty\\}$. Since it holds that
$\mathcal{B}(\bar{\mathbb{R}})|_{\mathbb{R}}=\mathcal{B}(\mathbb{R})$, any
$\mathcal{B}(\mathbb{R})$-measurable function can be identified with a
$\mathcal{B}(\bar{\mathbb{R}})$-measurable function (Klenke, 2013, p. 38,
Corollary 1.87). Finally, our theoretical results will rely on the notions of
probability kernels, Polish spaces, and exponential families, whose
definitions are recalled in Appendix A.

###  2.1 Main Assumptions and Definitions

We assume that, inside the class of functions that is encoded analytically in
the model, we are given a distribution over loss-functions with a specific
structure, which is modelled by a random variable:

###### Assumption 2.1

We are given a probability space
$\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}$ and, for some fixed
$N\in\mathbb{N}$, we are given random variables
$\mathbf{\mathfrak{S}},\mathbf{\mathfrak{S}}_{1},...,\mathbf{\mathfrak{S}}_{N}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\to\Theta$,
where $\Theta$ is a Polish space. Further, we are given a non-negative and
measurable loss-function $\ell:\mathbb{R}^{n}\times\Theta\to[0,+\infty]$.

###### Remark 2.2

The assumption of a Polish space is not restrictive, yet suffices to ensure
the existence of regular versions of the conditional probability distribution,
that is, for the disintegration of a joint distribution into a marginal and a
corresponding probability kernel (see Kallenberg, 2021, Thm. 8.5).

Then, ideally, we would like to find a solution to each realization of the
random objective:

Find $x^{\prime}:\Theta\to\mathbb{R}^{n}$, s.t.
$x^{\prime}(\mathbf{\mathfrak{S}})\in\operatorname*{arg\,min}_{x\in\mathbb{R}^{n}}\ell(x,\mathbf{\mathfrak{S}})\quad
a.s.$ (1)

However, we will only solve a relaxed version of (1) and provide
_generalization bounds_ for the _average performance_. Implicitly, we assume
that we are given several realizations
$\ell(\cdot,\theta_{1}),...,\ell(\cdot,\theta_{N})$ with similar structure, as
$\theta_{1},...,\theta_{N}\overset{iid}{\sim}\mathbb{P}_{\mathbf{\mathfrak{S}}}$.
This is close to the “standard” learning problem, where an algorithm gets
trained on a data set and then applied to unseen data. For this, we need to
have a data set:

###### Definition 2.3

The measurable function
$\mathbf{\mathfrak{D}}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\longrightarrow\Theta^{N},\quad\omega\mapsto\left(\mathbf{\mathfrak{S}}_{1},...,\mathbf{\mathfrak{S}}_{N}\right)(\omega)$
is called a _data set_. Here, $\Theta^{N}$ is called the _data-space_. If the
induced distribution
$\mathbb{P}_{\mathbf{\mathfrak{D}}}=\mathbb{P}_{(\mathbf{\mathfrak{S}}_{1},...,\mathbf{\mathfrak{S}}_{N})}$
factorizes into the product of the marginals, that is, if
$\mathbb{P}_{\mathbf{\mathfrak{D}}}=\bigotimes_{i=1}^{N}\mathbb{P}_{\mathbf{\mathfrak{S}}_{i}}$,
it is called _independent_ and if, additionally, it satisfies
$\mathbb{P}_{\mathbf{\mathfrak{D}}}=\bigotimes_{i=1}^{N}\mathbb{P}_{\mathbf{\mathfrak{S}}}$,
it is called _i.i.d._

PAC-Bayesian generalization bounds involve a so-called posterior distribution,
which usually is a “data-dependent distribution”. Often, this term is left
unspecified. Yet, as also pointed out by Rivasplata et al. (2020), this is an
instance of a probability kernel (also called a “stochastic”- or “Markov
kernel”). Another commonly used name is “regular conditional probability”,
following the definition of a regular conditional distribution through
probability kernels (Catoni, 2004; Alquier, 2008). In our opinion, the notion
of a probability kernel is the least ambiguous one, leading to the following
definition:

###### Definition 2.4

Let $\mathbf{\mathfrak{D}}$ be a data set with data-space $\Theta^{N}$, and
let ${X}$ be a measurable space. A probability kernel from $\Theta^{N}$ to
${X}$ is called a _data-dependent distribution_ on ${X}$.

For solving problem (1), for every realization $\theta$ of
$\mathbf{\mathfrak{S}}$, we apply an optimization algorithm $\mathcal{A}$ to
$\ell(\cdot,\theta)$. For this, we consider a similar setting as London
(2017), that is, randomized algorithms are considered as deterministic
algorithms with randomized hyperparameters:

###### Definition 2.5

Let ${H}$ be a Polish space and $n\in\mathbb{N}$. A measurable function

$\mathcal{A}:{H}\times\mathbb{R}^{n}\times\Theta\longrightarrow\mathbb{R}^{n},\quad(\alpha,x^{(0)},\theta)\mapsto\mathcal{A}(\alpha,x^{(0)},\theta)\,,$

is called a _parametric algorithm_ and ${H}$ is called the _hyperparameter
space_ of $\mathcal{A}$. $\mathbb{R}^{n}$ is the space of the optimization
variable, $\Theta$ the space of the parameters of the problem instance, and
${H}$ the space of the hyperparameters of the algorithm.

###### Remark 2.6

$\mathcal{A}$ corresponds to the _whole_ algorithm, that is, its output is the
final iterate of the optimization procedure, which has to be stopped after
some iteration $k\in\mathbb{N}$, and could be contrasted with the notion of an
_estimator_ in statistics. However, one can also model a single update-step,
and recover $\mathcal{A}$ as k-fold concatenation of this update-step.

Here, learning $\mathcal{A}$ refers to finding a distribution $\mathbb{Q}$ on
${H}$ based on its performance on a data set $\mathbf{\mathfrak{D}}$. In the
PAC-Bayesian approach to learning and generalization, one needs a reference
distribution, which can (and should) encode prior knowledge about suitable
choices of hyperparameters, called the _prior_ :

###### Assumption 2.7

We are given a parametric optimization algorithm $\mathcal{A}$ with Polish
hyperparameter space ${H}$, and a (prior) distribution
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$ on ${H}$ that is induced by a random
variable
$\mathbf{\mathfrak{H}}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\longrightarrow{H}$,
which is independent of the data set $\mathbf{\mathfrak{D}}$ and
$\mathbf{\mathfrak{S}}$. Further, the initialization
$x^{(0)}\in\mathbb{R}^{n}$ is given and fixed.

###### Notation 2.8

We will refer to the random variable $\mathbf{\mathfrak{S}}$ as the
_parameters_ of the loss-function $\ell$ and to $\mathbf{\mathfrak{H}}$ as the
_hyperparameters_ of the algorithm $\mathcal{A}$. Further, to simplify the
notation, we also use the short-hand
$\ell(\alpha,\theta):=\ell(\mathcal{A}(\alpha,\theta),\theta)$. Furthermore,
if not needed explicitly, $x^{(0)}$ and
$\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}$ will not be mentioned anymore
in the following.

Now, we define the risk of a randomized parametric optimization algorithm as
usual:

###### Definition 2.9

Suppose $\mathbf{\mathfrak{S}}$ and $\ell$ satisfy Assumption 2.1, and
$\mathcal{A}$ satisfies Assumption 2.7. Further, let
$\mathbf{\mathfrak{D}}=(\mathbf{\mathfrak{S}}_{1},...,\mathbf{\mathfrak{S}}_{N})$
be an i.i.d. data set with distribution
$\bigotimes_{i=1}^{N}\mathbb{P}_{\mathbf{\mathfrak{S}}}$ that is independent
of $\mathbf{\mathfrak{S}}$. The _risk_ of $\mathcal{A}$ is defined as the
measurable function:

$\displaystyle\mathcal{R}:{H}\longrightarrow[0,+\infty],\
\alpha\mapsto\mathbb{E}\left[\ell\left(\mathcal{A}(\alpha,\mathbf{\mathfrak{S}}),\mathbf{\mathfrak{S}}\right)\right]=\mathbb{E}[\ell(\alpha,\mathbf{\mathfrak{S}})]=\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)]\,.$

Similarly, the _empirical risk_ of $\mathcal{A}$ on $\mathbf{\mathfrak{D}}$ is
defined as the measurable function:

$\displaystyle\hat{\mathcal{R}}:{H}\times\Theta^{N}\longrightarrow[0,+\infty],\
(\alpha,\mathbf{\mathfrak{D}})\mapsto\hat{\mathcal{R}}(\alpha,\mathbf{\mathfrak{D}})=\frac{1}{N}\sum_{i=1}^{N}\ell\left(\alpha,\mathbf{\mathfrak{S}}_{i}\right)\,.$

The following theory is based on exponential families (see Definition A.4).
This is a class of distributions on some space ${X}$, which is mainly build
upon three parts: a function $\eta$ on some index set $\Lambda$, a function
$T$ on ${X}$, and a distribution $\mu$ on ${X}$. Since we want to have a
distribution over hyperparameters $\alpha\in{H}$, we will choose ${X}={H}$.
However, then, the exponential family
$(\mathbb{Q}_{\lambda})_{\lambda\in\Lambda}$ is not data-dependent. Therefore,
we will define _data-dependent exponential families_ by extending the domain
of $T$ from ${X}$ to ${X}\times\Theta^{N}$.

###### Definition 2.10

Let $\mathbf{\mathfrak{D}}$ be a data set with data-space $\Theta^{N}$, and
let $\Lambda$ be a non-empty index set. A family of probability measures
$(\mathbb{Q}_{\lambda})_{\lambda\in\Lambda}$ on a measurable space ${X}$ is
called an _data-dependent exponential family_ (in $\eta$ and $T$), if there is
a dominating probability measure $\mu$, that is,
$(\mathbb{Q}_{\lambda})_{\lambda\in\Lambda}\subset\mathcal{P}(\mu)$, functions
$\eta:\Lambda\to\mathbb{R}^{k}$, $A:\Lambda\to(0,+\infty)$, and measurable
functions $T:{X}\times\Theta^{N}\to\mathbb{R}^{k}$, $h:{X}\to(0,+\infty)$,
such that we have
$\mathbb{Q}_{\lambda}(d)=hA(\lambda)\exp\left(\langle\eta(\lambda),T(\cdot,d)\rangle\right)\cdot\mu$
for every $\lambda\in\Lambda$, $d\in\Theta^{N}$, that is,
$\mathbb{Q}_{\lambda}(d,\mathsf{B})=\int_{\mathsf{B}}h(x)A(\lambda)\exp(\langle\eta(\lambda),T(x,d)\rangle)\
\mu(dx)$, $\mathsf{B}\in\mathcal{B}({X})$.

###### Remark 2.11

If $\eta$ describes a lower-dimensional manifold in $\mathbb{R}^{k}$,
$(\mathbb{Q}_{\lambda})_{\lambda\in\Lambda}$ is called a _curved_ exponential
family (Efron, 1975), whose properties might differ from the ones for _linear_
exponential families ($\eta(\lambda)=\lambda\in\mathbb{R}^{k}$), for example,
convexity of the map $\lambda\mapsto A(\lambda)$.

We introduce data-dependency through $T$, since it strongly affects the shape
of the distribution and, contrary to $\eta$, is defined on the underlying
space ${X}$. Since we want to learn a distribution over hyperparameters
$\alpha\in{H}$, we make the following assumption:

###### Assumption 2.12

On the hyperparameter space ${H}$, we are given a data-dependent exponential
family $(\mathbb{Q}_{\lambda})_{\lambda\in\Lambda}$ in $\eta$ and $T$ with
dominating probability measure $\mu=\mathbb{P}_{\mathbf{\mathfrak{H}}}$, such
that the map $\alpha\mapsto
h(\alpha)\exp(\langle\eta(\lambda),T(\alpha,d)\rangle)$ is non-trivial and
integrable w.r.t. $\mathbb{P}_{\mathbf{\mathfrak{H}}}$ for every
$\lambda\in\Lambda$, $d\in\Theta^{N}$, that is,
$\mathbb{E}_{\mathbf{\mathfrak{H}}}\left[h\exp(\langle\eta(\lambda),T(\cdot,d)\rangle)\right]\in(0,\infty)$.

###### Remark 2.13

  * (i)

In PAC-Bayes, the dominating measure $\mathbb{P}_{\mathbf{\mathfrak{H}}}$ is
usually referred to as _prior_ and every distribution
$\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})$ is referred to
as a _posterior_. This deviates from the standard definitions of prior and
posterior in Bayesian statistics.

  * (ii)

The integrability assumption is slightly restrictive, as it affects the choice
of $h,\eta$ and $T$. However, in Section 4 we will _construct_ $\eta$ and $T$
such that this holds anyway.

The last integral in Assumption 2.12 will be of great interest in the
following. Here, we will use a similar notation as in Barndorff–Nielsen (2014)
and denote

$\begin{split}&c(\lambda,d):=\int_{{H}}h(\alpha)\exp(\langle\eta(\lambda),T(\alpha,d)\rangle)\
\mathbb{P}_{\mathbf{\mathfrak{H}}}(d\alpha)=\mathbb{E}_{\mathbf{\mathfrak{H}}}\left\\{h\exp(\langle\eta(\lambda),T(\cdot,d)\rangle)\right\\}\,,\\\
&\kappa(\lambda,d):=\log\left(c(\lambda,d)\right)=\log\left(\mathbb{E}_{\mathbf{\mathfrak{H}}}\left\\{h\exp\left(\langle\eta(\lambda),T(\cdot,d)\rangle\right)\right\\}\right),\end{split}$
(2)

or for short, $\kappa=\log(c)$. With this notation, it holds that
$A(\lambda)=c(\lambda)^{-1}$. Through this, as shown in Lemma B.1, every
member of the data-dependent exponential family is indeed a data-dependent
distribution on ${H}$ in the sense of Definition 2.4. Finally, we will
restrict the space $\Lambda$ to a compact set. This is needed to get a uniform
bound in $\lambda$, and such ideas appeared in the literature (see Langford
and Caruana, 2001; Catoni, 2007; Alquier, 2021).

###### Assumption 2.14

$\Lambda$ is a compact set with finite covering
$\mathcal{O}:=\\{\mathsf{O}_{1},...,\mathsf{O}_{K}\\}$, that is,
$\Lambda\subset\bigcup_{i=1}^{K}\mathsf{O}_{i}$, such that there is a constant
$C_{\mathcal{O}}$, which, for every $d\in\Theta^{N}$, allows for the bound
$\max_{i=1,...,K}\
\sup_{\lambda,\lambda^{\prime}\in\mathsf{O}_{i}}\kappa(\lambda,d)-\kappa(\lambda^{\prime},d)\leq
C_{\mathcal{O}}$.

###### Remark 2.15

The existence of a general finite covering for $\Lambda$ is a consequence of
the compactness. However, the existence of the constant $C_{\mathcal{O}}$ is
non-trivial. It does hold, for example, if $\Lambda$ is a finite set
($K=|\Lambda|$, $C_{\mathcal{O}}=0$), or, if $\left(\Lambda,\rho\right)$ is a
compact metric space and $\kappa$ is Lipschitz-continuous in $\lambda$
(uniformly in $d$) with Lipschitz constant L, such that
$C_{\mathcal{O}}=L\cdot\max_{i=1,...,K}\mathrm{diam}\ \mathsf{O}_{i}$, where
the diameter of a set $\mathsf{A}$ is given by $\mathrm{diam}\
\mathsf{A}=\sup_{x,y\in\mathsf{A}}\rho(x,y)$.

##  3 General PAC-Bayesian Theorem

In this section we prove the following general PAC-Bayesian bound for data-
dependent exponential families, which then can be specialized into a
generalization bound of the learned parametric optimization algorithm
$\mathcal{A}$. It is based on the following two lemmas, whose proofs can be
found in Appendix C and D, respectively. The first lemma is a form of the
Donsker–Varadhan variational formulation and yields uniformity in the
distributions $\mathbb{Q}$:

###### Lemma 3.1

Suppose that Assumption 2.12 holds and define $\kappa$ as in (2). Then for
every $\lambda\in\Lambda$ and $d\in\Theta^{N}$:

$\kappa(\lambda,d)=\sup_{\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})}\mathbb{Q}\left[\langle\eta(\lambda),T(\cdot,d)\rangle+\log(h)\right]-D_{\rm{KL}}(\mathbb{Q}\
\|\ \mathbb{P}_{\mathbf{\mathfrak{H}}})\,.$

Furthermore, for every $\lambda\in\Lambda$, the supremum is attained at
$\mathbb{Q}_{\lambda}(d)$.

Similarly, the second lemma yields uniformity in $\lambda\in\Lambda$. For
this, we need to be able to cover $\Lambda$, and to control
$\lambda\mapsto\kappa(\lambda,d)$ for every $d\in\Theta^{N}$, as specified in
Assumption 2.14.

###### Lemma 3.2

Suppose that Assumption 2.14 holds and assume that
$\mathbb{P}\\{\kappa(\lambda,\mathbf{\mathfrak{D}})>s\\}\leq\exp(-s)$ for all
$s\in\mathbb{R}$ and $\lambda\in\Lambda$. Then
$\mathbb{P}\\{\sup_{\lambda\in\Lambda}\kappa(\lambda,\mathbf{\mathfrak{D}})\leq\log(K/\epsilon)+C_{\mathcal{O}}\\}\geq
1-\epsilon$.

Taken together, these two lemmas allow us to state the PAC-Bayesian theorem in
its general form:

###### Theorem 3.3

Suppose that Assumptions 2.12 and 2.14 hold, and assume that
$\mathbb{E}_{\mathbf{\mathfrak{D}}}[c(\lambda,\cdot)]\leq 1$ for all
$\lambda\in\Lambda$. Then, it holds that:

$\displaystyle\mathbb{P}\Bigl{\\{}$ $\displaystyle\forall\lambda\in\Lambda,\
\forall\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}}):$
$\displaystyle\mathbb{Q}\left[\langle\eta(\lambda),T(\cdot,d)\rangle+\log(h)\right]|_{d=\mathbf{\mathfrak{D}}}\leq
D_{\rm{KL}}(\mathbb{Q}\ \|\
\mathbb{P}_{\mathbf{\mathfrak{H}}})+\log\left(\frac{K}{\varepsilon}\right)+C_{\mathcal{O}}\Bigr{\\}}\geq
1-\varepsilon\,.$

Proof  For every $\lambda\in\Lambda$, $c(\lambda,\mathbf{\mathfrak{D}})$ is a
non-negative random variable, and $x\mapsto\exp(x)$ is monotonically
increasing. Thus, Markov’s inequality yields for $s\in\mathbb{R}$:

$\mathbb{P}\left\\{c(\lambda,\mathbf{\mathfrak{D}})>\exp(s)\right\\}\leq\frac{\mathbb{E}[c(\lambda,\mathbf{\mathfrak{D}})]}{\exp(s)}\leq\exp(-s)\,.$

Thus,
$\mathbb{P}\left\\{\kappa(\lambda,\mathbf{\mathfrak{D}})>s\right\\}\leq\exp(-s)$,
because $c(\lambda,\mathbf{\mathfrak{D}})>\exp(s)$ is equivalent to
$\kappa(\lambda,\mathbf{\mathfrak{D}})=\log(c(\lambda,\mathbf{\mathfrak{D}}))>s$.
Hence, Lemma 3.2 is applicable and gives:

$\mathbb{P}\left\\{\sup_{\lambda\in\Lambda}\kappa(\lambda,\mathbf{\mathfrak{D}})\leq\log\left(\frac{K}{\varepsilon}\right)+C_{\mathcal{O}}\right\\}\geq
1-\varepsilon\,.$

Using Lemma 3.1 gives:

$\mathbb{P}\left\\{\sup_{\lambda\in\Lambda}\sup_{\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})}\mathbb{Q}[\langle\eta(\lambda),T(\cdot,d)\rangle+\log(h)]|_{d=\mathbf{\mathfrak{D}}}-D_{\rm{KL}}(\mathbb{Q}\
\|\
\mathbb{P}_{\mathbf{\mathfrak{H}}})\leq\log\Bigl{(}\frac{K}{\varepsilon}\Bigr{)}+C_{\mathcal{O}}\right\\}\geq
1-\varepsilon\,.$

Simply rearranging and reformulating yields the result.

###### Remark 3.4

  * (i)

Note that the statement is still true for a data-dependent prior
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$: Denoting the additional data set by
$\mathbf{\mathfrak{D}}^{\prime}$, one needs to assume that
$\mathbb{E}\left[c(\lambda,(\mathbf{\mathfrak{D}},\mathbf{\mathfrak{D}}^{\prime}))\right]\leq
1$. Also Lemma 3.1 still applies “pointwise” (with the appropriate notational
adjustments). Intuitively, this is due to the fact that the PAC-bound is a
statement about relative, and not absolute values.

  * (ii)

In Section 4 we provide sufficient conditions s.t.
$\mathbb{E}[c(\lambda,\mathbf{\mathfrak{D}})]\leq 1$ holds for all
$\lambda>0$.

  * (iii)

$\log(K)$ bears the intrinsic dimension of $\Lambda$ and thus, in full
generality, might be large. For our generalization bound, however, it has only
a minor influence, since there $\Lambda\subset\mathbb{R}$, and the the
empirical risk is typically much larger.

For the rest of the paper, we will have $h\equiv 1$, such that $\log(h)\equiv
0$. The following corollary, whose proof is given in Appendix E, shows how to
transform Theorem 3.3 into a high-probability bound on the risk.

######  Corollary 3.5 (PAC-Bayesian Generalization Bound)

Denote the coordinates of $T$ and $\eta$ by
$T=(T^{(1)},T^{(2)},...,T^{(k)})=:(T^{(1)},T^{(r)})$ and
$\eta=(\eta^{(1)},\eta^{(2)},...,\eta^{(k)})=:(\eta^{(1)},\eta^{(r)})$. If
$T^{(1)}=\mathcal{R}-\hat{\mathcal{R}}$ and $\eta^{(1)}>0$, the following are
equivalent for any $\lambda\in\Lambda$, $d\in\Theta^{N}$,
$\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})$:

  * (i)

$\mathbb{Q}\left[\langle\eta(\lambda),T(\cdot,d)\rangle\right]\leq
D_{\rm{KL}}(\mathbb{Q}\ \|\
\mathbb{P}_{\mathbf{\mathfrak{H}}})+\log\left(\frac{K}{\varepsilon}\right)+C_{\mathcal{O}}$,

  * (ii)

$\mathbb{Q}[\mathcal{R}]\leq\mathbb{Q}[\hat{\mathcal{R}}(\cdot,d)]+\frac{1}{\eta^{(1)}(\lambda)}\left(D_{\rm{KL}}(\mathbb{Q}\
\|\
\mathbb{P}_{\mathbf{\mathfrak{H}}})+\log\left(\frac{K}{\varepsilon}\right)+C_{\mathcal{O}}-\mathbb{Q}\left[\langle\eta^{(r)}(\lambda),T^{(r)}(\cdot,d)\rangle\right]\right)\,.$

In particular, if Theorem 3.3 applies, we can replace (i) with (ii).

##  4 Learning-to-Optimize with Guarantees

Here, we consider properties of optimization algorithms, that assert the
necessary condition $\mathbb{E}[c(\lambda,\mathbf{\mathfrak{D}})]\leq 1$ for
all $\lambda\in\Lambda$ to employ the PAC-Bayesian bound from Section 3.

###  4.1 Guaranteed Convergence

In the next theorem, the additional assumption on $\mathcal{A}$ is sufficient
to ensure the assumptions of Theorem 3.3. Essentially, it requires the loss of
the algorithm’s output to be bounded. Yet, as shown in Section 4.2, it is too
restrictive to learn hyperparameters that allow for a significant acceleration
compared to the standard choices from a worst-case analysis.

###### Theorem 4.1

Suppose that $\mathbf{\mathfrak{S}}$ and $\ell$ satisfy Assumption 2.1, and
suppose that $\mathcal{A}$ satisfies Assumption 2.7. Further, assume that
there is a constant $C\geq 0$ and a measurable function
$\rho:{H}\longrightarrow[0,\infty)$, such that for every $\alpha\in{H}$ it
holds that $\ell(\alpha,\cdot)\leq C\rho(\alpha)\ell(x^{(0)},\cdot)$
$\mathbb{P}_{\mathbf{\mathfrak{S}}}$-a.s. Furthermore, let
$\mathbf{\mathfrak{D}}$ be a corresponding i.i.d. data set of size
$N\in\mathbb{N}$. Finally, assume that
$\mathbb{E}\bigl{[}\ell(x^{(0)},\mathbf{\mathfrak{S}})^{2}\bigr{]}<\infty$,
and define $\eta:(0,\infty)\longrightarrow\mathbb{R}^{2}$ and
$T:{H}\times\Theta^{N}\longrightarrow\mathbb{R}^{2}$ through:

$\displaystyle\eta(\lambda):=\left(\lambda,\ -\frac{\lambda^{2}}{2}\right)\,,\
\ \ T(\alpha,d):=\left(\mathcal{R}(\alpha)-\hat{\mathcal{R}}(\alpha,d),\
\frac{\rho^{2}(\alpha)C^{2}}{N}\mathbb{E}\left[\ell(x^{(0)},\mathbf{\mathfrak{S}})^{2}\right]\right)\,.$

Then it holds that $\mathbb{E}[c(\lambda,\mathbf{\mathfrak{D}})]\leq 1$ for
all $\lambda>0$.

Proof  Since $\mathbf{\mathfrak{H}}$ and $\mathbf{\mathfrak{D}}$ are
independent, their joint distribution is given by the product measure
$\mathbb{P}_{\mathbf{\mathfrak{D}}}\otimes\mathbb{P}_{\mathbf{\mathfrak{H}}}$.
Further, $\exp\geq 0$, such that Fubini’s theorem is applicable and allows to
do the integration iteratively, that is:

$\mathbb{E}\left[\exp\left(\lambda(\mathcal{R}(\mathbf{\mathfrak{H}})-\hat{\mathcal{R}}(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{D}}))\right)\right]=\mathbb{E}\left\\{\mathbb{E}\left[\exp\left(\lambda(\mathcal{R}(\alpha)-\hat{\mathcal{R}}(\alpha,\mathbf{\mathfrak{D}}))\right)\right]\Bigl{|}_{\alpha=\mathbf{\mathfrak{H}}}\right\\}\,.$

Hence, first consider the inner integral for a fixed $\alpha\in{H}$. Then, by
definition and the i.i.d. assumption one gets:

$\displaystyle\mathbb{E}\left[\exp\left(\lambda(\mathcal{R}(\alpha)-\hat{\mathcal{R}}(\alpha,\mathbf{\mathfrak{D}}))\right)\right]$
$\displaystyle=\mathbb{E}\left\\{\exp\left(-\frac{\lambda}{N}\sum_{i=1}^{N}\left(\ell(\alpha,\mathbf{\mathfrak{S}}_{i})-\mathbb{E}[\ell(\alpha,\mathbf{\mathfrak{S}})]\right)\right)\right\\}$
$\displaystyle=\prod_{i=1}^{N}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left\\{\exp\left(-\frac{\lambda}{N}\left(\ell(\alpha,\cdot)-\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)]\right)\right)\right\\}\,.$

The loss-function is non-negative and ,by assumption, $\mathcal{A}$ can be
bounded $\mathbb{P}_{\mathbf{\mathfrak{S}}}$-a.s. Thus, for every
$\alpha\in{H}$, $\ell(\alpha,\mathbf{\mathfrak{S}})$ is a non-negative random
variable with finite second-moment, as
$\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)^{2}]\leq
C^{2}\rho(\alpha)^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[\ell(x^{(0)},\cdot)^{2}\right]<\infty$.
Thus, by Lemma B.2, the boundedness assumption on $\mathcal{A}$, and the
monotonicity of the exponential function, one gets the bound:

$\displaystyle\mathbb{E}_{\mathbf{\mathfrak{S}}}\left\\{\exp\left(-\frac{\lambda}{N}\left(\ell(\alpha,\cdot)-\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)]\right)\right)\right\\}$
$\displaystyle\leq\exp\Bigl{(}\frac{\lambda^{2}}{2N^{2}}\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)^{2}]\Bigr{)}$
$\displaystyle\leq\exp\Bigl{(}\frac{\lambda^{2}}{2N^{2}}C^{2}\rho(\alpha)^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(x^{(0)},\cdot)^{2}]\Bigr{)}\,.$

Therefore we have the following bound:

$\displaystyle\mathbb{E}\left[\exp\left(\lambda(\mathcal{R}(\alpha)-\hat{\mathcal{R}}(\alpha,\mathbf{\mathfrak{D}}))\right)\right]$
$\displaystyle\leq\prod_{i=1}^{N}\exp\Bigl{(}\frac{\lambda^{2}}{2N^{2}}C^{2}\rho(\alpha)^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[\ell(x^{(0)},\cdot)^{2}\right]\Bigr{)}$
$\displaystyle=\exp\Bigl{(}\frac{\lambda^{2}}{2N}C^{2}\rho(\alpha)^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[\ell(x^{(0)},\cdot)^{2}\right]\Bigr{)}\,.$

This can be rearranged into
$\mathbb{E}\left\\{\exp\left(\lambda(\mathcal{R}(\alpha)-\hat{\mathcal{R}}(\alpha,\mathbf{\mathfrak{D}}))-\frac{\lambda^{2}}{2}\frac{C^{2}}{N}\rho(\alpha)^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(x^{(0)},\cdot)^{2}]\right)\right\\}\leq
1$, as the right-hand side does not depend on $\mathbf{\mathfrak{D}}$. Since
$\mathbf{\mathfrak{H}}$ and $\mathbf{\mathfrak{D}}$ are independent, and
$\alpha\in{H}$ was arbitrary, this inequality does hold
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s. Therefore, one directly gets the
bound
$\mathbb{E}\left\\{\exp\left(\lambda(\mathcal{R}(\mathbf{\mathfrak{H}})-\hat{\mathcal{R}}(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{D}}))-\frac{\lambda^{2}}{2}\frac{C^{2}}{N}\rho(\mathbf{\mathfrak{H}})^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(x^{(0)},\cdot)^{2}]\right)\right\\}\leq
1$. Now, again by Fubini’s theorem, one can also switch the order of
integration to get:

$\mathbb{E}\left\\{\mathbb{E}\left[\exp\left(\lambda(\mathcal{R}(\mathbf{\mathfrak{H}})-\hat{\mathcal{R}}(\mathbf{\mathfrak{H}},d))-\frac{\lambda^{2}}{2}\frac{C^{2}}{N}\rho(\mathbf{\mathfrak{H}})^{2}\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(x^{(0)},\cdot)^{2}]\right)\right]\Bigl{|}_{d=\mathbf{\mathfrak{D}}}\right\\}\leq
1\,.$

Inserting the definition of $\eta$ and $T$ gives
$\mathbb{E}\left\\{\mathbb{E}\left[\exp\left(\langle\eta(\lambda),T(\mathbf{\mathfrak{H}},d)\rangle\right)\right]|_{d=\mathbf{\mathfrak{D}}}\right\\}\leq
1$. Here, the inner term is the same as
$\mathbb{E}\left[\exp\left(\langle\eta(\lambda),T(\mathbf{\mathfrak{H}},d)\rangle\right)\right]=\int_{H}\exp\left(\langle\eta(\lambda),T(\alpha,d))\rangle\right)\
\mathbb{P}_{\mathbf{\mathfrak{H}}}(d\alpha)=c(\lambda,d)$. Hence, this is the
same as $\mathbb{E}[c(\lambda,\mathbf{\mathfrak{D}})]\leq 1$.

###### Remark 4.2

The argument still works for a data-dependent prior, if the corresponding data
sets $\mathbf{\mathfrak{D}}^{\prime}$ and $\mathbf{\mathfrak{D}}$ are
independent: While interchanging the integration w.r.t.
$\mathbf{\mathfrak{D}}^{\prime}$ and $\mathbf{\mathfrak{H}}$ is not allowed,
an interchange w.r.t. $\mathbf{\mathfrak{H}}$ and $\mathbf{\mathfrak{D}}$ is
still valid (under the integral), that is, for a function $f$ it would hold
$\mathbb{E}\left[f(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{D}},\mathbf{\mathfrak{D}}^{\prime})\right]=\mathbb{E}_{\mathbf{\mathfrak{D}}^{\prime}}\left\\{\mathbb{E}_{\mathbf{\mathfrak{H}}|\mathbf{\mathfrak{D}}^{\prime}=d^{\prime}}\left[\mathbb{E}_{\mathbf{\mathfrak{D}}}\left[f(\alpha,\cdot,d^{\prime})\right]|_{\alpha=\mathbf{\mathfrak{H}}}\right]\right\\}=\mathbb{E}_{\mathbf{\mathfrak{D}}^{\prime}}\left\\{\mathbb{E}_{\mathbf{\mathfrak{D}}}\left[\mathbb{E}_{\mathbf{\mathfrak{H}}|\mathbf{\mathfrak{D}}^{\prime}=d^{\prime}}\left[f(\cdot,d,d^{\prime})\right]|_{d=\mathbf{\mathfrak{D}}}\right]\right\\}$,
and the inner term is $\leq 1$ in any case. The premise here is that the
boundedness assumption on $\mathcal{A}$ holds.

###  4.2 Conditioning on Convergence

Most of the time, the previous approach learns hyperparameters that ensure
convergence, since the boundedness assumption on $\mathcal{A}$ implicitly
requires weak theoretical worst-case estimates _almost surely_. If, for
example, the considered class of functions
$(\ell(\cdot,\theta))_{\theta\in\Theta}$ is that of general quadratic
functions, the convergence behaviour is accurately represented by analytic
quantities from a worst-case analysis. Thereby, the boundedness prevents
“aggressive” step-size parameters that lie outside the worst-case convergence
regime, as they would lead to a diverging behaviour, which increases the
incurred empirical risk dramatically. Thus, to motivate the upcoming
discussion, consider the following thought-experiment:

###### Example 4.3

Consider $\ell(x,\theta):=\frac{\theta}{2}x^{2}$ and assume that the chosen
algorithm is gradient descent, that is
$x^{(k+1)}=x^{(k)}-\tau\ell^{\prime}(x^{(k)},\theta)$. For a given $\theta$,
the optimal step-size is $\tau=\frac{1}{\theta}$, which gives convergence in
one step. However, if $\theta$ is given by samples from the distribution
$\mathbb{P}_{\mathbf{\mathfrak{S}}}=0.99\delta_{1}+0.01\delta_{100}$, a worst-
case analysis would suggest to take $\tau_{w}=\frac{1}{100}$. In this case, we
would have an algorithm that converges in a single step for 1% of the problem
instances, while having a linear convergence rate of $(\frac{99}{100})^{k}$
for the other 99%. Another choice is to take $\tau_{d}=1$, which leads to an
algorithm that does converge in a single step for 99% of the problem
instances, but diverges in 1% of the cases. By restricting to the 99% of the
cases where convergence does occur, the overall difference in speed is
drastic.

Hence, in this section, a different approach is taken: We actually allow for
divergence, if it only occurs in rare cases with a controllable probability,
that is, “almost surely” is relaxed to “with a sufficiently large
probability”. Essentially, we only consider the loss for all those
hyperparameters, where the loss is bounded by a certain constant, as well as
the probability for that to occur. Then, in Section 4.2.1, we develop a
technique that allows the user to actually control this probability. Clearly,
a stronger guarantee trades for convergence speed.

###### Definition 4.4

Given a measurable function $g:\Theta\to\mathbb{R}$, the (parametric)
_sublevel set_ $\mathsf{L_{g}}\subset{H}\times\Theta$ is defined as
$\mathsf{L_{g}}=\\{(\alpha,\theta)\in{H}\times\Theta\ :\
\ell(\alpha,\theta)\leq g(\theta)\\}$. The sections of $\mathsf{L_{g}}$ for
fixed $\alpha\in{H}$ will be denoted by $\mathsf{L_{g,\alpha}}$.

###### Remark 4.5

The concrete choice of $g$ will influence the PAC-Bayesian bound. Again, for
proving the result, the essential property is that the loss after applying
$\mathcal{A}$ can be bounded.

In Lemma B.3 we show that $\mathsf{L_{g}}$ is indeed a measurable set. This is
not obvious, as the loss function and the algorithm are composed in a non-
standard way. This result further implies that the sections
$\mathsf{L_{g,\alpha}}$ are measurable, too. Since ${H}$ and $\Theta$ are
Polish spaces, the product ${H}\times\Theta$ is again Polish. Hence, there
exists a regular version of the conditional probability of
$\mathbf{\mathfrak{S}}$, given $\mathbf{\mathfrak{H}}$, that is, a kernel
${H}\to\Theta,\
(\alpha,\mathsf{B})\mapsto\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{B}]$.
By Witting (2013, Thm. 1.122, p.124), this determines a regular version of the
conditional probability of $(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})$,
given $\mathbf{\mathfrak{H}}$, through ${H}\to{H}\times\Theta,\
(\alpha,\mathsf{B})\mapsto\mathbb{P}_{(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{B}]:=\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{B}_{\alpha}]$.
Here, for any measurable function $f:{H}\times\Theta\to\mathbb{R}$ for which
$\mathbb{E}[f(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})]$ exists, it holds
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s.:

$\mathbb{E}\left\\{f(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})\ |\
\mathbf{\mathfrak{H}}=\alpha\right\\}=\int_{\Theta}f(\alpha,\theta)\
\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}(d\theta)\,.$
(3)

Thus, we get $\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s. the equality
$\mathbb{P}\\{(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})\in\mathsf{B}\ |\
\mathbf{\mathfrak{H}}=\alpha\\}=\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{B}_{\alpha}]$.
In particular, this applies to the sublevel set $\mathsf{L_{g}}$, and the map
$\alpha\mapsto\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{L_{g,\alpha}}]$
is measurable.

###### Definition 4.6

Let $\mathsf{L_{g}}$ be a parametric sublevel set. Define the _sublevel
probability_ as the measurable function $\alpha\mapsto
p(\alpha):=\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{L_{g,\alpha}}]$.

###### Lemma 4.7

Suppose Assumption 2.7 holds, and let
$\mathbb{P}_{\mathbf{\mathfrak{S}}}[\mathsf{L_{g,\alpha}}]>0$ for every
$\alpha\in{H}$. Then we have $\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s.:

  * (i)

$p(\alpha)=\mathbb{P}_{\mathbf{\mathfrak{S}}}[\mathsf{L_{g,\alpha}}]$,

  * (ii)

$\mathbb{E}\left\\{\ell(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})\cdot\mathds{1}_{\mathsf{L_{g}}}(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})\
|\
\mathbf{\mathfrak{H}}=\alpha\right\\}=\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)\mathds{1}_{\mathsf{L_{g,\alpha}}}]=p(\alpha)\cdot\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)\
|\ \mathsf{L_{g,\alpha}}]$.

Proof  By (3) and the independence of $\mathbf{\mathfrak{S}}$ and
$\mathbf{\mathfrak{H}}$, we have
$\mathbb{E}\left\\{\ell(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})\mathds{1}_{\mathsf{L_{g}}}(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{S}})\
|\
\mathbf{\mathfrak{H}}=\alpha\right\\}=\int_{\Theta}\ell(\alpha,\theta)\mathds{1}_{\mathsf{L_{g,\alpha}}}(\theta)\
\mathbb{P}_{\mathbf{\mathfrak{S}}}(d\theta)=\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)\mathds{1}_{\mathsf{L_{g,\alpha}}}]$
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s., which shows the first equality of
(ii). Since $\mathbb{P}_{\mathbf{\mathfrak{S}}}[\mathsf{L_{g,\alpha}}]>0$, the
elementary conditional expectation is defined as
$\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)\ |\
\mathsf{L_{g,\alpha}}]=\frac{\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)\mathds{1}_{\mathsf{L_{g,\alpha}}}]}{\mathbb{P}_{\mathbf{\mathfrak{S}}}[\mathsf{L_{g,\alpha}}]}$.
Again by independence we have $\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s. the
equality
$\mathbb{P}_{\mathbf{\mathfrak{S}}}[\mathsf{L_{g,\alpha}}]=\mathbb{P}_{\mathbf{\mathfrak{S}}|\mathbf{\mathfrak{H}}=\alpha}[\mathsf{L_{g,\alpha}}]=p(\alpha)$,
which shows (i) and the second equality of (ii).

This construction allows us to give a more fine-grained analysis of the
algorithm, as it allows to trade the boundedness assumption for the sublevel
probability. This basically extends a worst-case analysis, which would
correspond to an uniform upper bound. Motivated by Lemma 4.7, to actually
analyze the algorithm in this setting, we define the _sublevel risk_ and its
empirical counterpart as the expect loss _conditioned_ on the sublevel set:

###### Definition 4.8

Let $\mathsf{L_{g}}$ be a parametric sublevel set. Then the _sublevel risk_
$\mathcal{R}_{g}:{H}\to[0,+\infty]$ is defined as the conditional expectation
of the loss given $\mathsf{L_{g,\alpha}}$:

$\displaystyle\alpha\mapsto\mathcal{R}_{g}(\alpha):=\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[\ell(\alpha,\cdot)\
|\
\mathsf{L_{g,\alpha}}\right]=\begin{cases}\frac{1}{p(\alpha)}\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)\mathds{1}_{\mathsf{L_{g,\alpha}}}],&\text{if
$p(\alpha)>0$}\,;\\\ 0,&\text{otherwise}\,.\end{cases}$

Given a data set
$\mathbf{\mathfrak{D}}=(\mathbf{\mathfrak{S}}_{1},...,\mathbf{\mathfrak{S}}_{N})$,
the _empirical sublevel risk_
$\hat{\mathcal{R}}_{g}:{H}\times\Theta^{N}\to[0,+\infty]$ is defined as
$(\alpha,\mathbf{\mathfrak{D}})\mapsto\hat{\mathcal{R}}_{g}(\alpha,\mathbf{\mathfrak{D}}):=\frac{1}{p(\alpha)}\frac{1}{N}\sum_{i=1}^{N}\mathds{1}_{\mathsf{L_{g,\alpha}}}(\mathbf{\mathfrak{S}}_{i})\ell(\alpha,\mathbf{\mathfrak{S}}_{i})$.

The following theorem is a direct generalization of Theorem 4.1. Especially,
note that the additional assumption on $\mathcal{A}$ is not needed anymore.

###### Theorem 4.9

Suppose that $\mathbf{\mathfrak{S}}$ and $\ell$ satisfy Assumption 2.1, and
suppose that $\mathcal{A}$ satisfies Assumption 2.7. Further, let
$\mathbf{\mathfrak{D}}$ be a corresponding i.i.d. data set of size
$N\in\mathbb{N}$, and let $\mathsf{L_{g}}$ be a parametric sublevel set with
sublevel probability $p$. Assume that
$\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{p>0\\}=1$ and
$\mathbb{E}_{\mathbf{\mathfrak{S}}}[g^{2}]<\infty$. Define
$\eta:(0,\infty)\longrightarrow\mathbb{R}^{2}$ and
$T:{H}\times\Theta^{N}\longrightarrow\mathbb{R}^{2}$ as

$\displaystyle\eta(\lambda):=\left(\lambda,\ -\frac{\lambda^{2}}{2}\right)\,,\
\ \
T(\alpha,d):=\left(\mathcal{R}_{g}(\alpha)-\hat{\mathcal{R}}_{g}(\alpha,d),\
\frac{1}{p(\alpha)^{2}N}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[g^{2}\mathds{1}_{\mathsf{L_{g,\alpha}}}\right]\right)\,.$

Then, for all $\lambda>0$, it holds that
$\mathbb{E}[c(\lambda,\mathbf{\mathfrak{D}})]\leq 1$.

Proof  The proof is very similar to the proof of Theorem 4.1 and basically
uses the same reasoning. Let
$\ell_{g}(\alpha,\theta):=\mathds{1}_{\mathsf{L_{g,\alpha}}}(\theta)\ell(\alpha,\theta)$.
Since $\mathbf{\mathfrak{H}}$ and $\mathbf{\mathfrak{D}}$ are independent, and
$\exp\geq 0$, one gets from Fubini’s theorem:

$\mathbb{E}\left[\exp(\lambda(\mathcal{R}_{g}(\mathbf{\mathfrak{H}})-\hat{\mathcal{R}}_{g}(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{D}})))\right]=\mathbb{E}\left\\{\mathbb{E}\left[\exp(\lambda(\mathcal{R}_{g}(\alpha)-\hat{\mathcal{R}}_{g}(\alpha,\mathbf{\mathfrak{D}})))\right]\Bigl{|}_{\alpha=\mathbf{\mathfrak{H}}}\right\\}\,.$

Thus, first consider a fixed $\alpha\in{H}$ with $p(\alpha)>0$. Then, by
definition and the i.i.d. assumption, it holds that:

$\displaystyle\mathbb{E}\left[\exp(\lambda(\mathcal{R}_{g}(\alpha)-\hat{\mathcal{R}}_{g}(\alpha,\mathbf{\mathfrak{D}})))\right]$
$\displaystyle=\mathbb{E}\left\\{\exp\left(-\frac{\lambda}{Np(\alpha)}\sum_{i=1}^{N}\left(\ell_{g}(\alpha,\mathbf{\mathfrak{S}}_{i})-\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell_{g}(\alpha,\cdot)]\right)\right)\right\\}$
$\displaystyle=\prod_{i=1}^{N}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left\\{\exp\left(-\frac{\lambda}{Np(\alpha)}\left(\ell_{g}(\alpha,\cdot)-\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell_{g}(\alpha,\cdot)]\right)\right)\right\\}\,.$

$\ell_{g}(\alpha,\cdot)$ is non-negative, and by definition of the parametric
sublevel set has a finite second-moment, that is
$\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell_{g}(\alpha,\cdot)^{2}]\leq\mathbb{E}_{\mathbf{\mathfrak{S}}}[g^{2}\mathds{1}_{\mathsf{L_{g,\alpha}}}]<\infty$.
Hence, by Lemma B.2 we have the inequality
$\mathbb{E}_{\mathbf{\mathfrak{S}}}\left\\{\exp\left(-\frac{\lambda}{Np(\alpha)}\left(\ell_{g}(\alpha,\cdot)-\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell_{g}(\alpha,\cdot)]\right)\right)\right\\}\leq\exp\left(\frac{\lambda^{2}}{2N^{2}p(\alpha)^{2}}\mathbb{E}\left[\ell_{g}(\alpha,\cdot)^{2}\right]\right)$.
Thus:

$\mathbb{E}\left[\exp(\lambda(\mathcal{R}_{g}(\alpha,\cdot)-\hat{\mathcal{R}}_{g}(\alpha,\mathbf{\mathfrak{D}})))\right]\leq\exp\left(\frac{\lambda^{2}}{2Np(\alpha)^{2}}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[g^{2}\mathds{1}_{\mathsf{L_{g,\alpha}}}\right]\right)\,.$

This can be rearranged into
$\mathbb{E}\left\\{\exp\left(\lambda(\mathcal{R}_{g}(\alpha)-\hat{\mathcal{R}}_{g}(\alpha,\mathbf{\mathfrak{D}}))-\frac{\lambda^{2}}{2Np(\alpha)^{2}}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[g^{2}\mathds{1}_{\mathsf{L_{g,\alpha}}}\right]\right)\right\\}\leq
1$, since the right-hand side is independent of $\mathbf{\mathfrak{D}}$. As
this holds for any $\alpha$ with $p(\alpha)>0$, which in turn does hold
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$-a.s., we get
$\mathbb{E}\left\\{\exp\left(\lambda(\mathcal{R}_{g}(\mathbf{\mathfrak{H}})-\hat{\mathcal{R}}_{g}(\mathbf{\mathfrak{H}},\mathbf{\mathfrak{D}}))-\frac{\lambda^{2}}{2Np(\mathbf{\mathfrak{H}})^{2}}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[g^{2}\mathds{1}_{\mathsf{L_{g,\alpha}}}\right]|_{\alpha=\mathbf{\mathfrak{H}}}\right)\right\\}\leq
1$. Changing the order of integration with Fubini’s theorem, we get:

$\mathbb{E}\left\\{\mathbb{E}\left[\exp\left(\lambda(\mathcal{R}_{g}(\mathbf{\mathfrak{H}})-\hat{\mathcal{R}}_{g}(\mathbf{\mathfrak{H}},d))-\frac{\lambda^{2}}{2Np(\mathbf{\mathfrak{H}})^{2}}\mathbb{E}_{\mathbf{\mathfrak{S}}}\left[g^{2}\mathds{1}_{\mathsf{L_{g,\alpha}}}\right]|_{\alpha=\mathbf{\mathfrak{H}}}\right)\right]\Bigl{|}_{d=\mathbf{\mathfrak{D}}}\right\\}\leq
1\,.$

Using the definition of $\eta$ and $T$, this is the same as
$\mathbb{E}\left\\{\mathbb{E}\left[\exp(\langle\eta(\lambda),T(\mathbf{\mathfrak{H}},d)\rangle)\right]|_{d=\mathbf{\mathfrak{D}}}\right\\}\leq
1\,.$ The inner term can be rewritten as
$\mathbb{E}\left[\exp(\langle\eta(\lambda),T(\mathbf{\mathfrak{H}},d)\rangle)\right]=\int_{H}\exp(\langle\eta(\lambda),T(\alpha,d)\rangle)\
\mathbb{P}_{\mathbf{\mathfrak{H}}}(d\alpha)=c(\lambda,d)$. Hence, in total we
get $\mathbb{E}\left[c(\lambda,\mathbf{\mathfrak{D}})\right]\leq 1$.

###### Remark 4.10

Since $\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{p>0\\}=1$, under the prior the
algorithm should be able to “reach” the sublevel set. This is a constraint on
the _support_ (see Definition A.5) of $\mathbb{P}_{\mathbf{\mathfrak{H}}}$,
which is not satisfied without further ado. Section 5 provides a construction
for achieving this.

####  4.2.1 Implementing the Non-divergence – Speed Trade-Off

In Section 4.2, care has to be taken in the choice of the prior
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$: Just minimizing the upper bound as much
as possible can lead to a neglection of a high sublevel probability, that is,
the algorithm is especially fast on a small subset of the parameters, while it
diverges for the rest, because the term $\frac{1}{p(\alpha)}$ might not
compensate for the smaller sublevel risk. Thus, if a certain sublevel
probability $\varepsilon_{\rm{conv}}\in[0,1]$ has to be ensured, one has to
enforce it. We propose to use absolute continuity:

###### Lemma 4.11

Let $\varepsilon_{\rm{conv}}\in[0,1]$ and assume that
$p(\mathbf{\mathfrak{H}})\geq\varepsilon_{\rm{conv}}$ a.s. Then, for every
$\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})$ we have
$\mathbb{Q}\\{p<\varepsilon_{\rm{conv}}\\}=0$.

Proof  Since
$\\{p\geq\varepsilon_{\rm{conv}}\\}=\\{p<\varepsilon_{\rm{conv}}\\}^{c}$, we
have $\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{p<\varepsilon_{\rm{conv}}\\}=0$.
Thus, the result follows by definition of absolute continuity.
Though the proof is trivial, this lemma has a very important consequence,
which we want to stress: If one can guarantee that a required property is
satisfied for the prior, it will be _preserved_ during the PAC-Bayesian
learning process, that is, if the prior only puts mass on hyperparameters that
ensure a certain sublevel probability, the posterior will do the same. We will
enforce this constraint in the learning procedure, which is discussed next.

##  5 Learning Procedure

$\alpha$$f(\alpha)$ 1) Find initialization by imitation.
$\alpha$$f(\alpha)$$supp(\mathbb{P}_{\mathbf{\mathfrak{H}}})$ 2) Locating the
prior. $\alpha$$f(\alpha)$ 3) Construct the prior. $\alpha$$f(\alpha)$ 4)
Compute the posterior. Figure 2: Learning procedure: 1) Imitation learning. 2)
Probabilistically constrained stochastic empirical risk minimization. 3)
Construct prior through sampling. 4) Compute posterior by performing the PAC-
Bayesian learning step.

This section deals with the learning procedure, that is, how the abstract
framework discussed in Sections 3 and 4 is actually implemented. Hence, this
is also the beginning of the second part of the paper, which is of a more
practical kind and less theoretical. The resulting learning procedure is
visualized in Figure 2 and consists of four steps:

  * (i)

Find an initialization that is trainable in a stable way: This is due to the
fact that one might include, for example, a neural network in the update step
of $\mathcal{A}$, which, if initialized randomly, might predict points at the
very beginning, that are so far off that one encounters numerical
instabilities. Hence, we first train the algorithm to “follow” another
algorithm $\mathcal{A}^{\prime}$, just for the purpose of non-divergence.

  * (ii)

Locate the prior by finding a point $\alpha^{(0)}\in{H}$ that satisfies the
constraint in Subsection 4.2 and yields a good performance. For this, we
perform a constrained version of stochastic empirical risk minimization with a
new, specifically designed loss function.

  * (iii)

Starting from $\alpha^{(0)}$, construct the prior distribution by running a
constrained version of a sampling algorithm.

  * (iv)

Compute the posterior distribution by finding the optimal
$\lambda^{*}\in\Lambda$ and performing the reweighting with the closed form
for $\mathbb{Q}_{\lambda^{*}}$.

Hence, first, in Subsection 5.1 we will identify the optimal posterior
$\mathbb{Q}^{*}$ in the abstract setting. Second, in Section 5.2, we realize
the constraining procedure, which allows to keep an optimization or sampling
procedure inside an abstractly defined set $\mathsf{A}$. Lastly, we put things
together in Section 5.3 and finalize the learning procedure by concretizing
each single step: Section 5.3.1 deals with the pre-computation phase in (i),
Section 5.3.2 deals with the initialization in (ii), Section 5.3.3 provides
the concrete choice of prior distribution in (iii), and Section 5.3.4 finally
yields the posterior distribution in (iv).

###  5.1 Minimization of the PAC-Bound

Because learning is phrased as minimizing the PAC-Bayesian upper-bound, this
is the first thing to do. Hence, in this subsection we consider $\eta$, $T$
and item $(ii)$ from Corollary 3.5. We seek for $\lambda\in\Lambda$ and
$\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})$ that minimize
the upper-bound in Corollary 3.5 $(ii)$, that is, we want to solve:

$\inf_{\lambda\in\Lambda}\
\inf_{\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})}\
\mathbb{Q}[\hat{\mathcal{R}}(\cdot,d)]+\frac{\left(D_{\rm{KL}}(\mathbb{Q}\ \|\
\mathbb{P}_{\mathbf{\mathfrak{H}}})+\log\left(\frac{K}{\varepsilon}\right)+C_{\mathcal{O}}-\mathbb{Q}\left[\langle\eta^{(r)}(\lambda),T^{(r)}(\cdot,d)\rangle\right]\right)}{\eta^{(1)}(\lambda)}\,.$

By factoring out $-\frac{1}{\eta^{(1)}(\lambda)}$ again, this is actually the
same as:

$\displaystyle\inf_{\lambda\in\Lambda}\
-\frac{1}{\eta^{(1)}(\lambda)}\left(\sup_{\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})}\mathbb{Q}[\langle\eta(\lambda),\tilde{T}(\cdot,d)\rangle]-D_{\rm{KL}}(\mathbb{Q}\
\|\
\mathbb{P}_{\mathbf{\mathfrak{H}}})-\log\left(\frac{K}{\varepsilon}\right)-C_{\mathcal{O}}\right),$

where $\tilde{T}(\alpha,d):=\bigl{(}-\hat{\mathcal{R}}(\alpha,d),\
T^{(r)}(\alpha,d)\bigr{)}$. Since $\log(K/\varepsilon)+C_{\mathcal{O}}$ is a
constant, Lemma 3.1 shows that the term inside the brackets is given by
$\tilde{\kappa}(\lambda,d)-\log(K/\varepsilon)-C_{\mathcal{O}}$, where
$\tilde{\kappa}$ corresponds to the exponential family
$(\tilde{\mathbb{Q}}_{\lambda})_{\lambda\in\Lambda}$ built upon $\tilde{T}$
and $\eta$ (with $h\equiv 1$). Furthermore, it shows that the (in this sense)
optimal posterior distribution
$\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})$ is given by the
corresponding member of the data-dependent exponential family
$\tilde{\mathbb{Q}}_{\lambda}(d)\propto\exp(\langle\eta(\lambda),\tilde{T}(\cdot,d)\rangle)\cdot\mathbb{P}_{\mathbf{\mathfrak{H}}}$,
usually called the Gibbs posterior (Alquier, 2021). By denoting
$F(\lambda,d):=-\frac{1}{\eta^{(1)}(\lambda)}(\tilde{\kappa}(\lambda,d)-\log(K/\varepsilon)-C_{\mathcal{O}})$,
one is left with solving the following problem:

$\inf_{\lambda\in\Lambda}\ F(\lambda,d),$ (4)

which for $\Lambda\subset\mathbb{R}$ is one-dimensional. Based on Theorem 4.1
and Theorem 4.9, we have to restrict to $\Lambda\subset(0,+\infty)$, such that
the solution to (4) can be seen as an approximation to the global minimum
$\inf_{\lambda>0}F(\lambda,d)$. For the latter one, one can show that the
solution set lies in a compact interval
$[\Lambda_{\rm{min}},\Lambda_{\rm{max}}]$, with $\Lambda_{\rm{min}}>0$, since
$F(\lambda,d)\to\infty$ as $\lambda\to 0$ or $\lambda\to\infty$. Under our
assumptions, $F(\cdot,d)$ is continuously differentiable. Hence, since
$\Lambda$ is compact, $F(\cdot,d)$ is Lipschitz-continuous on $\Lambda$ and
the minimum in (4) is attained. For a finite
$\Lambda\subset[\Lambda_{\rm{min}},\Lambda_{\rm{max}}]$, the optimization
reduces to grid search. For $\Lambda=[\Lambda_{\rm{min}},\Lambda_{\rm{max}}]$,
we employ grid search as initialization for gradient-based optimization.
However, we still might get stuck in a (close to optimal) local minimum. The
computational bottleneck is given by evaluating
$\lambda\mapsto\tilde{\kappa}(\lambda,d)$. In Sections 5.3.3 and 5.3.4 we will
ensure that this is cheap.

###  5.2 Sampling under Probabilistic Constraints

$\mathsf{C}$${Y}$${X}$$x$$\mathsf{C}_{x}$$x^{\prime}$$\mathsf{C}_{x^{\prime}}$$\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]$${X}$$p_{l}$$p_{u}$$1$$0$$\mathrm{supp}(\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}})$
Figure 3: Construction of $\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$: On the
left, the set $\mathsf{C}\subset{X}\times{Y}$ and two of its sections
$\mathsf{C}_{x},\mathsf{C}_{x^{\prime}}\subset{Y}$ are visualized. On the
right, the function $p(x)=\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]$,
the interval $[p_{l},p_{u}]$, and the resulting support
$\mathrm{supp}(\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}})$ of
$\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$ are visualized. Note that,
contrary to the visualization here, $p$ can actually be _highly
discontinuous_.

In this section, we describe a methodology that allows for sampling from a
distribution that is _probabilistically constrained_ in the following sense:
We are given two independent random variables
$\mathbf{\mathfrak{X}}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\to{X}$,
$\mathbf{\mathfrak{Y}}:\bigl{(}\Omega,\mathcal{F},\mathbb{P}\bigr{)}\to{Y}$
taking values in the Polish spaces ${X}$ and ${Y}$, with joint and marginal
distributions $\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}$,
$\mathbb{P}_{\mathbf{\mathfrak{X}}}$ and $\mathbb{P}_{\mathbf{\mathfrak{Y}}}$,
respectively. Further, we consider a measurable set
$\mathsf{C}\subset{X}\times{Y}$, and we want to generate samples
$\mathbf{\mathfrak{X}}=x\in{X}$, such that the probability of
$(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})$ lying in $\mathsf{C}$, given
$\mathbf{\mathfrak{X}}=x$, takes values in a certain interval:

$\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})|\mathbf{\mathfrak{X}}=x}[\mathsf{C}]=\mathbb{P}_{\mathbf{\mathfrak{Y}}|\mathbf{\mathfrak{X}}=x}[\mathsf{C}_{x}]\in[p_{l},p_{u}]\subset[0,1]\,.$

By independence of $\mathbf{\mathfrak{X}}$ and $\mathbf{\mathfrak{Y}}$, this
is $\mathbb{P}_{\mathbf{\mathfrak{X}}}$-almost surely the same as
$\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]\in[p_{l},p_{u}]$, and we
will use the later formulation from now on. While having an abstract
description of $\mathsf{C}$, we do not have any _explicit_ information about
its geometrical or topological properties, let alone access to a distance or a
projection onto it, and the marginal $\mathbb{P}_{\mathbf{\mathfrak{Y}}}$ can
only be accessed via samples
$\mathbf{\mathfrak{Y}}_{1},...,\mathbf{\mathfrak{Y}}_{n}$, $n\in\mathbb{N}$.
However, since we are able to evaluate, for every $x\in{X}$ and
$n\in\mathbb{N}$, whether $\mathbf{\mathfrak{Y}}_{n}\in\mathsf{C}_{x}$, this
allows to define a function $p:{X}\to[0,1]$, given by:

$\displaystyle\mathbb{P}_{(\mathbf{\mathfrak{X}},\mathbf{\mathfrak{Y}})}[\mathsf{C}]$
$\displaystyle=(\mathbb{P}_{\mathbf{\mathfrak{X}}}\otimes\mathbb{P}_{\mathbf{\mathfrak{Y}}})[\mathsf{C}]=\int_{X}\int_{Y}\mathds{1}_{\mathsf{C}}(x,y)\
\mathbb{P}_{\mathbf{\mathfrak{Y}}}(dy)\
\mathbb{P}_{\mathbf{\mathfrak{X}}}(dx)=\int_{X}\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]\
\mathbb{P}_{\mathbf{\mathfrak{X}}}(dx)$ $\displaystyle=:\int_{X}p(x)\
\mathbb{P}_{\mathbf{\mathfrak{X}}}(dx)\,.$

Then $p$ is indeed a measurable function, taking values in $[0,1]$. Thus, for
$p_{l},p_{u}\in[0,1]$ with $p_{l}<p_{u}$, we can define a measurable set
$\mathsf{A}:=\bigl{\\{}x\in{X}\ :\
\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]\in[p_{l},p_{u}]\bigr{\\}}$,
which yields a new measure $\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$ on
${X}$ by restricting to $\mathsf{A}$, that is, for a measurable set
$\mathsf{B}\subset{X}$ it holds:

$\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}[\mathsf{B}]:=\left((\mathds{1}_{[p_{l},p_{u}]}\circ
p)\cdot\mathbb{P}_{\mathbf{\mathfrak{X}}}\right)[\mathsf{B}]=\left(\mathds{1}_{\mathsf{A}}\cdot\mathbb{P}_{\mathbf{\mathfrak{X}}}\right)[\mathsf{B}]=\mathbb{P}_{\mathbf{\mathfrak{X}}}[\mathsf{A}\cap\mathsf{B}]\,.$

Therefore, as stated before, we have the following goal:

Goal: Sample from $\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$, that is, get
$\mathbf{\mathfrak{X}}_{1},...,\mathbf{\mathfrak{X}}_{K}\sim\mathbb{P}_{\mathbf{\mathfrak{X}}}$,
such that
$\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]|_{x=\mathbf{\mathfrak{X}}_{i}}\in[p_{l},p_{u}]$.

This construction is depicted in Figure 3: The left figure visualizes the
sections $\\{\mathsf{C}_{x}\\}_{x\in{X}}$ of the set $\mathsf{C}$, while the
right figure shows the corresponding construction of the support of
$\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$. Here, we assume that the imposed
constraint is realizable.

###### Assumption 5.1

The measure $\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$ has a non-empty
support.

###### Remark 5.2

Assumption 5.1 includes two things: (i) There have to be points $x\in{X}$
(hyperparameters) for which the property in $\mathsf{C}$ does hold, that is,
it has to be _realizable_. (ii) The support of
$\mathbb{P}_{\mathbf{\mathfrak{X}}}$ (prior) has to have a non-empty
intersection with this region.

###### Example 5.3

Consider the random variables $\mathbf{\mathfrak{S}}$ and
$\mathbf{\mathfrak{H}}$ from Section 4. By Lemma 4.11 we want to have
$p(\mathbf{\mathfrak{H}})\in[\varepsilon_{\rm{conv}},1]$, where the sublevel
probability is given as
$p(\alpha)=\mathbb{P}_{\mathbf{\mathfrak{S}}}[\mathsf{L}_{\alpha}]$ (Lemma
4.7), and the sublevel set $\mathsf{L}\subset{H}\times\Theta$ is measurable by
Lemma B.3. Thus, this corresponds to the identification ${X}={H}$,
${Y}=\Theta$, and $p_{l}=\varepsilon_{\rm{conv}}$, $p_{u}=1$.

####  5.2.1 Incorporation into a Sampling Procedure

By definition it holds that
$\mathbb{P}_{\mathbf{\mathfrak{X}}}=\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$
on $\mathsf{A}$. Hence, the only distinction between samples from
$\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$ and samples from
$\mathbb{P}_{\mathbf{\mathfrak{X}}}$ is the restriction to the set
$\mathsf{A}$. Therefore, since many sampling algorithms sample from
unnormalized measures anyways, it suffices to be able to sample from
$\mathbb{P}_{\mathbf{\mathfrak{X}}}$, if the restriction to $\mathsf{A}$ can
be satisfied in other ways. Thus, we have to integrate this constraint into a
sampling procedure for $\mathbb{P}_{\mathbf{\mathfrak{X}}}$. Because we do not
have any geometrical or topological information about the set $\mathsf{C}$, we
have to resort to statistical information: We have access to i.i.d. samples
$\mathbf{\mathfrak{Y}}_{1},...,\mathbf{\mathfrak{Y}}_{n}\sim\mathbb{P}_{\mathbf{\mathfrak{Y}}}$,
and, for a given $x\in{X}$, we are able to evaluate the i.i.d. Bernoulli-
distributed random variables
$\mathbf{\mathfrak{I}}_{n}:=\mathds{1}\\{\mathbf{\mathfrak{Y}}_{n}\in\mathsf{C}_{x}\\}$,
$n\in\mathbb{N}$. These have the parameter
$\mathbb{P}\\{\mathbf{\mathfrak{I}}_{n}=1\\}=\mathbb{P}\\{\mathbf{\mathfrak{Y}}_{n}\in\mathsf{C}_{x}\\}=\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x}]=p(x)$.
Thus, by estimating $p(x)$ with an estimator $\hat{p}(x)$, we can approximate
the constraint $\mathsf{A}$ with $\hat{\mathsf{A}}$:

$\mathsf{A}=\\{x\in{X}\ :\ p(x)\in[p_{l},p_{u}]\\}\approx\\{x\in{X}\ :\
\hat{p}(x)\in[p_{l},p_{u}]\\}=:\hat{\mathsf{A}}\,.$

Given a sample
$\mathbf{\mathfrak{X}}_{i}\sim\mathbb{P}_{\mathbf{\mathfrak{X}}}$, we resort
to a simple accept-reject mechanism as in _Metropolis-Hastings_ -type
algorithms (Robert and Casella, 2004), to decide whether
$\mathbf{\mathfrak{X}}_{i}\in\mathsf{A}$, that is, whether
$\mathbf{\mathfrak{X}}_{i}$ can actually be regarded as a sample from
$\tilde{\mathbb{P}}_{\mathbf{\mathfrak{X}}}$.

###### Remark 5.4

We believe that this is the only reasonable choice here, which keeps an
iterative algorithm _inside_ $\hat{\mathsf{A}}$. However, it does not provide
a way _into_ $\hat{\mathsf{A}}$, let alone $\mathsf{A}$.

Algorithm 1 Iterative estimation of the probability $p$

$p_{l},p_{u},\varepsilon\in[0,1]$.

$a,b\leftarrow 1,1$ $\triangleright$ Initialize with uninformative prior.

while $Q_{a,b}(p_{u})-Q_{a,b}(p_{l})\geq\varepsilon$ do $\triangleright$
$Q_{a,b}$ is the quantile function for $\mathrm{Beta}(a,b)$.

Draw $\mathbf{\mathfrak{I}}\sim\mathrm{Ber}(p)$

$a\leftarrow a+\mathbf{\mathfrak{I}}$ and $b\leftarrow
b+(1-\mathbf{\mathfrak{I}})$

end while

We estimate $p(x)$ in a Bayesian way, as it allows us to balance accuracy
against computational complexity through uncertainty-quantification, which we
use as a stopping criterion: We place a Beta-prior
$\mathrm{Beta}(a^{(0)},b^{(0)})$ over the interval $[0,1]$. As we do not have
prior knowledge, and the map $x\mapsto p(x)$ can be discontinuous222Consider
learning the step-size parameter $\alpha>0$ for gradient descent on quadratic
functions with largest eigenvalue $L$: The algorithm converges for
$\alpha<\frac{2}{L}$ ($p(\alpha)=1$) and diverges for $\alpha>\frac{2}{L}$
($p(\alpha)=0$) ., we use a noninformative prior (Berger, 1985, Ch. 3.3), that
is, $a_{0}=b_{0}=1$. Since the Beta distribution is the conjugate prior for
the Bernoulli distribution (Berger, 1985, p.130), that is, the posterior is
again a Beta-distribution, after observing a sample
$\mathbf{\mathfrak{I}}_{k+1}$, the parameters $a^{(k)},b^{(k)}$ get updated
as:

$a^{(k+1)}=a^{(k)}+\mathbf{\mathfrak{I}}_{k+1},\qquad
b^{(k+1)}=b^{(k)}+(1-\mathbf{\mathfrak{I}}_{k+1})\,.$

Hence, we can perform the estimation iteratively: We only draw a new sample
$\mathbf{\mathfrak{I}}_{n+1}$ as long as
$Q^{(n)}(p_{u})-Q^{(n)}(p_{l})\geq\varepsilon$, where $Q^{(n)}$ denotes the
quantile-function of the distribution $\mathrm{Beta}(a^{(n)},b^{(n)})$, and
$p_{u},p_{l},\varepsilon\in[0,1]$ are parameters that specify the accuracy of
the estimation. Finally, one can use the posterior mean
$\frac{a^{(n)}}{a^{(n)}+b^{(n)}}$ or posterior mode
$\frac{a^{(n)}-1}{a^{(n)}+b^{(n)}-2}$ (provided $a^{(n)},b^{(n)}>1$) as point
estimate $\hat{p}_{x}$. This procedure is summarized in Algorithm 1 and
depicted in Figure 4.

###### Remark 5.5

By adjusting $p_{l},p_{u}$ or $\varepsilon$, one can balance between accuracy
and computational complexity. Yet, the number of iterations needed also
depends crucially on the true probability: For $p\approx 0$ or $p\approx 1$,
the uncertainty decreases significantly faster than for $p\approx 0.5$.

Figure 4: Iterative estimation of $p(x)$: The black line shows the density
$f_{a^{(n)},b^{(n)}}$ of $\mathrm{Beta}(a^{(n)},b^{(n)})$ after having
observed $\mathbf{\mathfrak{I}}_{1},...,\mathbf{\mathfrak{I}}_{n}$. The red
dotted line indicates the true probability $p$, which we are trying to
estimate, and the blue dashed lines indicate the lower and upper quantiles
$q_{l},q_{u}$ corresponding to $p_{l}$ and $p_{u}$. The procedure stops as
soon as $q_{u}-q_{l}<\varepsilon$, which is indicated by the double-headed
arrow. Here, we use $q_{l}=0.05,q_{u}=0.95,\varepsilon=0.15$ and $p=0.8$.

####  5.2.2 Broader Context

Different, yet conceptually similar ideas have appeared in the context of
cutting the computational cost of Bayesian Markov Chain Monte Carlo algorithms
through subsampling: Korattikara et al. (2014) introduce sequential hypothesis
tests to reach the binary accept-reject decision in the Metropolis-Hastings
algorithm on subsamples of the data set. Similarly, Bardenet et al. (2014)
estimate the accept-reject step on a random subset of the data, in a way that
guarantees to coincide with the true accept-reject step with a user-specified
probability. Maclaurin and Adams (2014) introduce an auxiliary binary variable
$z_{n}\in\\{0,1\\}$, which allows to query only a subset of the data for the
computation of the exact likelihood. And Quiroz et al. (2018) use a
subsampling and bias-correction strategy to speed-up the sampling procedure.
Here, Bardenet et al. (2017) provides a summary and discussion of different
approximations and their biases (and errors). We leave the corresponding
analysis for our proposed approximation to future work.

####  5.2.3 Choice of the Sampling Procedure

Often, the hyperparameters $\alpha\in{H}$ of the learned algorithm
$\mathcal{A}$ are high-dimensional. Thus, we resort to using _stochastic
gradient Langevin dynamics_ (Welling and Teh, 2011) as the underlying sampling
algorithm, which gets constrained to the set $\hat{\mathsf{A}}$ by use of the
previously described procedure. This is summarized in Algorithm 2. However, if
it fits the application, other sampling algorithms can be used, too. The
computational overhead of the additional estimation depends on the cost of
evaluating $\mathds{1}\\{\mathbf{\mathfrak{Y}}_{n}\in\mathsf{C}_{x}\\}$. In
our case it is rather expensive, as every sample $\mathbf{\mathfrak{I}}_{n}$
requires to run the algorithm $\mathcal{A}$. However, this is to be expected,
as the “prediction” of an optimization algorithm corresponds to approximating
the solution of a minimization problem.

###### Remark 5.6

Algorithm 2 requires to start in the set $\hat{\mathsf{A}}$. If such a point
is not known, one can still run the algorithm and just “start” the accept-
reject mechanism as soon as one has found a point $x\in\hat{\mathsf{A}}$.
However, it is not guaranteed that such a point will actually be found.

Algorithm 2 Probabilistically constrained sampling

$p_{l},p_{u}\in[0,1]$, $n_{\rm{max}}\in\mathbb{N}$,
$x_{0}\in\hat{\mathsf{A}}$.

$n\leftarrow 0$ and $x\leftarrow x_{0}$

while $n\leq n_{\rm{max}}$ do

1) Draw a proposal $x^{\prime}$ with SGLD starting from $x$.

2) Estimate
$p(x^{\prime})=\mathbb{P}_{\mathbf{\mathfrak{Y}}}[\mathsf{C}_{x^{\prime}}]$ by
$\hat{p}(x^{\prime})$ with Algorithm 1.

if $\hat{p}(x^{\prime})\in[p_{l},p_{u}]$ then

$x\leftarrow x^{\prime}$

else

Reject $x^{\prime}$.

end if

end while

###  5.3 Putting Things Together

In the abstract learning procedure described in Sections 3 and 4 everything
relies on the existence of a prior distribution
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$ that a) gives a good performance of the
algorithm in terms of a small risk, and b) ensures the user-specified sublevel
probability. Both a) and b) have to be enforced through the _construction_ of
the prior $\mathbb{P}_{\mathbf{\mathfrak{H}}}$, which is why often a _data-
dependent_ prior is used (Parrado-Hernández et al., 2012; Lever et al., 2013;
Dziugaite and Roy, 2017, 2018; Dziugaite et al., 2021). Note that a) and b)
are qualitatively different: a) asks for a certain location and concentration
of $\mathbb{P}_{\mathbf{\mathfrak{H}}}$, while b) puts hard constraints on its
support, that is, it _excludes_ many frequently used distributions with
unconstrained support. Therefore, a significant part of the learning procedure
deals with the _construction_ of this prior only. Since the prior will be
data-dependent, yet has to be independent of the data set that is used in the
PAC-Bayesian learning procedure in Section 5.3.4, we have to split the data
set $\mathbf{\mathfrak{D}}$ into independent parts
$\mathbf{\mathfrak{D}}_{\rm{prior}}$, $\mathbf{\mathfrak{D}}_{\rm{val}}$,
$\mathbf{\mathfrak{D}}_{\rm{train}}$ and $\mathbf{\mathfrak{D}}_{\rm{test}}$,
where $\mathbf{\mathfrak{D}}_{\rm{prior}}$ and
$\mathbf{\mathfrak{D}}_{\rm{val}}$ are used for the construction of the prior
distribution, $\mathbf{\mathfrak{D}}_{\rm{train}}$ is used for the PAC-
Bayesian learning step, and $\mathbf{\mathfrak{D}}_{\rm{test}}$ is the test
set to evaluate the performance, that is, it is only needed for the
experiments. Nevertheless, for notational simplicity, we will use the generic
$\mathbf{\mathfrak{D}}$, implicitly assuming the above partitioning. As the
data set is _fixed_ now, it will be described by the corresponding
_realization_ $D=\\{\theta_{1},...,\theta_{N}\\}$ of the random variable
$\mathbf{\mathfrak{D}}=\\{\mathbf{\mathfrak{S}}_{1},...,\mathbf{\mathfrak{S}}_{N}\\}$.

###### Remark 5.7

Through the choice of the sampling algorithm, the concrete learning procedure
described here mainly applies to the case ${H}=\mathbb{R}^{d}$,
$d\in\mathbb{N}$. Nevertheless, the general methodology is still applicable to
other Polish spaces, if this choice can be adjusted accordingly.

####  5.3.1 Finding a Trainable Initialization

To increase numerical stability, we start with “imitation learning” (Chen et
al., 2020a), that is, the algorithm $\mathcal{A}$ should “follow” another
algorithm $\mathcal{A}^{\prime}$, for example, gradient descent. For this, we
minimize the mean squared error between the iterates of the two algorithms:
Given a starting point $x^{(0)}\in\mathbb{R}^{n}$, an iteration number
$s\in\mathbb{N}$, and a parameter $\theta\in\Theta$, denote the first $s$
iterates of $\mathcal{A}(\alpha,\theta,x^{(0)})$ by
$x^{(1)},...,x^{(s)}\in\mathbb{R}^{n}$ and the ones of
$\mathcal{A}^{\prime}(x^{(0)},\theta)$ by
$y^{(1)},...,y^{(s)}\in\mathbb{R}^{n}$. Then, define the loss as the mean
squared error over these iterations:

$\ell_{\rm{init}}(\alpha,\theta,x^{(0)},s):=\frac{1}{s}\sum_{k=1}^{s}\|x^{(k)}-y^{(k)}\|_{2}^{2}\,.$

In each iteration, that is, each prediction of tuples
$(x^{(1)},y^{(1)}),...,(x^{(s)},y^{(s)})$, the parameters $\theta$, $x^{(0)}$
and $s$ are randomized as described in Section 5.3.2. We postpone the
discussion of this step, as it is an integral part of the learning step
described there, while being less important here. We terminate the process as
soon as the average loss over $n_{\rm{init}}\in\mathbb{N}$ iterations reaches
a certain level of accuracy $\varepsilon>0$, that is, as soon as
$\frac{1}{n_{\rm{init}}}\sum_{i=1}^{n_{\rm{init}}}\ell_{\rm{init}}(\alpha_{i},\theta_{i},x^{(0)}_{i},s)<\varepsilon$.
The choice of $n_{\rm{init}}$ and $\varepsilon$ depends on the problem at
hand. Usually, rough estimates suffice ($\varepsilon\leq 10$), as the purpose
is to prevent divergence, and not actual imitation of $\mathcal{A}^{\prime}$.
The procedure is summarized in Algorithm 3.

Algorithm 3 Procedure to find an initialization

Data set $D_{\rm{prior}}$, $x^{(0)}\in\mathbb{R}^{n}$,
$s,n_{\rm{init}}\in\mathbb{N}$ and $\varepsilon>0$.

$m\leftarrow+\infty$ and sample $\theta\sim\mathit{U}_{D_{\rm{prior}}}$.

while $\frac{1}{n_{\rm{init}}}m>\varepsilon$ do

$m\leftarrow 0$

for $i=1,...,n_{\rm{init}}$ do

1) Compute $(x^{(1)},y^{(1)}),...,(x^{(s)},y^{(s)})$ with
$\mathcal{A}(\alpha,\theta,x^{(0)})$ and
$\mathcal{A}^{\prime}(\theta,x^{(0)})$, resp.

2) Compute
$\ell_{\rm{init}}(\alpha,\theta,x^{(0)},s)=\frac{1}{s}\sum_{k=1}^{s}\|x^{(k)}-y^{(k)}\|_{2}^{2}$.

3) Update $m\leftarrow m+\ell_{\rm{init}}(\alpha,\theta,x^{(0)},s)$

4) Update $\alpha$ by backpropagation and Adam. $\triangleright$ Other
algorithms possible.

5) Update $\theta$, $x^{(0)}$ and $s$ based on Section 5.3.2.

end for

end while

####  5.3.2 Locating the Prior

Algorithm 4 Procedure to locate the prior

Data sets $D_{\rm{prior}},D_{\rm{val}}$, numbers
$n_{\rm{max}},n,s\in\mathbb{N}$ with $s\leq n$, initialization $x^{(0)}$ and
thresholds $p_{l},p_{u}\in[0,1]$ with $p_{l}<p_{u}$.

………………………………………………………………………………………………………………………….

Set $x\leftarrow x^{(0)}$ and $b\leftarrow\mathbf{false}$, and sample
$\theta\sim\mathrm{U}_{D_{\rm{prior}}}$ $\triangleright$ $b$ = Point inside
constraint?

………………………………………………………………………………………………………………………….

for $i=1,...,n_{\rm{max}}$ do $\triangleright$ Other stopping criteria
possible.

……………………………………………………………………………………………………………………..

1.a) Compute $x^{(1)},...,x^{(s)}$ with $\mathcal{A}(\alpha,\theta,x^{(0)})$.

1.b) Compute
$\ell_{\rm{train}}(\alpha,\theta,x^{(0)},s)=\sum_{i=1}^{s}\mathds{1}_{\\{\ell(x^{(i-1)},\theta)>0\\}}\frac{\ell(x^{(i)},\theta)}{\ell(x^{(i-1)},\theta)}$.

1.c) Construct a proposal $\tilde{\alpha}$ by using backpropagation and Adam.

……………………………………………………………………………………………………………………..

2) Estimate $p(\tilde{\alpha})$ by $\hat{p}(\tilde{\alpha})$ with Algorithm 1
on $D_{\rm{val}}$.

if $\hat{p}(\tilde{\alpha})\in[p_{l},p_{u}]$ then $\triangleright$ If point
inside constraint, just update.

$\alpha\leftarrow\tilde{\alpha}$ and $b\leftarrow\mathbf{true}$

else$\triangleright$ If not…

if $b=\mathbf{true}$ then $\triangleright$ …reject moving outside constraint.

Reject $\tilde{\alpha}$, set $x^{(0)}\leftarrow x$, sample
$\theta\sim\mathrm{U}_{D_{\rm{prior}}}$, and continue with 1).

else$\triangleright$ …accept, if constraint has not been found yet.

$\alpha\leftarrow\tilde{\alpha}$

end if

end if

……………………………………………………………………………………………………………………..

3) Draw $\mathbf{\mathfrak{r}}\sim\mathrm{Ber}(\frac{s}{n})$.

if $\mathbf{\mathfrak{r}}=0$ then

$x^{(0)}\leftarrow x^{(s)}$

else $\mathbf{\mathfrak{r}}=1$

$x^{(0)}\leftarrow x$ and sample $\theta\sim\mathrm{U}_{D_{\rm{prior}}}$

end if

……………………………………………………………………………………………………………………..

end for

Empirically, the performance of the learned algorithm is significantly
improved by the following two design choices, namely randomizing the
trajectory length and using the ratio of successive losses. The motivation is
to prevent overfitting and to learn a scale-independent contraction of the
loss:
1) Ratio of Losses: Since the algorithm $\mathcal{A}(\alpha,\cdot)$ should
achieve a small risk
$\mathcal{R}(\alpha)=\mathbb{E}_{\mathbf{\mathfrak{S}}}[\ell(\alpha,\cdot)]$,
yet minimizing $\mathcal{R}$ is intractable, the canonical loss function to be
minimized is
$\hat{\mathcal{R}}(\alpha,D_{N})=\frac{1}{N}\sum_{i=1}^{N}\ell(\alpha,\theta_{i})$.
As ${H}$ might be high-dimensional, we resort to stochastic empirical risk
minimization, that is, in each iteration the observed loss would be of the
form $\ell(\alpha,\theta_{i})$. While this kind of loss was used extensively
before, for learning-to-optimize it has a strong disadvantage: The overall
outcome gets penalized only after the application of the _whole_ algorithm,
that is, after $n_{\rm{max}}$ iterations of the update-step. Thus, it does not
take the trajectory into account. Further, often it is hard to minimize and
does not lead to the desired performance. To circumvent this, Andrychowicz et
al. (2016) proposed to use the sum over the incurred losses, that is, to
consider
$\tilde{\ell}_{\rm{train}}(\alpha,\theta,x^{(0)}):=\sum_{i=1}^{n}\ell(x^{(i)},\theta)$.
Again, this formulation has a decisive flaw: Under most objectives, if the
algorithm performs reasonably well, the loss at the beginning is _several
orders_ of magnitude larger than the loss at the end, that is,
$\ell(x^{(1)},\theta)\gg\ell(x^{(n)},\theta)$. Hence,
$\tilde{\ell}_{\rm{train}}$ mainly penalizes the loss at the beginning,
leading to an algorithm that minimizes the loss very fast in early iterations,
but slows down _a lot_ in later iterations. This is due to
$\tilde{\ell}_{\rm{train}}$ being _scale-sensitive_. Additionally, the
incurred loss might vary strongly with the initialization $x^{(0)}$ alone,
thereby introducing ambiguity into the incurred losses. We propose to use the
_ratio of losses_ of successive iterates:

$\ell_{\rm{train}}(\alpha,\theta,x^{(0)},s):=\sum_{i=1}^{s}\mathds{1}_{\\{\ell(x^{(i-1)},\theta)>0\\}}\frac{\ell(x^{(i)},\theta)}{\ell(x^{(i-1)},\theta)}\,,\qquad
s\in\mathbb{N},\ s\leq n\,.$

This has several advantages: First, the loss is not scale-sensitive anymore,
such that it favours hyperparameters that yield a good performance in each
iteration. Second, there is no ambiguity in the observed loss through the
initialization, as the only criterion is a strong contraction of the loss
(instead of a small loss). Intuitively, the algorithm $\mathcal{A}$ gets
trained to have a fast _convergence rate_. Third, the incurred losses do not
vary too much, which empirically makes it easier to choose hyperparameters of
the learning procedure. However, it also has a disadvantage: If the function
values do indeed _converge_ in a setting where the optimal loss is strictly
greater than zero, this gets fully penalized, as then
$\frac{\ell(x^{(i)},\theta)}{\ell(x^{(i-1)},\theta)}\equiv 1$. For now, we do
not know how to get rid off this problem (apart from just stopping the
iterations in case of convergence) while keeping the advantages, yet we
(successfully) combat it with a decreasing step-size during training.
2) Randomized Trajectory Length: Training $\mathcal{A}$ with fixed
initialization $x^{(0)}$ and fixed trajectory length $n\in\mathbb{N}$, leads
to overfitting, that is, applying $\mathcal{A}$ at another starting point
$\tilde{x}^{(0)}$, or applying it for more than $n$ iterations typically leads
to divergence. To get rid off this, we propose the following randomization:
Given $\mathbb{N}\ni s\leq n$, set $p:=\frac{s}{n}$ and $y:=x^{(0)}$.

  * 0)

Sample a parameter $\theta$ uniformly at random from $D_{N}$.

  * 1)

Compute $x^{(1)},...,x^{(s)}$ with $\mathcal{A}(\alpha,\theta,y)$ and the loss
$\ell_{\rm{train}}(\alpha,\theta,y,s)$, and update $\alpha$.

  * 2)

Sample $\mathbf{\mathfrak{r}}^{(k)}\sim\mathrm{Ber}(p)$. If
$\mathbf{\mathfrak{r}}^{(k)}=0$, set $y:=x^{(s)}$ and go to step 1). If
$\mathbf{\mathfrak{r}}^{(k)}=1$, set $y:=x^{(0)}$ and go to step 0).

The random variable $\mathbf{\mathfrak{r}}^{(k)}$ decides whether the
algorithm gets restarted from $x^{(0)}$ with a new parameter $\tilde{\theta}$,
or if one continuous the current trajectory. To update $\alpha$, we use
backpropagation and Adam. The choice $p=\frac{s}{n}$ ensures that the expected
trajectory-length equals $n$: Define $\tau:=\inf\\{k\in\mathbb{N}\ :\
\mathbf{\mathfrak{r}}^{(k)}=1\\}$. Then, $\tau\sim\mathrm{Geo}(p)$ is a
geometrically distributed with expectation $\mathbb{E}[\tau]=\frac{1}{p}$.
Therefore, the actual length $\mathbf{\mathfrak{L}}=s\cdot\tau$ of the
trajectory is (basically) a geometrically distributed random variable with
$\mathbb{E}[\mathbf{\mathfrak{L}}]=s\mathbb{E}[\tau]=s\frac{1}{p}=s\frac{n}{s}=n$.

###### Remark 5.8

  * (i)

Similarly to Andrychowicz et al. (2016), during training we omit the
computation of second order derivatives. Additionally, and surprisingly, it
usually suffices to consider single iterates, that is $s=1$, where the
iterates are treated _independently_ in the computational graph. That amounts
to learning an update step that is agnostic to the recurrent nature of the
optimization algorithm and just learns to adapt to the local geometry of the
loss function _along the iterations_.

  * (ii)

Splitting the trajectory into several parts has been used before, especially
in model-free approaches based on recurrent neural networks (Andrychowicz et
al., 2016; Chen et al., 2017; Metz et al., 2019), where the main motivation
are training instabilities and memory bottlenecks. Similarly, randomizing the
trajectory length did appear in Wichrowska et al. (2017). Yet, to our
knowledge, our specific randomization is new, theoretically motivated, and
solves the problem of generalizing to more iterations quite effectively.

Empirically, by using these design choices we are able to train $\mathcal{A}$
in a way such that it performs equally well on several orders of magnitude of
the loss function, and that it can be used for more iterations than it was
actually trained for. The procedure is summarized in Algorithm 4 and consists
of three steps:

  * 1)

Given $\theta,x^{(0)}$ and $s\in\mathbb{N}$, compute a new proposal
$\tilde{\alpha}$ by performing one step of Adam on
$\ell_{\rm{train}}(\alpha,\theta,x^{(0)},s)$, corresponding to $s$ iterations
of $\mathcal{A}$.

  * 2)

Check for the constraint $p(\tilde{\alpha})\in[p_{l},p_{u}]$ by estimating it
with Algorithm 1. If the constraint is satisfied (or one did not find any
point inside the constraint so far), accept $\tilde{\alpha}$. Otherwise,
reject it and sample a new problem instance
$\tilde{\theta}\sim\mathrm{U}_{D_{N}}$.

  * 3)

If $\tilde{\alpha}$ got accepted, update $\theta$ and $x^{(0)}$ by randomizing
the length of the trajectory.

####  5.3.3 Constructing the Prior

Algorithm 5 Procedure to construct the prior

Data sets $D_{\rm{prior}}$ (sampling) and $D_{\rm{val}}$ (constraint),
$n_{s}\in\mathbb{N}$ and
$\alpha\in\mathrm{supp}({\tilde{\mathbb{P}}_{\mathbf{\mathfrak{H}}}})$.

1) Starting from $\alpha$, run Algorithm 2 (with $\ell_{train}$, and
randomized trajectory length) to get the points
$\alpha_{1},...,\alpha_{n_{s}}\in{H}$.

2) Evaluate $\varphi_{\rm{prior}}$ on $\\{\alpha_{1},...,\alpha_{n_{s}}\\}$ by
evaluating $\hat{\mathcal{R}}_{g}$ corresponding to $D_{\rm{prior}}$.

3) Compute $\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{\alpha_{j}\\}$, that is,
$\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{\alpha_{j}\\}=\sigma\left(\varphi_{\rm{prior}}(\alpha_{1}),...,\varphi_{\rm{prior}}(\alpha_{n_{s}})\right)_{j}\,.$

Besides the performance and the sublevel guarantees, the only assumption on
the prior $\mathbb{P}_{\mathbf{\mathfrak{H}}}$ is its independence of
$\mathbf{\mathfrak{D}}_{\rm{train}}$. Further, by Lemma 3.1 the functional
form of the posterior is fully specified, namely it is of the form:

$\mathbb{Q}_{\lambda}(d)\propto\exp(\varphi_{\lambda}(\cdot,d))\cdot\mathbb{P}_{\mathbf{\mathfrak{H}}},\qquad\lambda\in\Lambda\,,$
(5)

where the potential is given by
$\varphi_{\lambda}(\alpha,d)=\langle\eta(\lambda),\tilde{T}(\alpha,d)\rangle$.
Hence, for mathematical convenience, we will construct
$\mathbb{P}_{\mathbf{\mathfrak{H}}}$ by approximating the distribution
$\mathbb{P}^{\prime}$ given by

$\mathbb{P}^{\prime}\propto\exp\left(-\hat{\mathcal{R}}_{g,\rm{prior}}-\iota_{[p_{l},p_{u}]}\circ
p\right)\cdot\mu\,,$

where $\mu$ denotes another measure on ${H}$, which allows to sample from
$\mathbb{P}^{\prime}$ (possibly unnormalized). In our experiments it holds
${H}=\mathbb{R}^{d}$ and we choose $\mu=\lambda^{d}$, where $\lambda^{d}$ is
the $d$-dimensional Lebesgue measure. For sampling, we use (in principal) the
_Stochastic Gradient Langevin Dynamics_ algorithm. However, for ease of
implementation and its practicability, we use the output of the
backpropagation algorithm as proxy for the (sub)gradient, which is an element
of a so-called _conservative set-valued field_ (Bolte and Pauwels, 2021),
which can be shown (under mild assumptions) to be $\lambda^{d}$-almost
everywhere the same as the true gradient (Bianchi et al., 2022). Since we have
to resort to a sampling algorithm to get points
$\alpha_{1},...,\alpha_{n_{s}}\in{H}$, $n_{s}\in\mathbb{N}$, we define the
prior distribution directly as a discrete distribution, that is
$\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{\alpha\\}:=\frac{1}{Z}\sum_{i=1}^{n_{s}}a_{i}\delta_{\alpha_{i}}\\{\alpha\\}$.
Thus, $\mathbb{P}_{\mathbf{\mathfrak{H}}}$ is the suitable normalized discrete
measure on ${H}$ corresponding to $\alpha_{1},...,\alpha_{n}$, where the
coefficients and the normalization constant are given by
$Z=\sum_{i=1}^{n_{s}}a_{i}$ with
$a_{i}=\exp\left(-\hat{\mathcal{R}}_{g,\rm{prior}}(\alpha_{i})-\iota_{[p_{l},p_{u}]}\left(\hat{p}(\alpha_{i})\right)\right)$.
When $\alpha_{1},...,\alpha_{n_{s}}\in{H}$ are given, the corresponding
probabilities can equivalently be expressed with the so-called _softmax_
function
$\sigma(x_{1},...,x_{n})_{j}=\frac{\exp(x_{j})}{\sum_{i=1}^{n_{s}}\exp(x_{i})}$:
Define the potential $\varphi_{\rm{prior}}:{H}\to\mathbb{R}$ as
$\varphi_{\rm{prior}}(\alpha)=-\hat{\mathcal{R}}_{g,\rm{prior}}(\alpha)-\iota_{[p_{-},p_{+}]}\left(\hat{p}(\alpha)\right)$
to get:

$\displaystyle\mathbb{P}_{\mathbf{\mathfrak{H}}}\\{\alpha_{j}\\}$
$\displaystyle=\frac{\exp\left(\varphi_{\rm{prior}}(\alpha_{j})\right)}{\sum_{i=1}^{n_{s}}\exp\left(\varphi_{\rm{prior}}(\alpha_{i})\right)}=\sigma\left(\varphi_{\rm{prior}}(\alpha_{1}),...,\varphi_{\rm{prior}}(\alpha_{n_{s}})\right)_{j}\,.$

Hence, the prior distribution is given by evaluating the softmax-function,
where the potentials $\varphi_{\rm{prior}}$ have to be computed _only once_
for every $\alpha_{i}$, $i=1,...,n_{s}$. This is summarized in Algorithm 5.
Note that for sampling $\alpha_{1},...,\alpha_{n_{s}}$, the loss function
$\ell_{train}$ is used again.

###### Remark 5.9

As one would approximate the intractable integrals with Monte-Carlo estimates
anyway, choosing a discrete measure is less restrictive than it seems, and it
has the additional advantage of allowing for _exact_ instead of approximate
inference.

####  5.3.4 Computing the Posterior

Algorithm 6 Procedure to construct the posterior

Points $\\{\alpha_{1},...,\alpha_{n_{s}}\\}$, values
$\\{\varphi_{\rm{prior}}(\alpha_{1}),...,\varphi_{\rm{prior}}(\alpha_{n_{s}})\\}$,
data set $D_{\rm{train}}$.

1) Evaluate $\tilde{T}(\alpha_{i},d)$, $i=1,...,n_{s}$, on $D_{\rm{train}}$.

2) Setup
$\\{\varphi_{\lambda}(\alpha_{1},d),...,\varphi_{\lambda}(\alpha_{n_{s}},d)\\}$
as functions in $\lambda$.

3) Solve $\lambda^{*}\in\operatorname*{arg\,min}_{\lambda\in\Lambda}\
F(\lambda,d)$. $\triangleright$ $F(\lambda^{*},d)$ is the predicted PAC-bound.

4) Compute
$\mathbb{Q}_{\lambda^{*}}(d,\\{\alpha_{j}\\})=\sigma\left(\varphi_{\lambda^{*}}(\alpha_{1},d),...,\varphi_{\lambda^{*}}(\alpha_{n_{s}},d)\right)_{j}$,
$j=1,...,n_{s}$.

5) Optional: Choose
$\alpha^{*}=\operatorname*{arg\,max}_{i=1,...,n_{s}}\mathbb{Q}_{\lambda^{*}}(d,\\{\alpha_{i}\\})$
as final point-estimate.

Given a discrete prior $\mathbb{P}_{\mathbf{\mathfrak{H}}}$, every posterior
$\mathbb{Q}\in\mathcal{P}(\mathbb{P}_{\mathbf{\mathfrak{H}}})$ is also
discrete with the same support $\\{\alpha_{1},...,\alpha_{n_{s}}\\}$. Then, by
the closed-form solution (5) for the supremum, for every $\lambda\in\Lambda$
the optimal posterior $\mathbb{Q}_{\lambda}(d)$ is given by:

$\displaystyle\mathbb{Q}_{\lambda}(d,\\{\alpha_{j}\\})$
$\displaystyle=\frac{\exp\left(\langle\eta(\lambda),\tilde{T}(\alpha_{j},d)\rangle+\varphi_{\rm{prior}}(\alpha_{j})\right)}{\sum_{i=1}^{n_{s}}\exp\left(\langle\eta(\lambda),\tilde{T}(\alpha_{i},d)\rangle+\varphi_{\rm{prior}}(\alpha_{j})\right)}=\sigma\left(\varphi_{\lambda}(\alpha_{1},d),...,\varphi_{\lambda}(\alpha_{n},d)\right)_{j}\,,$

with the potential $\varphi_{\lambda}:{H}\times\Theta^{N}\to\mathbb{R},\
(\alpha,d)\mapsto\langle\eta(\lambda),\tilde{T}(\alpha,d)\rangle+\varphi_{\rm{prior}}(\alpha)$.
Thus, to get the distribution $\mathbb{Q}_{\lambda}(d)$ as a function of
$\lambda$, one has to compute $\tilde{T}(\alpha_{i},d)$ _only once_ for every
$i=1,...,n_{s}$, such that it can be evaluated with the softmax function.
Hence, the only missing ingredient is the optimal $\lambda^{*}\in\Lambda$,
which is found as described in Section 5.1. After evaluating the potentials
$\varphi_{\lambda}(\cdot,d)$, which has to be done anyways, evaluating
$\tilde{\kappa}(\cdot,d)$ in $\lambda$ is cheap. The process is summarized in
Algorithm 6.

##  6 Experiments

We consider the _smooth and strongly convex_ problem of minimizing quadratic
functions with varying strong convexity and smoothness constants, a _high-
dimensional_ image processing problem, the _non-smooth_ Lasso problem, and the
_non-smooth and non-convex_ problem of training a neural network. We use the
following training procedure in all experiments:
$N=N_{\rm{prior}}+N_{\rm{train}}+N_{\rm{val}}+N_{\rm{test}}$ denotes the total
number of datapoints, and we use $N_{\rm{prior}}=...=N_{\rm{test}}=250$.
(Sub)Gradients are defined by the output of backpropagation as it is
implemented in PyTorch (Paszke et al., 2019), and we use
$g(\theta):=a\ell(x^{(0)},\theta)^{b}$, $a,b>0$, to define the sublevel set
$\mathsf{L_{g}}$. In Algorithm 1, we use $p_{l}=0.95$, $p_{u}=1.0$,
$q_{l}=0.01$, $q_{u}=0.99$, and $\varepsilon=0.075$. Thus, the algorithm
should reach $\mathsf{L_{g}}$ in at least 95% of the cases, and for the
estimation of the sublevel probability it should concentrate 98% of the mass
within a distance of 0.075. In Algorithm 2, we use _stochastic gradient
Langevin dynamics_ to draw $10^{2}$ samples, where we decay the step-size
starting from $\alpha=10^{-6}$. In Algorithm 3, we use Adam with an initial
step-size of $10^{-3}$, which gets reduced by a factor of 0.5 every 200
iterations, until an accuracy of $\varepsilon=10^{-2}$ is reached, or for at
most $n_{\rm{init}}=10^{3}$ iterations. In Algorithm 4, we use Adam with an
initial step-size of $10^{-4}$, which gets reduced by a factor of 0.5 every
$2\cdot 10^{4}$ iterations, for a total of $n_{\rm{max}}=2\cdot 10^{5}$
iterations. We use a trajectory length of $s=1$, that is, only single points,
and update the constraint only every $2\cdot 10^{4}$ iterations (with a reset
to previous iterates, if we have left the set $\hat{\mathsf{A}}$). In
Algorithm 6, we use a finite $\Lambda$ with $|\Lambda|=75\cdot 10^{3}$, and an
accuracy (of the PAC-bound) of $\varepsilon=0.05$. As we contrast the learned
algorithm to first-order methods, in each iteration $\mathcal{A}$ has access
to iterates, (sub)gradients, and function values, and the update is solely
based on these. Here, we perform preprocessing: The (sub)gradient is split
into its norm $\|\nabla\ell(x^{(k)})\|$ and the corresponding unit vector
$d^{(k)}_{1}$. Further, the norm is transformed to
$n^{(k)}_{1}:=\log(1+\|\nabla\ell(x^{(k)})\|)$ to be less scale-sensitive. The
iterates $x^{(k)},x^{(k-1)}$ are combined into the momentum term
$m^{(k)}:=x^{(k)}-x^{(k-1)}$, which also is split into the unit vector
$d^{(k)}_{2}$ and the transformed norm $n^{(k)}_{2}$. In the evaluation, we
will always show the _loss over the iterations_ in the upper left plot, the
_performance in terms of computation time_ in the upper right plot, the _loss
histogram with predicted PAC-bound_ in the lower left plot, and the final
_estimate for the sublevel probability_ , that is, whether the
probabilistically constrained optimization/sampling procedure did work
correctly, in the lower right plot. Finally, note that we always show the
performance of a single sample $\alpha^{*}$, and not the mean performance.

###### Remark 6.1

  * (i)

We always use the output of the backpropagation algorithm instead of exact
(sub-)gradients, that is, the learned algorithms _do not rely on smoothness_.

  * (ii)

We use 100 samples only, as they are very costly: To evaluate the potentials
$\varphi_{\rm{prior}}$ and $\varphi_{\lambda}(\cdot,d)$ on a single sample
$\alpha\in{H}$, one has to compute all losses $\ell(\alpha,\theta_{i})$,
$i=1,...,N_{\rm{prior}}+N_{\rm{train}}$, that is, “solving”
$N_{\rm{prior}}+N_{\rm{train}}$ optimization problems.

###  6.1 Quadratics

$d^{(k)}_{1}$$d^{(k)}_{2}$$d^{(k)}_{1}\odot
d^{(k)}_{2}$Conv2d(3,16,1,bias=F)ReLUConv2d(16,16,1,bias=F)Conv2d(16,16,1,bias=F)ReLUConv2d(16,16,1,bias=F)Conv2d(16,16,1,bias=F)ReLUConv2d(16,1,1,bias=F)$d^{(k)}$$n^{(k)}_{1}$$n^{(k)}_{2}$Linear(2,8,bias=F)ReLULinear(8,8,bias=F)Linear(8,8,bias=F)ReLULinear(8,8,bias=F)Linear(8,8,bias=F)ReLULinear(8,1,bias=F)$s^{(k)}$$x^{(k)}$$x^{(k+1)}:=x^{(k)}+s^{(k)}\cdot
d^{(k)}$ Figure 5: Update step of $\mathcal{A}$ for quadratic problems: The
directions $d^{(k)}_{1}$, $d^{(k)}_{2}$ and $d^{(k)}_{1}\odot d^{(k)}_{2}$ are
inserted as different channels into the Conv2d-block, which performs $1\times
1$ “convolutions”, that is, the algorithm acts coordinate-wise on the input.
The scales $n^{(1)}$, $n^{(2)}$ get transformed separately by the fully-
connected block.

As first problem we consider strongly convex quadratic functions with varying
strong convexity, varying smoothness and varying right-hand side, that is,
each optimization problem is of the form:

$\min_{x\in\mathbb{R}^{n}}\frac{1}{2}\|Ax-b\|^{2}\,,\qquad
A\in\mathbb{R}^{n\times n},\ b\in\mathbb{R}^{n}\,.$

Thus, the parameters are given by
$\theta=(A,b)\in\mathbb{R}^{n^{2}+n}=:\Theta$, while the optimization variable
is $x\in\mathbb{R}^{n}$, where we use $n=200$.

##### Construction of the Parameters.

To control the strong-convexity and smoothness of $\ell$, we specify intervals
$[m_{-},m_{+}],[L_{-},L_{+}]\subset(0,+\infty)$, and sample
$m_{1},...,m_{N}\overset{iid}{\sim}\mathrm{U}_{[m_{-},m_{+}]}$,
$L_{1},...,L_{N}\overset{iid}{\sim}\mathrm{U}_{[L_{-},L_{+}]}$. Then, the
matrices $A_{j}$, $j=1,...,N$, are created as _diagonal matrices_ with entries
$a_{ii}^{j}=\sqrt{m_{j}}+i\cdot\frac{\sqrt{L_{j}}-\sqrt{m_{j}}}{n}$,
$i=1,...,n$, that is, we linearly interpolate from $\sqrt{m_{j}}$ to
$\sqrt{L_{j}}$. Hence, the matrix $A_{j}^{T}A_{j}$ has smallest and largest
eigenvalue $m_{j}$ and $L_{j}$, respectively. To change the relative position
between the ellipsoid of the quadratic and the initialization, we randomize
the right-hand side by sampling
$b_{1},...,b_{N}\overset{iid}{\sim}\mathcal{N}(\mu,\Sigma)$, where we create
$\mu$ and $\Sigma=C^{T}C$ by sampling
$\mu_{i},C_{i,k}\overset{iid}{\sim}\mathrm{U}_{[-5,5]}$, $i,k=1,...,n$.

Figure 6: Upper left: The dashed lines represent the mean losses, while the
dotted lines represent the median losses over the iterations $x^{(i)}$,
$i=0,...,700$. The shaded regions represent the 10th to 90th percentile (80 %
of the test data). The learned algorithm $\mathcal{A}$ is shown in blue, while
HBF is shown in red. Upper right: The different lines indicate the cumulative
computation time the two algorithms need to solve all the test problems up to
a certain accuracy (in function values) measured by
$\ell(x^{(i)},\theta)<\varepsilon$. However, note that, since the learned
algorithm is not guaranteed to stop by itself, both algorithms are run for
maximally $n_{\rm{max}}=1\mathrm{e}4$ iterations. Lower left: Loss histogram
after $n_{\rm{train}}=350$ iterations with predicted PAC-bound. Lower right:
Estimated sublevel probability. The orange dashed line shows the point
estimate, while the orange solid line shows the $\mathrm{Beta}$-posterior
arising from a non-informative prior. The gray dotted lines indicate the
constraints $p_{l},p_{u}$ and show the feasible region as green line.

##### Baseline.

Assuming that it is not feasible to compute $L_{\theta}$ and $m_{\theta}$ for
each problem instance $\ell_{\theta}$ separately, the given class of functions
is $L_{+}$-smooth and $m_{-}$-strongly convex. Hence, as baseline we use
_heavy-ball with friction_ (HBF) (Polyak, 1964), whose update is given by
$x^{(k+1)}=x^{(k)}-\tau\nabla f(x^{(k)})+\beta\left(x^{(k)}-x^{(k-1)}\right)$,
where the optimal worst-case convergence rate is attained for
$\tau=\Bigl{(}\frac{2}{\sqrt{L_{+}}+\sqrt{\mu_{-}}}\Bigr{)}^{2},\beta=\Bigl{(}\frac{\sqrt{L_{+}}-\sqrt{\mu_{-}}}{\sqrt{L_{+}}+\sqrt{\mu_{-}}}\Bigr{)}^{2}$
(Nesterov, 2018).

##### Algorithm.

The algorithmic update of the learned algorithm $\mathcal{A}$ is visualized in
Figure 5 and consists of two blocks: The update-block combines the gradient
direction $d^{(k)}_{1}$, the momentum direction $d^{(k)}_{2}$, and their
“interaction” $d^{(k)}_{1}\odot d^{(k)}_{2}$ into the new update-direction
$d^{(k)}$, while the other block computes a step-size based on the
corresponding logarithmically transformed norms $n^{(1)}$ and $n^{(2)}$.
Interestingly, the theoretically redundant concatenation of two linear layers
(1x1 convolutional layers, resp.) does seem to increase training stability.

##### Results.

The upper left plot of Figure 6 shows that the learned algorithm outperforms
HBF by orders of magnitude and, despite being trained for $n_{\rm{train}}=350$
iterations, one can use it for at least 700 iterations in most cases ($\geq
90\%$). However, the mean indicates that there are single instances for which
instabilities occur, and, by comparing it to the median, one observes that the
mean is far from being representative of the “typical” performance. Further,
the algorithm performs well on very different orders of magnitude, ranging
from $1\mathrm{e}5$ to $1\mathrm{e}{-20}$. The upper right plot confirms that
also in terms of computation time the learned algorithm is way faster than
heavy-ball, and the lower left plot shows that the predicted PAC-bound is
quite tight. Lastly, the lower right plot shows that the algorithm will not
converge in all of the cases, yet satisfies the specified constraints $p_{l}$
and $p_{u}$.

###  6.2 Image Processing

We consider (gray-scale) _image denoising/deblurring_ with a regularizer given
as smooth approximation to the $L_{1}$-norm of the image derivative, that is,
problems of the form:

$\min_{x\in\mathbb{R}^{n}}\frac{1}{2}\|Ax-b\|^{2}+\lambda\sum_{i,j=1}^{n}\sqrt{(D_{h}x)_{i,j}^{2}+(D_{w}x)_{i,j}^{2}+\varepsilon^{2}}\qquad\lambda\in\mathbb{R},\
A,D_{h},D_{w}\in\mathbb{R}^{n\times n},\ b\in\mathbb{R}^{n}\,.$

The matrix $A$ describes the “blurring” of the image, while $D_{h}$ and
$D_{w}$ are the discrete image derivatives in h- and w-direction,
respectively, which are used to penalize local changes in the image. We use
images of height $N_{h}=250$ and width $N_{w}=\texttt{int}(0.75\cdot
N_{h})=187$. Thus, the dimension $n$ of the optimization space is given by
$n=46750$. Further, as parameters $\theta$ we use the observed image and the
regularization parameter, that is,
$\theta=(b,\lambda)\in\mathbb{R}^{n+1}=:\Theta$.

$d^{(k)}_{1}$$d^{(k)}_{2}$$d^{(k)}_{1}\odot
d^{(k)}_{2}$Conv2d(3,16,3,bias=F)ReLUConv2d(16,16,3,bias=F)ReLUConv2d(16,1,3,bias=F)$d^{(k)}$$n^{(k)}_{1}$$n^{(k)}_{2}$Linear(2,8,bias=F)ReLULinear(8,8,bias=F)ReLULinear(8,1,bias=F)$s^{(k)}$$x^{(k)}$$x^{(k+1)}:=x^{(k)}+s^{(k)}\cdot
d^{(k)}$ Figure 7: Update step of $\mathcal{A}$ for the image-processing
problems: The directions $d^{(k)}_{1}$, $d^{(k)}_{2}$ and $d^{(k)}_{1}\odot
d^{(k)}_{2}$ are inserted as different channels (in the shape of the image)
into the Conv2d-block, which performs a $3\times 3$-convolution. The scales
$n^{(1)}$, $n^{(2)}$ get transformed separately by the fully-connected block.

##### Construction of the Parameters.

Throughout, we use $\varepsilon=0.01$. For computational efficiency, the
matrices $A,D_{h},D_{w}$ are implemented through the convolution of the image
$x$ with a corresponding kernel (with reflective boundary conditions). For
$A$, we use a standard $(5\times 5)$-Gaussian kernel, while $D_{h}$ and
$D_{w}$ are given through the kernels:

$k_{h}=\begin{pmatrix}0&0&0\\\ 0&-1&0\\\
0&1&0\end{pmatrix}\in\mathbb{R}^{3\times 3}\qquad\text{and}\qquad
k_{w}=\begin{pmatrix}0&0&0\\\ 0&-1&1\\\
0&0&0\end{pmatrix}\in\mathbb{R}^{3\times 3}\,.$

Additionally, after blurring an image with $A$, we add centered Gaussian noise
$\varepsilon_{i,j}$ with standard deviation $\sigma=\frac{25}{256}$ to each
pixel, that is, $b_{i,j}=(Ax^{*})_{i,j}+\varepsilon_{i,j}$ with
$\varepsilon_{i,j}\overset{iid}{\sim}\mathcal{N}\left(0,\sigma\right)$,
$i=1,...,N_{h}$, $j=1,...,N_{w}$. The regularization parameters
$\lambda_{i}\in\mathbb{R}$, $i=1,...,N$, are given by sampling uniformly, that
is, $\lambda_{i}\overset{iid}{\sim}\mathrm{U}_{[\lambda_{-},\lambda_{+}]}$,
where we use $\lambda_{-}=0.05$ and $\lambda_{+}=0.5$.

##### Baseline.

Since the problem is smooth and convex, yet not strongly convex, the baseline
algorithm is given by the _accelerated gradient descent_ algorithm due to
Nesterov (1983). Its update is given by first computing
$y^{(k+1)}=x^{(k)}+\frac{t_{k}-1}{t_{k+1}}(x^{(k)}-x^{(k-1)})$ followed by
setting $x^{(k+1)}=y^{(k)}-\tau\nabla f(y^{(k+1)})$. We use the optimal
choices $t_{k+1}=\frac{1}{2}\left(1+\sqrt{1+4t_{k}^{2}}\right)$ and
$\tau=\frac{1}{L}$. Here, the smoothness constant $L$ is given by the largest
eigenvalue of $A^{T}A+\frac{\lambda}{\varepsilon}D^{T}D$, where
$D\in\mathbb{R}^{2n\times n}$ is given by “stacking” $D_{h}$ and $D_{w}$, that
is, $D=\begin{pmatrix}D_{h}&D_{w}\end{pmatrix}^{T}$.

Figure 8: Upper left: Dashed lines represent mean losses, while dotted lines
show median losses over the iterations $x^{(i)}$, $i=0,...,50$. The shaded