Ghahramani Julyan Arbel

Riccardo Corradin

Micha L Lewandowski

Comment on Article by

We propose a simulation study to emphasise the difference between Variation of Information and Binder's loss functions in terms of number of clusters estimated by means of (1) the use of the MCMC output only and (2) a "greedy" method.

Wade and Ghahramani's paper is a very neat contribution to Bayesian cluster analysis in at least two respects: (i) by formalizing cluster credible coverage via Hasse diagrams, and (ii) by recasting the problem in a decision theory framework, with tangible improvements brought by the Variation of Information (VI) loss function [START_REF] Meilȃ | Comparing clusterings-an information based distance[END_REF]) over [START_REF] Binder | Bayesian cluster analysis[END_REF][START_REF] Dahl | Model-based clustering for expression data via a dirichlet process mixture model. Bayesian inference for gene expression and proteomics[END_REF].

We propose a simulation study implementing two algorithms provided by Wade and Ghahramani's package mcclust.ext for finding the argument minimizing the posterior expected loss: (1) the draw algorithm, which restricts the minimization problem to the MCMC output, and (2) the greedy algorithm, which is more reliable as it also scans the neighbouring clusters of the MCMC output, but with a larger computational cost. While increasing the sample size, we point out the radically different behavior of the number of clusters estimated under VI and Binder, especially with the greedy algorithm.

Our simulation study is based on the same data generation as in the first example of Section 6.1 in [START_REF] Wade | Bayesian cluster analysis: Point estimation and credible balls[END_REF]: a mixture of four Gaussian distributions equally weighted with means (±2, ±2) and identity covariance matrix. We estimated the model using a marginal approach provided by BNPmix 1 R package. We synthesised the output with mcclust.ext package. 2 The Dirichlet process mixture model was estimated with mass parameter fixed to 1, and by specifying an independent base measure on locations and scales, with a 0-vector prior mean for the location component and an identity matrix prior mean for the scale component (25 000 iterations with 5 000 burnin period). We considered four different sample sizes n = {20, 40, 100, 300}.

The results are shown in Figure 1. With the draw algorithm, the cluster estimates under both losses are quite close in terms of number of clusters. In contrast, the greedy algorithm leads to cluster estimates obtained via Binder's loss function with excessive size, while that obtained via VI remains coherent with the number of components of the model (four).

Similarly to the authors' finding, ours' indicates that Binder's loss function exhibits an undesirable property of overestimating the number of clusters (Miller andHarrison, 2013, 2014). Variation of Information tends to lessen this problem. As alluded to by the * Univ. Grenoble Alpes, Inria, CNRS, LJK, 38000 Grenoble, France. julyan.arbel@inria.fr; michal.lewandowski@inria.fr † DISMEQ, University of Milano Bicocca, 20126 Milano MI, Italy. riccardo.corradin@unimib.it 1 Package available at https://github.com/rcorradin/BNPmix, can be installed via devtools.

2 Code of the simulation study available at https://github.com/rcorradin/WGdiscussion.

imsart-ba ver. 2014/10/16 file: main.tex date: March 22, 2018

Comment on Article by Wade and Ghahramani

authors, a theoretical study of the asymptotic behavior of the VI estimator would be very timely. Especially in light of the recent contribution by [START_REF] Rajkowski | Analysis of MAP in CRP Normal-Normal model[END_REF] about the asymptotic behavior of the cluster estimator under the 0 -1 loss (MAP estimator).


Julyan Arbel
email: julyan.arbel@inria.fr

󰂏

Clustering Milky Way's Globulars: a Bayesian Nonparametric Approach

This chapter presents a Bayesian nonparametric approach to clustering, which is particularly relevant when the number of components in the clustering is unknown. The approach is illustrated with the Milky Way's globulars, that are clouds of stars orbiting in our galaxy. Clustering globulars is key for better understanding the Milky Way's history. We define the Dirichlet process and illustrate some alternative definitions such as the Chinese restaurant process, the Pólya Urn, the Ewens sampling formula, the stick-breaking representation through some simple R code. The Dirichlet process mixture model is presented, as well as the R package BNPmix implementing Markov chain Monte Carlo sampling. Inference for the clustering is done with the variation of information loss function.

R requirements

The code used during the presentation of the Stat4Astro summer school is available at the url: https://github.com/jarbel/Stat4Astro-Autrans. Additionally, the code used to generate the plots of this chapter is displayed in the text. This requires the following R packages: ggplot2, hexbin, viridis, gridExtra, ggpubr, rgl for graphical tools, reshape2 for operations on similarity matrices, mcclust and mcclust.ext for clustering estimation. needed_packages <-c("ggplot2", "hexbin", "viridis", "gridExtra", "ggpubr", "rgl", "reshape2", "dplyr", "mcclust") new_packages <-needed_packages[ !(needed_packages %in% installed.packages()[, "Package"])] if (length(new_packages))

install.packages(new_packages) lapply(needed_packages, require, character.only = TRUE)

install.packages("mcclust.ext_1.0.tar.gz", repos = NULL, type = "source") file.remove("mcclust.ext_1.0.tar.gz") library("mcclust.ext")

Introduction and motivation

Globulars1 are sets of stars orbiting some galactic center. The globular data we are considering here was studied in the 2015 Edition of the Stat4Astro school by [START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF] who used phylogenetic classification. The data are available on GitHub and can be downloaded as follows: spectra <-read.csv( "https://github.com/jarbel/Stat4Astro-Autrans/blob/ master/Talk_Arbel/bnp_code/data/GC4c_groups.dat", sep="")

It lists a total of dim(spectra) [START_REF] Aldous | Exchangeability and related topics[END_REF]=54 globulars for which dim(spectra) [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF]=7 variables are available2 : GC stands for the globular identifier; logTe is the logarithm of the maximum effective temperature on the horizontal branch; FeH denotes the metallicity; MV is the absolute V magnitude, which relates to both the brightness and the mass of the globular; Age of the globular; Grp4c and Grp3c are the phylogenetic classifications of [START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF] obtained by using respectively the four variables logTe, FeH, MV, Age and the three variables logTe, FeH, MV.

## [START_REF] Aldous | Exchangeability and related topics[END_REF] 54 7 ## [START_REF] Aldous | Exchangeability and related topics[END_REF] "GC" "logTe" "FeH" "MV" "Age" "Grp4c" "Grp3c" names(spectra) ## [START_REF] Aldous | Exchangeability and related topics[END_REF] "GC" "logTe" "FeH" "MV" "Age" "Grp4c" "Grp3c"

By using the additional spatial coordinates available on the Wikipedia list 2 , we can obtain a spatial scatterplot of the globulars. In Figure 1, the globulars present in the study (that is for which we have measurements for the above mentionned variables) are depicted in purple, the others in green

The two clusterings Grp4c and Grp3c are of respective size 3 and 4, and the cluster sizes are obtained with the table command as follows: table(spectra$Grp3c) ## ## 1 2 3 4 ## 17 7 21 9 11 18 From now on, we will focus on the three variables logTe, FeH, MV only. A 2-D representation for each pair of variables, as well as a side density plot for each variable can be obtained by: col <-plasma(1) alpha <-.8 spVF <-ggscatter(spectra, x = "MV", y = "FeH", color = col, size = 3, alpha = alpha)+ border() spLF <-ggscatter(spectra, x = "logTe", y = "FeH", color = col, size = 3, alpha = alpha)+ border() spLV <-ggscatter(spectra, x = "logTe", y = "MV", color = col, size = 3, alpha = alpha) + border() Vplot <-ggdensity(spectra, "MV", fill = col) Lplot <-ggdensity(spectra, "logTe", fill = col) Fplot <-ggdensity(spectra, "FeH", fill = col) + rotate() Vplot <-Vplot + clean_theme() Lplot <-Lplot + clean_theme() Fplot <-Fplot + clean_theme() ggarrange(Vplot, Lplot, NULL, spVF, spLF, Fplot, NULL, spLV, NULL, ncol = 3, nrow = 3, align = "hv", widths = c(2, 2, 1), heights = c(1, 2, 2))

table(spectra$Grp4c) ## ## 1 2 3 ## 25

MV

These scatter plots hardly identify any clustering structure amongst the globulars. However, astrophysicists expect several populations of globulars identified by similar time of formation, and similar chemical and physical conditions. The granularity of the populations of globulars may depend on the size of the considered sample: a small number of observations would likely lead to little discriminating power, while a large sample size would provide more evidence for identifying more distinct clusters. This setting where the number of clusters might grow with the sample size is well suited to a Bayesian nonparametric approach to clustering, which this chapter aims at introducing.

The rest of the chapter is organised as follows. Model-based clustering and the Dirichlet process are introduced in Section 3 and Section 4. We conclude with an illustration to the Globular dataset.

Model-based clustering

Mixture models are creating flexible models starting from simple ones. For instance, combining two unimodal densities p 1 and p 2 into πp 1 + (1 -π)p 2 , for π ∈ (0, 1) can create a bimodal distribution. The aim of mixtures is to increase the modelling capacities by combining simple distributions into flexible distributions which might display multimodality, skewness, etc. A mathematical definition of a mixture density is a convex combination of densities. Each density can be interpreted as a sub-population. Observations are associated to sub-populations through latent (un-observed) variables called allocation variables, which play a major role in devising sampling algorithms.

Consider a parametric family of distributions with densities p(•|φ), φ element of some parameter space Θ, with respect to a common measure. Broadly speaking, mixtures operate in a discrete way (finite or infinite countable number of classes) or in a continuous way (infinite uncountable number of classes). As an example of continuous mixtures, a Student t-distribution can be written as a mixture of Gaussian kernels with fixed mean and variance averaged over an inverse gamma distribution (playing the role of a mixing distribution). We will instead focus here on discrete mixtures.

Going back to the simple mixture of two densities, the mixture density takes on the form

πp 1 (•|φ 1 ) + (1 -π)p 2 (•|φ 2 )
, where the unknown parameters are (π, φ 1 , φ 2 ). Modelling data x = (x 1 , . . . , x N ) with such a density can be thought as the following two-step procedure, for i ∈ {1, . . . , N }:

1. toss a coin, which samples from one of the two classes θ i ∈ {φ 1 , φ 2 }. In other words, the class identified with parameter φ 1 has probability π, while the φ 2 class has the complement probability 1 -π;

2. sample the observation from the corresponding density x i ∼ p(x|θ i ).

The coin-tossing step can be equally thought of as a throw of dice whose number of facets equals the number of components in the mixture, and each facet has a probability given by the component weight. This mixing distribution takes the mathematical form of a convex combination πδ φ1 + (1 -π)δ φ2 of Dirac masses at φ 1 and φ 2 , each of which identifies with a class of the mixture.

We speak about nonparametric (discrete) mixtures when the number of classes is (countable) infinite. In which case the dice has an infinite number of facets. We index the classes by the positive integers, and denote the parameters by φ k and weights by π k . Then those weights must sum up to one, and the mixing measure can be written as

G = ∞ 󰁛 k=1 π k δ φ k ,
which is a probability measure. Sampling from the mixture distribution

p G (•) = ∞ 󰁛 k=1 π k p(•|φ k )
can again be done in the sequential way, for i ∈ {1, . . . , N }:

