Tr ∇ *s*(*z*) . **0** −*P*(*θ*) *P*(*θ*) *A*(*θ*) / *ρ*(*z*) − *η*∇(*φ*(*z*)). **0 0 0** *P*(*θ*)2**Σ(***θ***)** / *ρ*(*z*) = Tr ∇ *s*(*z*) . **0** −*P*(*θ*) *P*(*θ*) *A*(*θ*) / *ρ*(*z*) − ∇(*φ*(*z*)). **0 0 0** *ηP*(*θ*)2**Σ(***θ***)** / *ρ*(*z*) = Tr ∇ *s*(*z*) . **0** −*P*(*θ*) *P*(*θ*) *A*(*θ*) / *ρ*(*z*) − ∇(*φ*(*z*)). **0 0 0** *A***(***θ***)** / *ρ*(*z*) = Tr ∇ *s*(*z*) − ∇(*φ*(*z*)). **0** −*P*(*θ*) *P*(*θ*) *A*(*θ*) / *ρ*(*z*) = 0

and then ∇*φ*(*z*) = *s*(*z*) proving Theorem 2.

#### **Appendix B. I-SGD Method Proofs and Details**

*Appendix B.1. Proof of Corollary 1*

The requirement *C***(***θ***)** 0 ∀*θ*, ensures that the injected noise covariance is valid. The composite noise matrix is equal to **Σ(***θ***)** = **Λ**. Since ∇**Σ(***θ***)** = ∇**Λ** = **0** and *ηP***(***θ***)** = **Λ**−1 by construction, then Theorem 1 is satisfied.

## *Appendix B.2. Proof of Optimality of* **Λ**

Our design choice is to select *λ*(*p*) = *β*(*p*). By the assumptions, the matrix *B***(***θ***)** is diagonal, and consequently *C***(***θ***)** = **Λ** − *B***(***θ***)** is diagonal as well. The preconditioner **Λ** must be chosen to satisfy the positive semidefinite constraint, i.e., *C***(***θ***)***ii* ≥ 0 ∀*i*, ∀*θ*. Equivalently, we must satisfy *λ*(*p*) − *bj*(*θ*) ≥ 0 ∀*j* ∈ *Ip*, ∀*p*, ∀*θ*, where *Ip* is the set of indices of parameters belonging to *pth* layer. By assumption 3, i.e., *β*(*p*) = ∑*k*∈*Ip bk*(*θ*), it is easy to show that *bj*(*θ*), *j* ∈ *Ip*, is upper bounded as *bj*(*θ*) ≤ *β*(*p*). To satisfy the positive semidefinite requirement in all cases the minimum valid set of *λ*(*p*) is then determined as *λ*(*p*) = *β*(*p*).

## *Appendix B.3. Algorithmic Details*