1. throw the infinite dice G: θ i ∼ G;
2. sample the observation from the corresponding density x i ∼ p(x|θ i ).

The mixture density p G has an infinite number of parameters: the vector of weights (π 1 , π 2 , . . .) which is an element of the infinite simplex, and the vector of class parameters (φ 1 , φ 2 , . . .). A Bayesian approach to the problem requires endowing the parameters with some prior distribution. The Dirichlet process is the most prominent example of such a prior distribution. It is parametrised by the base measure, denoted by G 0 , and the precision parameter, a positive scalar denoted by α. Figure 2 displays plate representations of a Bayesian and a non Bayesian (or frequentist) approach to the mixture model. Such a plate representation makes it clear that a Bayesian approach add a hierarchical layer (that of the prior) to the non Bayesian set-up.

Choosing the number of components to include in a finite mixture is often a delicate question, as one rarely knows it a priori. There are mainly three strategies: i) to fit several finite mixture models with a range of plausible values for the number of components; then choose the model which maximises some information criterion; ii) consider the number of components as Bayesian parameter, which is to say endow it with a prior distribution on positive integers, such as a Poisson distribution. Such a model is termed a mixture of finite mixtures (MFM) by [START_REF] Miller | Mixture models with a prior on the number of components[END_REF]; or iii) let this number be a priori infinite (or as large as needed) and select it or estimate it a posteriori. We focus on this last option in the next section.

Bayesian nonparametrics around the Dirichlet process

Most of the models described in the previous chapters of this book are parametric: they can be described by a finite and fixed number of parameters. This number of parameters is independent of the dataset. Very convenient models in a Bayesian setting include conjugate models, where the posterior distribution has the same form as the prior distribution. For example, normal prior and normal likelihood, or beta prior and binomial likelihood, or gamma prior and Poisson likelihood. In contrast, we are dealing in this chapter with nonparametric models. Let us put it straightaway, the nonparametric saying is not the most fortunate, as nonparametric models do have parameters, many of them. [START_REF] Bernardo | Bayesian theory[END_REF] define a Bayesian nonparametric model as a probability model with infinitely many parameters, also referred to by [START_REF] Müller | Bayesian nonparametric inference-why and how[END_REF] as a model with massively many parameters. I think there are three ways this large number of parameters can be thought of: the number of parameters is infinite, or it is random, or it is growing with the data sample size.

It goes without saying that parametric models are easier to handle than their nonparametric counterpart, that are computationally and analytically more challenging. It is also true that interpreting a small and fixed number of parameters is likely to be easier than for, say, an infinite number of them. However, the computational and analytical extra burden is the price to pay for flexibility. Nonparametric models are less prone to misspecification than parametric models, which require a strong belief in the particular structure they imply. Without such a belief in the parametric assumption, a model might not be reliable. By nature, nonparametric models are well suited to study curves in general, so typical applications include density estimation and regression estimation. Clustering is also an example which has historically played a central role within Bayesian nonparametrics.

In this section we provide basic facts about the Dirichlet process, precisely different possible representations: the definition via the finite dimensional marginal distribution in Definition 4.2, the Chinese restaurant process in Proposition 4.3, the posterior distribution in Theorem 4.1, the Pólya Urn in Proposition 4.2 and the stick-breaking representation in Theorem 4.2.

Definition

The Dirichlet process [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF] plays a central role in Bayesian nonparametrics. A Dirichlet process can be viewed as a random variable where the variable is a probability measure. It has two parameters: the base measure, denoted by G 0 , and the precision parameter, a positive scalar denoted by α.

[9] defines the Dirichlet process by its finite marginal distributions: ie, how does the random measure spread its mass in the sample space? The answer is: as a Dirichlet distribution (which, in passing, explains the name). Recall that the Dirichlet distribution generalises the Beta distribution to any dimension k ≥ 2. The Dirichlet distribution in R k is restricted to the unit simplex, with density parametrised by k positive scalars α 1 , . . . , α k and proportional to

x α1-1 1 • • • x α k -1 k . Definition 4.1 (Dirichlet distribution).
A Dirichlet distribution on a simplex ∆ K is a probability distribution with parameters α i > 0 and a density function

f (x 1 , . . . , x K ; α 1 , . . . , α K ) = 1 B(α) K 󰁜 i=1 x αi-1 i .
It is common to refer to Dirichlet distribution as Dir(x 1 , . . . , x k ). Let us remark that Remark 4.1. The Dirichlet distribution is conjugate prior for the multinomial distribution.

Consider a finite partition of the space, denoted by (A 1 , . . . , A k ). The mass allocated by a Dirichlet process G to each region A j is a random variable G(A j ), thus giving rise to a random vector G(A 1 ), . . . , G(A k ) for the whole partition. Then G follows a Dirichlet process with parameters α and G 0 if the random vector G(A 1 ), . . . , G(A k ) is distributed as a Dirichlet distribution with parameters αG 0 (A 1 ), . . . , αG 0 (A k ). Definition 4.2 (Dirichlet process, [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF]). A random probability measure G follows a Dirichlet process with parameters α and G 0 on some space if for any finite partition (A 1 , . . . , A k ) of the space,

(G(A 1 ), . . . , G(A k )) ∼ Dir(αG 0 (A 1 ), . . . , αG 0 (A k )).
Note that this definition entails the strong result that such finite dimensional distributions consistently define a stochastic process. Figure 3 shows different realisations of a Dirichlet process in R 2 , with a standard Gaussian base measure G 0 and precision parameter equal to one, and with different partitions. Of course, not all of the sample space R 2 can be represented, but most of the mass is captured in the represented part of the space.

Properties

From Definition 4.2, we see that any measure set A of the space receives mass according to the following Beta distribution

P (A) ∼ Beta(αG 0 (A), α(1 -G 0 (A)). (1)
We have the following moments.

E[G(A)] = G 0 (A) Var[G(A)] = G 0 (A)(1 -G 0 (A)) 1 + α Cov[G(A), G(B)] = G 0 (A ∩ B) -G 0 (A)G 0 (B) 1 + α
Proof. We will make use of Equation [START_REF] Aldous | Exchangeability and related topics[END_REF]. From this we obtain

E(P (A)) = αP 0 (A) α(P 0 (A) + 1 -P 0 (A)) = P 0 (A)
and

Var(P (A)) = α 2 P 0 (A)(1 -P 0 (A)) α 2 (α + 1) .
We derive the covariance term in two cases, firstly taking into consideration the one with A ∩ B = ∅. In that case the sample space is partitioned into A, B and (A∪B) c , the complementary set of A ∪ B, which is equal to

(A ∪ B) c = A c ∩ B c
. Therefore we may write a joint probability vector 󰀃 P (A), P (B),

P (A c ∩ B c ) 󰀄 ∼ Dir 󰀃 αP 0 (A), αP 0 (B), αP 0 (A c ∩ B c ) 󰀄 ,
and hence Cov(P (A), P (B)) = -P 0 (A)P 0 (B)/(1 + α). In the more general case one may decompose

A = (A ∩ B) ∪ (A ∩ B c ) B = (B ∩ A) ∪ (B ∩ A c ),
so that

Cov(P (A), P (B)) = Cov(P (A ∩ B) + P (A ∩ B c ), P (B ∩ A) + P (B ∩ A c ))
and so forth using the linearity of covariance.

In passing, this shows that any two measurable parts of the space which are non intersecting receive mass from the Dirichlet process with negative correlation. Which makes sense since the total measure is constrained to be a probability measure, so more mass in some part of the space means less in the rest of the space.

Another central property is the conjugacy of the Dirichlet process: if some data are sampled from a Dirichlet process G, then the posterior distribution of G conditional on the data is still a Dirichlet process. Theorem 4.1 (Dirichlet process posterior distribution, [START_REF] Ferguson | A Bayesian analysis of some nonparametric problems[END_REF]). Let data X = (X 1 , . . . , X n ) be distributed according to the model

P ∼ DP(αG 0 ) X 1 , . . . , X n |G iid ∼ G.
Then the posterior distribution of G is given by

G|X 1 , . . . , X n ∼ DP 󰀣 αP 0 + n 󰁛 i=1 δ i 󰀤 .
Proof. This posterior can be obtained by remarking that for any finite measurable partition (A 1 , . . . , A k ), the posterior distribution of P (A 1 ), . . . , P (A k ) depends on the observations only via their cell counts (this comes from the tail-free property of the DP). Denote N j = #{1 ≤ i ≤ n : x i ∈ A j }, i.e. the number of observations in each cell of the partition of (A 1 , . . . , A k ).

(A k )), N |P ∼ Multinom k (n, α),
and hence we obtain the prior of the form

p(α) ∝ α αP0(A1)-1 1 • • • α αP0(A k )-1 k , while sampling model is p(N |α) ∝ α N1 1 • • • α N k k .
This results in the posterior of form

p(α|N ) ∝ α αP0(A1)+N1-1 1 • • • α αP0(A k )+N k -1 k = Dir k 󰀃 αP 0 (A 1 ) + N 1 , . . . , αP 0 (A k ) + N k 󰀄 .
Note that the notational convention adopted for the parameters of the Dirichlet process is interesting here: it provides the simple 'product' parameter αP 0 + 󰁓 n i=1 δ i , which can also be decoupled into a posterior precision parameter α n and a posterior base measure G n given by

α n = α + n, G n = α α + n P 0 + 1 α + n n 󰁛 i=1 δ Xi .
This posterior conjugacy, in turn, provides a simple form of the predictive distribution, that is the distribution of a new observation conditional on the data. This predictive is referred to as the Pólya urn, or Blackwell-MacQueen urn. Proposition 4.2 (Predictive distribution, Pólya urn, [START_REF] Blackwell | Ferguson distributions via Pólya urn schemes[END_REF]). In the model of Proposition 4.1, the predictive distribution for a new observation X n+1 is given by

X n+1 |X 1 , . . . , X n ∼ α α + n P 0 + 1 α + n n 󰁛 i=1 δ i .
Proof. This property is a result of taking the expected value of the posterior given in Proposition 4.1.

Such a predictive distribution induces ties in the observations with positive probability. More precisely, n observations sampled from a DP induce a partition of the integers 1, . . . , n. The distribution of this random partition is called the Chinese restaurant process 3 . This culinary metaphor describes the random partition induced by the DP as follows. Customers join a populated table with probability n j /(α + n), where n j denotes the number of clients already sitting around the table or sit at new table with probability α/(α + n). Proposition 4.3 (Chinese restaurant process, [START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF]). A random sample X 1:n from a DP with precision parameter α induces a partition of {1, . . . , n} into k sets of sizes n 1 , . . . , n k with probability

p(n 1 , . . . , n k ) = p({n 1 , . . . , n k }) = α k Γ(α) Γ(α + n) k 󰁜 j=1 Γ(n j ).
Proof. We will use the Pólya urn schema slightly changed by using n 1 , . . . , n k

P(X n+1 |X 1:n ) = α α + n P 0 + 1 α + n k 󰁛 j=1 n j δ X * j .
By exchangeability, the distribution of {n 1 , . . . , n k } does not depend on the order of the observations. Let's compute p(n 1 , . . . , p k ) as the probability of one draw where the first table consists of first n 1 observations etc.