In this section, we provide further details about the practical implementation of the proposed scheme. At any (discrete) time instant a minibatch version of the gradient is computed that is distributed, according to the hypotheses of the main paper, as *g*(*θ*) ∼ *N*(∇ *f*(*θ*), 2*b*(*θ*). Since we assumed that the second-order moment is a good approximation of the variance, we can estimate *b*(*θ*) as 1 2 (*g*(*θ*) *g*(*θ*)). In practice, we found that the following running average estimation procedure to be the most robust

$$b(\theta) \leftarrow \mu b(\theta) + (1 - \mu)\frac{1}{2}(\mathbf{g}(\theta) \odot \mathbf{g}(\theta))\tag{A1}$$

where *μ* ∈ (0, 1]. In all experiments we considered *μ* = 0.5

After a warmup period, the various *λ*(*p*) , layer per layer, are estimated as *λ*(*p*) = ∑*k*∈*Ip bk*(*θ*) and kept constant until the end. The estimation procedure continues during sampling phase,

as the quantity *λ*(*p*) − *b*(*θ*) is necessary at every step. As the learning rate is derived as 2 *λ*(*p*) , we found that the usage of second-order moments instead of variances, and in certain cases temperature scaling, kept the simulated trajectories more stable.

#### **Appendix C. Methodology**

We hereafter present additional implementation details.

#### *Appendix C.1. Regression Tasks, with Simple Models*

For this set of experiments we considered , the BASELINE is obtained by running the ADAM optimizer for 20,000 steps with learning rate 0.01 and default parameters. At test time we use 100 samples to estimate the predictive posterior distribution, using Equation (3), for the *sampling* methods (I-SGD,SGLD,SGHMC,VSGD), with a keep-every value equal to 1000. The I-SGD and VSGD sampling methods are started from the BASELINE. For I-SGD we selected temperature 0.01, while for SGHMC and SGLD we do performed experiments for temperatures 1 and 0.01. We modified the implementation of VSGD as the original implementation produced unstable learning rates (as noticed also in [9]). A simple and effective solution we implement that we kept throughout the experimental campaigns is to divide the learning rate by the number of parameters (thus performing variational inference on a tempered version of the posterior). For SGLD the learning rate decay is the one suggested in [2], with initial and finial learning rate equal to 10−6 and 10−8 respectively. For MCD we collected 1000 samples with standard dropout rate of 0.5. All our experiments use 10-splits. The considered batch size is 64 for all methods.

#### *Appendix C.2. Classification Task,* CONVNET

For the LENET-5 on MNIST experiment, we do consider also the SWAG algorithm. At test time we use 30 samples for all methods. Baselines are again trained using ADAM optimizer for 20,000 steps with learning rate 0.01 and default parameters. For I-SGD and SGHMC we collected samples for the different temperatures of 1 and 0.01. SGLD has initial and final learning rates of 10−3 and 10−5. For all the sampling methods we do collect 100 samples with a keep-every of 10,000 steps. SWAG results are obtained by collecting the statistics over 300 epochs using ADAM optimizer and decreasing the learning rate every epoch in accordance with the original paper schedule [9]. DROP results are obtained by training the networks with SGD, with learning rate 0.005 and momentum 0.5. The number of collected samples for this method is 1000. The batch size for all the methods is 128.

As explained in the main text, we performed an ablation study on the considered baselines. In Table A1 we do report the results for the additional variants obtained by early stopping (10,000 iterations instead of 20,000) BASELINE S, to ablate overfitting, and BASELINE L, by training for 30,000 iterations. Finally, we include the best performing BASELINE R, obtained starting from BASELINE, reducing the learning rate by a factor of 10 and training for 10,000 more iterations.

**Table A1.** Baselines comparison for classification on MNIST dataset.

#### *Appendix C.3. Classification Task, Deeper Models*

We here report details for the RESNET-18 on CIFAR10 experiments. The BASELINE is obtained with ADAM optimizer with learning rate 0.01 decreased by a factor of 10 every 50 epochs for a total of 200 epochs and weight decay of 0.05. For this set of experiments no

temperature scaling was required. We could not find good hyperparameters for the SGLD scheme. Concerning I-SGD, SGHMC and VSGD the keep-every value is chosen as 10,000 and the number of collected samples is 30. For SWAG we used the default parameters described in [9]. Notice that for SWAG we performed the following ablation study: we trained the networks considering as loss function the joint log-likelihood and included or not the suggested weight decay of the original work [9]. From a purely Bayesian perspective no weight decay should be considered to be the information is implicit in the prior; however, we found that without the extra decay SWAG was not able to obtain competitive results. As underlined in Section 4.1, not necessarily a better posterior approximation translates into better empirical results.

## *Appendix C.4. Definition of the Metrics*

For regression datasets, we consider RMSE and MNLL. Consider a single datapoint *Ui* = (*xi*, *yi* ), with *xi* the input of the model and *yi* the true corresponding output. The output of the model, for a single sample of parameters *θj*, is *y*ˆ *θj* (*xi*). RMSE is defined as

1 *N N* ∑ *i*=1 ||*yi* − *μ*(*xi*)||2, where *μ*(*xi*) is the empirical mean 1 *NMC NMC* ∑ *j*=1 *y*ˆ *θj* (*xi*). MNLL is defined instead as ( 1 *N N* ∑ 1 2 log(2*πσ*2 *i* ) + 1 2 ||*yi*−*μ*(*xi*)||2 *σ*2 , where *σ*2 *i* is the empirical variance.

*i*=1 For classification datasets, we consider ACC,MNLL and entropy. Consider a single datapoint *Ui* = (*xi*, *yi*), with *xi* the input of the model and *yi* the true corresponding label. The output of the model, for a single sample of parameters *θj*, is the *Ncl* vector **p***θj* (*xi*). The

*i*

averaged probability vector for a single sample is **p**(*xi*) = 1 *NMC NMC* ∑ *i*=1 **p***θj* (*xi*).ACC is defined as 1 *N N* ∑ *i*=1 1(arg max **p**(*xi*) = *yi*). MNLL is computed as 1 *N N* ∑ *i*=1 log **p***yi* (*xi*) . Entropy, as stated in the main text, is instead computed according to 1 *N* ∑ *Ncl* ∑ **p***k*(*xi*)log(**p***k*(*xi*)) .

*N*

*i*=1

*k*=1

## **References**

## *Article* **PAC-Bayes Unleashed: Generalisation Bounds with Unbounded Losses**

**Maxime Haddouche 1, Benjamin Guedj 2,3,\*, Omar Rivasplata 2 and John Shawe-Taylor 2**

**Abstract:** We present new PAC-Bayesian generalisation bounds for learning problems with unbounded loss functions. This extends the relevance and applicability of the PAC-Bayes learning framework, where most of the existing literature focuses on supervised learning problems with a bounded loss function (typically assumed to take values in the interval [0;1]). In order to relax this classical assumption, we propose to allow the range of the loss to depend on each predictor. This relaxation is captured by our new notion of *HYPothesis-dependent rangE* (HYPE). Based on this, we derive a novel PAC-Bayesian generalisation bound for unbounded loss functions, and we instantiate it on a linear regression problem. To make our theory usable by the largest audience possible, we include discussions on actual computation, practicality and limitations of our assumptions.

**Keywords:** statistical learning theory; PAC-Bayes; generalisation bounds

## **Citation:** Haddouche, M.; Guedj, B.; Rivasplata, O.; Shawe-Taylor, J. PAC-Bayes Unleashed: Generalisation Bounds with Unbounded Losses. *Entropy* **2021**, *23*, 1330. https://doi.org/10.3390/ e23101330

Academic Editor: Boris Ryabko

Received: 22 August 2021 Accepted: 25 September 2021 Published: 12 October 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

## **1. Introduction**

Since its emergence in the late 1990s, the PAC-Bayes theory (see the seminal works of [1–3], the recent survey by [4] and work by [5]) has been a powerful tool to obtain generalisation bounds and to derive efficient learning algorithms. Generalisation bounds are helpful for understanding how a learning algorithm may perform on future similar batches of data. While the classical generalization bounds typically address the performance of individual predictors from a given hypothesis class, PAC-Bayes bounds typically address a randomized predictor defined by a distribution over the hypothesis class.

PAC-Bayes bounds were originally meant for binary classification problems [6–8], but the literature now includes many contributions involving any bounded loss function (without loss of generality, with values in [0; 1]), not just the binary loss. Our goal is to provide new PAC-Bayes bounds that are valid for unbounded loss functions, and thus extend the usability of PAC-Bayes to a much larger class of learning problems. To do so, we reformulate the general PAC-Bayes theorem of [9] and use it as basic building block to derive our new PAC-Bayes bound.

Some ways to circumvent the bounded range assumption on the losses have been explored in the recent literature. For instance, one approach consists of assuming a tail decay rate on the loss, such as sub-gaussian or sub-exponential tails [10,11]; however, this approach requires the knowledge of additional parameters. Some other works have also looked into the analysis for heavy-tailed losses, e.g., ref. [12] proposed a polynomial moment-dependent bound with *f*-divergences, while [13] devised an exponential bound that assumes the second (uncentered) moment of the loss is bounded by a constant (with a truncated risk estimator, as recalled in Section 4 below). A somewhat related approach was explored by [14], who do not assume boundedness of the loss, but instead control higher-order moments of the generalization gap through the Efron-Stein variance proxy. See also [5].

99

We investigate a different route here. We introduce the *HYPothesis-dependent rangE* (HYPE) condition, which means that the loss is upper-bounded by a term that depends on the chosen predictor (but does not depend on the data). Thus, effectively, the loss may have an arbitrarily large range. The HYPE condition allows us to derive an upper bound on the exponential moment of a suitably chosen functional, which, combined with the general PAC-Bayes theorem, leads to our new PAC-Bayes bound. To illustrate it, we instantiate the new bound on a linear regression problem, which additionally serves the purpose of illustrating that our HYPE condition is easy to verify in practice, given an explicit formulation of the loss function. In particular, we shall see in the linear regression setting that a mere use of the triangle inequality is enough to check the HYPE condition. The technical assumptions on which our results are based are comparable to those of the classical PAC-Bayes bounds; we state them in full detail, with discussions, for the sake of clarity and to make our work accessible.

**Our contributions are twofold.** (i) We propose PAC-Bayesian bounds holding with unbounded loss functions, therefore overcoming a limitation of the mainstream PAC-Bayesian literature for which a bounded loss is usually assumed. (ii) We analyse the bound, its implications, limitations of our assumptions, and their usability by practitioners. We hope this will extend the PAC-Bayes framework into a widely usable tool for a significantly wider range of problems, such as unbounded regression or reinforcement learning problems with unbounded rewards.

**Outline.** Section 2 introduces our notation and definition of the HYPE condition and provides a general PAC-Bayesian bound, which is valid for any learning problem complying with a mild assumption. For the sake of completeness, we present how our approach (designed for the unbounded case) behaves in the bounded case (Section 3). This section is not the core of our work, but rather serves as a safety check and particularises our bound to more classical PAC-Bayesian assumptions. We also provide numerical experiments. Section 4 introduces the notion of *softening functions* and particularises Section 2's PAC-Bayesian bound. In particular, we make explicit all terms in the right-hand side. Section 5.1 extends our results to linear regression (which has been studied from the perspective of PAC-Bayes in the literature, most recently by [15]). We also experimentally illustrate the behaviour of our bound. Finally, Section 6 presents, in detail, related works and Section 7 contains all proofs of the original claims we make in the paper.

#### **2. Framework and Preliminary Results**

The learning problem is specified by three variables (H, Z, -) consisting of a set H of predictors, the data space Z, and a loss function -: H×Z→ R+.

For a given positive integer *m*, we consider size-*m* datasets. The space of all possible datasets of this fixed size is S = Z*m*; an arbitrary element of this space is *s* = (*z*1, ... , *zm*). We denote *S* as a random dataset: *S* = (*Z*1, ... , *Zm*) where the random data points *Zi* are independent and sampled from the same distribution *μ* over Z. We call *μ* the data-generating distribution. The assumption that the *Zi*'s are *independent and identically distributed* is typically called the i.i.d. data assumption. It means that the random sample *S* (of size *m*) has distribution *μ*⊗*m* which is the product of *m* copies of *μ*.

For any predictor *h* ∈ H, we define the *empirical risk* of *h* over a sample *s*, denoted *Rs*(*h*), and the *theoretical risk* of *h*, denoted *R*(*h*), as:

$$R\_{\delta}(h) = \frac{1}{m} \sum\_{i=1}^{m} \ell(h, z\_i) \qquad \text{and} \qquad \mathcal{R}(h) = \mathbb{E}\_{\mu}[\ell(h, Z)]$$

respectively, where E*μ*[-(*h*, *Z*)] denotes the expectation with respect to *Z* ∼ *μ*. Finally, we define the *risk gap* Δ*s*(*h*) = *R*(*h*) − *Rs*(*h*) for any *h* ∈ H and *s* ∈ S. Often, Δ*s*(*h*) is referred to as the generalisation gap.

Notice that for a random dataset *S*, the empirical risk *RS*(*h*) is random, with expected value E*μ*⊗*m* [*RS*(*h*)] = *R*(*h*), where E*μ*⊗*m* the expectation under the distribution of the random sample *S*.

In general, E*μ*[·] denotes an expectation under the distribution *μ*. When we want to emphasize the role of the random variable *Z* ∼ *μ* we write E*Z*[·] or E*Z*∼*μ*[·] instead of E*μ*[·]. We use a similar convention for expectations related to any other distributions and random quantities. We now introduce the key concept to our analysis.

**Definition 1.** (HYPE). *A loss function* - : H×Z → R+ *is said to satisfy the hypothesisdependent range (*HYPE*) condition if there exists a function K* : H → R+\{0} *such that* sup*z*∈Z -(*h*, *z*) ≤ *K*(*h*) *for every predictor h. We then say that is* HYPE(*K*) *compliant.*

Let M+ 1 (H) be the set of probability distributions on H. We assume that all considered probability measures on H are defined on a fixed *σ*-algebra over H, while the notation M+ 1 (H) hides the *σ*-algebra, for simplicity. For *P*, *P* ∈ M+ 1 (H), the notation *P P* indicates that *P* is absolutely continuous with respect to *P* (i.e., *P* (*A*) = 0 if *P*(*A*) = 0 for measurable *A* ⊂ H). We write *P* ∼ *P* to indicate that *P P* and *P P* , i.e., these two distributions are absolutely continuous with respect to each other.

We now recall a result from Germain et al. [9]. Note that while implicit in many PAC-Bayes works (including theirs), we make it explicit that both the prior *P* and the posterior *Q* must be absolutely continuous with respect to each other. We discuss this restriction below.

**Theorem 1.** (Adapted from [9], Theorem 2.1.) *For any P* ∈ M+ 1 (H) *with no dependency on data, for any function F* : R+ × R+ → R*, define the exponential moment:*

$$\chi := \mathbb{E}\_S \mathbb{E}\_{h \sim P} \left[ \mathfrak{e}^{F(\mathbb{R}\_S(h), \mathbb{R}(h))} \right] \dots$$

*If F is convex, then for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over random samples S, simultaneously for all Q* ∈ M+ 1 (H) *such that Q* ∼ *P we have:*

$$F(\mathbb{E}\_{h \sim \mathbb{Q}}[R\_S(h)], \mathbb{E}\_{h \sim \mathbb{Q}}[R(h)]) \le \text{KL}(\mathbb{Q}||P) + \log\left(\frac{\chi}{\delta}\right).$$

The proof is deferred to Section 7.1. Note that the proof in [9] requires that *P Q*, although it is not explicitly stated; we highlight this in our own proof. While *Q P* is classical and necessary for the KL(*Q*||*P*) to be meaningful, *P Q* appears to be more restrictive. In particular, we have to choose *Q* such that it has the exact same support as *P* (e.g., choosing a Gaussian and a truncated Gaussian is not possible). However, we can still apply our theorem when *P* and *Q* belong to the same parametric family of distributions, e.g., both 'full-support' Gaussian or Laplace distributions, but these are just two examples and there are many others.

Note that Alquier et al. [10] (Theorem 4.1) adapted a result from Catoni [8], which only requires *Q P*. This comes at the expense of what Alquier et al. [10] (Definition 2.3) called a *Hoeffding's assumption*, which means that the exponential moment *χ* is assumed to be bounded by a function depending only on the hyperparameters (such as the dataset size *m* or parameters given by Hoeffding's assumption). Our analysis does not require this assumption, which might prove restrictive in practice.

Theorem 1 may be seen as a basis to recover many classical PAC-Bayesian bounds. For instance, *F*(*x*, *y*) = 2*m*(*x* − *y*)2, recovers McAllester's bound as recalled in [4] (Theorem 1). To get a usable bound, the outstanding task is to bound the exponential moment *χ*. Note that a previous attempt has been made in [11], as described in Section 6.1 below. Furthermore, under the assumption that the distribution *P* has no dependency on the data, we may swap the order of integration in the exponential moment thanks to Fubini-Tonelli's theorem and the positiveness of the exponential:

$$\chi = \mathbb{E}\_{h \sim P} \mathbb{E}\_S \left[ \mathfrak{e}^{F(\mathcal{R}\_S(h), \mathcal{R}(h))} \right] \dots$$

This is the starting point for the way that the exponential moment was handled in several works in the PAC-Bayes literature. Essentially, for a fixed *h*, one may upper-bound the innermost expectation (with respect to *S*) using standard exponential moment inequalities.

In this work, we will use Theorem 1 with *F*(*x*, *y*) = *mαD*(*x*, *y*), where *α* > 0, and *D* : R+ × R+ → R is a convex function. In this case, the high-probability inequality of the theorem takes the form:

$$\begin{split} \left| D\left( \mathbb{E}\_{\mathbb{H}\sim Q}[R\_S(h)], \mathbb{E}\_{\mathbb{H}\sim Q}[R(h)] \right) \right| &\leq \\ & \frac{1}{m^d} \left( \mathrm{KL}(Q||P) + \log \left( \frac{1}{\delta} \mathbb{E}\_{\mathbb{H}\sim P} \mathbb{E}\_S \, e^{m^d D\left(R\_S(h), R(h)\right)} \right) \right) . \end{split} \tag{1}$$

Our goal is to control E*S emαD*(*RS*(*h*),*R*(*h*)) for a fixed *h*, when *D*(*x*, *y*) = *y* − *x*. This will readily give us control on the exponential moment *χ*. To do so, we propose the following theorem:

**Theorem 2.** *Let h* ∈ H *be a fixed predictor and α* ∈ R*. If the loss function is* HYPE(*K*) *compliant, then for* Δ*S*(*h*) = *R*(*h*) − *RS*(*h*) *we have:*

$$\mathbb{E}\_{\mathcal{S}}\left[e^{m^{\alpha}\Delta\_{\mathcal{S}}(h)}\right] \le \exp\left(\frac{K(h)^2}{2m^{1-2\alpha}}\right).$$

**Proof.** Let *h* ∈ H. Then:

$$\begin{aligned} \mathbb{E}\_{\mathbb{S}}\left[e^{\mathsf{w}^{\mathsf{a}}\Delta\_{\mathbb{S}}(h)}\right] &= \mathbb{E}\left[\exp\left(m^{\mathsf{a}-1}\sum\_{i=1}^{m}\left(l(h,Z\_{i})-R(h)\right)\right)\right] \\ &= \mathbb{E}\left[\prod\_{i=1}^{m}\exp\left(m^{\mathsf{a}-1}\left(\ell(h,Z\_{i})-R(h)\right)\right)\right] \\ &= \prod\_{i=1}^{m}\mathbb{E}\left[\exp\left(m^{\mathsf{a}-1}\left(\ell(h,Z\_{i})-R(h)\right)\right)\right]. \end{aligned}$$

We now apply Hoeffding's lemma, for any *i* ∈ {1..*m*}, the random (in *Zi*) variable -(*h*, *Zi*) − *R*(*h*) is centered, taking values in [−*K*(*h*); *K*(*h*)], so that:

$$\mathbb{E}\left[\exp\left(m^{\mathfrak{a}-1}(\ell(h,Z\_i)-R(h))\right)\right] \le \exp\left(m^{2\mathfrak{a}-2}\frac{4K(h)^2}{8}\right).$$

and finally:

$$\mathbb{E}\_S\left[e^{m^a \Lambda\_S(h)}\right] \le \prod\_{i=1}^m \exp\left(m^{2a-2} \frac{4K(h)^2}{8}\right) = \exp\left(\frac{K(h)^2}{2m^{1-2a}}\right).$$

The strength of this result lies in the fact that *K*(*h*)2 *m*1−2*α* , is a decreasing factor in *m*, when *α* ≤ 1/2, and more generally, one can control how fast the exponential moment will explode when *m* grows by the choice of the hyperparameter *α*.

For convenient cross-referencing, we state the following rewriting of Theorem 1.

**Theorem 3.** *Let the loss be* HYPE(*K*) *compliant. For any P* ∈ M+ 1 (H) *with no data dependency, for any α* ∈ R *and for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over size-m random samples S, simultaneously for all Q such that Q* ∼ *P we have:*

$$\mathbb{E}\_{h\sim Q}[R(h)] \le \mathbb{E}\_{h\sim Q}[R\_S(h)] + \frac{1}{m^a} \left( \text{KL}(Q||P) + \log \frac{\mathbb{E}\_{h\sim P}\left[\exp\left(\frac{K(h)^2}{2m^{1-2a}}\right)\right]}{\delta} \right).$$

**Proof.** We first apply Theorem 1 with *F*(*x*, *y*) = *mα*(*y* − *x*). More precisely, we use Equation (1) with *D*(*x*, *y*) = *y* − *x*. We then conclude with Theorem 2.

#### **3. Safety Check: The Bounded Loss Case**

*3.1. Theoretical Results*

At this stage, the reader might wonder whether this new approach allows for the recovery of known results in the bounded case: the answer is yes.

In this section, we study the case where is bounded by some constant *C* ∈ R+ \ {0}. In other words, we consider the case that sup*h* sup*z* -(*h*, *z*) ≤ *C*. We provide a bound, valid for any choice of "priors" *P* and "posteriors" *Q* such that *P* ∼ *Q*, which is an immediate corollary of Theorem 3.

**Proposition 1.** *Let be* HYPE(*K*) *compliant, with K*(*h*) = *C constant, and let α* ∈ R*. Let P* ∈ M+ 1 (H) *be a distribution with no data dependency. Then, for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over random m-samples S, simultaneously for all Q* ∈ M+ 1 (H) *such that Q* ∼ *P we have:*

$$\mathbb{E}\_{h \sim Q}[R(h)] \le \mathbb{E}\_{h \sim Q}[R\_S(h)] + \frac{\mathsf{KL}(Q||P) + \log(1/\delta)}{m^a} + \frac{\mathsf{C}^2}{2m^{1-a}}.$$

**Remark 1.** *We provide Proposition 1 to evaluate the robustness of our approach. For instance, by comparing it with the PAC-Bayesian bound found in Germain et al. [11]. This discussion can be found in Section 6.1, where the bound from Germain et al. [11] is presented in detail.*

**Remark 2.** *At first glance, a naive remark: in order to control the rate of convergence of all the terms of the bound in Proposition 1 (as is often the case in classical PAC-Bayesian bounds), then the only case of interest is in fact α* = 1 2 *. However, one could notice that the factor C*2 *is not optimisable, while the KL is. In this way, if it appears that C*2 *is too big, in practice, one wants to have the ability to attenuate its influence as much as possible and this may lead us to consider α* < 1/2*. The following lemma answers this question.*

**Lemma 1.** *For any given K*1 > 0*, the function fK*1 (*α*) := *K*1 *mα* + *C*2 *m*1−*α reaches its minimum at*

$$\alpha\_0 = \frac{1}{2} + \frac{1}{2\log(m)}\log\left(\frac{2K\_1}{C^2}\right).$$

**Proof.** The explicit calculus of the *f K*1 and the resolution of *f K*1 (*α*) = 0 provides the result.

**Remark 3.** *Lemma 1 indicates that with a fixed "prior" P and "posterior" Q, taking K*1 = KL(*Q*||*P*) + log(1/*δ*)*, gives the optimised value of the bound in Proposition 1. We numerically show in Section 3.2 (first experiment there) that optimising α leads to significantly better results.*

Now the only remaining question is how to optimise the KL divergence. To do so, we may need to fix an "informed prior" to minimise the KL divergence with an interesting posterior. This idea has been studied by [16,17] and, more recently, by Mhammedi et al. [18], Rivasplata et al. [5], among others. We will adapt it to our problem in the simplest way.

We now introduce some additional notation. For a sample *s* = (*z*1, ... , *zm*) and *k* ∈ {1..*m*}, we define *s*≤*k* := {*z*1, ... , *zk*} and *s*>*k* := {*zk*+1, ... , *zm*}. Then, similarly, for a random sample *S*, we have the splits *S*≤*k* and *S*>*k*.

**Proposition 2.** *Let be* HYPE(*K*) *compliant, with constant K*(*h*) = *C, and α*1, *α*2 ∈ R*. Consider any "priors" P*1 ∈ M+ 1 (H) *(possibly dependent on S*>*m*/2*) and P*2 ∈ M+ 1 (H) *(possibly dependent* *on S*≤*m*/2*). Then, for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over random size-m samples S, simultaneously for all Q* ∈ M+ 1 (H) *such that Q* ∼ *P*1 *and Q* ∼ *P*2 *we have:*

$$\begin{split} \mathbb{E}\_{h \sim Q}[R(h)] &\leq \mathbb{E}\_{h \sim Q}[R\_S(h)] + \frac{1}{2} \left( \frac{\mathrm{KL}(Q || P\_1) + \log(2/\delta)}{(m/2)^{a\_1}} + \frac{\mathcal{C}^2}{2(m/2)^{1-a\_1}} \right) \\ &+ \frac{1}{2} \left( \frac{\mathrm{KL}(Q || P\_2) + \log(2/\delta)}{(m/2)^{a\_2}} + \frac{\mathcal{C}^2}{2(m/2)^{1-a\_2}} \right). \end{split}$$

**Proof.** Let *P*1, *P*2, *Q* be as stated in Proposition 2. We first notice that by using Proposition 1 on the two halves of the sample, we obtain, with a probability of at least 1 − *δ*/2:

$$\mathbb{E}\_{\hbar \sim Q}[R(\hbar)] \le \mathbb{E}\_{\hbar \sim Q} \left[ \frac{1}{m/2} \sum\_{i=1}^{m/2} \ell(\hbar, Z\_i) \right] + \frac{\mathrm{KL}(Q || P\_1) + \log(2/\delta)}{(m/2)^{a\_1}} + \frac{\mathcal{C}^2}{2(m/2)^{1-a\_1}}$$

and also with probability at least 1 − *δ*/2:

$$\mathbb{E}\_{h\sim Q}[\mathbb{R}(h)] \le \mathbb{E}\_{h\sim Q} \left[ \frac{1}{m/2} \sum\_{i=1}^{m/2} \ell(h, Z\_{m/2+i}) \right] + \frac{\mathrm{KL}(Q || P\_2) + \log(2/\delta)}{(m/2)^{a\_2}} + \frac{\mathcal{C}^2}{2(m/2)^{1-a\_2}}.$$

Hence, with a probability of at least 1 − *δ*, both inequalities hold, and the result follows by adding them and dividing by 2.

**Remark 4.** *One can notice that the main difference between Proposition 2 and Proposition 1 lies in the implicit PAC-Bayesian paradigm that our priors must not depend on the data. With this last proposition, we implicitly allow P*1 *to depend on S*>*m*/2 *and P*2 *on S*≤*m*/2*, which can in practice lead to far more accurate priors. We numerically show this fact in Section 3.2's second experiment. Note that this idea is not new and has been studied, for instance, in [19] for the specific case of SVMs.*

#### *3.2. Numerical Experiments*

Our experimental framework has been inspired by the work of [18].

**Settings.** We generate synthetic data for classification, and we are using the 0–1 loss. The data space is Z = X ×Y = R*d* × {0, 1} with *d* ∈ N. The set of predictors H is parameterised with *d*-dimensional 'weight' vectors: H = {*hw* : X→Y| *w* ∈ R*d*}. For simplicity, we identify *hw* with *w* and we also identify the space H, with the weight space W = R*d*. For *z* = (*x*, *y*) ∈ Z and *w* ∈ W, we define the loss as -(*w*, *z*) := |1 *φ*(*wx*) > 1/2 − *y*|, where *φ*(*r*) = 1 1+*e*−*r* . We want to learn an optimised predictor given a dataset *S* = (*Zi*)*i*=1..*m* where *Zi* = (*Xi*,*Yi*). To do so, we use *regularised logistic regression* and compute:

$$\psi(S) := \arg\min\_{w \in \mathcal{W}} \lambda \frac{||w||^2}{2} - \frac{1}{m} \sum\_{i=1}^{m} y\_i \log \left( \phi(w^\top \mathbf{x}\_i) \right) + (1 - y\_i) \log \left( 1 - \phi(w^\top \mathbf{x}\_i) \right) \tag{2}$$

where *λ* is a fixed regularisation parameter.

We also restrict the probability distributions (over W = R*d*), considered for this learning problem. We consider the Gaussian distribution N (*w*, *σ*2 *Id*) with centre *w* ∈ R*d* and diagonal covariance *σ*2 *Id* ∈ R*d*×*d* with *σ*2 > 0.

**Parameters.** We set *δ* = 0.05, *λ* = 0.01. We approximately solve Equation (2) by using the minimize function of the optimisation module in Python, with the Powell method. To approximate gaussian expectations, we use Monte-Carlo sampling.

**Synthetic data.** We generate synthetic data for *d* = 10 according to the following process: for a fixed sample size *m*, we draw *X*1, ..., *Xm* under the multivariate Gaussian distribution N (0, *Id*) and for each *i* we compute the label if *Xi* as: *Yi* = 1{*φ*(*w*∗*xi*) > 1/2} where *w*∗ is the vector formed by the *d* first digits of the number *π*.

**Normalisation trick.** Given the predictors shape, we notice that for any *w* ∈ W:

$$\mathbb{1}\{\phi(w^\top x) > 1/2\} = 1 \quad \Leftrightarrow \quad \frac{1}{1 + \exp(-w^\top x)} > \frac{1}{2} \quad \Leftrightarrow \quad w^\top x < 0.1$$

Thus, the value of the prediction is exclusively determined by the sign of the inner product, and this quantity is definitely not influenced by the norm of the vector. Then, for any sample *S*, we call the **normalisation trick** the fact of considering *w*ˆ(*S*)/||*w*ˆ(*S*)|| instead of *w*ˆ(*S*) in our calculations. This process will not deteriorate the quality of the prediction and will considerably enhance the value of the KL divergence.

#### 3.2.1. First experiment

Our goal here is to highlight the point discussed in Remark 2, e.g., the influence of the parameter *α* in Proposition 1. We arbitrarily fix *σ*2 0 = 1/2, and define our *naive prior* as *P*0 = N (0, *σ*2 0 *Id*). For a fixed dataset *S*, we define our posterior as *P*(*S*) := N (ˆ *h*(*S*), *σ*2 *Id*), with *σ*2 ∈ {1/2, ... , 1/2*J* } (for *J* = log2(*m*)) such that it is minimising the bound among candidates. We computed two curves: first, Proposition 1 with *α* = 1/2 second, Proposition 1 again with *α* equals to the value proposed in Lemma 1. Notice that to compute this last bound, we first optimised our choice of posterior with *α* = 1/2 and then optimised *α*, to be consistent with Lemma 1. Indeed, we proved this lemma by assuming that the KL divergence was already fixed, hence our optimisation process is in two steps. Note that we chose to apply the normalisation trick here, we then obtained the left curve of Figure 1.

**Discussion.** From this curve, we formulate several remarks. First, we remark on this specific case, our theorem provides a tight result in practice (with an error rate lesser than 10% for the bound with optimised alpha). Second, we can now confirm that choosing an optimised *α* leads to a tighter bound. In further studies, it will be relevant to adjust *α* with regards to the different terms of our bound instead of looking for an identical convergence rate for all terms.

#### 3.2.2. Second Experiment

We now study Proposition 2 to see if an informed prior effectively provides a tighter bound than a naive one. We will use the notations introduced in Proposition 2. For a dataset *S*, we define *w*1(*S*) = *w*(*S*>*m*/2) as the vector resulting from the optimisation of Equation (2) on *S*>*m*/2. Similarly, we define *w*2(*S*) := *w*(*S*≤*m*/2). We arbitrarily fix *σ*2 0 = 1/2, and define our *informed priors* as: *P*1 = N (*w*1(*S*), *σ*2 0 *Id*) and *P*2 = N (*w*2(*S*), *σ*2 0 *Id*). Finally, we define our posterior as *P*(*S*) := N (*w*ˆ(*S*), *σ*2 *Id*), with *σ*2 ∈ {1/2, ..., 1/2*J* } (for *J* = log2(*m*)) with *σ*2 optimising the bound among the same candidate than the first experiment. We computed two curves: first, Proposition 1 with *α* optimised accordingly to Lemma 1 secondly, Proposition 2 with *α*1, *α*2 optimised as well, and informed priors as defined above. We chose to not apply the normalisation trick here, we then obtained the right curve of Figure 1.

**Discussion.** It is clear, that with this framework, having an informed prior is a powerful tool to enhance the quality of our bound. Notice that we voluntarily chose to not apply the normalisation trick here. The reason is that this trick appears to be too powerful in practice, and applying it leads to counterproductive results; to highlight our point: the bound without informed prior would be tighter than the one with informed prior. Furthermore, this trick is linked to the specific structure of our problem and is not valid for any classification problem. Thus, the idea of providing informed priors remains an interesting tool for most cases.

**Figure 1.** Above, result of the first experiment which highlight the importance of optimising *α*. Below, result of the second experiment which show how effective an informed prior is.

#### **4. PAC Bayesian Bounds with Smoothed Estimator**

We now move on to control the right-hand side term in Theorem 3 when *K* is not constant. A first step is to consider a transformed estimate of the risk, inspired by the truncated estimator from [20], also used in [21], and more recently in [13]. The following is inspired by the results of [13], which we summarise in Section 6.

The idea is to modify the estimator *RS*(*h*) for any *h* by introducing a threshold *t* and a function *ψ* which will attenuate the influence of the empirical losses (-(*h*, *Zi*))*i*=1..*m* that exceed *t*.

**Definition 2.** *ψ*-risks. *For every t* > 0*, ψ* : R+ → R+*, for any h* ∈ H*, we define the empirical ψ-risk RS*,*ψ*,*t and the theoretical ψ-risk Rψ*,*t as follows:*

$$R\_{S, \psi, t}(h) := \frac{t}{m} \sum\_{i=1}^{m} \psi\left(\frac{\ell(h, Z\_i)}{t}\right) \quad \text{and} \quad R\_{\psi, t}(h) = \mathbb{E}\_{\mu} \left[t \, \psi\left(\frac{\ell(h, Z)}{t}\right)\right]$$

*where Z* ∼ *μ. Notice that* E*S* 7 *RS*,*ψ*,*t*(*h*) 8 = *Rψ*,*t*(*h*)*.*

We now focus on what we call *softening functions*, i.e., functions that will temper high values of the loss function -.

**Definition 3.** (Softening function). *We say that ψ* : R+ → R+ *is a softening function if:*

*We let* F *denote the set of all softening functions.*

**Remark 5.** *Notice that those three assumptions ensure that ψ is continuous at* 1*. For instance, the functions f* : *x* → *x*1{*x* ≤ 1} + 1{*x* > 1} *and g* : *x* → *x*1{*x* ≤ 1} + (2 √*x* − 1)1{*x* > 1} *are in* F*. In Section 6 we compare these softening functions and those used by Holland [13].*

Using *ψ* ∈ F, for a fixed threshold *t* > 0, the softened loss function *tψ* -(*h*,*z*) *t* verifies for any *h* ∈ H, *z* ∈ Z:

$$t\,\psi\left(\frac{\ell(h,z)}{t}\right) \le t\,\psi\left(\frac{K(h)}{t}\right).$$

because *ψ* is non-decreasing. In this way, the exponential moment in Theorem 3 can be far more controllable. The trade-off lies in the fact that softening - (instead of taking directly -) will deteriorate our ability to distinguish between two bad predictions when both of them are greater than *t*. For instance, if we choose *ψ* ∈ F such as *ψ* = 1 on [1; +∞) and *t* > 0, if *ψ*(-(*h*, *z*)/*t*) = 1 for a certain pair (*h*, *z*), then we cannot tell how far -(*h*, *z*) is from *t* and we only can affirm that -(*h*, *z*) ≥ *t*.

We now move on to the following lemma, which controls the shortfall between E*h*∼*Q*[*R*(*h*)] and E*h*∼*Q*[*Rψ*,*t*(*h*)] for all *Q* ∈ M+ 1 (H), for a given *ψ* and *t* > 0. To do that, we assume that *K* admits a finite moment under any posterior distribution:

$$\forall \mathcal{Q} \in \mathcal{M}\_1^+(\mathcal{H}), \ \mathbb{E}\_{\hbar \sim \mathcal{Q}}[\mathcal{K}(\hbar)] \prec +\infty. \tag{3}$$

For instance, in the case of H identified with a weight space W = R*N*, and if *K* is polynomial in ||*w*|| (where ||.|| denotes the Euclidean norm), then this assumption holds if we consider Gaussian priors and posteriors.

**Lemma 2.** *Assume that Equation (3) holds, and let ψ* ∈ F*, Q* ∈ M+ 1 (H), *t* > 0*. We have:*

$$\mathbb{E}\_{h\sim Q}[\mathcal{R}(h)] \le \mathbb{E}\_{h\sim Q}[\mathcal{R}\_{\emptyset,t}(h)] + \mathbb{E}\_{h\sim Q}[\mathcal{K}(h)\mathbbm{1}\{\mathcal{K}(h)\geq t\}].$$

**Proof.** Let *ψ* ∈ F, *Q* ∈ M+ 1 (H), *t* > 0. We have, for *h* ∈ H :

$$\begin{aligned} &R(h) - R\_{\psi,t}(h) \\ &= \mathbb{E}\_{Z \sim \mu} \left[ \ell(h, Z) - t\psi\left(\frac{\ell(h, Z)}{t}\right) \right], \end{aligned}$$

and using that ∀*x* ∈ [0, 1], *ψ*(*x*) = *x*,

$$=\mathbb{E}\_{Z\sim\mu}\left[\left(\ell(h,Z)-t\psi\left(\frac{\ell(h,Z)}{t}\right)\right)\mathbb{1}\{\ell(h,Z)\geq t\}\right]$$

while using that -(*h*, *z*) ≤ *K*(*h*),

$$=\mathbb{E}\_{Z\sim\mu}\left[\left(\ell(h,Z)-t\psi\left(\frac{\ell(h,Z)}{t}\right)\right)\mathbb{1}\{\ell(h,Z)\geq t\}\mathbb{1}\{K(h)\geq t\}\right]$$

and continuing:

≤ E*Z*∼*μ*[-(*h*, *Z*)1{-(*h*, *Z*) ≥ *t*}]1{*K*(*h*) ≥ *t*} (*ψ* ≥ 0) ≤ *K*(*h*)P*Z*∼*μ*{-(*h*, *Z*) ≥ *t*}1{*K*(*h*) ≥ *t*} ( -(*h*, *Z*) ≤ *K*(*h*))

Finally, by crudely bounding the probability by 1, we get:

$$\mathcal{R}(h) \le \mathcal{R}\_{\emptyset, t}(h) + \mathcal{K}(h)\mathbf{1}\{\mathcal{K}(h) \ge t\}.$$

Hence the result by integrating over H with respect to *Q*.

Finally we present the following theorem, which provides a PAC-Bayesian inequality bounding the theoretical risk by the empirical *ψ*-risk for *ψ* ∈ F.

**Theorem 4.** *Let be* HYPE(*K*) *compliant, and assume K satisfies Equation (3). Then for any P* ∈ M+ 1 (H) *with no data dependency, for any α* ∈ R*, for any ψ* ∈ F *and for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over size-m random samples S, simultaneously for all Q such that Q* ∼ *P we have:*

$$\begin{split} \mathbb{E}\_{h \sim Q} [\mathcal{R}(h)] &\leq \mathbb{E}\_{h \sim Q} \left[ \mathbb{R}\_{S, \psi, t}(h) \right] + \mathbb{E}\_{h \sim Q} [K(h) \mathbb{1} \{ K(h) \geq t \}] \\ &+ \frac{\mathrm{KL}(Q || P) + \log \left( \frac{1}{\delta} \right)}{m^{a}} \\ &+ \frac{1}{m^{a}} \log \left( \mathbb{E}\_{h \sim P} \left[ \exp \left( \frac{t^{2}}{2m^{1-2a}} \psi \left( \frac{K(h)}{t} \right)^{2} \right) \right] \right) . \end{split}$$

.

**Proof.** Let *ψ* ∈ F, we define the *ψ*-loss:

$$\ell\_2(h, z) = t\psi\left(\frac{\ell(h, z)}{t}\right).$$

Since *ψ* is non decreasing, we have for all (*h*, *z*) ∈H×Z:

$$\ell\_2(h, z) \le t \psi\left(\frac{K(h)}{t}\right) := K\_2(h).$$

Thus, we apply Theorem 3 to the learning problem defined with -2: for any *α* and *δ* ∈ (0, 1), with probability at least 1 − *δ* over size-*m* random samples *S*, simultaneously for all *Q* such that *Q* ∼ *P* we have:

E*h*∼*Q* 7 *Rψ*,*t*(*h*) 8 ≤ E*h*∼*Q* 7 *RS*,*ψ*,*t*(*h*) 8 + KL(*Q*||*P*) + log1 *δ mα* + 1 *mα* log E*h*∼*P* . exp *K*2(*h*)2 2*m*1−2*α* /.

We then add E*h*∼*Q*[*K*(*h*)1{*K*(*h*) ≥ *t*}] on both sides of the latter inequality and apply Lemma 2.

**Remark 6.** *Notice that the function ψ* : *x* → *x*1{*x* ≤ 1} + 1{*x* > 1} *is such that for any given prior P we have* E*h*∼*P* . exp *t* 2 2*m*1−2*α ψ K*(*h*) *t* 2 / < +∞*. So the exponential moment can be controlled with a good choice of ψ. Thus the strength of Theorem 4 is to provide a PAC-Bayesian bound valid for any set of posterior measures verifying Equation (3). The choice of ψ minimising the bound is still an open problem.*

#### **5. The Linear Regression Problem**

#### *5.1. Theoretical Result*

We now focus on the celebrated linear regression problem and see how our theory translates to that particular learning problem. We assume that the data is a size-*m* random sample *S* = (*Zi*)*i*=1..*m* where the *Zi* are i.i.d. drawn from the distribution *μ*, and *Zi* = (*Xi*,*Yi*) with *Xi* ∈ R*N*, *Yi* ∈ R.

Our goal here is to find the most accurate predictor *hw* (with *w* ∈ R*N*), with respect to the loss function -(*hw*, *z*) = |*w*, *x* − *y*|, where *z* = (*x*, *y*). We will make the following mild assumption: there exists *B*, *C* ∈ R+\{0} such that for all *z* = (*x*, *y*) drawn under *μ*:

$$||x|| \le B \quad \text{and} \quad |y| \le C$$

where ||.|| is the norm associated to the classical inner product of R*N*. Under this assumption we note that for all *z* = (*x*, *y*) drawn according to *μ*, we have:

$$\ell(\mathfrak{h}\_w, z) = |\langle w, x \rangle - y| \le |\langle w, x \rangle| + |y| \le ||w||. ||x|| + |y| \le B||w|| + C.$$

Thus we define *K*(*hw*) = *B*||*w*|| + *C* for *w* ∈ R*N*. If we first restrict ourselves to the framework of Section 2, we want to use Theorem 3 and doing so, our goal is to bound *ξ* := E*w*∼*P* exp *K*(*w*)2 2*m*1−2*α* . The shape of *K* invites us to consider a Gaussian prior. Indeed, we notice that if *P* = N (0, *σ*2*IN*) with 0 < *σ*2 < *m*1−2*α B*2 , then *ξ* < +∞. Notice that we cannot take just any Gaussian prior, however with a small *α*, the condition 0 < *σ*2 < *m*1−2*α B*2 may become quite loose. Thus, we have the following:

**Theorem 5.** *Let α* ∈ R *and N* ≥ 6*. Assume that the loss is* HYPE(*K*) *compliant with K*(*h*) = *B*||*h*|| + *C, with B* > 0, *C* ≥ 0*. For a prior distribution, consider any Gaussian P* = N (0, *σ*2*IN*) *with σ*2 = *t m*1−2*α B*2 *,* 0 < *t* < 1*. Then, for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over size-m random samples S, simultaneously for all Q* ∈ M+ 1 (H) *such that P* ∼ *Q we have:*

$$\begin{split} \mathbb{E}\_{h \sim Q}[R(h)] &\leq \mathbb{E}\_{h \sim Q}[R\_S(h)] + \frac{\mathrm{KL}(Q || P) + \log(2/\delta)}{m^a} + \frac{\mathcal{C}^2}{2m^{1-a}} \left( 1 + f(t)^{-1} \right) \\ &+ \frac{N}{m^a} \left( \log \left( 1 + \left( \frac{\mathcal{C}}{\sqrt{2f(t)m^{1-2a}}} \right) \right) + \log \left( \frac{1}{\sqrt{1-t}} \right) \right) \end{split}$$

*where f*(*t*) = 1−*t t .*

The proof is deferred to Section 7.2. To compare our result with those found in the literature, we can fix *α* = 1/2. Doing so, we lose the dependency in *m* for the choice of the variance of the prior (which now only depends on *B*), but we recover the classic decreasing factor 1/√*m*.

**Remark 7.** *Notice that for now we did not use Section 4, even if we could (because K is polynomial in* ||*w*|| *and we consider Gaussian priors and posteriors, so Equation (3) is satisfied). Doing so, we obtained a bound which appears to depend linearly on the dimension N. In practice, N may be too big, and in this case, introducing an adapted softening function ψ (one can think for instance of ψ*(*x*) = *x*1{*x* ≤ 1} + 1{*x* > 1}*) is a powerful tool to attenuate the weight of the exponential moment. This also extends the class of authorised Gaussian priors by avoidance, to stick with a variance σ*2 = *t m*1−2*α B*2 *,* 0 < *t* < 1*.*

#### *5.2. Numerical Experiment*

#### 5.2.1. Setting

In this section we apply Theorem 5 on a concrete linear regression problem. The situation is as follows: we want to approximate the function *f*(*x*) = '*w*∗, *x* , where *w*∗ ∈ R*d*. We assume that W = [−*c*, *c*] *d* so that *w*∗ lies in an hypercube centred at 0 of half-side *c* > 0, i.e., the set {(*wi*)*i*=1,...,*d* | ∀*i*, |*wi*| ≤ *c*}. Doing so we have ||*w*∗|| ≤ *c* √*d*.

Furthermore, we assume that input data are drawn inside a hypercube of half-side *e* > 0, i.e., X = [−*e*,*e*] *d*. Doing so we have for any data *x*, ||*x*|| ≤ *e* √*d*.

For any data *x* ∈ R*d*, we define *y* = *f*(*x*). As before, we identify the hypothesis set H with the weight space W = R*d*. As described in Section 5.1, we set -(*hw*, *x*, *y*) = |*w*, *x* − *y*|. We then remark that for any (*w*, *x*, *y*):

$$\begin{aligned} \ell(h\_{w^\*} \mathbf{x}, y) &\leq |\langle w, \mathbf{x} \rangle| + |y| \leq ||w|| ||\mathbf{x}|| + |\sqrt{\langle w^\*, \mathbf{x} \rangle}| \\ &\leq \varepsilon \sqrt{d} ||w|| + \sqrt{||w^\*|| |. ||\mathbf{x}||} \leq \varepsilon \sqrt{d} ||w|| + \sqrt{c \sqrt{d} \, \boldsymbol{\omega} \sqrt{d}} \\ &\leq \varepsilon \sqrt{d} ||w|| + \sqrt{c \underline{d}} .\end{aligned}$$

Then we can define *B* = *e* √*d* and *C* = √*cde* to apply Theorem 5. We restrict (as before) the class of distributions over W to be *d*-dimensional Gaussians:

$$\left\{ \mathcal{N}(w, \sigma^2 I\_d) \mid w \in \mathcal{H}, \sigma^2 \in \mathbb{R}^+ \right\},$$

which is the set of candidate distributions for this learning problem. Recall that in practice, given a fixed *α* ∈ R, we are only allowed to consider priors such that their variance *σ*2 ∈ 0; *m*1−2*α B*2 . We want to learn an optimised predictor (posterior) given a random dataset *S* = ((*Xi*,*Yi*))*i*=1,...,*m*. To do so, we consider synthetic data.

#### 5.2.2. Synthetic Data

We draw *w*∗ under a Gaussian (with mean 0 and standard deviation equal to 5) truncated to the hypercube centered at 0 of the half-side *c* > 0. We generate synthetic data according to the following process: for a fixed sample size *m*, we draw *X*1, ... , *Xm* under a Gaussian (with mean 0 and standard deviation equal to 5) truncated to the hypercube centered at 0 of the half-side *e* > 0.

#### 5.2.3. Experiment

First, we fix *c* = *e* = 10. Our goal here is to obtain a generalisation bound on our problem. We fix arbitrarily, for a fixed *α* ∈ R, *t*0 = 1/2 and *σ*2 0 = *t*0 *m*1−2*α B*2 and we define our *naive prior* as *P*0 = N (0, *σ*2 0 *Id*). For a given dataset *S*, we define our posterior as *Q*(*S*) := N (*w*ˆ(*S*), *σ*2 *Id*), with *σ*2 ∈ {*σ*2 0/2, ..., *σ*2 0/2*J* } (*J* = log2(*m*)), such that it is minimising the bound among candidates. Note that all the previously defined parameters are dependent on *α*, which is why we choose *α* ∈ {*i*/step | 0 ≤ *i* ≤ step} for step a fixed integer (in practice step = 8 or 16) and we take the value of *α* minimising the bound among the candidates as well. Figure 2 contains two figures, one with *d* = 10, the other with *d* = 50. On each figure are computed the right-hand side term in Theorem 5 with an optimised *α* for each step.

#### 5.2.4. Discussion

To the the best of our knowledge, this is the first attempt to numerically compute PAC-Bayes bounds for unbounded problems, making it impossible to compare to other results. We stress, however, that obtaining numerical values for the bound without assuming a bounded loss is a significant first step. Furthermore, we consider a rather hard problem: *f* is not linear, so we cannot rely on a linear approximation fitting perfectly data, and the larger the dimension, the larger the error, as illustrated by Figure 2. Thus, for any posterior *Q*, the quantity E*h*∼*Q*[*R*(*h*)] is potentially large in practice and our bound might not be tight. Finally, notice that optimising *α* (instead of taking *α* = 1/2 to recover a classic convergence rate) leads to a significantly better bound. A numerical example of this assertion is presented in Section 3.2. We aim to conduct further studies to consider the convergence rate as an hyperparameter to optimise, rather than selecting the same rate for all terms in the bound.

**Figure 2.** Evaluation of the right hand side in Theorem 5 with *d* = 10 and *d* = 50.

## **6. Existing Work**

*6.1. Germain et al., 2016*

In Germain et al. [11] (Section 4), a PAC-Bayesian bound has been provided for all *sub-gamma* losses with a variance *t* 2 and scale parameter *c* > 0, under a data distribution *μ* and a prior *P*, i.e. losses such that for every *λ* ∈ 0, 1 *c* the following is satisfied:

$$\log\left(\frac{1}{\delta}\mathbb{E}\_{h\sim P}\mathbb{E}\_{\mathcal{S}}\,e^{\lambda\left(R(h)-R\_{\mathcal{S}}(h)\right)}\right) \leq \frac{t^2}{c^2}(-\log(1-c\lambda)-\lambda c) \leq \frac{\lambda^2 t^2}{2(1-c\lambda)}.$$

Note that a sub-gamma loss (with regards to *μ* and *P*) is potentially unbounded. Germain et al. then propose the following PAC-Bayesian bound:

**Theorem 6.** *Ref. [11]. If the loss is sub-gamma with a variance t* 2 *and scale parameter c, under the data distribution μ and a fixed prior P* ∈ H*, then for any δ* ∈ [0; 1]*, with probability* 1 − *δ over size-m random samples, simultaneously for all Q P we have:*

$$\mathbb{E}\_{h \sim Q}[R(h)] \le \mathbb{E}\_{h \sim Q}[R\_S(h)] + \frac{\text{KL}(Q||P) + \log(1/\delta)}{m} + \frac{t^2}{2(1-c)}$$

.

Theorem 6 will be quoted several times in this paper given that it is a concrete PAC Bayesian bound provided with the will to overcome the constraint of a bounded loss. It is also one of the only one found in the literature.

Can we apply this theorem to the bounded case? The answer is yes: we remark that thanks to Hoeffding's lemma, if is bounded by *C* > 0, then for any *h* ∈ H it holds that *RS*(*h*) − *R*(*h*) ∈ [−*C*, *C*] almost surely. So, ∀*λ* ∈ R, logE*z*∼*μ eλ*(*R*(*h*)−*RS*(*h*) ≤ *λ*2*C*2 2 . Therefore, for any prior *P*, we have:

$$\log \mathbb{E}\_{h \sim P} \mathbb{E}\_{z \sim \mu} \left[ \mathcal{e}^{\lambda(R(h) - R\_S(h))} \right] \le \frac{\lambda^2 C^2}{2}.$$

Thus, is sub-gamma with variance *C*2 and scale parameter 0. Then, Theorem 6 can be applied with *t* 2 = *C*2, *c* = 0.

**Comparison with Proposition 1.** We remark that by taking *K* = *C* and *α* = 1 in Proposition 1, we are recovering Theorem 6. However, our approach allows us to say that if we can obtain a more precise form of *K* such that ∀*h* ∈ H, *K*(*h*) ≤ *C* and *K* is non-constant, Theorem 3, will ensure that:

$$\frac{1}{m^{\alpha}}\log\left(\mathbb{E}\_{h\sim P}\left[\exp\left(\frac{K(h)^2}{2m^{1-2\alpha}}\right)\right]\right) \le \frac{C^2}{2m^{1-\alpha}}.$$

Thus, having precise information on the behavior of the loss function -, with regards to the predictor *h*, allows us to obtain a tighter control of the exponential moment, and hence a tighter bound.

**Remark 8.** *We can see that Theorem 6 cannot control the factor C*2/2*. However, Ref. [11] remarked on this apparent weakness and partially corrected this issue [11] (Section 4, Equations (13) and (14)). Indeed, they proposed to balance the influence of m between the different terms of the PAC-Bayes bound by providing the same convergence rate in* 1/√*m to all terms.*

*We can then see Proposition 1 as a proper generalisation of Germain et al. [11] (Section 4, Equations (13) and (14)). Indeed, our bound exhibits properly the influence of the parameter α. Thus, we understand (and Lemma 1 proves it) that the choice of α deserves a study in itself in the way it is now a parameter of our optimisation problem. This fact has already been highlighted in Alquier et al. [10] (Theorem 4.1) (where λ* := *mα).*

## *6.2. Holland, 2019*

In [13], Holland proposed a PAC Bayesian inequality with unbounded loss. For that, he introduced a function *ψ* verifying a few specific conditions, different to those used in Section 4 to define our set of softening functions. Indeed, he considered a function *ψ* such that:

$$-\log\left(1-u+\frac{u^2}{b}\right) \le \psi(u) \le \log\left(1+u+\frac{u^2}{b}\right).\tag{4}$$

We remark that, as Holland did, we supposed that our softening functions are nondecreasing. We chose softening functions to be equal to the identity function (*x* → *x*) on [0, 1], which is quite restrictive. However, we are imposing softening functions to be lesser than the identity on [1, +∞); whereas, Holland supposed *ψ* to be bounded and satisfy Equation (4). A concrete example of such a function *ψ*, lies in the piecewise polynomial function of Catoni and Giulini [21], defined by:

$$\psi(u) = \begin{cases} -2\sqrt{2}/3 & \text{if } u \le -\sqrt{2} \\ u - u^3/6 & \text{if } u \in [-2\sqrt{2}/3, 2\sqrt{2}/3] \\ 2\sqrt{2}/3 & \text{otherwise} \end{cases}$$

As in Section 4, we are considering the *ψ*-empirical risk *RS*,*ψ*,*t* for any *t* > 0. Holland provided his theorem given the fact the following assumptions are realised:

**Theorem 7.** *Ref. [13]. Let P be a prior distribution on model* H*. Let the three assumptions listed above hold. Setting t* 2 = *mM*2/(2 log(*δ*−1))*, then for any δ* ∈ [0; 1]*, with probability of at least* 1 − *δ over the random draw of the size-m sample S, simultaneously for all Q it holds that:*

$$\begin{split} \mathbb{E}\_{h \sim Q}[R(h)] &\leq \mathbb{E}\_{h \sim Q} \left[ R\_{S, \psi, t}(h) \right] + \frac{1}{\sqrt{m}} \left( \mathrm{KL}(Q || P) + \frac{1}{2} \log \left( \frac{8 \pi M\_2}{\delta^2} \right) - 1 \right) \\ &+ \frac{1}{\sqrt{m}} \nu^\*(\mathcal{H}) + O\left(\frac{1}{m}\right) \end{split}$$

*where:*

$$\nu^\*(\mathcal{H}) := \frac{\mathbb{E}\_{\hbar \sim P} \Big[ \exp \left( \sqrt{m} (R(\hbar) - R\_{S, \psi, t}(\hbar)) \right) \Big]}{\mathbb{E}\_{\hbar \sim P} \Big[ \exp \left( R(\hbar) - R\_{S, \psi, t}(\hbar) \right) \Big]}.$$

## **7. Proofs**

*7.1. Proof of Theorem 1*

**Proof.** Let *F* : R+ × R+ → R be a convex function, *P* a fixed prior, and *δ* ∈ [0, 1]. Since E*h*∼*P eF*(*RS*(*h*),*R*(*h*)) is a nonnegative random variable, we know that, by Markov's inequality, for any *h* ∈ H :

$$\mathbb{P}\left(\mathbb{E}\_{\mathrm{h}\sim P}\left[e^{F(\mathcal{R}\_{\mathrm{S}}(h),\mathcal{R}(h))}\right] > \frac{1}{\delta}\mathbb{E}\_{\mathrm{S}}\,\mathbb{E}\_{\mathrm{h}\sim P}\left[e^{F(\mathcal{R}\_{\mathrm{S}}(h),\mathcal{R}(h))}\right]\right) \leq \delta.$$

So with probability of at least 1 − *δ*, we have:

$$\mathbb{E}\_{\mathsf{h}\sim P}\left[e^{F\left(\mathcal{R}\_{\mathcal{S}}(h),\mathcal{R}(h)\right)}\right] \stackrel{\sim}{\leq} \frac{1}{\delta} \mathbb{E}\_{\mathcal{S}}\,\mathbb{E}\_{\mathsf{h}\sim P}\left[e^{F\left(\mathcal{R}\_{\mathcal{S}}(h),\mathcal{R}(h)\right)}\right] = \frac{\mathcal{X}}{\delta}.$$

Applying the log function on each side of this inequality gives us with probability of at least 1 − *δ* over samples *S*:

$$\log\left(\mathbb{E}\_{h\sim P}\left[e^{F(R\_S(h),R(h))}\right]\right) \le \log\left(\frac{\chi}{\delta}\right).$$

We now rename *A* := log E*h*∼*P eF*(*RS*(*h*),*R*(*h*)) .

Furthermore, if we denote by *dQ dP* the Radon-Nikodym derivative of *Q* with respect to *P* when *Q P*, we then have, for all *Q* such that *Q* ∼ *P*:

$$\begin{split} A &= \log \left( \mathbb{E}\_{h \sim Q} \left[ \frac{dP}{dQ} e^{F(R\_S(h), R(h))} \right] \right) \\ &= \log \left( \mathbb{E}\_{h \sim Q} \left[ \left( \frac{dQ}{dP} \right)^{-1} e^{F(R\_S(h), R(h))} \right] \right) \end{split} \tag{4}$$

and by concavity of log and Jensen's inequality,

$$\begin{aligned} &\geq -\mathbb{E}\_{h\sim Q} \Big[ \log \left( \frac{dQ}{dP} \right) \Big] + \mathbb{E}\_{h\sim Q} [F(R\_S(h), R(h))] \\ &= -\operatorname{KL}(Q||P) + \mathbb{E}\_{h\sim Q} [F(R\_S(h), R(h))] \end{aligned}$$

while by convexity of *F* with Jensen's inequality,

$$0 \geq -\mathsf{KL}(Q||P) + F(\mathbb{E}\_{h \sim Q}[\mathcal{R}\_S(h)], \mathbb{E}\_{h \sim Q}[R(h)]).$$

Hence, for *Q* such that *Q* ∼ *P*,

$$F\left(\mathbb{E}\_{h\sim Q}[R\_S(h)], \mathbb{E}\_{h\sim Q}[R(h)]\right) \le \text{KL}(Q||P) + A.$$

So with probability 1 − *δ*, for *Q* such that *Q* ∼ *P*,

$$F\left(\mathbb{E}\_{h\sim\mathcal{Q}}[R\_S(h)], \mathbb{E}\_{h\sim\mathcal{Q}}[R(h)]\right) \leq \mathsf{KL}(Q||P) + \log\left(\frac{\mathcal{X}}{\delta}\right).$$

This completes the proof of Theorem 1.

## *7.2. Proof of Theorem 5*

We first provide a technical property. Recall that:

$$\xi = \mathbb{E}\_{h \sim P} \left[ \exp \left( \frac{K(h)^2}{2m^{1-2\alpha}} \right) \right].$$

**Proposition 3.** *Let α* ∈ R*. Suppose the loss is* HYPE(*K*) *compliant with K*(*h*) = *B*||*h*|| + *C, with B* > 0*, C* ≥ 0*. Then, for any Gaussian prior P* = N (0, *σ*2*IN*) *with σ*2 = *t m*1−2*α B*2 *,* 0 < *t* < 1 *and N* ≥ 6 *we have:*

$$\xi \le 2 \exp\left(\frac{\mathbb{C}^2}{2m^{1-2\alpha}f(t)}(1+f(t))\right) \frac{1}{\left(\sqrt{1-t}\right)^N} \left(1+\left(\frac{\mathbb{C}}{\sqrt{2f(t)m^{1-2\alpha}}}\right)\right)^{N-1}$$

*with f*(*t*) = 1−*t t .*

**Proof.** We recall that *σ*2 = *t m*1−2*α B*2 . By expliciting the expectation and *K*(*h*) we thus obtain:

$$\begin{split} \mathcal{E}\_{\mathbb{S}} &= \left(\frac{1}{\sqrt{2\pi\sigma^{2}}}\right)^{N} \int\_{\mathbb{H}\in\mathbb{R}^{N}} \exp\left(\frac{(\mathcal{B}||h||+\mathcal{C})^{2}}{2m^{1-2\alpha}} - \frac{||h||^{2}\mathcal{B}^{2}}{2tm^{1-2\alpha}}\right) dh \\ &= \left(\frac{1}{\sqrt{2\pi\sigma^{2}}}\right)^{N} \int\_{\mathbb{H}\in\mathbb{R}^{N}} \exp\left(-\frac{1}{2m^{1-2\alpha}} \left(f(t)B^{2}||h||^{2} - 2BC||h|| - \mathcal{C}^{2}\right)\right) dh \\ &= \left(\frac{1}{\sqrt{2\pi\sigma^{2}}}\right)^{N} \int\_{\mathbb{H}\in\mathbb{R}^{N}} \exp\left(-\frac{B^{2}f(t)}{2m^{1-2\alpha}} \left(||h||^{2} - \frac{2C||h||}{Bf(t)} - \frac{\mathcal{C}^{2}}{B^{2}f(t)}\right)\right) dh \\ &= \exp\left(\frac{\mathcal{C}^{2}}{2m^{1-2\alpha}f(t)}(1+f(t))\right) \frac{1}{(\sqrt{2\pi\sigma^{2}})^{N}} \int\_{\mathbb{H}\in\mathbb{R}^{N}} \exp\left(-\frac{B^{2}f(t)}{2m^{1-2\alpha}} \left(||h|| - \frac{\mathcal{C}}{Bf(t)}\right)^{2}\right) dh. \end{split}$$

We will use the spherical coordinates in *N*-dimensional Euclidean space given in [22]:

$$
\varphi: (h\_1, \dots, h\_N) \to (r, \varphi\_1, \dots, \varphi\_{N-1}),
$$

where especially *r* = ||*h*|| and also the Jacobian of *φ* is given by:

$$d^N V = r^{N-1} \prod\_{k=1}^{N-2} \sin^k(\varphi\_{N-1-k}) = r^{N-1} d\_{S^{N-1}} V.$$

Let us also precise that as given in Blumenson [22] (page 66), we have that the surface of the sphere of radius 1 in *N*-dimensional space is:

$$\int\_{\mathcal{P}\_1,\dots,\mathcal{P}\_{N-1}} d\_{S^{N-1}} V \, d\boldsymbol{\varrho}\_1 \dots d\boldsymbol{\varrho}\_{N-1} = \frac{2\sqrt{\pi}^N}{\Gamma\left(\frac{N}{2}\right)}$$

where Γ is the Gamma function defined as:

$$
\Gamma(\mathbf{x}) = \int\_0^{+\infty} t^{x-1} e^{-t} dt \quad \text{for } \ge -1.
$$

Then, if we set:

$$A := \int\_{h \in \mathbb{R}^N} \exp\left(-\frac{B^2 f(t)}{2m^{1-2\alpha}} \left(||h|| - \frac{\mathbb{C}}{B f(t)}\right)^2\right) dh$$

we obtain by a change of variable:

$$\begin{split} A &= \int\_{r,\varphi\_{1},\ldots,\varphi\_{N-1}} \exp\left(-\frac{B^{2}f(t)}{2m^{1-2\alpha}} \left(r-\frac{\mathbb{C}}{Bf(t)}\right)^{2}\right) d^{N}V \, dr d\varphi\_{1} \ldots d\varphi\_{N-1} \\ &= \left(\frac{2\sqrt{\pi}^{N}}{\Gamma\left(\frac{N}{2}\right)}\right) \int\_{r=0}^{+\infty} \exp\left(-\frac{B^{2}f(t)}{2m^{1-2\alpha}} \left(r-\frac{\mathbb{C}}{Bf(t)}\right)^{2}\right) r^{N-1} dr \\ &= \left(\frac{2\sqrt{\pi}^{N}}{\Gamma\left(\frac{N}{2}\right)}\right) \int\_{r=-\frac{C}{2f(t)}}^{+\infty} \left(r+\frac{\mathbb{C}}{Bf(t)}\right)^{N-1} \exp\left(-\frac{B^{2}f(t)}{2m^{1-2\alpha}}r^{2}\right) dr \\ &= \left(\frac{2\sqrt{\pi}^{N}}{\Gamma\left(\frac{N}{2}\right)}\right) \sum\_{k=0}^{N-1} \binom{N-1}{k} \left(\frac{\mathbb{C}}{Bf(t)}\right)^{N-k-1} \int\_{r=-\frac{C}{2f(t)}}^{+\infty} r^{k} \exp\left(-\frac{B^{2}f(t)}{2m^{1-2\alpha}}r^{2}\right) dr. \end{split}$$

We fix a random variable *X* such that:

$$X \sim \mathcal{N}\left(0, \frac{m^{1-2a}}{B^2(f(t))}\right).$$

We then have for any *k* positive integer, if *k* is even:

$$\begin{split} \int\_{r=-\frac{C}{Bf(t)}}^{+\infty} r^k \exp\left(-\frac{B^2f(t)}{2m^{1-2\alpha}}r^2\right) dr &\leq \int\_{r=-\infty}^{+\infty} r^k \exp\left(-\frac{B^2f(t)}{2m^{1-2\alpha}}r^2\right) dr\\ &\leq \sqrt{2\pi \frac{m^{1-2\alpha}}{B^2f(t)}} \mathbb{E}[|X|^k]. \end{split}$$

And if *k* is odd:

$$\begin{split} \int\_{r=-\frac{C}{Bf(t)}}^{+\infty} r^k \exp\left(-\frac{B^2f(t)}{2m^{1-2a}}r^2\right) dr &\leq \int\_{r=0}^{+\infty} r^k \exp\left(-\frac{B^2f(t)}{2m^{1-2a}}r^2\right) dr\\ &\leq \sqrt{2\pi \frac{m^{1-2a}}{B^2f(t)}} \mathbb{E}[|X|^k \mathbb{1}(X \geq 0)] \\ &\leq \sqrt{2\pi \frac{m^{1-2a}}{B^2f(t)}} \mathbb{E}[|X|^k]. \end{split}$$

So we have:

$$A \le \left(\frac{2\sqrt{\pi}^N}{\Gamma\left(\frac{N}{2}\right)}\right) \sum\_{k=0}^{N-1} \binom{N-1}{k} \left(\frac{\mathcal{C}}{Bf(t)}\right)^{N-k-1} \sqrt{2\pi \frac{m^{1-2\alpha}}{B^2f(t)}} \mathbb{E}[|X|^k].$$

As precised in [23], we have for any *k*:

$$\mathbb{E}[|X|^k] = \left(\sqrt{\frac{m^{1-2\alpha}}{B^2 f(t)}}\right)^k 2^{k/2} \frac{\Gamma\left(\frac{k+1}{2}\right)}{\sqrt{\pi}}.$$

So finally:

$$A \le 2\sqrt{\pi}^N \sum\_{k=0}^{N-1} \binom{N-1}{k} \left(\frac{\mathbb{C}}{Bf(t)}\right)^{N-k-1} \left(\sqrt{\frac{2m^{1-2a}}{B^2f(t)}}\right)^{k+1} \frac{\Gamma\left(\frac{k+1}{2}\right)}{\Gamma\left(\frac{N}{2}\right)}.$$

**Lemma 3.** *If N* ≥ 6*, then:*

$$\max\_{k=0\ldots N-1} \frac{\Gamma\left(\frac{k+1}{2}\right)}{\Gamma\left(\frac{N}{2}\right)} = 1.$$

**Proof.** As precised in the introduction of Srinivasan and Zvengrowski [24], Gauss [25] (page 147) proved that on the interval [*x*0, +∞) where *x*0 ∈ [1.46, 1.47], Γ is a monotonic increasing function. So, for *N* − 1 ≥ *k* ≥ 2, Γ( *k*+1 2 ) ≤ Γ( *N* 2 ). And because Γ(1/2) = √*π*, Γ(1) = 1, we have:

$$\max\_{k=0\ldots N-1} \frac{\Gamma\left(\frac{k+1}{2}\right)}{\Gamma\left(\frac{N}{2}\right)} = \max\left(\frac{\sqrt{\pi}}{\Gamma\left(\frac{N}{2}\right)}, \frac{\Gamma\left(\frac{N-1+1}{2}\right)}{\Gamma\left(\frac{N}{2}\right)}\right) = \max\left(\frac{\sqrt{\pi}}{\Gamma\left(\frac{N}{2}\right)}, 1\right)$$

Because *N* ≥ 6, and Γ is monotone and increasing on [3; +∞], we have Γ(*N*/2) ≥ Γ(3) ≥ √*π*. Hence the result.

Using Lemma 3 allows us to write:

$$A \le 2\sqrt{\pi}^N \sum\_{k=0}^{N-1} \binom{N-1}{k} \left(\frac{\mathbb{C}}{Bf(t)}\right)^{N-k-1} \left(\sqrt{\frac{2m^{1-2\alpha}}{B^2f(t)}}\right)^{k+1}.$$

We recall that *σ*2 = *t m*1−2*α B*2 and *f*(*t*) = 1−*t t* . Then we can write:

$$A \le 2\sqrt{\pi}^N \sum\_{k=0}^{N-1} \binom{N-1}{k} \left(\frac{\mathbb{C}}{Bf(t)}\right)^{N-k-1} \left(\sqrt{\frac{2\sigma^2}{1-t}}\right)^{k+1}.$$

We now conclude with the final bound on *ξ*:

*ξ* ≤ exp *C*2 2*m*1−2*α f*(*t*)(1 + *f*(*t*)) 1 ( √ 2*πσ*2)*N A* ≤ exp *C*2 2*m*1−2*α f*(*t*)(1 + *f*(*t*)) 1 ( √ 2*πσ*2)*N* 2 √*πN* ∑*N*−1 *k*=0 ( *N*−1 *k* ) *C B f*(*t*) *N*−*k*−1 & 2*σ*2 1−*t k*+1 ≤ 2 exp *C*2 2*m*1−2*α f*(*t*)(1 + *f*(*t*)) ∑*N*−1 *k*=0 ( *N*−1 *k* ) *C B f*(*t*) *N*−*k*−1& 1 1−*t k*+1 & *B*2 2*tm*1−2*α N*−*k*−1 ≤ 2 exp *C*2 2*m*1−2*α f*(*t*)(1 + *f*(*t*)) ∑*N*−1 *k*=0 ( *N*−1 *k* ) *C* √*t* (1−*t*) √ 2*m*1−2*α N*−*k*−1& 1 1−*t k*+1 ≤ 2 exp *C*2 2*m*1−2*α f*(*t*) (1+*f*(*t*)) ( √1−*t*) *N* ∑*N*−1 *k*=0 ( *N*−1 *k* ) √ *C* 2 *f*(*t*)*m*1−2*α N*−*k*−1 ≤ 2 exp *C*2 2*m*1−2*α f*(*t*) (1+*f*(*t*)) ( √1−*t*) *N* 1 + √ *C* 2 *f*(*t*)*m*1−2*α N*−1 .

This completes the proof of Proposition 3.

**Proof of Theorem 5.** We combine Theorem 3 with Proposition 3. We also upper-bound *N* − 1 by *N*.

**Author Contributions:** Conceptualization, M.H., B.G. and J.S.-T.; Formal analysis, M.H., B.G. and O.R.; Project administration, B.G.; Supervision, B.G.; Writing—original draft, M.H., B.G. and O.R.; Writing—review and editing, M.H., B.G., O.R. and J.S.-T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported in part by the U.S. Army Research Laboratory and the U. S. Army Research Office, and by the U.K. Ministry of Defence and the U.K. Engineering and Physical Sciences Research Council (EPSRC) under grant number EP/R013616/1. BG acknowledges partial support from the French National Agency for Research, grants ANR-18-CE40-0016-01 and ANR-18-CE23- 0015-02.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**

## *Article* **Differentiable PAC–Bayes Objectives with Partially Aggregated Neural Networks**

**Felix Biggs 1 and Benjamin Guedj 1,2,\***

**\*** Correspondence<EMAIL_ADDRESS>
**Abstract:** We make two related contributions motivated by the challenge of training stochastic neural networks, particularly in a PAC–Bayesian setting: (1) we show how averaging over an ensemble of stochastic neural networks enables a new class of partially-aggregated estimators, proving that these lead to unbiased lower-variance output and gradient estimators; (2) we reformulate a PAC–Bayesian bound for signed-*output* networks to derive in combination with the above a directly optimisable, differentiable objective and a generalisation guarantee, without using a surrogate loss or loosening the bound. We show empirically that this leads to competitive generalisation guarantees and compares favourably to other methods for training such networks. Finally, we note that the above leads to a simpler PAC–Bayesian training scheme for sign-*activation* networks than previous work.

**Keywords:** statistical learning theory; PAC–Bayes theory; deep learning

## **1. Introduction**

**Citation:** Biggs, F.; Guedj, B. Differentiable PAC–Bayes Objectives with Partially Aggregated Neural Networks. *Entropy* **2021**, *23*, 1280. https://doi.org/10.3390/e23101280

Academic Editor: Udo Von Toussaint

Received: 22 August 2021 Accepted: 27 September 2021 Published: 29 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The use of stochastic neural networks has become widespread in the PAC–Bayesian and Bayesian deep learning [1] literature as a way to quantify predictive uncertainty and obtain generalisation bounds. PAC–Bayesian theorems generally bound the expected loss of *randomised* estimators, so it has proven easier to obtain non-vacuous numerical guarantees on generalisation in such networks.

However, we observe that when training these in the PAC–Bayesian setting, the objective used is generally somewhat divorced from the bound on misclassification loss itself, often because non-differentiability leads to difficulties with direct optimisation. For example, Langford and Caruana [2], Zhou et al. [3], and Dziugaite and Roy [4] all initially train non-stochastic networks before using them as the mode of a distribution, with variance chosen, respectively, through a computationally-expensive sensitivity analysis, as a proportion of weight norms, or by optimising an objective with both a surrogate loss function and a different dependence on the Kullback–Leibler (KL) divergence from their bound.

In exploring methods to circumvent this gap, we also note that PAC–Bayesian bounds can often be straightforwardly adapted to aggregates or averages of estimators, leading directly to analytic and differentiable objective functions (for example, [5]). Unfortunately, averages over deep stochastic networks are usually intractable or, if possible, very costly (as found by [6]).

Motivated by these observations, our main contribution is to obtain a compromise by defining new and general "*partially-aggregated*" Monte Carlo estimators for the average output and gradients of deep stochastic networks (Section 3), with the direct optimisation of PAC–Bayesian bounds in mind. Although our main focus here is on the use of this estimator in a PAC–Bayesian application, we emphasise that the technique applies generally to stochastic networks and thus has links to other variance-reduction techniques for training them, such as the pathwise estimator used in the context of neural networks by [7] amongst

many others or Flipout [8]; indeed, it can be used in combination with these techniques. We provide proofs (Section 4) that this application leads to lower variances than a Monte Carlo forward pass and lower variance final-layer gradients than REINFORCE [9].

A further contribution of ours is a first application of this general estimator to nondifferentiable "signed-output" networks (with a final output ∈ {−1, +1} and arbitrarily complex other structure, see Section 4). As well as reducing variances as stated above, a small amount of additional structure in combination with partial-aggregation enables us to extend the pathwise estimator to the other layers, which usually requires a fully differentiable network and eases training by reducing the variance of gradient estimates.

We adapt a binary classification bound (Section 5) from Catoni [10] to these networks, yielding straightforward and directly differentiable objectives when used in combination with aggregation. Closing this gap between objectives and bounds leads to improved theoretical properties.

Further, since most of the existing PAC–Bayes bounds for neural networks have a heavy dependency on the distance from initialisation of the obtained solution, we would intuitively expect these lower variances to lead to faster convergence and tighter bounds (from finding low-error solutions nearer to the initialisation). We indeed observe this experimentally, showing that training PAC–Bayesian objectives in combination with partial aggregation leads to competitive experimental generalisation guarantees (Section 6), and improves upon naive Monte Carlo and REINFORCE.

As a useful corollary, this application also leads us to a similar but simpler PAC– Bayesian training method for sign-activation neural networks than Letarte et al. [6], which successfully aggregated networks with all sign activation functions ∈ {+1, −1} and a non-standard tree structure, but incurred an exponential KL divergence penalty and a heavy computational cost (so that in practice they often resorted to a Monte Carlo estimate). Further, the lower variance of our obtained estimator predictions enables us to use the Gibbs estimator directly (where we draw a single sample function for every new example), leading to a modified bound on the misclassification loss which is a factor of two tighter without a significant performance penalty.

We discuss further and outline future work in Section 7.

#### **2. Background**

We begin here by setting out our notation and the requisite background.

Generally, we consider parameterised functions, { *fθ* : X → Y|*θ* ∈ Θ ⊂ R*N*}, in a specific form, choosing X ⊂ R*d*0 and an arbitrary output space Y which could be for example {−1, +1} or R. We wish to find functions minimizing the out-of-sample risk, *R*(*f*) = E(*x*,*y*)∼D-(*f*(*x*), *y*), for some loss function -, for example the 0-1 misclassification loss for classification, -0−1(*y*, *y* ) = **1**{*y* -= *y* }, or the binary linear loss, lin(*y*, *y* ) = 1 2 (1 − *yy* ), with Y = {+1, −1}. These must be chosen based on an i.i.d. sample *S* = {(*xi*, *yi*)}*m i*=1 ∼ D*m* from the data distribution D, using the surrogate of in-sample empirical risk, *RS*(*f*) = 1 *m* ∑*m i*=1 -(*f*(*xi*), *yi*). We denote the expected and empirical risks under the misclassification and linear losses, respectively *R*0−1, *R*lin, *R*0−1 *S* and *R*lin *S* .

In this paper, we consider learning a distribution (PAC–Bayesian posterior), *Q*, over the parameters *θ*. PAC–Bayesian theorems then provide bounds on the expected generalization risk of randomised classifiers, where every prediction is made using a newly sampled function from our posterior, *fθ*, *θ* ∼ *Q*.

We also consider averaging the above to obtain *Q*-aggregated prediction functions,

$$F\_Q(\mathfrak{x}) := \mathbb{E}\_{\theta \sim Q} f\_\theta(\mathfrak{x}). \tag{1}$$

In the case of a convex loss function, Jensen's inequality lower bounds the risk of the randomised function by its *Q*-aggregate: -(*FQ*(*x*), *y*) ≤ E*f*∼*Q*-(*f*(*x*), *y*). The equality is achieved by the linear loss, a fact we will exploit to obtain an easier PAC–Bayesian optimisation objective in Section 5.

#### *2.1. Analytic Q-Aggregates for Signed Linear Functions*

*Q*-aggregate predictors are analytically tractable for "signed-output" functions (here the sign function and "signed" functions have outputs ∈ {+1, −1}, as the terminology "binary", used sometimes in the literature, suggests to us too strongly an output ∈ {0, 1}) of the form *fw*(*x*) = sign(*w* · *x*) under a normal distribution, *Q*(*w*) = *N*(*μ*,I), as specifically considered in a PAC–Bayesian context for binary classification by [5], obtaining an differentiable objective similar to the SVM. Provided *x* -= **0**:

$$F\_Q(\mathbf{x}) := \mathbb{E}\_{\mathbf{w} \sim N(\mu, \mathbb{I})} \operatorname{sign}(\mathbf{w} \cdot \mathbf{x}) = \operatorname{erf}\left(\frac{\mu \cdot \mathbf{x}}{\sqrt{2} ||\mathbf{x}||}\right). \tag{2}$$

In Section 4, we will consider aggregating signed output (*f*(*x*) ∈ {+1, −1}) functions of a more general form.

#### *2.2. Monte Carlo Estimators for More Complex Q-Aggregates*

The framework of *Q*-aggregates can be extended to less tractable cases (for example, with *fθ* a randomised or a "Bayesian" neural network, see, e.g., [1]) through a simple and unbiased Monte Carlo approximation:

$$F\_Q(\mathbf{x}) = \mathbb{E}\_{\theta \sim Q} f\_\theta(\mathbf{x}) \approx \frac{1}{T} \sum\_{t=1}^T f\_{\theta^t}(\mathbf{x}) := \mathbb{P}\_Q(\mathbf{x}).\tag{3}$$

If we go on to parameterize our posterior *Q* by *φ* ∈ Φ ⊂ R*N* as *Qφ* and wish to obtain gradients without a closed form for *FQφ* (*x*) = E*θ*∼*Qφ fθ* (*x*), there are two possibilities. One is REINFORCE [9], which requires only a differentiable density function, *qφ*(*θ*) and makes a Monte Carlo approximation to the left hand side of the identity ∇*φ*E*θ*∼*qφ fθ* (*x*) = E*θ*∼*qφ* [ *fθ* (*x*)∇*φ* log *qφ*(*θ*)].

The other is the pathwise estimator, which additionally requires that *fθ* (*x*) be differentiable w.r.t. *θ*, and that the probability distribution chosen has a standardization function, *Sφ*, which removes the *φ* dependence, turning a parameterised *qφ* into a nonparameterised distribution *p*: for example, *Sμ*,*σ*(*X*)=(*X* − *μ*)/*σ* to transform a general normal distribution into a standard normal. If this exists, the right hand side of ∇*φ*E*θ*∼*qφ fθ* (*x*) = E∼*p*∇*φ f S*−1 *φ* () (*x*) generally yields lower-variance estimates than REIN-FORCE (see for a modern survey [11]).

The variance introduced by REINFORCE can make it difficult to train neural networks when the pathwise estimator is not available, for example when non-differentiable activation functions are used. Below we find a compromise between the analytically closed form of (2) and the above estimator that enables us to make differentiable certain classes of network and extend the pathwise estimator where otherwise it could not be used. Through this we are able to stably train a new class of network.

## *2.3. PAC–Bayesian Approach*

We use PAC–Bayes in this paper to obtain generalisation guarantees and theoreticallymotivated training methods. The primary bound utilised is based on the following theorem, valid for a loss taking values in [0, 1]:

**Theorem 1** ([10], Theorem 1.2.6)**.** *Given probability measure P on hypothesis space* F *and α* > 1*, for all Q on* F *with probability at least* 1 − *δ over S* ∼ D*m,*

$$\mathbb{E}\_{f \sim Q} \mathbb{R}(f) \le \inf\_{\lambda > 1} \Phi\_{\lambda/m}^{-1} \left[ \mathbb{E}\_{f \sim Q} R\_S(f) + \frac{a}{\lambda} \Delta \right]$$

$$with \ \Phi\_{\gamma}^{-1}(t) = \frac{1 - \exp(-\gamma t)}{1 - \exp(-\gamma)} \, and \ \Delta = \mathrm{KL}(Q|P) - \log \delta + 2 \log \left( \frac{\log a^2 \lambda}{\log a} \right).$$

This slightly opaque formulation (used previously by [3]) gives essentially identical results when KL/*m* is large to the better-known "small-kl" PAC–Bayes bounds originated by Langford and Seeger [12], Seeger et al. [13]. It is chosen because it leads to objectives that are *linear* in the empirical loss and KL divergence, like

$$\mathbb{E}\_{f \sim Q} R\_S(f) + \frac{\mathbb{KL}(Q|P)}{\lambda}. \tag{4}$$

This objective is minimised by a Gibbs distribution and is closely related to the evidence lower bound (ELBO) usually optimised by Bayesian Neural Networks [1]. Such a connection has been noted throughout the PAC–Bayesian literature; we refer the reader to [14] or [15] for a formalised treatment.

#### **3. The Partial Aggregation Estimator**

Here we outline our main contribution: a reformulation of *Q*-aggregation for neural networks leading to different, lower-variance, Monte Carlo estimators for their outputs and gradients. These estimators apply to networks with a dense final layer, and arbitrary stochastic other structure (for example convolutions, residual layers or a non-feedforward structure). Specifically, they take the form

$$f\_{\theta}(\mathbf{x}) = A(\mathbf{z}\mathbf{w}\cdot\boldsymbol{\eta}\_{\theta^{\neg \mathbf{w}}}(\mathbf{x})) \tag{5}$$

with *θ* = vec(*w*, *θ¬w*) ∈ Θ ⊂ R*D*, *w* ∈ R*d*, and *θ¬w* ∈ Θ*¬w* ⊂ R*D*−*d* the parameter set excluding *w*, for the non-final layers of the network. These non-final layers are included in *ηθ¬w* : X→A*d* ⊆ R*d* and the final activation is *A* : R → Y. For simplicity we have used a one-dimensional output but we note that the formulation and below derivations trivially extend to a vector-valued function with elementwise activations. We require the distribution over parameters to factorise like *Q*(*θ*) = *Qw*(*w*)*Q¬w*(*θ¬w*), which is consistent with the literature.

We recover a similar functional form to that considered in Section 2.1 by rewriting the function as *A*(*w* · *a*) with *a* ∈ A*d* the randomised hidden-layer activations. The "aggregated" activation function on the final layer, which we define as *I*(*a*) := - *A*(*w* · *a*) d*Qw*(*w*), may then be analytically tractable. For example, with *w* ∼ *N*(*μ*,I) and a sign final activation, we recall (2) where *I*(*a*) = erf *μ*·*a* √2 *a* .

Using these definitions we can write the *Q*-aggregate in terms of the conditional distribution on the activations, *a*, which takes the form *Q*˜ *¬w*(*a*|*x*) := (*η*(·)(*x*)) ◦ *Q¬w*, (i.e., the distribution of *ηθ¬w* (*x*)|*x*, with *θ¬w* ∼ *Q¬w*). The *Q*-aggregate can then be stated as

$$\begin{split} F\_{\mathbb{Q}}(\mathbf{x}) &:= \mathbb{E}\_{\theta \sim \mathcal{Q}}[f\_{\theta}(\mathbf{x})] \\ &= \int\_{\theta^{\neg w}} \left[ \int\_{\mathbb{R}^{d}} A(w \cdot \eta\_{\theta^{\neg w}}(\mathbf{x})) \, \mathrm{d}Q^{w}(w) \right] \mathrm{d}Q^{\neg w}(\theta^{\neg w}) \\ &= \int\_{\theta^{\neg w}} I(\eta\_{\theta^{\neg w}}(\mathbf{x})) \, \mathrm{d}Q^{\neg w}(\theta^{\neg w}) \\ &= \int\_{\mathcal{A}^{d}} I(\mathbf{a}) \, \mathrm{d}\{ (\eta\_{(\cdot)}(\mathbf{x})) \diamond Q^{\neg w} \}(\mathbf{a}) \\ &=: \int\_{\mathcal{A}^{d}} I(\mathbf{a}) \, \mathrm{d}Q^{\neg w}(\mathbf{a}|\mathbf{x}). \end{split}$$

In most cases, the final integral cannot be calculated exactly or involves a large summation, so we resort to a Monte Carlo estimate, for each *x* drawing *T* samples of the randomised activations, {*at* }*T t*=1 ∼ *Q*˜ *¬w*(*a*|*x*) to obtain the "partially-aggregated" estimator

$$F\_Q(\mathbf{x}) = \int\_{\mathcal{A}^d} I(\mathbf{a}) \, \mathbf{d} \tilde{Q}^{\sim \text{av}}(\mathbf{a}|\mathbf{x}) \approx \frac{1}{T} \sum\_{t=1}^T I(\mathbf{a}^t) = \mathbb{P}\_Q^\*(\mathbf{x}). \tag{6}$$

This is quite similar to the original estimator from (3), but in fact the aggregation of the final layer may significantly reduce the variance of the outputs and also make better gradient estimates possible, as we will show below.

#### *3.1. Reduced Variance Estimates*

**Proposition 1.** Lower variance outputs: *For a neural network as defined by Equation* (5) *and the unbiased Q-aggregation estimators defined by Equations* (3) *and* (6)*,*

$$\mathbb{V}\_Q[\mathcal{F}\_Q^\*(\mathfrak{x})] \le \mathbb{V}\_Q[\mathcal{F}\_Q(\mathfrak{x})].$$

**Proof.** Treating *a* as a random variable, always conditioned on *x*, we have

$$\begin{split} & \mathbb{V}\_{Q}[\hat{\mathbb{P}}\_{Q}(\mathbf{x})] - \mathbb{V}\_{Q}[\hat{\mathbb{P}}\_{Q}^{\star}(\mathbf{x})] = \mathbb{E}\_{Q} |\hat{\mathbb{P}}\_{Q}(\mathbf{x})|^{2} - \mathbb{E}\_{Q} |\hat{\mathbb{P}}\_{Q}^{\star}(\mathbf{x})|^{2} \\ &= \frac{1}{T} \mathbb{E}\_{\mathfrak{a}|\mathbf{x}} \left[ \mathbb{E}\_{\mathfrak{w}} |A(\mathfrak{w} \cdot \mathfrak{a})|^{2} - |\mathbb{E}\_{\mathfrak{w}} A(\mathfrak{w} \cdot \mathfrak{a})|^{2} \right] \\ &= \frac{1}{T} \mathbb{E}\_{\mathfrak{a}|\mathbf{x}} [\mathbb{V}\_{\mathfrak{w}}[A(\mathfrak{w} \cdot \mathfrak{a})]] \ge 0. \end{split}$$

$$\mathbb{0}$$

From the above we see that the aggregate outputs estimated through partial-aggregation have lower variances. Next, we consider the two unbiased gradient estimators for the distribution over final-layer weights, *w*, arising from partial-aggregation or REINFORCE (as would be used, for example, where the final layer is non-differentiable). Assuming *Qw* has a density, *qφ*(*θw*), parameterised by *φ*, these use forward samples of {*wt* , *θ¬w*,(*t*)}*T t*=1 as:

$$\begin{aligned} \hat{\mathbf{G}}(\boldsymbol{\mathfrak{x}}) &:= \frac{1}{T} \sum\_{t=1}^{T} A \left( \boldsymbol{\mathfrak{w}}^{t} \cdot \boldsymbol{\mathfrak{y}}^{t} \right) \nabla\_{\boldsymbol{\Phi}} \log q\_{\boldsymbol{\Phi}}(\boldsymbol{\mathfrak{w}}^{t}), \\ \hat{\mathbf{G}}^{\*}(\boldsymbol{\mathfrak{x}}) &:= \frac{1}{T} \sum\_{t=1}^{T} \nabla\_{\boldsymbol{\Phi}} I\_{\boldsymbol{\Phi}\boldsymbol{\Phi}}(\boldsymbol{\mathfrak{y}}^{t}). \end{aligned}$$

**Proposition 2.** Lower variance gradients: *Under the conditions of Proposition 1 and the above definitions,*

$$\operatorname{Cov}\_{\mathbb{Q}}[\hat{\mathbf{G}}^\*(\mathfrak{x})] \preceq \operatorname{Cov}\_{\mathbb{Q}}[\hat{\mathbf{G}}(\mathfrak{x})],$$

*where A* " *B* ⇐⇒ *B* − *A is positive semi-definite. Thus, for all u* -= **0***,* V[*G*ˆ ∗(*x*) · *u*] ≤ V[*G*ˆ(*x*) · *u*]*.*

**Proof.** Writing *v* := ∇*φ* log *qφ*(*w*) and using the unbiasedness of the estimators,

$$\begin{split} & \text{Cov}\_{Q}[\hat{\mathbf{G}}\_{Q}(\mathbf{x})] - \text{Cov}\_{Q}[\hat{\mathbf{G}}\_{Q}^{\star}(\mathbf{x})] \\ &= \mathbb{E}\_{Q}[\hat{\mathbf{G}}\_{Q}(\mathbf{x})\hat{\mathbf{G}}\_{Q}(\mathbf{x})^{T}] - \mathbb{E}\_{Q}[\hat{\mathbf{G}}\_{Q}^{\star}(\mathbf{x})\hat{\mathbf{G}}\_{Q}^{\star}(\mathbf{x})^{T}] \\ &= \frac{1}{T} \mathbb{E}\_{\mathbf{a}|\mathbf{x}} \left[ \mathbb{E}\_{\mathbf{w}}[A(\mathbf{w}\cdot\mathbf{a})^{2}\mathbf{w}^{\star}] - \nabla\_{\Phi}I\_{\eta\_{\phi}}(\boldsymbol{\eta}^{t}) \left(\nabla\_{\Phi}I\_{\eta\_{\phi}}(\boldsymbol{\eta}^{t})^{T}\right) \right] \\ &= \frac{1}{T} \mathbb{E}\_{\mathbf{a}|\mathbf{x}} \left[ \text{Cov}\_{\mathbf{w}}[A(\mathbf{w}\cdot\mathbf{a})\nabla\_{\Phi}\log q\_{\phi}(\mathbf{w})] \right] \succeq 0 \end{split}$$

where in the final line we have used that ∇*φ Iqφ* (*ηt* ) = ∇*φ*E*w*[*A*(*w* · *ηt* )] = E*w*[*A*(*w* · *ηt* )*v*].

## *3.2. Single Hidden Layer*

For clarity (and to introduce notation to be used in Section 4.2) we will briefly consider the case of a neural network with one hidden layer, *fθ* (*x*) = *A*2(*w*2 · *A*1(*W*1*x*)). The randomised parameters are *θ* = vec(*w*2, *W*1), *W*1 ∈ R*d*1×*d*0 , *w*2 ∈ R*d*1 and the elementwise activations are *A*1 : R*d*1 → A*d*1 1 ⊆ R*d*1 and *A*2 : R → Y. We choose the distribution *Q*(*θ*) =: *Q*2(*w*2)*Q*1(*W*1) to factorise over the layers. This is identical to the above and sets *ηW*1 (*x*) = *A*1(*W*1*x*).

Sampling *a* is straightforward if sampling *W*1 is. Further, if the final layer aggregate is differentiable, and so is the hidden layer activation *A*1, we may be able to use the lowervariance pathwise gradient estimator for gradients with respect to *Q*1. We note that this may be possible even if the activation *A*2 is not differentiable, as in Section 4, where we extend the pathwise estimator where we could not otherwise use it.

Computationally, we may implement the above by analytically finding the distribution on the "pre-activations" *W*1*x* (trivial for the normal distribution) before sampling this and passing through the activation. With the pathwise estimator this is known as the "local reparameterization trick" [16], which can lead to considerable computational savings on parallel minibatches compared to direct hierarchical sampling, *a* = *A*1(*W*1*x*) with *W*1 ∼ *Q*1. We will utilise this in all our reparameterizable dense networks, and a variation on it to save computation when using REINFORCE in Sections 4.2 and 6.

#### **4. Aggregating Signed-Output Networks**

Here we consider a first practical application of the aggregation estimator to stochastic neural networks with a final dense sign-activated layer. We have seen above that this partial aggregation leads to better-behaved training objectives and lower-variance gradient estimates across arbitrary other network structure, It may also allow use of pathwise gradients for the other layers, which would not be possible otherwise due to the nondifferentiability of the final layer.

Specifically, these networks take the form of Equation (5) with the final layer activation a sign function and weights drawn from a unit variance normal distribution, *Qw*(*w*) = N (*μ*,I). The aggregate *I*(*a*) is given by Equation (2). Normally-distributed weights are chosen because of the simple analytic forms for the aggregate (reminiscent of the tanh activation occasionally used for neural networks) and KL divergence (effectively an *L*2 regularisation penalty); we note however that closed forms are available for other commonlyused distributions such as the Laplace.

Using Equations (3) and (6) with independent samples {(*wt* , *θ¬w*,(*t*))}*T t*=1 ∼ *Q* and *ηt* := *ηθ¬w*,(*t*)(*x*) leads to the two unbiased estimators for the output (henceforth assuming the technical condition P*η*|*x*{*η* = **0**} = 0 that allows aggregation to be well-defined).

$$\mathcal{F}\_Q(\mathbf{x}) := \frac{1}{T} \sum\_{t=1}^T \text{sign}(\mathbf{w}^t \cdot \boldsymbol{\eta}^t) \tag{7}$$

$$\mathcal{F}\_Q^\*(\mathbf{x}) := \frac{1}{T} \sum\_{t=1}^T \text{erf}\left(\frac{\boldsymbol{\mu} \cdot \boldsymbol{\eta}^t}{\sqrt{2}||\boldsymbol{\eta}^t||}\right). \tag{8}$$

It follows immediately from Propositions 1 and 2 that the latter and the associated gradient estimators have lower variances than the former or the REINFORCE gradient estimates (which we would otherwise have to use due to the non-differentiability of the final layer).

#### *4.1. Lower Variance Estimates of Aggregated Sign-Output Networks*

We clarify the situation with the lower variance estimates further below. In particular, we find that the reduction in variance from using the partial-aggregation estimator is controlled by the norm *μ* , so that for small *μ* (as could be expected early in training) the difference can be large, while as *μ* grows, the difference in variance is controlled and we could reasonably revert to the Monte Carlo (or Gibbs) estimator. Note also that as *FQ*(*x*) → ±1 (as expected after training), both variances disappear.

We also show that a stricter condition than Proposition 2 holds on the variances of the aggregated gradients here, and thus that the non-aggregated gradients are noisier in all cases than the aggregate.

**Proposition 3.** *With the definitions given by Equation* (7)*, for all x* ∈ R*d*0 *, T* ∈ N*, and Q with normally-distributed final layer,*

$$0 \le \mathbb{V}\_{\mathcal{Q}}[\hat{F}\_{\mathcal{Q}}(\mathfrak{x})] - \mathbb{V}\_{\mathcal{Q}}[\hat{F}\_{\mathcal{Q}}^\*(\mathfrak{x})] \le \frac{1}{T} \left(1 - \left| \operatorname{erf} \left( \frac{||\mu||}{\sqrt{2}} \right) \right|^2 \right).$$

**Proof.** The left identity follows directly from Proposition 1. We also have

$$\begin{split} \mathbb{V}\_{\mathbb{Q}}[\mathbb{f}\_{\mathbb{Q}}(\mathbf{x})] - \mathbb{V}\_{\mathbb{Q}}[\mathbb{f}\_{\mathbb{Q}}^{\*}(\mathbf{x})] &= \frac{1}{T} \mathbb{E}\_{\mathbf{a}|\mathbf{x}}[\mathbb{V}\_{\mathbf{w}}[\text{sign}(\mathbf{w} \cdot \mathbf{a})]] \\ &= \frac{1}{T} \left( 1 - \mathbb{E}\_{\mathbf{a}|\mathbf{x}} \left| \text{erf} \left( \frac{\boldsymbol{\mu} \cdot \boldsymbol{\eta}}{\sqrt{2} ||\boldsymbol{\eta}||} \right) \right|^{2} \right) \\ &\leq \frac{1}{T} \left( 1 - \left| \text{erf} \left( \frac{||\boldsymbol{\mu}||}{\sqrt{2}} \right) \right|^{2} \right) . \end{split}$$

**Proposition 4.** *Under the conditions of Proposition 3,*

$$\text{Cov}[\hat{\mathbf{G}}^\*(\mathbf{x})] \preceq \text{Cov}[\hat{\mathbf{G}}(\mathbf{x})] + \frac{1 - 2/\pi}{T} \mathbb{I}.$$

*Thus, for all u with u* = 1*,*

$$\mathbb{V}[\mathbf{\hat{G}}^\*(\mathbf{x}) \cdot \boldsymbol{\mu}] \le \mathbb{V}[\mathbf{\hat{G}}(\mathbf{x}) \cdot \boldsymbol{\mu}] + \frac{1 - 2/\pi}{T} \cdot \boldsymbol{\mu}$$

**Proof.** It is straightforward to show that

$$\begin{aligned} \boldsymbol{\hat{G}}(\mathbf{x}) &:= \frac{1}{T} \sum\_{t=1}^{T} \text{sign}(\boldsymbol{w}^{t} \cdot \boldsymbol{\eta}^{t}) (\boldsymbol{\mu} - \boldsymbol{w}^{t}) \\ \boldsymbol{\hat{G}}^{\*}(\mathbf{x}) &:= \frac{1}{T} \sum\_{t=1}^{T} \frac{\boldsymbol{\eta}^{t}}{||\boldsymbol{\eta}^{t}||} \sqrt{\frac{2}{\pi}} \exp\left[ -\frac{1}{2} \left( \frac{\boldsymbol{\mu} \cdot \boldsymbol{\eta}^{t}}{||\boldsymbol{\eta}^{t}||} \right)^{2} \right] \\ \boldsymbol{\text{Cov}}[\boldsymbol{\hat{G}}(\mathbf{x})] &= \frac{1}{T} \left( \mathbb{I} - \boldsymbol{G} \boldsymbol{G}^{T} \right) \\ \boldsymbol{\text{Cov}}[\boldsymbol{\hat{G}}^{\*}(\mathbf{x})] &= \frac{1}{T} \left( \mathbb{E} \left[ \frac{\boldsymbol{\eta} \boldsymbol{\eta}^{T}}{||\boldsymbol{\eta}||^{2}} \frac{2}{\pi} \boldsymbol{\varepsilon}^{-\left(\frac{\boldsymbol{\eta} \cdot \boldsymbol{\eta}}{\|\boldsymbol{\eta}\|}\right)^{2}} \right] - \boldsymbol{G} \boldsymbol{G}^{T} \right) \end{aligned}$$

so for *u* -= **0**,

$$\begin{aligned} &\operatorname{Tu}^{\operatorname{T}}\left(\operatorname{Cov}[\hat{\mathbf{G}}(\mathbf{x})]-\operatorname{Cov}[\hat{\mathbf{G}}^{\*}(\mathbf{x})]\right)\mathfrak{u} \\ &=||\mathfrak{u}||^{2}-\frac{2}{\pi t}\mathbb{E}\left[\frac{|\mathbf{u}\cdot\boldsymbol{\eta}|^{2}}{||\boldsymbol{\eta}||^{2}}e^{-\left(\frac{\mathbf{p}\cdot\boldsymbol{\eta}}{\|\boldsymbol{\eta}\|}\right)^{2}}\right] \geq||\mathfrak{u}||^{2}\left(1-\frac{2}{\pi}\right)>0. \end{aligned}$$

Above we have brought *u* inside the term with an expectation, which is then bounded using Cauchy–Schwarz on |*u* · *η*|/ *η* ≤ *u* , and *e*−|*t*| ≤ 1 for all *t* ∈ R.

#### *4.2. All Sign Activations*

Here we examine an important special case previously examined by Letarte et al. [6]: a feed-forward network with all sign activations and normal weights. This takes the form

$$f\_{\theta}(\mathbf{x}) = \text{sign}(\mathbf{w}\_{L} \cdot \text{sign}(W\_{L-1} \dots \text{sign}(W\_1 \mathbf{x}) \dots \mathbf{x})),$$

with *θ* := vec(*wL*, ... , *W*1) and *Wl* := [*wl*,1 ... *wl*,*dl* ] *T*; *l* ∈ {1, ..., *L*} are the layer indices. We choose unit-variance normal distributions on the weights, which factorise into *Ql*(*Wl*) = ∏*dl i*=1 *ql*,*i*(*wl*,*i*) with *ql*,*i* = N (*μl*,*i*, I*dl*−1 ). In the notation of Section 3, *ηθ¬w* (*x*) = sign(*WL*−1 . . . sign(*W*1*x*)...) is the final layer activation, which could easily be sampled by mapping *x* through the first *L* − 1 layers with draws from the weight distribution.

Instead, we go on to make an iterative replacement of the weight distributions on each layer by conditionals on the layer activations to obtain the summation

$$F\_Q(\mathbf{x}) = \sum\_{\substack{\mathbf{a}\_1 \in \{+1, -1\}^{d\_1}}} \tilde{Q}\_1(\mathbf{a}\_1|\mathbf{x}) \quad \times \quad \dots \tag{9}$$

$$\dots \quad \times \sum\_{\substack{\mathbf{a}\_{L-1} \in \{+1, -1\}^{d\_{L-1}}}} \tilde{Q}\_{L-1}(\mathbf{a}\_{L-1}|\mathbf{a}\_{L-2}) \text{ erf}\left(\frac{\mu\_L \cdot \mathbf{a}\_{L-1}}{\sqrt{2}||\mathbf{a}\_{L-1}||}\right). \tag{9}$$

The number of terms is exponential in the depth so we instead hierarchically sample the *al*. Like local reparameterisation, this leads to a considerable computational saving over sampling a separate weight matrix for every input. The conditionals can be found in closed form: we can factorise individual hidden units *Q*˜*l*(*al*|*al*−1) := ∏*dl i*=1 *q*˜*l*,*i*(*al*,*i*|*al*−1), and find their activation distributions (with *a*0 := *x* and *z* a dummy variable):

$$\begin{split} \tilde{q}\_{l,i}(a\_{l,i} = \pm 1 \mid a\_{l-1}) &= \int\_{0}^{\infty} \mathcal{N}\left(z; \pm \mu\_{l,i} \cdot a\_{l-1}, \lVert a\_{l-1} \rVert^{2}\right) dz \\ &= \frac{1}{2} \left[ 1 \pm \text{erf}\left(\frac{\mu\_{l,i} \cdot a\_{l-1}}{\sqrt{2} \lVert a\_{l-1} \rVert}\right) \right] .\end{split}$$

A marginalised REINFORCE-style gradient estimator for *conditional* distributions can then be used; this does not necessarily have better statistical properties but in combination with the above is much more computationally efficient. This idea of "conditional sampling" is inspired by the local reinforce trick. Using samples {(*at* 1 ... *at L*−1)}*T t*=1 ∼ *Q*˜,

$$\frac{\partial F\_Q(\mathbf{x})}{\partial \mu\_{l,i}} \approx \frac{1}{T} \sum\_{t=1}^T \text{erf}\left(\frac{\mu\_L \cdot a\_{L-1}^t}{\sqrt{2}||\mathbf{a}\_{L-1}^t||}\right) \frac{\partial}{\partial \mu\_{l,i}} \log \overline{q}\_{l,i}(a\_{l,i}^t|\mathbf{a}\_{l-1}^t). \tag{10}$$

This formulation along with Equation (9) resembles the PBGNet model of [6], but derived in a very different way. Indeed both are equivalent in the single-hidden-layer case, but with more layers PBGNet uses an unusual tree-structured network to make the individual activations independent and avoid an exponential computational dependency on the depth in Equation (9). This makes the above summation exactly calculable but is also still not efficient enough in practice, so they resort further to a Monte Carlo approximation: informally, this draws new samples for every layer *l* based on an average of those from the previous layer, *al*|{*a*(*t*) *l*−1}*T t*=1 ∼ 1 *T* ∑*T t*=1 *Q*˜(*al*|*a*(*t*) *l*−1).

This is all justified within the tree-structured framework but leads to an exponential KL penalty which—as hinted by Letarte et al. [6] and shown empirically in Section 6 makes PAC–Bayes bound optimisation strongly favour shallower such networks. Our formulation avoids this, is more general—applying to alternative network structures—and we believe it is significantly easier to understand and use in practice.

#### **5. PAC–Bayesian Objectives with Signed-Outputs**

We now move to obtain binary classifiers with guarantees for the expected misclassification error, *R*0−1, which we do by optimizing PAC–Bayesian bounds. Such bounds (as in Theorem 1) will usually involve the non-differentiable and non-convex misclassification loss -0−1. However, to train a neural network we need to replace this by a differentiable surrogate, as discussed in the introduction.

Here we adopt a different approach by using our signed-output networks, where since *f*(*x*) ∈ {+1, −1}, there is an exact equivalence between the linear and misclassification losses, -0−1(*f*(*x*), *y*) = lin(*f*(*x*), *y*), avoiding an extra factor of two from the inequality -0−1 ≤ 2lin.

Although we have only moved the non-differentiability into *f*, the form of a PAC–Bayesian bound and the linearity of the loss and expectation allow us to go further and aggregate,

$$\mathbb{E}\_{f \sim Q} \ell\_{0-1}(f(\mathbf{x}), y) = \mathbb{E}\_{f \sim Q} \ell\_{\text{lin}}(f(\mathbf{x}), y) = \ell\_{\text{lin}}(F\_{\mathbb{Q}}(\mathbf{x}), y) \tag{11}$$

which allows us to use the tools discussed in Section 4 to obtain lower-variance estimates and gradients. Below we prove a small result to show the utility of this:

**Proposition 5.** *Under the conditions of Proposition 3 and y* ∈ {+1, −1}*,*

$$\begin{aligned} \mathbb{V}\_Q[\ell\_{\text{lin}}(\mathcal{F}\_Q^\*(\mathbf{x}), y)] &\leq \mathbb{V}\_Q[\ell\_{\text{lin}}(\mathcal{F}\_Q(\mathbf{x}), y)] \\ &\leq \mathbb{V}\_{f \sim Q}[\ell\_{0-1}(f(\mathbf{x}), y)] = \frac{1}{4}(1 - |F\_Q(\mathbf{x})|^2). \end{aligned}$$

**Proof.**

$$\mathbb{E}\left[\mathbb{V}\_{\mathcal{Q}}[\ell\_{\text{lin}}(\mathbb{F}\_{\mathcal{Q}}(\mathbf{x}),y)]\right] = \mathbb{E}\_{\mathcal{Q}}\left|\frac{1}{2}(y\mathbb{F}\_{\mathcal{Q}}(\mathbf{x}) - y\mathbb{F}\_{\mathcal{Q}}(\mathbf{x}))\right|^2 = \frac{1}{4}\mathbb{V}\_{\mathcal{Q}}[\mathbb{F}\_{\mathcal{Q}}(\mathbf{x})]^2$$

and a similar result for *F*ˆ∗ *Q*. *f* = *F*ˆ *Q* if *T* = 1 and lin(*f*(*x*), *y*) = -0−1(*f*(*x*), *y*). The result then follows from this and Proposition 3.

Combining (11) with Theorem 1, we obtain a directly optimizable, differentiable bound on the misclassification loss without introducing the above-mentioned factor of 2.

**Theorem 2.** *Given P on θ and α* > 1*, for all Q on θ and λ* > 1 *simultaneously with probability at least* 1 − *δ over S* ∼ D*m,*

$$\mathbb{E}\_{\theta \sim Q} R^{0-1}(f\_{\theta}) \le \Phi\_{\lambda/m}^{-1} \left[ R\_S^{\text{lin}}(F\_Q) + \frac{\alpha}{\lambda} \Delta \right]$$

*with* Φ−1 *γ* (*t*) = 1−exp(−*γt*) 1−exp(−*γ*) *, fθ* : R*d* → {+1, −1}, *θ* ∈ Θ*, and* Δ = KL(*Q*|*P*) − log *δ* + 2 loglog *α*2*λ* log *α .*

Thus, for each *λ*, which can be held fixed ("**fix-***λ*") or simultaneously optimized throughout training for automatic regularisation tuning ("**optim-***λ*"), we obtain a gradient descent objective:

$$R\_S^{\rm lin}(\hat{F}\_Q^\*) + \frac{\mathbf{KL}(Q|P)}{\lambda}. \tag{12}$$

#### **6. Experiments**

All experiments (Table 1) run on "binary"-MNIST, dividing MNIST into two classes, of digits 0–4 and 5–9. Neural networks had three hidden layers with 100 units per layer and **sign**, sigmoid (**sgmd**) or **relu** activations, before a single-unit final layer with sign activation. *Q* was chosen as an isotropic, unit-variance normal distribution with initial means drawn from a truncated normal distribution of variance 0.05. The data-free prior *P* was fixed equal to the initial *Q*, as motivated by Dziugaite and Roy [4] (Section 5 and Appendix B).

**Table 1.** Average (from ten runs) binary-MNIST losses and bounds (*δ* = 0.05) for the best epoch and optimal hyperparameter settings of various algorithms. Hyperparameters and epochs were chosen by bound if available and non-vacuous, otherwise by training linear loss. Bold numbers indicate the best values and standard deviation is reported in italics.

The objectives **fix-***λ* and **optim-***λ* from Section 5 were used for batch-size 256 gradient descent with Adam [17] for 200 epochs. Every five epochs, the bound (for a minimising *λ*) was evaluated using the entire training set; the learning rate was then halved if the bound was unimproved from the previous two evaluations. The best hyperparameters were selected using the best bound achieved in these evaluations through a grid search of initial learning rates ∈ {0.1, 0.01, 0.001}, sample sizes *T* ∈ {1, 10, 50, 100}. Once these were selected training was repeated 10 times to obtain the values in Table 1.

*λ* in **optim-***λ* was optimised through Theorem 2 on alternate mini-batches with SGD and a fixed learning rate of 10−4 (whilst still using the objective (12) to avoid effectively scaling the learning rate with respect to empirical loss by the varying *λ*). After preliminary experiments in **fix-***λ*, we set *λ* = *m* = 60,000, the training set size, as is common in Bayesian deep learning.

We also report the values of three baselines: **reinforce**, which uses the fix-*λ* objective without partial-aggregation, forcing the use of REINFORCE gradients everywhere; **mlp**, an unregularised non-stochastic relu neural network with tanh output activation; and the PBGNet model (**pbg**) from Letarte et al. [6]. For the latter, a misclassification error bound obtained through -0−1 ≤ 2lin must be used as their test predictions were made through the sign of a prediction function ∈ [−1, +1], not ∈ {+1, −1}. Further, despite significant additional hyperparameter exploration, we were unable to train a three layer network through the PBGNet algorithm directly comparable to our method, likely because of the exponential KL penalty (in their Equation 17) within that framework; to enable comparison, we therefore allowed the number of hidden layers in this scenario to vary ∈ {1, 2, 3}. Other baseline tuning and setup was similar to the above, see the Appendix A for more details.

During evaluation **reinforce** draws a new set of weights for every test example, equivalent to the evaluation of the other models; but doing so during training, with multiple parallel samples, is prohibitively expensive. Two different approaches to straightforward, not partially-aggregated, gradient estimation for this case suggest themselves, arising from different approximations to the *Q*-expected loss of the minibatch, *B* ⊆ *S* (with data indices B). From the identities

$$\begin{aligned} \mathbb{E}\nabla\_{\boldsymbol{\theta}}\mathbb{E}\_{\theta \sim q\_{\boldsymbol{\theta}}}\mathbb{R}\_{B}(f\_{\boldsymbol{\theta}}) &= \mathbb{E}\_{\theta \sim q\_{\boldsymbol{\theta}}}\frac{1}{|B|}\sum\_{i \in B} \ell(f\_{\boldsymbol{\theta}}(\boldsymbol{x}\_{i}), y\_{i}) \nabla\_{\boldsymbol{\theta}}\log q\_{\boldsymbol{\theta}}(\boldsymbol{\theta}) \\ &= \frac{1}{|B|}\sum\_{i \in B} \mathbb{E}\_{\theta \sim q\_{\boldsymbol{\theta}}}\ell(f\_{\boldsymbol{\theta}}(\boldsymbol{x}\_{i}), y\_{i}) \nabla\_{\boldsymbol{\theta}}\log q\_{\boldsymbol{\theta}}(\boldsymbol{\theta}) \end{aligned}$$

we obtain two slightly different estimators for ∇*φ*E*θ*∼*qφRB*(*fθ* ):

$$\leq \frac{1}{T|B|} \sum\_{t=1}^{T} \sum\_{i \in B} \ell(f\_{\emptyset^{(t,i)}}(\mathfrak{x}\_i), \mathfrak{y}\_i) \nabla\_{\boldsymbol{\Phi}} \log q\_{\boldsymbol{\Phi}}(\boldsymbol{\theta}^{(t,i)})$$

$$\frac{1}{|T|B|}\sum\_{i\in B}\sum\_{t=1}^{T}\ell(f\_{\theta^t}(\mathbf{x}\_i),\underline{y}\_i)\nabla\_{\boldsymbol{\theta}}\log q\_{\boldsymbol{\theta}}(\boldsymbol{\theta}^t).$$

The first draws many more samples and has lower variance but is much slower computationally; even aside from the *O*(|*B*|) increase in computation, there is a slowdown as the optimised BLAS matrix routines cannot be used, and the very large matrices involved may not fit in memory (see for more information [16]).

Therefore, as is standard in the Bayesian Neural Network literature with the pathwise estimator, we use the latter formulation, which has a similar computational complexity to local-reparameterisation and our marginalised REINFORCE estimator (10). We should note though that in preliminary experiments, the alternate estimator did not appear to lead to improved results. This clarifies the advantages of marginalised sampling, which can lead to lower variance with a similar computational cost.

#### **7. Discussion**

The experiments demonstrate that partial-aggregation enables training of multi-layer non-differentiable neural networks in a PAC–Bayesian context, which is not possible with REINFORCE gradients and a multiple-hidden-layer PBGNet [6]. These obtained only vacuous bounds, and our misclassification bounds also improve those of a single-hidden-layer PBGNet.

Our experiments raise a couple of questions: firstly, why is it that lower variance estimates empirically lead to tighter bounds? We speculate that the faster convergence of SGD in this case takes us to a more "local" minimum of the objective, closer to our initialisation. Since most existing PAC–Bayes bounds for neural networks have a very strong dependence on this distance from initialisation through the KL term, this leads to tighter bounds. This distance could also be reduced through other methods we consider out-of-scope, such as the data-dependent bounds employed by Dziugaite and Roy [18] and Letarte et al. [6].

A second and harder question is asking why the non-stochastic mlp model obtains a lower overall error. The bound optimisation is empirically quite conservative, but does not necessarily lead to better generalisation; understanding this gap is a key question in the theory of deep learning.

In future work we will develop significant new tools to extend partial-aggregation to multi-class classification, and to improve test prediction bounds for sign(*F*ˆ *Q*(*x*)) with *T* > 1, as in PBGNet, which gave slightly improved predictive performance despite the inferior theoretical guarantees.

**Author Contributions:** Conceptualization, F.B. and B.G.; Formal analysis, F.B. and B.G.; Methodology, F.B. and B.G.; Project administration, B.G.; Software, F.B.; Writing—original draft, F.B.; Writing review and editing, F.B. and B.G. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This work is supported in part by the U.S. Army Research Laboratory and the U. S. Army Research Office, and by the U.K. Ministry of Defence and the U.K. Engineering and Physical Sciences Research Council (EPSRC) under grant number EP/R013616/1. BG acknowledges partial support from the French National Agency for Research, grants ANR-18-CE40-0016-01 and ANR-18-CE23-0015- 02. FB acknowledges the support from the Foundational Artificial Intelligence Centre for Doctoral Training at University College London.

**Data Availability Statement:** Not applicable

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Further Experimental Details**

*Appendix A.1. Aggregating Biases with the Sign Function*

We used a bias term in our network layers, leading to a simple extension of the above formulation, omitted in the main text for conciseness:

$$\mathbb{E}\_{w \sim \mathcal{N}(\mu, \Sigma), b \sim \mathcal{N}(\beta, \sigma^2)} \operatorname{sign}(w \cdot \mathbf{x} + b) = \operatorname{erf}\left(\frac{\mu \cdot \mathbf{x} + \beta}{\sqrt{2(\mathbf{x}^T \Sigma \mathbf{x} + \sigma^2)}}\right)$$

since *w* · *x* + *b* ∼ N (*μ* · *x* + *β*, *xT*Σ*x* + *σ*2) and

$$\begin{aligned} \mathbb{E}\_{z \sim \mathcal{N}(\boldsymbol{\alpha}, \boldsymbol{\beta}^2)} \operatorname{sign} z &= P(z \ge 0) - P(z < 0) \\ &= \left[ 1 - \Phi(-\boldsymbol{\alpha}/\beta) \right] - \Phi(-\boldsymbol{\alpha}/\beta) \\ &= 2\Phi(\boldsymbol{\alpha}/\beta) - 1 = \operatorname{erf}(\boldsymbol{\alpha}/\sqrt{2}\beta). \end{aligned}$$

The bias and weight co-variances were chosen to be diagonal with a scale of 1, which leads to some simplification in the above.

## *Appendix A.2. Dataset Details*

We used the MNIST dataset version 3.0.1, available online at http://yann.lecun.com/ exdb/mnist/ (accessed on 4 June 2021), which contains 60,000 training examples and 10,000 test examples, which were used without any further split, and rescaled to lie in the range [0, 1]. For the "binary"-MINST task, the labels +1 and −1 were assigned to digits in {5, 6, 7, 8, 9} and {0, 1, 2, 3, 4}, respectively, and images were scaled into the interval [0, 1].

#### *Appendix A.3. Hyperparameter Search for Baselines*

The baseline comparison values offered with our experiments were optimized similarly to the above, for completeness we report everything here.

The MLP model had three hidden ReLu layers of size 100 each trained with Adam, a learning rate ∈ {0.1, 0.01, 0.001} and a batch size of 256 for 100 epochs. Complete test and train evaluation was performed after every epoch, and in the absence of a bound, the model and epoch with lowest train linear loss was selected.

For PBGNet we choose the values of hyperparameters from within these values giving the least bound value. Note that, unlike in [6], we do not allow the hidden size to vary {∈ 10, 50, 100}, and we use the entire MNIST training set as we do not need a validation set. While attempting to train a three hidden layer network, we also searched through the hyperparameter settings with a batch size of 64 as in the original, but after this failed, we returned to the original batch size of 256 with Adam. All experiments were performed using the code from the original paper, available at https://github.com/gletarte/dichotomizeand-generalize (accessed on 4 June 2021).

Since we were unable to train a multiple-hidden-layer network through the PBGNet algorithm, for this model only we explored different numbers of hidden layers ∈ {1, 2, 3}.

#### *Appendix A.4. Final Hyperparameter Settings*

In Table A1 we report the hyperparameter settings used for the experiments in Table 1 after exploration. To save computation, hyperparameter settings that were not learning (defined as having a whole-train-set linear loss of > 0.45 after ten epochs) were terminated early. This was also done on the later evaluation runs, where in a few instances the fix-*λ* sigmoid network failed to train after ten epochs; to handle this we reset the network to obtain the main experimental results.

**Table A1.** Chosen hyperparameter settings and additional details for results in Table 1. Best hyperparameters were chosen by bound if available and non-vacuous, otherwise by best training linear loss through a grid search as described in Section 6 and Appendix A.3. Run times are rounded to nearest 5 min.

For clarity we repeat here the hyperparameter settings and search space:

## *Appendix A.5. Implementation and Runtime*

Experiments were implemented using Python and the TensorFlow library [19]. Reported approximate runtimes are for execution on a NVIDIA GeForce RTX 2080 Ti GPU.

#### **References**

**Dimitri Meunier 1 and Pierre Alquier 2,\***

**Abstract:** Online learning methods, similar to the online gradient algorithm (OGA) and exponentially weighted aggregation (EWA), often depend on tuning parameters that are difficult to set in practice. We consider an online meta-learning scenario, and we propose a meta-strategy to learn these parameters from past tasks. Our strategy is based on the minimization of a regret bound. It allows us to learn the initialization and the step size in OGA with guarantees. It also allows us to learn the prior or the learning rate in EWA. We provide a regret analysis of the strategy. It allows to identify settings where meta-learning indeed improves on learning each task in isolation.

**Keywords:** meta-learning; hyperparameters; priors; online learning; Bayesian inference; online optimization; gradient descent

## **1. Introduction**

In many applications of modern supervised learning, such as medical imaging or robotics, a large number of tasks is available but many of them are associated with a small amount of data. With few datapoints per task, learning them in isolation would give poor results. In this paper, we consider the problem of learning from a (large) sequence of regression or classification tasks with small sample size. By exploiting their similarities we seek to design algorithms that can utilize previous experience to rapidly learn new skills or adapt to new environments.

Inspired by human ingenuity in solving new problems by leveraging prior experience, *meta-learning* is a subfield of machine learning whose goal is to automatically adapt a learning mechanism from past experiences to rapidly learn new tasks with little available data. Since it "learns the learning mechanism" it is also referred to as *learning-to-learn* [1]. It is seen as a critical problem for the future of machine learning [2]. Numerous formulations exist for meta-learning and we focus on the problem of *online meta-learning* where the tasks arrive one at a time and the goal is to efficiently transfer information from the previous tasks to the new ones such that we learn the new tasks as efficiently as possible (this has also been refered to as *lifelong learning*). Each task is in turn processed *online*. To sum up, we have a stream of tasks and for each task a stream of observations.

In order to solve online tasks, diverse well-established strategies exist: perceptron, online gradient algorithm (OGA), online mirror descent, follow-the-regularized-leader, exponentially weighted aggregation (EWA, also refered to as *generalized Bayes* etc.) We refer the reader to [3–6] for introductions to these algorithms and to so-called regret bounds, that control their generalization errors. We refer to these algorithms as the *within-task* strategies. The big challenge is to design a meta-strategy that uses past experiences to adapt a within-task strategy to perform better on the next tasks.

In this paper, we propose a new meta-learning strategy. The main idea to learn the tuning parameters is to minimize its regret bound. We provide a meta-regret analysis for our strategy. We illustrate our results in the case where the within-task strategy is the online gradient algorithm, and exponentially weighted aggregation. In the case of OGA, the tuning parameters considered are the initialization and the gradient steps. For EWA,

**Citation:** Meunier, D.; Alquier, P. Meta-Strategy for Learning Tuning Parameters with Guarantees. *Entropy* **2021**, *23*, 1257. https://doi.org/ 10.3390/e23101257

Academic Editor: Gholamreza Anbarjafari

Received: 9 August 2021 Accepted: 23 September 2021 Published: 27 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

we consider either the learning rate, or the prior. In each case, we compare the regret incurred when learning the tasks in isolation to our meta-regret bound. This allows us to identify settings where meta-learning indeed improves on learning in isolation.

#### *1.1. Related Works*

Meta-learning is similar to multitask learning [7–9] in the sense that the learner faces many tasks to solve. However, in multitask learning, the learner is given a fixed number of tasks, and can learn the connections between these tasks. In meta-learning, the learner must prepare to face future tasks that are not given yet.