To proceed, let us use Pólya urn scheme: we denote nj = 󰁓 j i=1 n i and hence nk = n, the total number of observations. We can observe the following pattern: first ball opens new table, following n j -1 ones fill in that table and so forth. That quantity can be rewritten as

α k α(α + 1) . . . (α + n -1) k 󰁜 j=1 (n j -1)!,
where one can rewrite both terms using Gamma function Γ(x) = 󰁕 ∞ 0 u x-1 e -u du:

α(α + 1) . . . (α + n -1) = Γ(α + n) Γ(α) ,
and (n j -1)! = Γ(n j ).

Note that for ordered partitions we have

p(n 1 , . . . , n k ) = p(n 1 , . . . , n k ) k! .
The following lines of code sample observations from a Dirichlet process with a base measure in argument, which is also interpreted as the color distribution in the Pólya urn scheme. This is applied to sample 10 observations from a Dirichlet process with a Gaussian N (0, 1) base measure, which is also interpreted as the color distribution in the Pólya urn scheme, and precision parameter varying from 1, 10 to 100.

N_ball <-10 # with alpha = 1 polya_sample <-polya_urn_model(function() rnorm( 1 These experiments illustrate that large values of the mass parameter α tend to produce a larger number of distinct values in a sample of a given size. For instance, the probability that all n = 10 observations are distinct is equal to

P(X 1 , . . . , X n are pairwise distinct) = α n Γ(α) Γ(α + n) ,
which is approaching 1 as α grows, for n fixed.

In our case, this probability is respectively equal to 3e-07, 3e-02 and 0.6 for α equal to 1, 10 and 100, for n = 10.

Stick-breaking representation

The Dirichlet process is a discrete random probability measure which can be represented as a convex combination of infinitely many Dirac masses,

G = ∞ 󰁛 k=1 π k δ φ k .
The stick-breaking representation, due to [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF], provides a constructive way of building the weights (π k ) k of the Dirichlet process. This is done by sequentially breaking a stick of initial unit length, into pieces whose lengths correspond to the (π k ) k . More specifically, we require independent and identically distributed (iid) random variables V k ∼ Beta(1, α). The first weight π 1 corresponds to V 1 . This leaves a piece of length 1 -V 1 , which is broken at V 2 in order to define π 2 = V 2 (1 -V 1 ). And sequentially, the same procedure is applied to the remaining part, which equals (1 -V 1 )(1 -V 2 ) at this second step. It is easy to see that after k steps, one defines

π k = V k (1 -V 1 ) • • • (1 -V k-1
), and the remaining piece has length

(1 -V 1 ) • • • (1 -V k ).
The representation is completed by assuming iid draws from the base measure G 0 for the locations φ k , independent from the V k . Theorem 4.2 (Stick-breaking representation, [START_REF] Sethuraman | A constructive definition of Dirichlet priors[END_REF]).

Let V 1 , V 2 , . . . iid ∼ Beta(1, α) and φ 1 , φ 2 , . . . iid ∼ G 0 be independent random variables. Define π 1 = V 1 , π k = V k (1 -V 1 ) • • • (1 -V k-1 ), for any k ≥ 2. Then G = 󰁓 ∞ k=1 π k δ φ k ∼ DP(αG 0 ).
Proof. We provide a sketch of proof of this result in two steps. First, to show that the remaining stick length at step k, (1

-V 1 ) • • • (1 -V k ), converges to zero as k → ∞.
This ensures that the weights vector lives in the unit simplex, and in turn that the measure 󰁓 ∞ k=1 π k δ φ k is a probability measure. Second, to use the stick-breaking construction to show that the defined G satisfies the distributional equation

G d = V δ φ + (1 -V )G, (2)
where V ∼ Beta(1, α) and φ ∼ G, independently, whose only solution turns out to be the Dirichlet process, by properties of the Dirichlet distribution.

The following code implements the sampling of the first num_weights=50 weights of a DP. These are then plotted in the order they appear in the stick-breaking construction (left) as well as indexed by their corresponding location (right). Colors indicate the stick-breaking order of appearance. Note that the weights are not necessarily strictly decreasing, but only stochastically decreasing. This means they are decreasing in expectation, as can be easily checked:

E(π k ) = E(V k (1 -V 1 ) • • • (1 -V k-1 )) = EV k E(1 -V 1 ) • • • E(1 -V k-1 ) = 1 α + 1 󰀕 α α + 1 󰀖 k-1 .
stick_breaking_process = function(num_weights, alpha) { betas = rbeta(num_weights, 1, alpha) remaining_stick_lengths = c(1, cumprod(1 -betas))[1:num_weights] weights = remaining_stick_lengths * betas weights } num_weights <-50 draw_stick_breaking <-function(alpha) { labels <-1:num_weights locations <-rnorm(num_weights) SB_weights <-stick_breaking_process(num_weights, alpha) df <-data.frame(labels, locations, SB_weights) order_plot <ggplot(df, aes(labels, SB_weights, fill = as.factor(labels))) + geom_bar(stat = "identity") + theme(legend.position="none") location_plot <ggplot(df, aes(locations, SB_weights, fill = as.factor(labels))) + geom_bar(stat = "identity", width = .1) + theme(legend.position="none") grid.arrange(order_plot, location_plot, ncol = 2) } draw_stick_breaking(1) Combining the model-based clustering approach described in Section 3 with the DP, we are ready to work with nonparametric mixture models, where the Dirichlet process can be used as a prior distribution on the mixing probability measure [START_REF] Lo | On a class of Bayesian nonparametric estimates: I. Density estimates[END_REF].

The following code shows densities draws from DP mixture densities with N (0, 1) base measure, varying mass parameter in {1, 10, 100, 1000}, and centered Gaussian kernel with varying variance σ 2 in {0.2, 0.4, 1}. Enlarging α tends to produce flatter densities, while reducing σ 2 tends to produce more irregular ones.

u], N_xaxis)), PU) res <-apply(res, 1, mean) new_draw <-cbind(res, x_axis, alpha, sigma2[u]) result <-rbind(result, new_draw) } } result <-as.data.frame(result) names(result) <-c("density", "x", "alpha", "sigma2") DP_mixt <-qplot(data = result, y = density, x = x, geom = c("line", "area")) + facet_grid(alpha ~sigma2, labeller = label_bquote(rows = alpha == .(alpha), cols = sigma^2 == .(sigma2))) + aes(color = as.factor(alpha)) + theme(legend.position = "none") DP_mixt In this last section, we show how Dirichlet process mixtures can be applied to perform clustering of globulars of our galaxy.

σ 2 = 0.2 σ 2 = 0.4 σ 2 = 1 α = 1 α = 10 α = 100 α = 1000 -2 0 2 -2 0 2 -2 0

Markov chain Monte Carlo sampling

Inference in R for DPM is implemented in packages including DPpackage [START_REF] Jara | DPpackage: Bayesian non-and semi-parametric modelling in R[END_REF], BNPdensity [START_REF] Barrios | Modeling with normalized random measure mixture models[END_REF]. Here we present a recent package called BNPmix [START_REF] Arbel | Dirichlet process mixtures under affine transformations of the data[END_REF]. It uses C++ code with the Rcpp package, which makes the implementation efficient. The package is available on GitHub4 and can be installed with devtools as follows:

library(devtools) install_github("rcorradin/BNPmix")

This package uses d-dimensional location-scale Dirichlet process mixture of Gaussians. This means that it assumes a Gaussian kernel p(X|φ) = Φ d (X|µ, Σ), where the component specific parameter θ is composed of the mean µ and covariance matrix Σ. The base measure G 0 is chosen as an independent product of a Gaussian and an inverse-Wishart distributions

G 0 (dµ, dΣ; π) = N d (dµ; m 0 , B 0 ) × IW (dΣ; ν 0 , S 0 ), (3)
The hyperparameters (m 0 , B 0 , ν 0 , S 0 ) are typically expressed by an empirical Bayes approach. See [START_REF] Arbel | Dirichlet process mixtures under affine transformations of the data[END_REF] for details. Additionally, the precision parameter α is endowed with a gamma distribution with shape t 1 and scale t 2 .

Posterior inference is carried out by a Markov chain Monte Carlo (MCMC) algorithm. The function is called DPmixMulti, and it takes the following arguments:

• the data;

• a grid of the sample space, with the same dimension as that of the data, on which the density is to be evaluated;

• MCMC_param, a list of parameters for the MCMC including number of simulations nsim and burn-in nburn;

• starting_val, a list containing the initial values for the components;

• params, hyperparameters (m 0 , B 0 , ν 0 , S 0 ) for the base measure G 0 of the DPM and (t 1 , t 2 ) for the hyperprior on the precision parameter α. When omitted, those parameters are set in an empirical Bayes way by default.

A grid to evaluate the densities is defined as follows, spanning all data points, plus/minus 10% with respect to the extreme data points. The output MCMC_output of this function is a list of three objects:

• distribution the estimated (posterior mean) distribution evaluated on the grid.

• result_cluster a matrix, each row is an iteration, each column an observation, each entry is a latent component, thus this contains clusterings for each iteration of the chain.

• result_theta a vector, the value of DP precision parameter α over the iterations.

The posterior clustering information is enclosed in the second object listed above, that we can save in a matrix MCMC_clustering:

MCMC_clustering <-MCMC_output[[2]]

Clustering estimation

It should be noted that a DPM assumes a priori an infinite number of components in the mixture (see for instance the stick-breaking representation). As such, it is a misspecified model in essence when it comes to estimating a fixed clustering, that is a clustering which is not deemed to grow as the sampled data grows. However, in the spirit of the famous quote by George Box "All models are wrong but some are useful", misspecified models can be used as long as they provide insightful results.

Clustering with DPM can be done in different ways. Indeed, the MCMC output we have obtained so far essentially consists in many different clusterings of the data, among which one needs to choose a 'best' clustering for some decision rule.

A first approach consists in looking at the number of components, an estimator of which can be obtained with the mode a posterior (MAP) for instance. However, the posterior distribution of the number of components in a DPM turns out to be inconsistent under some model specifications [START_REF] Miller | A simple example of Dirichlet process mixture inconsistency for the number of components[END_REF]. Posterior consistency is a theoretical property of a posterior distribution: when more and more data are collected from some fixed data distribution, we say that the posterior distribution of some parameter is consistent if it converges to a point mass at the true fixed value of this parameter. The inconsistency property provided by [START_REF] Miller | A simple example of Dirichlet process mixture inconsistency for the number of components[END_REF] is as follows. Assume a DPM model with standard Gaussian base measure, precision parameter equal to one, and standard Gaussian kernel. Suppose data are generated from a standard Gaussian. This means that the data are sampled from a very simple mixture which admits a single component. [START_REF] Miller | A simple example of Dirichlet process mixture inconsistency for the number of components[END_REF] study the posterior distribution of the number of clusters in this situation, and they prove that it does not concentrate to a point mass at one, thus proving posterior inconsistency.

Let us plot the histogram of the number of clusters:

MCMC_number_cluster = apply(X = MCMC_clustering, MARGIN = 1, FUN = max) df_cluster <-data.frame(nb = MCMC_number_cluster) c <-ggplot(df_cluster, aes(nb,..density..)) c + geom_histogram(breaks=(0:max(df_cluster)), aes(fill=..density..) , col="white", alpha = .9 ) + scale_fill_viridis(option = "inferno") + labs(title="Histogram for number of clusters") + labs(x="Number of clusters", y="Density") The MAP estimator for the number of clusters is three, obtained as:

which.max(table(MCMC_number_cluster)) ## 3 ## 3
Note that the above histogram is also quite in agreement with an estimate of four clusters.

Stepping back, recall that a Bayes estimator is obtained from a formal decision theory rule: given a loss function, a Bayes estimator minimizes the posterior expected loss. For instance, with Euclidean parameter spaces,

• the L 2 , squared loss provides the posterior mean,

• the L 1 , absolute loss provides the posterior median,

• the 0 -1 loss provides the mode a posteriori (MAP).

We focus now on a loss function L on clusterings. The posterior expected loss of clustering c ′ , denoted by L(c ′ ), is obtained by averaging the loss with respect to posterior distribution, over the set of all partitions of the integers 1, . . . , n denoted by Several losses have been considered in the literature:

A n L(c ′ ) = 󰁛 c∈An L(c, c ′ )p(c|x),
• 0-1 loss [START_REF] Rajkowski | Analysis of MAP in CRP Normal-Normal model[END_REF],

• Binder loss [START_REF] Dahl | Model-based clustering for expression data via a Dirichlet process mixture model. Bayesian inference for gene expression and proteomics[END_REF],

• Variation of information (VI) [START_REF] Wade | Bayesian cluster analysis: Point estimation and credible balls[END_REF].

The 0-1 loss gives rise to the MAP estimator:

L 0-1 (c ′ ) = 󰁛 c∈An L 0-1 (c, c ′ )p(c|x) = 󰁛 c∈An, c∕ =c ′ p(c|x) = 1 -p(c ′ |x)
which is to say that the expected loss of c ′ is all the posterior mass except that of c ′ . So that it is easily minimized at the value c ′ which has maximum posterior weight:

ĉ = arg min c ′ ∈An L 0-1 (c ′ ) = arg max c ′ ∈An p(c ′ |x) := M AP.
Negative results by [START_REF] Rajkowski | Analysis of MAP in CRP Normal-Normal model[END_REF] show that the mode a posteriori (MAP) is inconsistent. Instead of the 0-1 loss, one can resort to the Variation of information (VI) which was devised by [START_REF] Meilă | Comparing clusterings-an information based distance[END_REF] for clustering comparison. It stems from information theory, and compares information (in terms of Shannon entropy H) in two clusterings with information shared between the two clusterings (I), see [START_REF] Meilă | Comparing clusterings-an information based distance[END_REF] Estimation under the Variation of information loss function was recently studied by [START_REF] Wade | Bayesian cluster analysis: Point estimation and credible balls[END_REF], with the interesting finding that this loss function tends to reduce the over-estimation of the number of cluster that is commonly obtained under Binder loss, for instance by [START_REF] Dahl | Model-based clustering for expression data via a Dirichlet process mixture model. Bayesian inference for gene expression and proteomics[END_REF]. See also [START_REF] Arbel | Discussion of "Bayesian Cluster Analysis: Point Estimation and Credible Balls[END_REF] for a further comparison between Binder and VI on varying sample sizes. We implement the Variation of information approach by using the mcclust.ext R package developed by [START_REF] Wade | Bayesian cluster analysis: Point estimation and credible balls[END_REF]. This requires to compute the posterior similarity matrix associated to the MCMC output, that is a matrix whose entries represent the posterior probability that two observations are clustered together (or rather, a Monte Carlo approximation of it).

posterior_similarity_matrix = comp.psm(MCMC_clustering)

It can be represented as follows, where darker colors depict a higher posterior probability of shared clustering. Of course, this optimization is only done approximately by the mcclust.ext R package, which scans the MCMC output and a neighbouring region of it in order to choose the best clustering.

clustering_Binder=minbinder.ext(posterior_similarity_matrix, method = "greedy") table(clustering_Binder$cl) ## ## 1 2 3 4 5 ## 16 6 29 2 1 clustering_VI=minVI(posterior_similarity_matrix, method = "greedy") table(clustering_VI$cl) ## ## 1 ## 54 Thus, the Binder loss produces an estimated clustering with 5 groups, whereas under the Variation of information, all the observations are gathered into a single cluster of size 54. The interpretation of these results is that under the conditions of the assumed Dirichlet process mixture model, data are not deemed heterogeneous enough to justify multiple clusters under VI loss, whereas they do under Binder loss. The estimated clusterings and the ones obtained by [START_REF] Fraix-Burnet | The environment of formation as a second parameter for globular cluster classification[END_REF] are represented below in a 3D-like fashion by using the rgl package function plot3d.

Figure 1 :

 1 Figure 1: Globular coordinates. Gobulars included in the present study are depicted in purple, the others in green

  (a) Frequentist setting: non random G. (b) Bayesian approach: G ∼ DP(αG0).

Figure 2 :

 2 Figure 2: Plate representations of a nonparametric mixture with mixing measure G.

Figure 3 :

 3 Figure 3: Dirichlet process realisations in R 2 , varying in rows. Different partitions of R 2 are used in columns. The base measure G 0 is a standard Gaussian and the precision parameter is equal to one.

Proposition 4 . 1 (

 41 Dirichlet process moments). If G ∼ DP(αG 0 ), then for any measurable sets A and B

  Then we have 󰀃 P (A 1 ), . . . , P (A k ) 󰀄󰀏 󰀏 X 1:n d = 󰀃 P (A 1 ), . . . , P (A k ) 󰀄󰀏 󰀏 N 1:k . Let us use the shorthand notation: α = (α 1 , . . . , α k ) = (P (A 1 ), . . . , P (A k )) and N = (N 1 , . . . , N k ). Then 󰀝 α ∼ Dir k (αP 0 (A 1 ), . . . , αP 0

  polya_urn_model <-function(base_measure, N_ball, alpha) { balls <-c() for (i in 1:N_ball) { if (runif(1) < alpha / (alpha + length(balls))) { # Add a new ball color. new_color <-base_measure() balls <-c(balls, new_color) } else { # Pick out a ball from the urn, and add back a # ball of the same color. ball <-balls[sample(1:length(balls), 1)] balls <-c(balls, ball) } } balls }

  seq(-3, 3, length = N_xaxis) result <-NULL for (alpha in alpha_vect) { PU <-polya_urn_model(function() rnorm(1), N_draws, alpha) for (u in 1:N_urns) { res <-mapply(function(mean) dnorm(x_axis, rep(mean, N_xaxis),

5

 5 Application to clustering of globulars of our galaxy

  data <-spectra[, c[START_REF] Antoniak | Mixtures of Dirichlet processes with applications to Bayesian nonparametric problems[END_REF][START_REF] Arbel | Discussion of "Bayesian Cluster Analysis: Point Estimation and Credible Balls[END_REF][START_REF] Arbel | Dirichlet process mixtures under affine transformations of the data[END_REF]] grid <-expand.grid(seq(range(data[,1])[1]-.1 * diff(range(data[,1])), range(data[,1])[2]+.1 * diff(range(data[,1])), length.out = 40), seq(range(data[,2])[1]-.1 * diff(range(data[,2])), range(data[,2])[2]+.1 * diff(range(data[,2])), length.out = 40), seq(range(data[,3])[1]-.1 * diff(range(data[,3])), range(data[,3])[2]+.1 * diff(range(data[,3])), length.out = 40)) We can now run the MCMC function DPmixMulti: MCMC_param <-list(nsim = 10^4, nburn = 5*10^3) MCMC_output <-DPmixMulti(data = as.matrix(data), grid = grid, MCMC_param = MCMC_param)

  and the decision is taken by choosing the best partition ĉ = arg min c ′ ∈An 󰁛 c∈An L(c, c ′ )p(c|x).

  for details:VI(c, 󰁥 c) = H(c) + H(󰁥 c) -2I(c, 󰁥 c).

  The next step consists in resolving the minimization problem ĉ = arg min c ′ ∈An 󰁛 c∈An VI(c, c ′ )p(c|x).

  subtitle = c("Binder loss, 5 groups", "VI loss, 1 group", "Fraix-Burnet (2009), 4 groups", "Fraix-Burnet (2009), 3 groups") clustering = cbind(clustering_Binder$cl, clustering_VI$cl, Grp3c, Grp4c) for(i in 1:4){ cl = clustering[, i] plot3d(logTe, FeH, MV, type="s", size=2, col=plasma(max(cl))[cl], box = FALSE, sub = subtitle[i]) }

Globulars are more commonly called globular clusters in the literature, though we shall prefer the phrasing 'globulars' to 'globular clusters' in order to avoid ambiguous terms like 'globular clusters clusters'.

A more comprehensive list of globulars can be accessed at https://en.wikipedia.org/wiki/List_of_ globular_clusters, though it contains information about magnitude and diameter only.

According to[START_REF] Aldous | Exchangeability and related topics[END_REF], the restaurant analogy is due to Jim Pitman and Lester Dubins.

At the url: https://github.com/rcorradin/BNPmix.

Acknowledgements

I would like to thank Michał Lewandowski for helping with typing parts of this chapter and with merging the Globular data by [10] with the Globular table that can be found on Wikipedia, see Footnote 2.


Mariia Vladimirova
email: mariia.vladimirova@inria.fr

Julyan Arbel
email: julyan.arbel@inria.fr

Pablo Mesejo
email: pmesejo@decsai.ugr.es

Bayesian neural networks increasingly sparsify their units with depth

We investigate deep Bayesian neural networks with Gaussian priors on the weights and ReLU-like nonlinearities, shedding light on novel sparsity-inducing mechanisms at the level of the units of the network, both pre-and post-nonlinearities.

The main thrust of the paper is to establish that the units prior distribution becomes increasingly heavy-tailed with depth. We show that first layer units are Gaussian, second layer units are sub-Exponential, and we introduce sub-Weibull distributions to characterize the deeper layers units. Bayesian neural networks with Gaussian priors are well known to induce the weight decay penalty on the weights. In contrast, our result indicates a more elaborate regularization scheme at the level of the units, ranging from convex penalties for the first two layers -weight decay for the first and Lasso for the second -to non convex penalties for deeper layers. Thus, despite weight decay does not allow for the weights to be set exactly to zero, sparse solutions tend to be selected for the units from the second layer onward. This result provides new theoretical insight on deep Bayesian neural networks, underpinning their natural shrinkage properties and practical potential.

Introduction

Neural networks (NNs) [START_REF] Bishop | Neural Networks for Pattern Recognition[END_REF], and their deep extensions [START_REF] Goodfellow | Deep learning[END_REF], have largely been used in many research areas such as image analysis [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF], signal processing [START_REF] Graves | Speech recognition with deep recurrent neural networks[END_REF], or reinforcement learning [START_REF] Silver | Mastering the game of go with deep neural networks and tree search[END_REF], just to name a few. These performances have greatly strengthened the line of research that aims at better understanding the driving mechanisms behind the effectiveness of deep neural networks. One important aspect of this analysis that has recently gained much attention is the study of distributional properties of the NNs through Bayesian inference.

Bayesian approaches investigate models by assuming a prior distribution on their parameters. Bayesian machine learning refers to extending standard machine learning approaches with posterior inference, a line of research pioneered by the works [START_REF] Neal | Bayesian training of backpropagation networks by the hybrid Monte Carlo method[END_REF]; [START_REF] Mackay | A practical Bayesian framework for backpropagation networks[END_REF] on Bayesian neural networks which now extends to a broad class of models, including Bayesian GAN [START_REF] Saatci | Bayesian GAN[END_REF]. See [START_REF] Polson | Deep learning: A Bayesian perspective[END_REF] for a review. The interest of the Bayesian approach to NNs is at least twofold. First, it offers a principled approach for modeling uncertainty of the training procedure, which is a limitation of standard NNs which only provide point estimates. A second main asset of Bayesian models is that they represent regularized versions of their classical counterparts. For instance, mode a posteriori (MAP) estimation of a Bayesian regression model with double exponential (Laplace) prior is equivalent to Lasso regression [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF], while a Gaussian prior leads to ridge regression. When it comes to neural networks, the regularization mechanism is also well appreciated in the literature, since neural networks traditionally suffer from overparameterization, resulting in overfitting.

Central in the field of regularization techniques is the weight decay penalty [START_REF] Krogh | A simple weight decay can improve generalization[END_REF], which is equivalent to MAP estimation of a Bayesian neural network with independent Gaussian priors on the weights. [START_REF] Srivastava | Dropout: A simple way to prevent neural networks from overfitting[END_REF] have suggested dropout as a regularization method in which neurons are randomly turned off. [START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning[END_REF] proved that the neural network trained with dropout is equivalent to a probabilistic model, i.e. a deep Gaussian process [START_REF] Damianou | Deep Gaussian processes[END_REF]. It leads to the consideration of such neural networks as Bayesian models.

This study is devoted to the investigation of hidden units prior distributions in Bayesian neural networks under assumption of independent Gaussian weights. We first describe a fully connected neural network architecture as illustrated in Figure 1. Given an input x ∈ R N , the -th hidden layer unit activations are defined as

g ( ) (x) = W ( ) h ( -1) (x), h ( ) (x) = φ(g ( ) (x)), (1)
where W ( ) is a weight matrix including the bias vector. A nonlinear activation function φ : R → R is applied element-wise, which is called nonlinearity, g ( ) = g ( ) (x) is a vector of pre-nonlinearities, and h ( ) = h ( ) (x) is a vector of post-nonlinearities. When we refer to either pre-or post-nonlinearities, we will use the notation U ( ) . Figure 1: Neural network architecture and characterization of the -layer units prior distribution as sub-Weibull distribution with tail parameter /2 (see Definition 4.1).

Contributions

In this paper, we extend the theoretical understanding of feedforward fully connected neural networks by studying prior distributions at the units level, under the assumption of independent and normally distributed weights. Our contributions are threefold:

(i) We define the notion of sub-Weibull distributions (Definition 4.1), which are characterized by tails lighter than (or equally light as) Weibull distributions; in the same way as sub-Gaussian or sub-Exponential distributions correspond to distributions with tails lighter than Gaussian and Exponential distributions, respectively. Sub-Weibull distributions are parameterized by a positive tail index θ and equivalent to sub-Gaussian for θ = 1/2 and sub-Exponential for θ = 1. We provide a moments characterization of the sub-Weibull property in Proposition 4.1.

(ii) As our main contribution, we prove in Theorem 3.1 that under some conditions on the nonlinear function φ, a Gaussian prior on the weights induces a sub-Weibull distribution on the units (both pre-and post-nonlinearities) with optimal tail parameter θ = /2 (see Figure 1). The condition on φ essentially imposes that φ strikes to +∞ or -∞ for large absolute values of the argument, as ReLU does. In the case of bounded support φ, like sigmoid or tanh, the units are bounded, making them de facto sub-Gaussian2 .

(iii) Lastly, we offer an interpretation of the main result from a sparsity-inducing viewpoint in Section 2. Heavy-tailed priors are known to induce a sparse model representation, such as the Lasso in a regression problem. We extrapolate this finding to the setting of a deep Bayesian neural network (BNN), showing that the units tend to be more sparsely represented as layers become deeper.

Related work

Studying the distributional behaviour of feedforward networks has been a fruitful avenue for understanding these models, as pioneered by the works of Radford Neal [START_REF] Neal | Bayesian training of backpropagation networks by the hybrid Monte Carlo method[END_REF][START_REF] Neal | Bayesian learning for neural networks[END_REF] and David MacKay [START_REF] Mackay | A practical Bayesian framework for backpropagation networks[END_REF]. The first results in the field addressed the limiting setting when the number of units per layer tends to infinity, also called the wide regime. [START_REF] Neal | Bayesian learning for neural networks[END_REF] proved that a single hidden layer neural network with normally distributed weights tends in distribution in the wide limit either to a Gaussian process [START_REF] Rasmussen | Gaussian Processes for Machine Learning[END_REF] or to an α-stable process, depending on how the prior variance on the weights is rescaled. In recent works, Matthews et al. (2018a), or its extended version [START_REF] Matthews | Gaussian process behaviour in wide deep neural networks[END_REF], and [START_REF] Lee | Deep neural networks as Gaussian processes[END_REF] extend the result of Neal to morethan-one-layer neural networks: when the number of hidden units grows to infinity, deep neural networks (DNNs) also tend in distribution to the Gaussian process, under the assumption of Gaussian weights for properly rescaled prior variances. For the rectified linear unit (ReLU) activation function, the Gaussian process covariance function is obtained analytically [START_REF] Cho | Kernel methods for deep learning[END_REF]. For other nonlinear activation functions, [START_REF] Lee | Deep neural networks as Gaussian processes[END_REF] use a numerical approximation algorithm.

Various distributional properties are also studied in neural networks regularization methods. The dropout technique [START_REF] Srivastava | Dropout: A simple way to prevent neural networks from overfitting[END_REF] was reinterpreted as a form of approximate Bayesian variational inference [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF][START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning[END_REF]. [START_REF] Gal | Dropout as a Bayesian approximation: Representing model uncertainty in deep learning[END_REF] built a connection between dropout and the Gaussian process. [START_REF] Kingma | Variational dropout and the local reparameterization trick[END_REF] proposed a way to interpret Gaussian dropout. They suggested variational dropout where each weight of a model has its individual dropout rate. Sparse variational dropout from [START_REF] Molchanov | Variational dropout sparsifies deep neural networks[END_REF] extends variational dropout to all possible values of dropout rates and leads to a sparse solution.

The approximate posterior is chosen to factorize either over rows or individual entries of the weight matrices. The prior usually factorizes in the same way.

The choice of the prior and its interaction with the approximating posterior family are studied in [START_REF] Hron | Variational Bayesian dropout: pitfalls and fixes[END_REF]. Performing dropout can be used as a Bayesian approximation. However, as noted by [START_REF] Duvenaud | Avoiding pathologies in very deep networks[END_REF], dropout has no regularization effect on infinitely-wide hidden layers.

The recent work by [START_REF] Bibi | Analytic expressions for probabilistic moments of pl-dnn with gaussian input[END_REF] provides the expression of the first two moments of the output units of a one layer neural network. Obtaining the moments is a first step to characterizing a whole distribution, however the methodology of [START_REF] Bibi | Analytic expressions for probabilistic moments of pl-dnn with gaussian input[END_REF] is limited to the first two moments and to one layer neural networks, while we address the problem in more generality for deep neural networks.

In the remainder of the paper, we present our three main contributions in reverse but perhaps more pedagogical order, that is starting with intuitions and interpretation (iii), then moving to theoretical results (ii) and ending up with the necessary statistical background and definitions (i). More specifically, Section 2 illustrates shrinkage and penalization techniques, providing an interpretation for our main contribution, Theorem 3.1, that is stated in Section 3. Section 4 introduces the sub-Weibull distribution family based on its tail behavior and proves a moment-based characterization. Conclusions and future work are reported in Section 5, while proofs and additional technical results are deferred to the Supplementary material.

2 Sparsity-inducing prior on the units

Short digest on penalized estimation

Our main theoretical contribution, Theorem 3.1, characterizes the marginal prior distribution of the network units as follows: when the depth increases, the distribution becomes more heavy-tailed, as will be precised in the next section. First, in this section, we provide an interpretation of the result in terms of sparsityinducing mechanism at the level of the units. To this aim, a short reminder about shrinkage methods and penalization is presented.

The shrinkage idea is probably best illustrated on the simple linear regression model, where the aim of shrinkage is to improve prediction accuracy by shrinking, or even putting exactly to zero, some coefficients in the regression. Under these circumstances, inference is also more interpretable since, by reducing the number of coefficients effectively used in the model, it is possible to grasp its salient features. Shrinking is performed by imposing a penalty on the size of the coefficients, which is equivalent to allowing for a given budget on their size. Denote the regression parameter by β ∈ R p , the regression sum-of-squares by R(β), and the penalty by λL(β), where L is some norm on R p and λ some positive tuning parameter. Then, the two formulations of the regularized problem min β∈R p R(β) + λL(β), and min β∈R p R(β) subject to L(β) ≤ t, are equivalent, with some one-to-one correspondence between λ and t, and are respectively termed the penalty and the constraint formulation. This latter formulation provides an interesting geometrical intuition of the shrinkage mechanism: the constraint L(β) ≤ t reads as imposing a total budget of t for the parameter size in terms of the norm L. If the ordinary least squares estimator βols lives in the L-ball with surface L(β) = t, then there is no effect on the estimation. In contrast, when βols is outside the ball, then the intersection of the lowest level curve of the sum-of-squares R(β) with the L-ball defines the penalized estimator.

The choice of the L norm has considerable effects on the problem, as can be sensed geometrically. Consider for instance L q norms, with q ≥ 0. For any q > 1, the associated L q norm is differentiable and contours have a round shape without sharp angles. In that case, the penalty effect is to shrink the β coefficients towards 0. The most well-known estimator falling in this class is the ridge regression obtained with q = 2, see Figure 2 top-left panel. In contrast, for any q ∈ (0, 1], the L q norm has some non differentiable points along the axis coordinates, see Figure 2 top-right and bottom panels. Such critical points are more likely to be hit by the level curves of the sum-of-squares R(β), thus setting exactly to zero some of the parameters. A very successful approach in this class is the Lasso obtained with q = 1. Note that the problem is computationally much easier in the convex situation which occurs only for q ≥ 1.

Figure 2: L 2/ -norm unit balls (in dimension 2) for layers = 1, 2, 3 and 10.

MAP on weights W is weight decay

These penalized methods have a simple Bayesian counterpart in the form of the mode a posteriori (MAP) estimator. In this context, the objective function R is the negative log-likelihood, while the penalty L is the negative log-prior. The objective function takes on the form of sum-of-squared errors for regression under Gaussian errors, and of cross-entropy for classification.

For neural networks, it is well-known that an independent Gaussian prior on the weights

π(W ) = L =1 i,j e -1 2 (W ( ) i,j ) 2
, is equivalent to the weight decay penalty, also known as ridge regression:

L(W ) = L =1 i,j (W ( ) i,j ) 2 = W 2 2 ,
where products and sums involving i and j above are over 1 ≤ i ≤ H -1 and 1 ≤ j ≤ H , H 0 and H L representing respectively the input and output dimensions.

MAP on units U induces sparsity

Now moving the point of view from weights to units leads to a radically different shrinkage effect. Let U ( ) m denote the m-th unit of the -th layer (either preor post-nonlinearity). We prove in Theorem 3.1 that conditional on the input x, a Gaussian prior on the weights translates into some prior on the units U

( ) m
that is marginally sub-Weibull with optimal tail index θ = /2. This means that the tails of U ( )

m satisfy P(|U ( ) m | ≥ u) ≤ exp -u 2/ /K 1 for all u ≥ 0, (2)
for some positive constant K 1 . The exponent of u in the exponential term above is optimal in the sense that Equation ( 2) is not satisfied with some parameter θ smaller than /2. Thus, the marginal density of U

( ) m on R is approximately proportional to π ( ) m (u) ≈ e -u 2/ /K1 .
The joint prior distribution for all the units U = (U ( ) m ) 1≤ ≤L,1≤m≤H can be expressed from all the marginal distributions by Sklar's representation theorem as

π(U ) = L =1 H m=1 π ( ) m (U ( ) m ) C(F (U )), (3)
where C represents the copula of U (which characterizes all the dependence between the units) while F denotes its cumulative distribution function. The penalty incurred by such a prior distribution is obtained as the negative log-prior,

L(U ) = - L =1 H m=1 log π ( ) m (U ( ) m ) -log C(F (U )), ≈ L =1 H m=1 |U ( ) m | 2/ -log C(F (U )), ≈ U (1) 2 2 + U (2) 1 1 + • • • + U (L) 2/L 2/L -log C(F (U )). ( 4
)
The first L terms in (4) indicate that some shrinkage operates at every layer of the network, with a penalty term that takes the form of the L 2/ norm. Thus, the deeper the layer, the stronger the sparsity at the level of the units, as summarized in Table 1.

Layer Penalty on W Penalty on U

1 W (1) 2 2 , L 2 U (1) 2 2 L 2 (weight decay) 2 W (2) 2 2 , L 2 U (2) L 1 (Lasso) W ( ) 2 2 , L 2 U ( ) 2/ 2/ L 2/

Bayesian neural networks have heavy-tailed deep units

The deep learning approach uses stochastic gradient descent and error back-propagation in order to fit the network parameters (W ( ) ) 1≤ ≤L , where iterates over all network layers. In the Bayesian approach, the parameters are random variables described by probability distributions.

Assumptions on neural network

We assume a prior distribution on the model parameters, that are the weights W . In particular, let all weights (including biases) be independent and have zero-mean normal distribution

W ( ) i,j ∼ N (0, σ 2 w ), (5)
for all 1 ≤ ≤ L, 1 ≤ i ≤ H -1 and 1 ≤ j ≤ H . Given some input x, such prior distribution induces by forward propagation 1 a prior distribution on the pre-nonlinearities and post-nonlinearities, whose tail properties are the focus of this section. To this aim, the nonlinearity φ is required to span at least half of the real line as follows. We introduce an extended version of the nonlinearity assumption from Matthews et al. (2018b):

Definition 3.1 (Extended envelope property for nonlinearities). A nonlinearity φ : R → R is said to obey the extended envelope property if there exist

c 1 , c 2 , d 1 , d 2 ≥ 0 such that the following inequalities hold |φ(u)| ≥ c 1 + d 1 |u| for all u ∈ R + or u ∈ R -, |φ(u)| ≤ c 2 + d 2 |u| for all u ∈ R. ( 6
)
The interpretation of this property is that φ must shoot to infinity at least in one direction (R + or R -, at least linearly (first line of ( 6)), and also at most linearly (second line of ( 6)). Of course, compactly supported nonlinearities such as sigmoid and tanh do not satisfy the extended envelope property but the majority of other nonlinearities do, including ReLU, ELU, PReLU, and SeLU.

Lemma 3.1. Let a nonlinearity φ : R → R obey the extended envelope property. Then for any symmetric random variable X the following asymptotic equivalence 3 holds

E[φ(X) k ] E[X k ], for k → ∞. ( 7
)
The proof is deferred to the Supplementary material.

Main theorem

This section postulates the rigorous result with a proof sketch. In Supplementary material one can find proofs of intermediate lemmas and a covariance theorem which states the non-negative covariance between post-nonlinearities.

Theorem 3.1 (Sub-Weibull units). Consider a feedforward Bayesian neural network with Gaussian priors (5) with nonlinearity φ satisfying the extended envelope condition of Definition 3.1. Then conditional on the input x, the marginal prior distribution induced by forward propagation (1) on any unit (pre-or postnonlinearity) of the -th hidden layer is sub-Weibull with optimal tail parameter θ = /2. That is for any 1 ≤ ≤ L, and for any

1 ≤ m ≤ H , U ( ) m ∼ subW( /2),
where a subW distribution is defined in Definition 4.1, and U

( )
m is either a pre-nonlinearity g ( )

m or a post- nonlinearity h ( ) m .
Proof. The idea is to prove by induction with respect to hidden layer depth that pre-and postnonlinearities satisfy the asymptotic moment equivalence

g ( ) k k /2 and h ( ) k k /2 .
The statement of the theorem then follows by the moment characterization of optimal sub-Weibull tail coefficient in Proposition 4.3.

According to Lemma 4.1, centering does not harm tails properties, then, for simplicity, we consider zero-mean distributions W ( ) i,j ∼ N (0, σ 2 w ). Base step: consider the distribution of the first hidden layer pre-nonlinearity g ( = 1). Since weights W m follow normal distribution and x is a feature vector, then each hidden unit W T m x follow also normal distribution g = W T m x ∼ N (0, σ 2 w x 2 ). Then, for normal zero-mean variable g, having variance σ 2 = σ 2 w x 2 , it holds the equality in sub-Gaussian property with variance proxy equals to normal distribution variance and from Lemma B.1:

g k √ k.
As activation function φ obeys extended envelope property, according to Lemma B.2, nonlinearity moments are asymptotic equivalent to symmetric variable moments φ(g

) k g k ∼ √ k.
It implies that first hidden layer post-nonlinearities h have sub-Gaussian distribution or sub-Weibull with tail parameter θ = 1/2 (Definition 4.1).

Inductive step: show that if the statement holds for -1, then it also holds for .

Suppose the post-nonlinearity of ( -1)-th hidden layer satisfies the moment condition. Hidden units satisfy the non-negative covariance theorem (Theorem C.1):

Cov h ( -1) s , h( -1) t ≥ 0, for any s, t ∈ N.

Let the number of hidden units in ( -1)-th layer equals to H. Then according to Lemma B.3, under assumption of zero-mean Gaussian weights, pre-nonlinearities of -th hidden layer g

( ) = H i=1 W ( -1) m,i h ( -1) i also sat- isfy the moment condition, but with θ = /2 g ( ) k k /2 .
Using Lemma B.2 from Supplementary material, one can show that post-nonlinearities h ( ) satisfy the same moment condition as pre-nonlinearities g ( ) . This finishes the proof.

We illustrate the result of Theorem 3.1 in Figure 3 which represents the first three hidden layers prenonlinearity marginal distributions (top panel). These densities are obtained as kernel density estimators from a sample of size 10 5 from the prior on the prenonlinearities, which is itself obtained by sampling 10 5 sets of weights W from the Gaussian prior (5) and forward propagation via (1). The three hidden layers of neural network have H 1 = 25, H 2 = 24 and H 3 = 4 hidden units, respectively. Being a linear combination involving symmetric weights W , pre-nonlinearities g are also symmetric, thus we visualize only their positive part. The input vector x ∈ R 50 is sampled from a standard normal distribution once for all at the start. The nonlinearity φ is the ReLU function. The prior distribution of post-nonlinearities has a Dirac mass at zero with a coefficient of 1/2 and they are no more symmetric. But the post-nonlinearity prior distribution tails remain the same as of pre-nonlinearities on R + , and is represented on the bottom panel of Fig-

Convolutional neural networks

Convolutional neural networks [START_REF] Goodfellow | Deep learning[END_REF] are a particular kind of neural network for processing data that has a known grid-like topology, which allows to encode certain properties into the architecture. These then make the forward function more efficient to implement and vastly reduce the amount of parameters in the neural network. Neurons in such networks are arranged in three dimensions: width, height and depth. There are three main types of layers that can be concatenated in CNN architectures: convolutional layer, pooling layer, and fully-connected layer (exactly as seen in standard Neural Networks). Convolutional layer computes dot products between a region in the inputs and its weights. Therefore, each region can be considered as a particular case of fully-connected layer. Pooling layer is performed to control overfitting and computations in deep architectures. The pooling layer operates independently on every depth slice of the input and residues it spatially.

The most commonly functions used in pooling layers are max pooling and average pooling. Proposition 3.1. The operations: 1. max pooling and 2. averaging do not modify the optimal tail parameter θ of sub-Weibull family. Consequently, the result of Theorem 3.1 carries over to Convolutional neural networks.

Proof. Let X i ∼ subW(θ) for 1 ≤ i ≤ N be units from one region where pooling operation is applied. Using Definition 4.1, for all x ≥ 0 and some constant K > 0 we have

P(|X i | ≥ x) ≤ exp -x 1/θ /K for all i.
1. Max pooling operation takes the maximum element in the region. Since X i , 1 ≤ i ≤ N are the elements in one region, we want to check if the tail of max 1≤i≤N X i obeys sub-Weibull property with optimal tail parameter equals to θ. Using probability properties, we get

P( max 1≤i≤N X i ≥ x) = P( ∪ 1≤i≤N X i ≥ x) ≤ N i=1 P(X i ≥ x).
Taking into account that the variables X i belong to one sub-Weibull family for all 1 ≤ i ≤ N , the distribution tail has the following form

P( max 1≤i≤N X i ≥ x) ≤ 2N exp -x 1/θ /K .
The constant 2 is appeared through distribution symmetry. Then, for K 1 = K/ log(2N ), the proposition for max pooling layer holds.

2. Apply averaging on the unit region:

P 1 N N i=1 X i ≥ x = ∩ 1≤i≤N P X i /N ≥ x .
Substituting the sub-Weibull property, we obtain the same upper bound

P 1 N N i=1 X i ≥ x ≤ exp -x 1/θ /K .
Summarizing and division by a constant does not influence the distribution tail, yielding the proposition result regarding the averaging operation.

Corollary 3.1. Consider a convolutional neural network containing convolutional, pooling and fully-connected layers under assumptions from Section 3.1. Then a unit of -th hidden layer has sub-Weibull distribution with optimal tail parameter θ = /2, where is the number of convolutional and fully-connected layers.

Proof. Proposition 3.1 implies that the pooling layer keeps the tail parameter. From discussion at the beginning of the section, the result of Theorem 3.1 is also applied to CNNs where the depth is considered as the number of convolutional and fully-connected layers.

Sub-Weibull distributions

Let X be a random variable. The following proposition states different equivalent distribution properties, such as tail decay and the growth of moments. The proof of this result shows how to transform one type of information about random variables into another.

Proposition 4.1 (Equivalent properties). Let X be a random variable. Then the following properties are equivalent; the parameters K i > 0 appearing in these properties differ from each other by at most an absolute constant factor.

1. The tails of X satisfy

P(|X| ≥ x) ≤ exp -x 1/θ /K 1 for all x ≥ 0.
2. The moments of X satisfy

X k = E[|X| k ] 1/k ≤ K 2 k θ for all k ≥ 1.
It is a short version of Proposition A.1 which states two additional equivalent properties in terms of the moment generating function of the random variable X. Proposition A.1 is deferred to and proved in Supplementary material.

Definition 4.1 (Sub-Weibull random variable). A random variable X that satisfies one of the equivalent properties of Proposition 4.1 is called a sub-Weibull random variable with the tail parameter θ, which is denoted by X ∼ subW(θ).

Informally, the tails of a subW(θ) distribution are dominated by (i.e. decay at least as fast as) the tails of a Weibull variable with the shape parameter equal to 1/θ [START_REF] Rinne | The Weibull distribution: a handbook[END_REF]. The larger tail parameter θ, the heavier the tails of the sub-Weibull distribution.

Sub-Gaussian and sub-Exponential variables, which are commonly used, are special cases of sub-Weibull random variables with tail parameter θ = 1/2 and θ = 1, respectively (see Table 2).

Distribution Tail Moments

Sub-Gaussian

F (x) ≤ e -λx 2 X k ≤ C √ k Sub-Exponential F (x) ≤ e -λx X k ≤ Ck Sub-Weibull F (x) ≤ e -λx 1/θ X k ≤ Ck θ
Table 2: Sub-Gaussian, sub-Exponential and sub-Weibull distributions comparison in terms of tail F (x) = P (X ≥ x) and moment condition, with λ and C some positive constants. The first two are a special case of the last with θ = 1/2 and θ = 1 respectively.

Proposition 4.2 (Inclusion). Let θ 1 and θ 2 such that 0 < θ 1 < θ 2 be tail proxy parameters for some sub-Weibull distributed variables. Then the following inclusion holds subW(θ 1 ) ⊂ subW(θ 2 ).

Proof. For X ∼ subW(θ 1 ), it holds that

X k ≤ K 2 k θ1 . Since k θ1 ≤ k θ2 for all k ≥ 1, this yields X k ≤ K 2 k θ2 , which by definition implies X subW(θ 2 ).
The following proposition is key in establishing that neural network units of layer are subW( /2), where /2 is optimal.

Proposition 4.3 (Optimal sub-Weibull tail coefficient and moment condition). Let θ > 0 and let X be a random variable satisfying the following asymptotic equivalence on moments3

X k k θ .
Then X is sub-Weibull distributed with optimal tail parameter θ, in the sense that for any θ < θ, X is not sub-Weibull with tail parameter θ .

Proof. Since X satisfies Condition 2 of Proposition 4.1, X ∼ subW(θ). Let θ < θ. Since X k k θ , there does not exist any constant K 2 such that X k ≤ K 2 k θ , so X is not sub-Weibull with tail proxy parameter θ .

It is typically assumed that the random variable X has zero mean. If this is not the case, we can always center X by subtracting the mean. Let us prove that variable centering does not change the tail parameter of sub-Weibull distribution it follows.

Lemma 4.1 (Centered variables). Centering does not harm tail properties. In particular, random variables X and (X -E[X]) belong to the same sub-Weibull family, i.e. with the same optimal tail proxy parameter.

The proof can be found in Supplementary material.

Conclusion and future work

Despite the ubiquity of deep learning throughout science, medicine and engineering, the underlying theory has not kept pace with applications for deep learning in general, and for neural networks in particular. In this paper, we have extended the state of knowledge on Bayesian neural networks by providing a characterization of the marginal prior distribution of the units. We proved that they are heavier-tailed as depth increases, and interpreted this result as a sparsity-inducing mechanism at the level of the units.

Since initialization and learning dynamics are key in modern machine learning in order to properly tune deep learning algorithms, a good implementation practice requires a proper understanding of the prior distribution at play and of the regularization it incurs.

We hope that our results will open avenues for further research. Firstly, Theorem 3.1 regards the marginal prior distribution of the units, while a full characterization of the joint distribution of all units U remains an open question. More specifically, a precise description of the copula defined in Equation (3) would provide valuable information about the dependence between the units, and also about the precise geometrical structure of the balls induced by that penalty. Secondly, the interpretation of our result (Section 2) is concerned with the mode a posteriori of the units, which is a point estimator. One of the benefits of the Bayesian approach to neural networks lies in its ability to provide a principled approach to uncertainty quantification, so that an interpretation of our result in terms of the full posterior distribution would be very appealing. Lastly, the practical potentialities of our results are many: to delve into Bayesian deep neural networks distributional properties and better comprehend their sparsifying mechanisms will contribute to design and understand regularization strategies to avoid overfitting and improve generalization. Future work will explore these directions.

A Equivalent sub-Weibull distribution properties

Proposition A.1 (Sub-Weibull distribution). Let X be a random variable. Then the following properties are equivalent; the parameters K i > 0 appearing in these properties differ from each other by at most an absolute constant factor.

1. The tails of X satisfy

P(|X| ≥ x) ≤ 2 exp -x 1/θ /K 1 for all x ≥ 0.
2. The moments of X satisfy

X k = E[|X| k ] 1/k ≤ K 2 k θ for all k ≥ 1. 3. The MGF of X 1/θ satisfies E exp λ 1/θ X 1/θ ≤ K 2 exp(K 1/θ 3 λ 1/θ )
for all λ such that |λ| ≤ 1 K3 .

4. The MGF of X 1/θ is bounded at some point, namely

E exp X 1/θ /K 4 ≤ 2.
Proof. 1 ⇒ 2. Assume property 1 holds. Applying the integral identity for |X| k , we obtain

E |X| k = ∞ 0 P |X| k > x dx = ∞ 0 P |X| > x 1/k dx ≤ ∞ 0 2 -x 1/(kθ) /K 1 dx = 2K kθ 1 kθ ∞ 0 e -u u kθ-1 du = 2K kθ 1 kθ Γ (kθ) ∼ K kθ 1 kθ (kθ -1) kθ-1 ∼ (K 1 kθ) kθ .
Taking the k-th root of the expression above yields property 2

X k (K 1 θ) θ k θ ≤ K 2 k θ , with K 2 = (K 1 θ) θ .
2 ⇒ 3. Assume property 2 holds. Recalling the Taylor series expansion of the exponential function, we obtain

E exp λ 1/θ X 1/θ = E 1 + ∞ k=1 (λ 1/θ |X| 1/θ ) k k! = 1 + ∞ k=1 λ k/θ E[|X| k/θ ] k! . Property 2 guarantees that E[|X| k ] ≤ K 2 k k/θ and E[|X| k/θ ] ≤ K 2 (k/θ) k for some K 2 . Stirling s approx- imation yields k! ≥ (k/e) k
. Substituting these two bounds, we get

E exp λ 1/θ X 1/θ ≤ ∞ k=1 λ k/θ K 2 (k/θ) k (k/e) k = ∞ k=0 K 2 (eλ 1/θ /θ) k = K 2 1 -eλ 1/θ /θ ,
provided that eλ 1/θ /θ < 1, in which case the geometric series above converges. To bound this quantity further, we can use the numeric inequality

1 1-x ≤ e 2x , which is valid for x ∈ [0, 1/2]. It follows that E exp λ 1/θ X 1/θ ≤ K 2 exp 2eλ 1/θ /θ
for all λ satisfying |λ| ≤ θ 2e θ . This yields property 3 with K 3 = (2e/θ) θ .

3 ⇒ 4. Assume property 3 holds. Take λ = 1/K 4 , where K 4 ≥ K 3 /(ln 2-ln K 2 ) θ . This yields property 4. 4 ⇒ 1. Assume property 4 holds. We may assume that K 4 = 1. Then, by Markov's inequality and property 3, we obtain

P |X| > x = P e |X| 1/θ > e x 1/θ ≤ E[e |X| 1/θ ] e x 1/θ ≤ 2e -x 1/θ /K1 .
This proves property 1 with K 1 = 1.

Remark A.1. The constant 2 that appears in some properties in Proposition A.1 does not have any special meaning. It is chosen for simplicity and can be replaced by other absolute constants.

B Intermediate lemmas

Introduce the definition of asymptotic equivalence between numeric sequences:

D > 0 such that d ≤ a k b k ≤ D, for all k ∈ N. (8)
Lemma B.1 (Gaussian moments). Let X be a normal random variable such that X ∼ N (0, σ 2 ), then the following asymptotic equivalence holds

X k √ k.
Proof. The moments of central normal absolute random variable |X| are equal to

E[|X| k ] = R |x| k p(x) dx = 2 ∞ 0 x k p(x) dx = 1 √ π σ k 2 k/2 Γ k + 1 2 . ( 9
)
We have the expression for the Gamma function

Γ(z) = 2π z z e z 1 + 1 12z + o 1 z . (10)
Substituting ( 10) into the central normal absolute moment (9), we obtain

E[|X| k ] = 1 √ π σ k 2 k/2 4π k + 1 k + 1 2e (k+1)/2 • 1 + 1 6(k + 1) + o 1 k = 2σ k √ 2e k + 1 e k/2 1 + 1 6(k + 1) + o 1 k .
Then the roots of absolute moments can be written in the form of

X k = σ e 1/(2k) k + 1 e 1 + 1 6(k + 1) + o 1 k 1/k = σ e √ k + 1 e 1/(2k) 1 + 1 6(k + 1)k + o 1 k 2 = σ e c k √ k + 1.
Here the coefficient c k denotes

c k = 1 e 1/(2k) 1 + 1 6(k + 1)k + o 1 k 2 → 1,
with k → ∞. Thus, asymptotic equivalence holds

X k √ k + 1 √ k.
Lemma B.2 (Nonlinearity moments). Let X be symmetric random variable, φ(x) be nonlinear function satisfying extended envelope property, then the following asymptotic equivalence holds φ(X) k X k .

Proof. According to extended envelope property,

E[φ(X) k ] E[X k ].
That means there exist constants d and D such that for all k ∈ N it holds

d ≤ E[φ(X) k ] E[|X| k ] ≤ D.
Observing that

d ≤ d 1/k ≤ φ(X) k X k ≤ D 1/k ≤ D ,
the bounding constants are d = min{1, d}, D = max{1, D}. It yields asymptotic equivalence

φ(X) k X k .
The lemma is proved.

Lemma B.3 (Multiplication moments). Let W and X be independent random variables such that W ∼ N (0, σ 2 ) and for some p > 0 it holds

X k k p . ( 11
)
Let W i be independent copies of W , and X i be copies of X, i = 1, . . . , H with non-negative covariance between moments of copies Cov X s i , X t j ≥ 0, for i = j, s, t ∈ N.

(12)

Then we have the following asymptotic equivalence

H i=1 W i X i k k p+1/2 . ( 13
)
Proof. Let us proof the statement, using mathematical induction.

Base case: show that the statement is true for H = 1. For independent variables W and X, we have

W X k = E[|W X| k ] 1/k = E[|W | k ]E[|X| k ] 1/k = W k X k . ( 14
)
Since the random variable W follows Gaussian distribution, then Lemma B.1 implies

W k √ k. ( 15
)
Substituting assumption (11) and weight norm asymptotic equivalence (15) into ( 14) leads to the desired asymptotic equivalence (13) in case of H = 1.

Inductive step: show that if for H = n -1 the statement holds, then for H = n it also holds.

Suppose for H = n -1 we have

n-1 i=1 W i X i k k p+1/2 . ( 16
)
Then, according to the covariance assumption ( 12), for

H = n we get n i=1 W i X i k k = n-1 i=1 W i X i + W n X n k k (17) ≥ k j=0 C j k n-1 i=1 W i X i j j W n X n k-j k-j . ( 18
)
Using the equivalence definition (Def. B.1), from the induction assumption ( 16) for all j = 0, . . . , k there exists absolute constant

d 1 > 0 such that n-1 i=1 W i X i j j ≥ d 1 j p+1/2 j . (19)
Recalling previous equivalence results in the base case, there exists constant m 2 > 0 such that

W n X n k-j k-j ≥ d 2 (k -j) p+1/2 k-j . ( 20
)
Substitute obtained bounds ( 19) and ( 20) into equation (17

) with denoted d = min{d 1 , d 2 }, obtain n i=1 W i X i k k ≥ d k k j=0 C j k j j (k -j) k-j p+1/2 = d k k k(p+1/2) k j=0 C j k j k j 1 - j k k-j p+1/2 . ( 21
)
Notice the lower bound of the following expression

k j=0 C j k j k j 1 - j k k-j p+1/2 ≥ k j=0 j k j 1 - j k k-j p+1/2 ≥ 2. ( 22
)
Substituting found lower bound ( 22) into ( 21), get

n i=1 W i X i k k ≥ 2 d k k k(p+1/2) > d k k k(p+1/2) . (23)
Now prove the upper bound. For random variables Y and Z the Holder's inequality holds

Y Z 1 = E [|Y Z|] ≤ E |Y | 2 E |Z| 2 1/2 = Y Z 2 Y Z 2 .
Holder's inequality leads to the inequality for

L k norm Y X k k ≤ Y k 2k Z k 2k . ( 24
)
Obtain the upper bound of

n i=1 W i X i k k from the norm property (24) for the random variables Y = n-1 i=1 W i X i k-j and Z = (W n X n ) j n i=1 W i X i k k = n-1 i=1 W i X i + W n X n k k (25) ≤ k j=0 C j k n-1 i=1 W i X i j 2j W n X n k-j 2(k-j) . ( 26
)
From the induction assumption ( 16) for all j = 0, . . . , k there exists absolute constant D 1 > 0 such that

n-1 i=1 W i X i j 2j ≤ D 1 (2j) p+1/2 j . ( 27
)
Recalling previous equivalence results in the base case, there exists constant D 2 > 0 such that

W n X n k-j 2(k-j) ≤ D 2 2(k -j) p+1/2 k-j . (28)
Substitute obtained bounds ( 27) and ( 28) into equation (25) with denoted

D = max{D 1 , D 2 }, obtain n i=1 W i X i k k ≤ D k k j=0 C j k 2j j 2(k -j) k-j p+1/2
.

Find an upper bound for 1 -j k k-j j k j p+1/2

.

Since expressions 1 -j k and j k are less than 1, then 1 -j k k-j j k j p+1/2

< 1 holds for all natural numbers p > 0. For the sum of binomial coefficients it holds the inequality

k j=0 C j k < 2 k . So the final upper bound is n i=1 W i X i k k ≤ 2 k D k (2k) k(p+1/2) . (29)
Hence, taking the k-th root of ( 23) and ( 29), we have upper and lower bounds which imply the equivalence for H = n and the truth of inductive step

d k p+1/2 ≤ n i=1 W i X i k ≤ D k p+1/2 ,
where d = d and D = 2 p+3/2 D. Since both the base case and the inductive step have been performed, by mathematical induction the equivalence holds for all

H ∈ N H i=1 W i X i k k p+1/2 .

C Covariance theorem

Theorem C.1 (Non-negative covariance between hidden units). Consider the deep neural network described in Section 3 with assumptions from Section 3.1. The covariance between hidden units of the same layer is non-negative. Moreover, for given -th hidden layer units h ( ) and h( ) , it holds Cov h ( ) s h( ) t ≥ 0, where s, t ∈ N.

For first hidden layer = 1 there is equality for all s and t.

Proof. Recall the covariance definition for random variables X and Y

Cov [X, Y ] = E[XY ] -E[X]E[Y ]. (30)
For Gaussian random variables X and Y , Cov [X, Y ] = 0 means X and Y are independent.

The proof is based on induction with respect to the hidden layer number.

In the proof let us make notation simplifications: w m = W m and w mi = W mi for all minH . If the index m is omitted, then w is some the vectors w m , w i is i-th element of the vector w m .

1. First hidden layer. Consider the first hidden layer units h (1) and h(1) . The covariance between units is equal to zero Cov h (1) , h(1) = Cov φ(g (1) ), φ( h(1) )

= Cov φ(w (1) x), φ( w(1) x) = 0, since the weights w (1) and w(1) are from N (0, σ 2 w ) and independent. Thus, the first hidden layer units are independent and its covariance (30) equals to 0. Moreover, since h (1) and h(1) are independent, then h (1) s and h(1) t are also independent and it holds Cov h (1) s h(1) t = 0, where s, t ∈ N.

2. Next hidden layers. Assume that ( -1)-th hidden layer has H -1 hidden units, where > 1. Then -th hidden layer pre-nonlinearity is equal to

g ( ) = H -1 i=1 w ( ) i h ( -1) i . (31)
We want to prove that the covariance (30) between -th hidden layer pre-nonlinearities is non-negative. Let us show firstly the idea of the proof in the case H -1 = 1 and then briefly the proof for any finite H -1 > 1, H -1 ∈ N.

2.1 One hidden unit. In the case H -1 = 1, the covariance (30) can be written in the form of

Cov w ( ) h ( -1) s , w( ) h ( -1) t (32) = E w ( ) s E w( ) t E h ( -1) s+t -E w ( ) s E w( ) t E h ( -1) s E h ( -1) t (33)
Regarding that the weighs are zero-mean distributed, its moments are equal to zero with odd order. When s (or t) is odd, the weight moment equals to zero E w ( ) s = 0 (or of order t) and both terms in equation (32) are equal to zero. Consider even moment orders s = 2s 1 and t = 2t 1 , where s 1 , t 1 ∈ N.

For obtaining the non-negative covariance, we need to prove that

E h ( -1) 2(s1+t1) ≥ E h ( -1) 2s1 E h ( -1) 2t1 .
Since a function f (x 1 , x 2 ) = x 1 x 2 is convex for x 1 ≥ 0 and x 2 ≥ 0, then, taking x 1 = h ( -1) 2s1 and x 2 = h ( -1) 2t1 , we have

E[f (x 1 , x 2 )] = E h ( -1) 2(t1+s1) (34)
and

f (E[x 1 ], E[x 2 ]) = E h ( -1) 2s1 E h ( -1) 2t1 .
(35) According to Jensen's inequality for convex function f , the expression ( 34) is more or equal to the expression ( 35) and the condition we need (34) is satisfied.

2.1. H hidden units. Now let us consider the covariance between pre-nonlinearities (31) for H -1 = H > 1. Raise the sum in the brackets to the power

H i=1 w ( ) i h ( -1) i s = = s n H =0 C n H s w ( ) H h ( -1) H n H H-1 i=1 w ( ) i h ( -1) i s-n H = • • • = = s n H =0 s-n H n H-1 =0
h ni+mi i , E [A n1:n H ] E [B m1:m H ] = = W n1:n H Wm1:m H • • E h s-H i=1 ni 1 H i=2 h ni i E h t-H i=1 mi 1 H i=2 h mi i ,
where W n1:n H Wm1:m H is the product of weights moments

W n1:n H Wm1:m H = = E w s-H i=1 ni 1 E wt-H i=1 mi 1 H i=2 E [w ni i ] E [ wmi i ] .
For W n1:n H Wm1:m H not equal to zero, all the powers must be even: 2s 1 = s-

and

f E[x 1 ], E[x 2 ] = E H i=1 h 2si i E H i=1 h 2ti i . (38)
According to Jensen's inequality for convex function f , the expression (37) is more or equal to the expression (38) and the condition we need (36) is satisfied.

3. Post-nonlinearities.

Let show the proof for the ReLU nonlinearity.

The distribution of the -th hidden layer prenonlinearity g ( ) is the sum of symmetric distributions, which are products of Gaussian variables w ( ) and nonnegative ReLU output, i.e. ( -1)-th hidden layer postnonlinearity h ( -1) . It leads that g ( ) follows symmetric distribution and the following inequality

That means if non-negative covariance is proven for pre-nonlinearities, for post-nonlinearities it is also nonnegative. Omit the proof for the other nonlinearities with extended envelope property, since instead of precise equation (39), the asymptotic equivalence for moments will be used for positive part and for negative part -precise expectation expression which depends on certain nonlinearity.

Figure 3 :

 3 Figure 3: Illustration of the first three layers hidden units marginal prior distributions.

  Definition B.1 (Asymptotic equivalence). Two sequences a k and b k are called asymptotic equivalent and denoted as a k b k if there exist constants d > 0 and

H

  i=1 n i , 2s 2 = n 2 , . . . , 2s H = n H , 2t 1 = t -H i=1 m i , 2t 2 = m 2 , . . . , 2t H = m H . Now we need to prove E f (x 1 , x 2 ) = x 1 x 2 is convex for x 1 ≥ 0 and x 2 ≥ 0, then, taking x 1 =

Table 1 :

 1 Comparison of Bayesian neural network shrinkage effect on weights W and units U .

  C n1:n H A n1:n H where C n1:n H = C n H

					s-H i=2 ni
			• • •	C n1:n H A n1:n H
					n1=0
	=
	n1:n H
					s	H-1 i=1 C ni s-H r=i+1 nr . And the
	same way for the second bracket
	H	w( ) i h ( -1) i	t	=
	i=1
		t	t-m H
	=
	m H =0	m H-1 =0

  C m1:m H B m1:m H . C m1:m H C n1:n H E [A n1:n H B m1:m H ] -n1:n H m1:m H C m1:m H C n1:n H E [A n1:n H ] E [B m1:m H ] .For covariance being non-negative it is enough to show that the differenceE [A n1:n H B m1:m H ] -E [A n1:n H ] E [B m1:m H ]is non-negative for all the numbers n i and m i . Consider some term with numbers n 1 , . . . , n H , m 1 , . . . , m H . Since the weights are Gaussian and independent, we have the following equation, omitting the superscript for simplicity,E [A n1:n H B m1:m H ] = = W n1:n H Wm1:m H •E h

		the form of
		Cov	 	H -1	w i h ( ) ( -1) i	s	,	H -1	w	( ) i h ( -1) i	t	  =
				i=1				i=1
		=
		n1:n H m1:m H
						i=1 (ni+mi) k-H	H
						1
											i=2
	t-H i=2 mi
	• • •	C m1:m H B m1:m H =
	m1=0
	=
	m1:m H
	So the covariance in our notations can be written in

A trivial version of our main result holds.

See Definition B.1 for the asymptotic equivalence definition in Supplementary material.

Acknowledgments

We would like to thank Stéphane Girard for fruitful discussions on Weibull-like distributions, and Cédric Févotte for pointing out the potential relationship of our heavy-tail result with sparsity-inducing priors.