# Diffusion Generative Modelling for Divide-and-Conquer MCMC

Connie Trojan
Department of Mathematics and Statistics
Lancaster University
Lancaster, LA1 4YF, UK
<EMAIL_ADDRESS>
&Paul Fearnhead
Department of Mathematics and Statistics
Lancaster University
Lancaster, LA1 4YF, UK
<EMAIL_ADDRESS>
&Christopher Nemeth
Department of Mathematics and Statistics
Lancaster University
Lancaster, LA1 4YF, UK
<EMAIL_ADDRESS>

###### Abstract

Divide-and-conquer MCMC is a strategy for parallelising Markov Chain Monte
Carlo sampling by running independent samplers on disjoint subsets of a
dataset and merging their output. An ongoing challenge in the literature is to
efficiently perform this merging without imposing distributional assumptions
on the posteriors. We propose using diffusion generative modelling to fit
density approximations to the subposterior distributions. This approach
outperforms existing methods on challenging merging problems, while its
computational cost scales more efficiently to high dimensional problems than
existing density estimation approaches.

##  1 Introduction

Markov chain Monte Carlo samplers are a common numerical tool for Bayesian
inference. However, each update step of a Metropolis-Hastings MCMC sampler
requires a calculation involving the full dataset to compute the acceptance
probability. This gives each iteration of the MCMC algorithm a linear
complexity in the number of datapoints, which is computationally impractical
for very large datasets, thus motivating the use of faster approximate
samplers.

There are two general approaches to MCMC to reduce the per-iteration cost (see
Bardenet et al., 2017, for an overview). One is based on approximate MCMC
algorithms that use a subset of the data at each iteration. Such methods
include stochastic gradient Langevin dynamics (Welling and Teh, 2011) or its
extensions (Nemeth and Fearnhead, 2021). The other is a divide-and-conquer
approach (Scott et al., 2016; Neiswanger et al., 2014), where the dataset is
partitioned before inference so that MCMC chains can be run in parallel on
each subset of the data and then merged. This latter approach has several
advantages over the former, as we can trivially run MCMC algorithms on data
subsets in an embarrassingly parallel fashion, across multiple CPUs, without
incurring a communication cost between CPUs. However, these methods are less
commonly used in practice due to the challenge of reliably merging the
information from the MCMC output on different data subsets. Whilst our work is
motivated by divide-and-conquer MCMC, similar challenges of merging
information from disjoint datasets arise in federated learning (see e.g. Li et
al., 2020), where the data is naturally partitioned into subsets that cannot
be merged due to privacy or communication constraints. The approach we present
to solve the challenge of parallel MCMC could also be used to perform Bayesian
inference in this setting as well.

#### Our contribution

In this paper, we present a new approach to posterior merging that leverages
the recent advances in diffusion generative modelling, leading to a method
that makes no distributional assumptions about the posteriors, such as
Gaussianity, yet is able to scale well to complex and high dimensional
posterior distributions where current methods for divide-and-conquer MCMC tend
to perform poorly.

##  2 Background on divide-and-conquer MCMC

Divide-and-conquer MCMC methods aim to generate samples from the posterior
distribution of parameter $\theta$, given samples from MCMC chains that are
conditioned on subsets of the full dataset, $Y$. In this setting, $Y$ is
partitioned into subsets $Y^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}$,
called shards, which are often divided between multiple machines. MCMC chains
are then run in parallel targeting the subposterior distributions
$p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}$ conditioned on each shard of
the data. If the prior distribution $p(\theta)$ in each subposterior is scaled
geometrically according to the number of shards $S$, then the subposteriors
can be multiplied to obtain the full posterior distribution:

$\displaystyle p^{\text{full}}(\theta)$
$\displaystyle:=p(\theta|Y)\propto\prod_{s}p(\theta)^{\frac{1}{S}}p(Y^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}|\theta)=\prod_{s}p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(\theta)\,.$
(1)

Whilst this is a simple analytic relationship between $p^{\text{full}}$ and
the $p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}$s, it is difficult to
produce samples from $p^{\text{full}}$ given samples from each
$p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}$. Current methods to do this
come in two categories:

#### Approximately Gaussian

Since the dataset is typically very large in this setting, a natural approach
is to appeal to the Bernstein-von Mises theorem and assume each subposterior
distribution can be approximated by a Gaussian distribution. Several
approximate methods in the literature are exact under this assumption, for
example, parametric Gaussian estimation (Neiswanger et al., 2014) or methods
based on transforming the subposterior samples to obtain approximate samples
from the full posterior distribution (e.g. Scott et al., 2016; Vyner et al.,
2023). These methods are computationally efficient since no further MCMC
sampling is required but are biased when the subposterior distributions are
non-Gaussian, e.g. when they are skewed or multi-modal. However, non-
Gaussianity is still quite common in the big data setting, since individual
shards may not have enough data to form an informative posterior distribution
for high-dimensional models.

#### Non-Gaussian approaches

Methods that are exact in the general setting do exist (see Chan et al.,
2023), but are much slower than approximate methods. Methods based on non-
parametric density estimation (for example, Neiswanger et al., 2014; Nemeth
and Sherlock, 2018) are asymptotically consistent in the number of MCMC
samples and do not require the subposterior distributions to be approximately
Gaussian. However, they struggle to scale to high-dimensional problems where a
large number of samples of $\theta$ is required to approximate the posterior
density since evaluating the approximation requires a computation over all of
these MCMC samples.

##  3 Methodology

Figure 1: Annealed diffusion sampling in the mixture of Gaussians example.
Full posterior in black.

Result: Density estimates $\hat{p}_{t}(\theta,t)$ interpolating between a
Gaussian and $\hat{p}^{\text{full}}(\theta)$

Require: Subposterior MCMC samples
$\\{\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{1:n_{s}}\\}_{s=1}^{S}$
and score function evaluations
$\\{\nabla^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{1:n_{s}}\\}_{s=1}^{S}$

for  _$s=1,\ldots,S$_ do

Affine transformation:

$\qquad\mu_{s}\leftarrow\text{mean}(\\{\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{1:n_{s}}\\})$;
$V_{s}\leftarrow\text{matsqrt}(\text{cov}(\\{\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{1:n_{s}}\\}))$

$\qquad\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i}\leftarrow
V_{s}^{-1}(\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i}-\mu_{s})\quad
i\in\\{1,\ldots,n_{s}\\}$

$\qquad\nabla^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i}\leftarrow
V_{s}\nabla^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i}\quad
i\in\\{1,\ldots,n_{s}\\}$

Train neural network parameters:

$\qquad\rho\leftarrow\operatorname*{argmin}_{\rho}\;\mathbb{E}_{t,\theta_{t},\theta_{0}}[L^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(t,\theta_{t},\theta_{0};\rho)]$,

where $t\sim U(0,1)$, $\theta_{t}\sim p_{t|0}(\theta_{t}|\theta_{0})$,
$\theta_{0}\sim\text{Categorical}(\\{\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{1:n_{s}}\\})$,
and $\qquad
L^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(t,\theta_{t},\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i};\rho)=\|-\nabla_{\theta_{t}}E^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(\theta_{t},t;\rho)-\kappa_{t}\nabla_{\theta_{t}}\log
p_{t|0}(\theta_{t}|\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i})-(1-\kappa_{t})m(t)^{-1}\nabla^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i}\|_{2}^{2}$

end for

$\hat{p}_{t}(\theta,t):=\exp\left(-\sum_{s=1}^{S}E^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(V_{s}^{-1}(\theta-\mu_{s}),t;\rho^{{\scriptscriptstyle(}s{\scriptscriptstyle)}})\right)$

Algorithm 1 Diffusion posterior approximation

#### Overview

Our proposed approach to this problem involves training a neural network to
approximate the unnormalised log-density of each subposterior distribution
using diffusion generative modelling. We can then use the sum of the
approximated log-densities in an MCMC sampler to sample from the full
posterior distribution. Each diffusion model learns a sequence of densities
interpolating between the subposterior distribution at time $t=0$, and a
Gaussian approximation at $t=1$. This results in a sequence of approximations
to the full posterior that interpolates between the parametric Gaussian
approximation of Neiswanger et al. (2014) and the final non-Gaussian posterior
approximation. Figure 1 illustrates the proposed diffusion posterior
approximation algorithm on the Gaussian mixture example in Section 5.2.
Starting from a reference distribution, e.g. a Gaussian distribution, and
implicitly smoothing the density estimates across time, the diffusion model
can successfully learn complex, high-dimensional distributions. Using the
sequence of densities in an annealed MCMC sampling procedure is also helpful
for sampling from distributions where ordinary MCMC sampling struggles.
Algorithm 1 gives an overview of the diffusion merging process. We retain the
score function evaluations
$\nabla^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i}:=\nabla
p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{i})$
obtained while running the subposterior MCMC samplers in order to use a
combination of the usual denoising score matching objective (Vincent, 2011)
and target score matching (De Bortoli et al., 2024), as described in Section
3.3.

###  3.1 Diffusion generative modelling

#### Stochastic differential equations

Diffusion-based generative models (e.g. Ho et al., 2020; Song et al., 2021),
work by defining a diffusion process that starts from a data distribution
$p_{0}$ and adds noise to the data until the distribution converges to a known
Gaussian prior distribution $p_{1}$. The time reversal of this process can be
used to generate new samples from the data distribution. The noising process
is a diffusion that is initialised at time $t=0$, which is the data
distribution, and defined by an SDE of the form:

$\displaystyle\mathrm{d}X_{t}=f(X_{t},t)\mathrm{d}t+g(t)\mathrm{d}W_{t}\,.$
(2)

Here, the drift term $f$ is usually linear in $X_{t}$ and controls the mean of
the process, while the diffusion term $g$ controls the rate at which Gaussian
noise is added. These are chosen so that the process converges to a Gaussian
distribution as $t\rightarrow\infty$ regardless of the form of the true data
density $p_{0}(x)$. By scaling the coefficients appropriately, it can
typically be assumed that $X_{t}$ is approximately distributed according to
the limiting distribution at time $t=1$. The most commonly used SDE in
generative modelling is the variance-preserving (VP) SDE (Song et al., 2021),
which is an inhomogenous Ornstein-Uhlenbeck process that converges to a
standard Gaussian distribution:

$\displaystyle\mathrm{d}X_{t}=-\frac{1}{2}\beta(t)X_{t}\,\mathrm{d}t+\sqrt{\beta(t)}\,\mathrm{d}W_{t}\,,\;\beta(t)>0\,.$
(3)

If $\beta$ is chosen to be a linearly increasing function of $t$, this can be
seen as a continuous version of the discrete sequence of noising kernels used
by Ho et al. (2020), who perturb the data $x_{0}$ by sampling
$x_{i}\sim\mathcal{N}(\sqrt{1-\beta_{i}}x_{i-1},\beta_{i})$. Song et al.
(2021) showed that the time reversal of an SDE of the form in (2) has the
following form, with time now running in reverse:

$\displaystyle\mathrm{d}X_{t}=\left[f(X_{t},t)-g(t)^{2}\,\nabla\log
p_{t}(X_{t})\right]\mathrm{d}t+g(t)\mathrm{d}\tilde{W}_{t}\,.$ (4)

where $p_{t}(x)$ is the density of the noised distribution at time $t$, i.e.
the marginal density of the SDE at time $t$ when initialised at $p_{0}$.

#### Score matching objective

In general, the density $p_{t}(x)$ is unknown, but we can approximate its
score function $\nabla\log p_{t}$ using a parameteric model $\psi(x,t;\rho)$.
The parameters $\rho$ of the function $\psi$ are estimated by minimising the
denoising score matching objective (Vincent, 2011)

$\displaystyle L_{DSM}(\rho,t)$
$\displaystyle=\mathbb{E}_{p_{0}(X_{0})p_{t|0}(X_{t}|X_{0})}\left[\|\psi(x_{t},t;\rho)-\nabla_{x_{t}}\log
p_{t|0}(x_{t}|x_{0})\|_{2}^{2}\right]\,.$ (5)

Note that this uses only the transition density $p_{t|0}$ of the diffusion
process, which is simple to calculate for linear SDEs since $X_{t}-X_{0}$ has
a Gaussian distribution whose parameters can be computed from the SDE
coefficients. The function $\psi$ is usually a single time-conditional neural
network fit over all values of $t$ in $(0,1]$, so that it implicitly smooths
score estimates across time. The full training objective is a weighted average
of $L_{DSM}(\rho,t)$ across time, with $t$ uniformly sampled on $(0,1]$.

#### Engery-based models

Since the score function of a distribution determines its density up to
normalising constant, we can also use diffusion modelling for unnormalised
density estimation, by parameterising the score function estimate as the
gradient of a density function. This idea was suggested by Salimans and Ho
(2021) as a way of ensuring that the score function approximation is a
conservative vector field. This is known as an energy-based parameterisation
because we model an energy function $E(x,t;\rho)$ and approximate the
unnormalised noised density $p_{t}$ by $\exp(-E(x,t;\rho))$. Salimans and Ho
(2021) proposed the parameterisation

$\displaystyle E(x,t;\rho)$
$\displaystyle=\frac{1}{2s(t)}||x-\psi(x,t;\rho)||^{2}_{2}\,,$ (6)

where $\psi(x,t;\rho):\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a neural
network and $s(t)^{2}$ is the variance of the noising kernel $p_{t|0}$. The
gradient $-\nabla_{x}E(x,t;\rho)$ is substituted into the usual score matching
objective in training, while $-E(x,t;\rho)$ models the log-density and can be
used in MCMC sampling.

###  3.2 Reparameterised stochastic differential equations

Approximating the target posterior distribution using diffusion models can be
challenging when the noise prior is significantly offset from the target
distribution. Approximations to the target distribution trained with denoising
score matching tend to poorly estimate the location and scale of the energy
function. This could be addressed by normalising the dataset to have mean zero
and unit covariance and will greatly improve the accuracy of the energy
estimates. This normalisation is helpful for this problem in several different
ways, easing neural network training and making the noise prior distributions
in the diffusion models closer to their target distributions.

#### Transformed SDE

For any SDE with a standard Gaussian limiting distribution, this preprocessing
is equivalent to using a modified SDE that converges to a Gaussian
approximation to the target distribution, rather than a standard Gaussian. If
$X_{t}$ evolves according to the linear SDE

$\displaystyle\mathrm{d}X_{t}$
$\displaystyle=f(t)X_{t}\mathrm{d}t+g(t)\mathrm{d}W_{t}\,,$ (7)

then by the Itô formula (Øksendal, 2000) the transformed process
$\tilde{X}_{t}=AX_{t}+b$ has SDE

$\displaystyle\mathrm{d}\tilde{X}_{t}$
$\displaystyle=f(t)(\tilde{X}_{t}-b)\mathrm{d}t+Ag(t)\mathrm{d}W_{t}\,,$ (8)

for $b\in\mathbb{R}^{d}$ and $A\in\mathbb{R}^{d\times d}$. By coupling, the
limiting distribution of this SDE is $N(b+m,AS^{2}A^{\top})$, where $m$ and
$S^{2}$ are the limiting mean and covariance of the original SDE. The score
functions for the densities $p_{t}$ and $\tilde{p}_{t}$ of $X_{t}$ and
$\tilde{X}_{t},$ respectively, are related as follows:

$\displaystyle\nabla\log\tilde{p}_{t}(x)$
$\displaystyle=\nabla\log\left\\{p_{t}(A^{-1}(x-b))|A^{-1}|\right\\}=\nabla\log
p_{t}(A^{-1}(x-b))\,.$ (9)

Thus, fitting a diffusion model to dataset $\tilde{X}$ using the transformed
SDE (8) is equivalent to fitting a diffusion model to the transformed dataset
$X=A^{-1}(\tilde{X}-b)$ using the original SDE (7) and transforming the
learned score function.

#### SDE for product targets

When composing diffusion models, choosing different transformations for each
component distribution is equivalent to using a different SDE for each
diffusion model. This is not an issue, since we can still derive a sequence of
noised densities that converges to the target product density as $t\rightarrow
0$, by transforming the learned density functions
$\hat{p}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{t}$:

$\displaystyle\hat{p}_{t}(x,t)$
$\displaystyle\propto\prod_{s}\hat{p}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{t}(A^{-1}_{s}(x-b_{s}),t)\,.$
(10)

In the divide-and-conquer setting, if the limiting distribution of SDE (7) is
a standard Gaussian, we propose choosing and $b_{s}$ and $A_{s}$ so that
$b_{s}$ is the sample mean of subposterior $s$ and $A_{s}A^{\top}_{s}$ is its
sample covariance. This makes the limiting distribution of SDE (8) a Gaussian
approximation to that subposterior. The product of these limiting
distributions, which is the noise prior $\hat{p}_{1}$ for the full posterior,
is then equal to the Gaussian approximation to the full posterior suggested by
Neiswanger et al. (2014) (see Appendix C). Since the noise prior in diffusion
models must be Gaussian, this choice is in a sense optimal for both the
subposteriors and the full posterior, as it makes the noise priors as similar
as possible to their respective target distributions. This means that our
sequence of densities $\hat{p}_{t}$ interpolates between a Gaussian
approximation to the full posterior and the learned non-Gaussian
approximation. We follow Vyner et al. (2023) in choosing $A$ to be the
symmetric positive-definite square root of the sample covariance matrix $V$,
$A=U\Lambda^{\frac{1}{2}}U^{\top}$, where $U$ and $\Lambda$ are the matrices
of eigenvectors and eigenvalues, respectively, in the eigendecomposition
$V=U\Lambda U^{\top}$.

This normalisation scheme makes it simpler to fit the neural network
approximation to the density across time, since it effectively standardises
the mean and variance of its inputs, which is known to improve fitting (Huang
et al., 2023; Karras et al., 2022). When the normalised dataset $X$ is used
with the variance preserving SDE, the mean and variance of $X_{t}$ are
invariant over time. The original SDE (7) has Gaussian marginal density
$p_{0|t}$ with mean $m(t)x_{0}$ and covariance matrix $S(t)^{2}$, so that
$X_{t}$ has mean $m(t)\mathbb{E}(X_{0})$ and variance $m(t)^{2}\,{\rm
Var}(X_{0})+S(t)^{2}$.

###  3.3 Alternative score matching objectives

#### Target score matching

Denoising score matching (DSM) often struggles to approximate the score
function at low noise levels, since the variance of its score estimates
explodes as $t\rightarrow 0$. De Bortoli et al. (2024) propose an alternative
objective called target score matching (TSM) that has lower variance near time
$t=0$, which can be used when it is possible to evaluate the unnoised log-
density function $p_{0}$. The target score matching loss proposed by De
Bortoli et al. (2024) is:

$\displaystyle L_{TSM}(\theta,t)$
$\displaystyle=\mathbb{E}_{X_{0},X_{t}}\left[\|\psi(x_{t},t;\theta)-m(t)^{-1}\nabla\log
p_{0}(x_{0})\|_{2}^{2}\right]\,,$ (11)

which is designed so that estimates of the score of $p_{t}$ at $x_{t}$ are
matched to a rescaling of the unnoised score of $p_{0}$ at $x_{0}$. The
variance of Monte Carlo estimates in $L_{TSM}$ is low near $t=0$, but
increases with $t$, exploding near $t=1$ for the variance preserving SDE,
where $m(t)$ tends to 0. As such, De Bortoli et al. (2024) suggest taking a
convex combination of the regression targets of DSM and TSM, weighted in
favour of TSM near $t=0$ and of DSM near $t=1$, yielding estimates of
$\nabla\log p_{t}(x)$ that are well behaved across time. Following their
suggestion, we minimise the objective function

$\displaystyle L(\rho,t)$
$\displaystyle=\mathbb{E}_{X_{0},X_{t}}\left[\|\psi(x_{t},t;\rho)-\kappa_{t}\nabla_{x_{t}}\log
p_{t|0}(x_{t}|x_{0})-(1-\kappa_{t})m(t)^{-1}\nabla\log
p_{0}(x_{0})\|_{2}^{2}\right],$ (12)

using a uniform weighting for $t$ and combination weights,

$\displaystyle\kappa_{t}=\frac{s(t)^{2}}{s(t)^{2}+m(t)^{2}\sigma^{2}_{\mathrm{data}}}\,,$
(13)

which is optimal when the target has distribution $p_{0}\sim
N(0,\sigma^{2}_{\mathrm{data}}I_{d})$. De Bortoli et al. (2024) show that
using this combined objective results in faster convergence to a mixture of
Gaussians target distribution than using either DSM or TSM, which shows slower
convergence than DSM since approximating the score well at times closer to 1
is important for SDE sampling.

#### Divide-and-conquer target score matching

In the case of gradient-based MCMC algorithms, such as the Hamiltonian Monte
Carlo algorithm, the score function on the subposterior is calculated are
stored as part of the MCMC training. This means that no additional
subposterior evaluations are needed for training the target score matching
objective. However, in order to use target score matching with normalised
data, we must also rescale these score evaluations by $A$. This is because the
training data for the neural network is the normalised dataset
$X=A^{-1}(\tilde{X}-b)$ while the score function evaluations we have are for
the density of the unnormalised $\tilde{X}$. We have
$X_{t}=m(t)A^{-1}\tilde{X}_{0}+W$, with $W\perp\\!\\!\\!\perp\tilde{X}_{0}$,
so the new regression target in TSM is
$m(t)^{-1}A\,\nabla\log\tilde{p}_{0}(\tilde{x}_{0})$, where
$\nabla\log\tilde{p}_{0}(\tilde{x}_{0})$ will have been computed and stored
whilst running the MCMC sampler on the subposterior distribution. Using a
normalised dataset in training means that its covariance is isotropic and
therefore $\sigma^{2}_{\mathrm{data}}=1$ in weighting (13).

###  3.4 Model structure

All of our examples use a residual MLP, in line with Du et al. (2023). Details
can be found in Appendix A.5. The final output of the neural network $\psi$
was used to parameterise the energy function as follows:

$\displaystyle E(x,t;\rho)$
$\displaystyle=-\frac{1}{2(m(t)^{2}+s(t)^{2})}||x-\psi(x,s(t);\rho)||^{2}_{2}\,,$
(14)

This is based on the parameterisation proposed by Salimans and Ho (2021)
described in Equation (6), with the output scaling adjusted to prevent the
energy function from degenerating as $t\rightarrow 0$. This choice was
inspired by the similarity of this parameterisation to the Gaussian density
function – the variance of $X_{t}$ will be $m(t)^{2}\,{\rm
Var}(X_{0})+s(t)^{2}I_{d}$, which is $(m(t)^{2}+s(t)^{2})I_{d}$ if the data is
normalised as suggested in Section 3.2. This means that a constant output of
$0$ will match the true energy function for a Gaussian target when a
normalised dataset is used. Regardless of the target distribution, $E$ will
tend to the energy function for the noise prior as $t\rightarrow\infty$ as
long as the output of $\psi$ tends to 0.

###  3.5 Sampling from compositions of diffusion models

In order to use diffusion modelling in divide-and-conquer MCMC, we need to be
able to combine different models to generate samples from the product of their
target distributions. A naive approach would be to add the component score
functions $\nabla\log p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{t}$
together to obtain the score of the product of the
$p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{t}$. At time $t=0,$ this is
exactly equal to the score of the target product distribution
$p^{\text{full}}_{t}$. However, this relationship does not hold for $t>0$
since we cannot interchange noising with multiplication of densities. Indeed,
simply substituting the score sum into the reverse SDE fails to generate
samples from the correct distribution. This fails even for Gaussian targets,
where for $t>0$, the score sum estimate corresponds to a Gaussian distribution
with a different mean and variance to $p^{\text{full}}_{t}$ (see Appendix D).

As we have an enery-based estimate for each noised sub-posterior, our solution
is to use an annealed MCMC sampling procedure (Geffner et al., 2023; Du et
al., 2023), since the sequence of densities, $\hat{p}_{t}$, obtained by
multiplying the noised component densities together interpolates smoothly
between a tractable noise prior and the target product distribution. We can
then transport samples from the new prior to the target by starting with a
sample of size $N$ from the prior at $t=1$ and then iteratively using a fixed
number of unadjusted MCMC updates to target $\hat{p}_{t}$ for a sequence of
predetermined timepoints $1>t_{2}>\ldots>t_{n}=0$ (see Appendix A.3). Using an
energy based parameterisation enables the use of a Metropolis-Hastings
adjusted sampler in this procedure, which has been shown to improve results
for compositional generation (Du et al., 2023).

Using the parameterisation described in Section 3.4 instead of the choice
given in Salimans and Ho (2021), we can obtain an estimate for the unnoised
target density at time $t=0$. We found in our experiments that in some cases
this could be used within an ordinary MCMC sampler to generate samples from
the full posterior with equivalent accuracy to the annealed sampling method.
For particularly challenging distributions, e.g. multimodal distributions, the
annealed sampling procedure had better mixing and allowed accurate recovery of
the mode weights.

##  4 Related work

#### Normalising flows

The most similar method in the literature is that of Mesquita et al. (2020),
who approximate the subposteriors with a discrete normalising flow, i.e. a
sequence of invertible neural network transformations mapping a reference
distribution onto the target. Since their work, diffusion-based modelling has
emerged as an alternative to normalising flows that is more efficient to train
and performs better on density estimation (Song et al., 2021). In addition,
our formulation enables the use of a wider class of non-invertible neural
network architectures and parameterises a density estimate that can be
evaluated without using the change-of-variables formula.

#### Neural density estimation

Score matching methods have previously been used in neural density estimation
by Saremi et al. (2018), who use denoising score matching to train a neural
network approximation to a kernel density estimate. Song and Ermon (2019) note
that this estimator can be difficult to use in MCMC sampling as fixing a small
bandwidth produces poor estimates in low density areas, causing poor mixing of
the MCMC sampler and a failure to recover mode weights in multimodal
distributions. This would cause issues in divide-and-conquer MCMC since the
subposterior distributions often have poor overlap with the full posterior, so
accuracy in low density areas is important. The use of diffusion models to
learn an energy function approximation was proposed by Salimans and Ho (2021)
and developed further by Du et al. (2023) in order to sample from compositions
of diffusion image generation models.

##  5 Experiments

The experiments in this section were chosen to demonstrate that diffusion
models can be used to accurately recover the full posterior distribution in
divide-and-conquer problems where other merging methods struggle, in
particular where the subposteriors have poor overlap or are significantly non-
Gaussian. We compare to the following methods in the literature:

  * •

Consensus Monte Carlo (CMC) (Scott et al., 2016), where a weighted average of
samples from each shard is taken.

  * •

Sub-posteriors with inflation, scaling and shifting (SwISS) (Vyner et al.,
2023), where subposteriors are transformed with an affine transformation to
approximate the full posterior.

  * •

Gaussian parametric density estimation (Neiswanger et al., 2014), which is
identical to the noise prior used in the diffusion approximation.

  * •

Semiparametric density estimation (Neiswanger et al., 2014), which is a
product of the Gaussian estimate and a nonparametric correction factor that
scales better to high dimensional problems than pure kernel density
estimation. This was implemented with the R package parallelMCMCcombine
(Miroshnikov and Conlon, 2014).

  * •

Gaussian process approximation to the log density (Nemeth and Sherlock, 2018),
using the GPJax python package (Pinder and Dodd, 2022).

The merging methods were compared numerically to samples generated from the
full posterior using three sample based discrepancy metrics: Mahalanobis
distance (Mah), integrated absolute distance (IAD), and mean absolute skew
deviation (Skew). Full experimental details as well as the numerical
comparison for the toy examples can be found in Appendix A. Code to reproduce
our experiments can be found at https://github.com/ctrojan/DiffusionDnC. We
report training times for the methods requiring an optimisation phase, and
sampling times for the methods requiring additional MCMC sampling.

The diffusion approximations used the variance preserving SDE. The neural
network architecture was the same for each experiment, with the exception of
the input and output layers which must have the same size as the dimension of
the target distribution. The number of training epochs was also chosen so that
the number of training updates was the same for each experiment. This was done
to highlight the fact that the diffusion models did not require additional
hyperparameter tuning to perform well across different experiments. The neural
network training on each shard can be done in parallel, so while training time
makes up the majority of the execution time, it does not depend on the number
of shards. This made the execution time of the diffusion merging algorithm
very similar across experiments, regardless of the complexity of the target
distribution, the number of shards, or the length of the MCMC chains sampled
from each subposterior.

Figure 2: Merged posterior contour plots for the toy logistic regression
example.

###  5.1 Toy logistic regression

Our first example is a synthetic logistic regression dataset, with a
1-dimensional covariate $x\sim N(0.5,1)$. The true value of the parameter of
interest is $\theta=(-3,-3)$, leading to a low positive rate of around 0.1.
1000 datapoints were generated and split across 15 shards, creating a
moderately challenging merging problem as the number of positive examples on
each shard varies considerably, so the subposteriors are both non-Gaussian and
very dissimilar. See Figure 2 for the posterior contour plots for each method
– note that the location is difficult to estimate without using density
estimation, likely because the subposteriors are skewed so assumptions of
Gaussianity do not hold. Only the Gaussian process and diffusion
approximations recovered the true posterior with reasonable accuracy.

###  5.2 Toy mixture of Gaussians

(a) Full posterior

(b) Subposteriors

Figure 3: Mixture of Gaussians posterior contour plots for $\theta_{1}$ and
$\theta_{2}$.

In this example, the data was drawn from a 1D mixture of 3 Gaussians,
$x\sim\frac{1}{3}N(\theta_{1},0.2)+\frac{1}{3}N(\theta_{2},0.2)+\frac{1}{3}N(\theta_{3},0.2)$
with $\theta=(0.4,0,-0.4)$. This gives a posterior distribution with 6 modes
since the likelihood is invariant to label switching of the $\theta_{i}$. The
generated dataset was of size 2000, and was split across 4 shards. See Figure
3 for an illustration of the full posterior and subposteriors – note that the
subposterior modes appear in differing locations, making it difficult to
recover the structure of the full posterior. See Figure 4 for a comparison of
the merged posterior contour plots for $\theta_{1}$ and $\theta_{2}$. The
diffusion approximation was the only method to accurately recover the full
posterior’s mode locations and weights.

Figure 4: Merged posterior contour plots for first two parameters in the
mixture of Gaussians example.

###  5.3 Robust linear regression on the combined cycle power plant dataset

Inspired by Chan et al. (2023), we consider a robust linear regression example
on the combined cycle power plant dataset (Tfekci and Kaya, 2014). This
dataset consists of 9568 hourly observations of 4 features as well as the net
hourly electrical output of the power plant, which is the regression target.
The model consists of a linear fit to the data with $t$-distributed errors to
increase robustness to outliers. We sample from the joint posterior
distribution of the regression coefficients $\beta\in\mathbb{R}^{5}$ and noise
scale $\sigma\in\mathbb{R}_{>0}$, fixing the degrees of freedom of the noise
distribution to 5. We split the data randomly across 8 shards, reporting
results in Table 1 as an average over 5 splits of the dataset, with standard
deviations in parentheses. In this example, the full marginal distribution of
$\sigma$ is very difficult to estimate from the subposteriors as is is highly
skewed and has poor overlap with the subposteriors. No method succeeded at
accurately recovering this marginal, but the SwISS and diffusion
approximations were closest to the true location and had the best accuracy
overall.

Table 1: Power plant discrepancies and average wall clock execution time. Method | Mah | IAD | Skew | Training | Sampling
---|---|---|---|---|---
Consensus | 6.25 (0.04) |  0.21 (0.02) |  0.06 (0.00) | - | -
SwISS |  4.14 (0.03) |  0.21 (0.02) |  0.06 (0.00) | - | -
Gaussian | 6.25 (0.04) |  0.21 (0.02) | 0.14 (0.00) | - | -
Semiparametric | 5.74 (0.89) | 0.90 (0.08) | 2.88 (1.38) | - | 255s
Gaussian process | 7.58 (0.07) | 0.28 (0.01) | 0.14 (0.01) | 82s | 1743s
Diffusion |  4.14 (0.04) |  0.21 (0.02) | 0.07 (0.01) | 100s | 5s

###  5.4 Logistic regression on the spambase dataset

We fit a logistic regression to the spambase dataset (Hopkins et al., 1999),
which consists of 4600 e-mails classified as spam or not spam, summarised by
57-dimensional vectors of word and character frequencies. The parameter space
is relatively high dimensional ($d=58$), creating a challenging merging
problem since the subposteriors vary in shape and have poor overlap, as well
as requiring a large number of samples to summarise each distribution. We
split the data randomly across 4 shards. Results are reported in Table 2 as an
average over 5 splits, with standard deviations in parentheses. The diffusion
merge significantly outperformed the other algorithms on Mahalanobis distance
and IAD, showing that the location and scale of the full posterior
distribution were recovered more accurately. However, the skew deviation was
higher than for the SwISS algorithm due to underestimation of the skew on the
subposteriors.

Table 2: Spambase discrepancies and average wall clock execution time. Method | Mah | IAD | Skew | Training | Sampling
---|---|---|---|---|---
Consensus | 6.00 (0.76) | 0.24 (0.02) | 0.26 (0.01) | - | -
SwISS | 7.04 (0.64) | 0.26 (0.02) |  0.22 (0.01) | - | -
Gaussian | 6.02 (0.77) | 0.24 (0.02) | 0.37 (0.00) | - | -
Semiparametric | 6.10 (0.35) | 0.29 (0.02) | 0.42 (0.04) | - | 304s
Gaussian process | 6.25 (0.96) | 0.29 (0.03) | 0.37 (0.00) | 69s | 1246s
Diffusion |  4.54 (0.79) |  0.17 (0.02) | 0.26 (0.01) | 149s | 4s

##  6 Discussion

In this paper we proposed the use of diffusion generative modelling to merge
MCMC samples generated in parallel from disjoint subsets of the full dataset.
The resulting method is embarrassingly parallel, with the exception of the
final merging step where a new MCMC run is performed using the diffusion
approximations. Our method outperformed existing merging algorithms in the
literature on complex and high dimensional posterior distributions. It is also
more computationally efficient on complex problems than existing density
estimation approaches. This is because the MCMC sampling stage is very
efficient – the density approximation is cheap to evaluate, with a cost that
is independent of the number of samples used to train it. The majority of the
computational cost comes from the training time of the neural networks, which
scales well to larger problems and can be done in parallel.

#### Limitations

The main limitation to our approach is that it is more computationally costly
than methods which do not require optimisation or further MCMC sampling, such
as consensus Monte Carlo and SwISS. In cases where the simplifying assumptions
made by these methods fail, however, it can be used to recover the full
posterior distribution with greater accuracy at a moderate computational cost.

## Acknowledgments and Disclosure of Funding

CT acknowledges the support of the EPSRC-funded EP/S022252/1 Centre for
Doctoral Training in Statistics and Operational Research in Partnership with
Industry (STOR-i); PF was supported by EPSRC grants EP/Y028783/1, EP/R034710/1
and EP/R018561/1; CN was supported by EPSRC grants EP/Y028783/1 and
EP/V022636/1.

## References

  * Bardenet et al. (2017) R. Bardenet, A. Doucet, and C. Holmes.  On Markov chain Monte Carlo methods for tall data.  _Journal of Machine Learning Research_ , 18(47):1–43, 2017.
  * Bradbury et al. (2018) J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang.  JAX: composable transformations of Python+NumPy programs, 2018.  URL http://github.com/google/jax.
  * Cabezas et al. (2024) A. Cabezas, A. Corenflos, J. Lao, and R. Louf.  BlackJAX: Composable Bayesian inference in JAX, 2024.  arXiv preprint arXiv:2402.10797.
  * Chan et al. (2023) R. S. Y. Chan, M. Pollock, A. M. Johansen, and G. O. Roberts.  Divide-and-conquer fusion.  _Journal of Machine Learning Research_ , 24(193):1–82, 2023.
  * De Bortoli et al. (2024) V. De Bortoli, M. Hutchinson, P. Wirnsberger, and A. Doucet.  Target score matching, 2024.  arXiv preprint arXiv:2402.08667.
  * DeepMind et al. (2020) DeepMind, I. Babuschkin, K. Baumli, A. Bell, S. Bhupatiraju, J. Bruce, P. Buchlovsky, D. Budden, T. Cai, A. Clark, I. Danihelka, A. Dedieu, C. Fantacci, J. Godwin, C. Jones, R. Hemsley, T. Hennigan, M. Hessel, S. Hou, S. Kapturowski, T. Keck, I. Kemaev, M. King, M. Kunesch, L. Martens, H. Merzic, V. Mikulik, T. Norman, G. Papamakarios, J. Quan, R. Ring, F. Ruiz, A. Sanchez, L. Sartran, R. Schneider, E. Sezener, S. Spencer, S. Srinivasan, M. Stanojević, W. Stokowiec, L. Wang, G. Zhou, and F. Viola.  The DeepMind JAX Ecosystem, 2020.  URL http://github.com/google-deepmind.
  * Du et al. (2023) Y. Du, C. Durkan, R. Strudel, J. B. Tenenbaum, S. Dieleman, R. Fergus, J. Sohl-Dickstein, A. Doucet, and W. S. Grathwohl.  Reduce, reuse, recycle: Compositional generation with energy-based diffusion models and MCMC.  In _Proceedings of the 40th International Conference on Machine Learning_ , volume 202 of _Proceedings of Machine Learning Research_ , pages 8489–8510. PMLR, 2023.
  * Geffner et al. (2023) T. Geffner, G. Papamakarios, and A. Mnih.  Compositional score modeling for simulation-based inference, 2023.  arXiv preprint arXiv:2209.14249v3.
  * Heek et al. (2023) J. Heek, A. Levskaya, A. Oliver, M. Ritter, B. Rondepierre, A. Steiner, and M. van Zee.  Flax: A neural network library and ecosystem for JAX, 2023.  URL http://github.com/google/flax.
  * Ho et al. (2020) J. Ho, A. Jain, and P. Abbeel.  Denoising diffusion probabilistic models.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 6840–6851. Curran Associates, Inc., 2020.
  * Hopkins et al. (1999) M. Hopkins, E. Reeber, G. Forman, and J. Suermondt.  Spambase.  UCI Machine Learning Repository, 1999.  DOI: https://doi.org/10.24432/C53G6X.
  * Huang et al. (2023) L. Huang, J. Qin, Y. Zhou, F. Zhu, L. Liu, and L. Shao.  Normalization techniques in training DNNs: Methodology, analysis and application.  _IEEE Transactions on Pattern Analysis & Machine Intelligence_, 45(08):10173–10196, 2023.
  * Karras et al. (2022) T. Karras, M. Aittala, T. Aila, and S. Laine.  Elucidating the design space of diffusion-based generative models.  _arXiv preprint 2206.00364_ , 2022.
  * (14) M. Kelly, R. Longjohn, and K. Nottingham.  The UCI machine learning repository.  https://archive.ics.uci.edu.
  * Kingma and Ba (2014) D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization, 2014.  arXiv preprint arXiv:1412.6980.
  * Li et al. (2020) T. Li, A. K. Sahu, A. Talwalkar, and V. Smith.  Federated learning: Challenges, methods, and future directions.  _IEEE Signal Processing Magazine_ , 37(3):50–60, 2020.
  * Mesquita et al. (2020) D. Mesquita, P. Blomstedt, and S. Kaski.  Embarrassingly parallel MCMC using deep invertible transformations.  In _Proceedings of The 35th Uncertainty in Artificial Intelligence Conference_ , volume 115 of _Proceedings of Machine Learning Research_ , pages 1244–1252. PMLR, 2020.
  * Miroshnikov and Conlon (2014) A. Miroshnikov and E. M. Conlon.  parallelMCMCcombine: An R package for Bayesian methods for big data and analytics.  _PLoS ONE_ , 9(9):e108425, 2014.
  * Neiswanger et al. (2014) W. Neiswanger, C. Wang, and E. P. Xing.  Asymptotically exact, embarrassingly parallel MCMC.  In _Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence_ , UAI’14, pages 623––632, Arlington, Virginia, USA, 2014. AUAI Press.
  * Nemeth and Fearnhead (2021) C. Nemeth and P. Fearnhead.  Stochastic gradient Markov chain Monte Carlo.  _Journal of the American Statistical Association_ , 116(533):433–450, 2021.
  * Nemeth and Sherlock (2018) C. Nemeth and C. Sherlock.  Merging MCMC subposteriors through Gaussian-process approximations.  _Bayesian Analysis_ , 13(2):507–530, 2018.
  * Øksendal (2000) B. Øksendal.  _Stochastic Differential Equations: An Introduction with Applications_.  Springer-Verlag, 5th edition, 2000.
  * Pinder and Dodd (2022) T. Pinder and D. Dodd.  GPJax: A Gaussian process framework in JAX.  _Journal of Open Source Software_ , 7(75):4455, 2022.
  * Ramachandran et al. (2017) P. Ramachandran, B. Zoph, and Q. V. Le.  Searching for activation functions, 2017.  arXiv preprint arXiv:1710.05941.
  * Salimans and Ho (2021) T. Salimans and J. Ho.  Should EBMs model the energy or the score?  In _Energy Based Models Workshop - ICLR 2021_ , 2021.
  * Saremi et al. (2018) S. Saremi, A. Mehrjou, B. Schölkopf, and A. Hyvärinen.  Deep energy estimator networks, 2018.  arXiv preprint arXiv:1805.08306.
  * Scott et al. (2016) S. L. Scott, A. W. Blocker, F. V. Bonassi, H. A. Chipman, E. I. George, and R. E. McCulloch.  Bayes and big data: the consensus Monte Carlo algorithm.  _International Journal of Management Science and Engineering Management_ , 11(2):78–88, 2016.
  * Song and Ermon (2019) Y. Song and S. Ermon.  Generative modeling by estimating gradients of the data distribution.  In _Advances in Neural Information Processing Systems_ , volume 32, pages 11895–11907. Curran Associates, Inc., 2019.
  * Song et al. (2021) Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, and B. Pool.  Score-based generative modeling through stochastic differential equations.  In _International Conference on Learning Representations_ , 2021.
  * Tfekci and Kaya (2014) P. Tfekci and H. Kaya.  Combined Cycle Power Plant.  UCI Machine Learning Repository, 2014.  DOI:https://doi.org/10.24432/C5002N.
  * Vincent (2011) P. Vincent.  A connection between score matching and denoising autoencoders.  _Neural Computation_ , 23(7):1661–1674, 2011.
  * Vyner et al. (2023) C. Vyner, C. Nemeth, and C. Sherlock.  SwISS: A scalable Markov chain Monte Carlo divide-and-conquer strategy.  _Stat_ , 12(1):e523, 2023.
  * Welling and Teh (2011) M. Welling and Y. W. Teh.  Bayesian learning via stochastic gradient Langevin dynamics.  In _Proceedings of the 28th International Conference on Machine Learning_ , ICML’11, pages 681–688. Omnipress, 2011.

##  Appendix A Experimental details

Experiments were run using Python 3 and the JAX package (Bradbury et al.,
2018). Real datasets were provided by the UCI machine learning repository
(Kelly et al., ) via the ucimlrepo Python package. Both the spambase and
combined cycle power plant datasets are released under a CC BY 4.0 licence.
Experiments were run on a Linux virtual machine with Ubuntu 22.04.1 LTS on
CPU, which was an Intel® Xeon(R) Gold 6248R CPU @ 3.00GHz × 4.

###  A.1 Discrepancy metrics

The merging methods were compared numerically to samples generated from the
full posterior using the following sample based discrepancy metrics:

#### Mahalanobis distance

$d_{Mah}=\sqrt{(\mu_{a}-\mu_{f})^{\top}\,V_{f}^{-1}\,(\mu_{a}-\mu_{f})}$,
where $\mu_{a}$ and $\mu_{f}$ are the means of the approximate and exact
samples respectively and $V_{f}$ is the sample covariance matrix of the exact
samples. This is useful for assessing how close the location of the
approximate samples is to that of the true distribution, accounting for
different scales in different directions.

#### Integrated absolute distance

$d_{IAD}=\frac{1}{2d}\sum_{i=1}^{d}\int_{\Theta_{i}}|\hat{\pi}_{a,i}(\theta_{i})-\hat{\pi}_{f,i}(\theta_{i})|\,\mathrm{d}\theta_{i}$
. Here, $\hat{\pi}_{.,i}$ is a kernel density estimate of the marginal density
of dimension $i$ based on the approximate ($\hat{\pi}_{a,i}$) or exact
($\hat{\pi}_{f,i}$) samples. Integration range was chosen by computing the
union of the 5-sigma intervals centered on the mean for both sets of samples,
which was sufficient to accurately obtain the number of decimal places
reported. This measures discrepancy on the 1D marginals.

#### Mean absolute skew deviation

$d_{\text{skew}}=\frac{1}{d}\sum_{i=1}^{d}|\hat{\gamma}_{f}^{(i)}-\hat{\gamma}_{a}^{(i)}|$,
where $\hat{\gamma}^{(i)}_{.}$ is a sample approximation of the skewness of
the 1D marginal density of component $i$. Here,
$\hat{\gamma}^{(i)}=\frac{1}{n}\sum_{j=1}^{n}[(\theta_{i}-\mu_{i})/\sigma_{i}]^{3}$,
$\mu_{i}$ and $\sigma_{i}$ being the sample mean and standard deviation of
component $i$. This is a measure of similarity in shape that is independent of
differences in location and scale between the distributions.

###  A.2 MCMC sampling

MCMC inference was carried out with the BlackJAX Python package (Cabezas et
al., 2024). The NUTS algorithm was used, with step sizes tuned to achieve an
acceptance rate close to 80%. In the toy examples, samplers were initialised
at an area of high posterior density and 10 000 samples were drawn from each
subposterior after a burn in of 100.

Since the evaluation time of the Gaussian process approximation scales
quadratically with the number of inducing points, we follow Nemeth and
Sherlock (2018) in thinning the MCMC chains before use. In each example the
subposterior chains were thinned to 1000 samples.

###  A.3 Diffusion sampling

For sampling from the diffusion approximations, we used the learned density
approximation at a fixed time $t=0$ in the NUTS algorithm, with the exception
of the mixture of Gaussians experiment where the annealed sampling procedure
of Du et al. (2023) was used. This algorithm is described in more detail in
Algorithm 2.

Result: Sample $\theta$ from the target distribution

Require: Sequence $p_{t}$ of target densities, MCMC update function and

its hyperparameters $\eta$, number of steps $n_{outer}$, $n_{inner}$

$t\leftarrow 1$

$\theta_{1}\sim p_{1}$

for  _$i\in\\{1,\ldots,n_{outer}\\}$_ do

$t\leftarrow t-1/n_{outer}$

$\theta_{t}\leftarrow\theta_{t-1}$

for  _$j\in\\{1,\ldots,n_{inner}\\}$_ do

$\theta_{t}\leftarrow\mbox{MCMC\\_update}(\theta_{t};p_{t},\eta)$

end for

end for

$\theta\leftarrow\theta_{t}$

Algorithm 2 Annealed MCMC sampling

###  A.4 Details for individual experiments

#### Toy logistic regression

In this example, the dataset consisted of a 1-dimensional covariate $x_{i}\sim
N(0.5,1)$ with binary class labels generated as
$y_{i}\sim\text{Bernoulli}(\text{sigmoid}(\beta_{1}+\beta_{2}x_{i}))$. The
true value of the parameter of interest was $\beta=(-3,-3)$. A normal prior of
$N(0,5)$ was placed over all parameters. An illustration of the subposteriors
can be found in Figure 5 \- note that they are very skewed and differ
substantially in location and scale. Most of them have poor overlap with the
full posterior (overlaid in blue) which is very concentrated around the true
parameter value. The numerical results of the algorithms are reported in Table
3.

Figure 5: Subposterior contours for the logistic regression with full posterior in blue. Table 3: Toy logistic regression discrepancies and wall clock time. Method | Mah | IAD | Skew | Training time | Sampling time
---|---|---|---|---|---
Consensus | 1.28 | 0.42 | 0.02 | - | -
SwISS | 2.57 | 0.75 | 0.03 | - | -
Gaussian | 1.28 | 0.41 | 0.17 | - | -
Semiparametric | 2.07 | 0.63 | 0.74 | - | 65s
Gaussian process | 0.08 | 0.04 | 0.22 | 36s | 2766s
Diffusion | 0.08 | 0.03 | 0.01 | 99s | 8s

#### Toy Gaussian mixture

A standard Normal prior placed over all parameters. To ensure good mixing
between modes when sampling from the full posterior and subposteriors, MCMC
samplers randomly permuted the $\theta_{i}$ at each step as in Neiswanger et
al. (2014), since this is a move that leaves the posterior density invariant.
Label switching was not used for sampling from the merged distributions in
order to give an accurate representation of the fitted density. Instead, the
annealed sampling procedure was used in the diffusion, and 10 independent
chains were run for the Gaussian process. The annealed sampling for the
diffusion model used HMC rather than NUTS since it allows the update steps to
be vectorised efficiently. The number of leapfrog steps was set at 3 and a
fixed step size was used across all times. 300 evenly spaced time points were
used, and 1 MCMC step was performed at each time point. Numerical results are
reported in Table 4.

Table 4: Mixture of Gaussians discrepancies and wall clock time. Method | Mah | IAD | Skew | Training time | Sampling time
---|---|---|---|---|---
Consensus | 0.14 | 0.53 | 0.16 | - | -
SwISS | 0.16 | 0.25 | 0.14 | - | -
Gaussian | 0.13 | 0.55 | 0.14 | - | -
Semiparametric | 0.14 | 0.52 | 0.12 | - | 24s
Gaussian process | 0.19 | 0.24 | 0.21 | 35s | 3412s
Diffusion | 0.11 | 0.04 | 0.12 | 98s | 24s

#### Power plant robust regression

Priors of $N(0,10^{2})$ were placed on the regression coefficients, and for
the noise standard deviation a chi-squared prior $\chi(1)$ with scale
parameter 10 was used. MCMC samplers were initialised at the mean of the prior
distribution. 50 000 samples were generated from each posterior after a burn
in of 100. For the merged distributions, the GP and diffusion samplers were
initialised at the mode of the Gaussian approximation to the full posterior.

#### Spambase logistic regression

Priors of $N(0,5)$ were placed on all parameters. The full posterior and
subposterior MCMC samplers were initialised at areas of high posterior density
by using 20 epochs of the ADAM optimiser with a learning rate of $10^{-2}$ to
approximately find the MAP parameter estimate. 50 000 samples were generated
from each posterior after a burn in of 100. For the merged distributions, the
GP and diffusion samplers were initialised at the mode of the Gaussian
approximation to the full posterior. Using the inverse of its covariance
matrix as the inverse mass matrix in NUTS sampling greatly reduced the number
of leapfrog steps required to sample from the full posterior approximation.

###  A.5 Neural network architecture

The neural network used was a residual MLP, implemented with the Flax Python
package (Heek et al., 2023). The noise standard deviation $s(t)$ was
concatenated to the input of each layer. In this architecture, each hidden
layer has the same dimension, with the exception of the output layer which has
the same dimension as $x$. After the first layer, the hidden layers are
organised in blocks of two with skip connections that add the input of the
block to the output of its hidden layers. The size of the network was kept the
same across experiments, with the exception of the output layer which must
have the same dimension as the target distribution. A neural network with 1
residual block and the hidden layer dimension 32 was sufficient for the
examples considered. The activation function for the hidden layers was
$\text{SiLU}(x)=x\,\text{sigmoid}(x)$, a smooth alternative to the ReLU
activation function (Ramachandran et al., 2017).

###  A.6 Training

The Adam optimiser (Kingma and Ba, 2014) as implemented by the Optax Python
package (DeepMind et al., 2020) was used for model training, with a batch size
of 32 and default hyperparameters. The Gaussian process parameters were fit
for 200 epochs. For the diffusion models, 500 epochs of training were used for
the toy examples, and 100 for the power plant and spambase examples.

##  Appendix B Scalability to very high-dimensional problems

Here, we report the results of an additional high-dimensional synthetic
logistic regression experiment to show that our method also scales well to
dimensions higher that those considered in the main paper. In this experiment,
the dataset consisted of 1000 realisations of a 99-dimensional covariate $x$
drawn from a standard Normal distribution, with binary class labels generated
as $y_{i}\sim\text{Bernoulli}(\text{sigmoid}(\beta_{0}+\beta^{\top}x_{i}))$.
The true value of $\beta_{0}$ was -5, and the remaining $\beta_{i}$ were drawn
from a standard Normal distribution. A prior of $N(0,5)$ was placed over all
parameters. The dataset was split across 4 shards and 50 000 samples were
drawn from all posterior distributions after a burn in of 100. Experimental
details were otherwise the same as in the other experiments.

This was repeated for 5 simulated datasets and the numerical results are
reported in Table 5 as an average over the 5 runs with standard deviations in
brackets. SwISS outperformed all other methods here as the subposteriors were
similarly shaped and had low skew. Note that with $d=100$ here, the non- and
semi- parametric density estimation approaches had significantly worse
performance than the other methods, while their execution time increased
significantly compared to the lower-dimensional experiments. In comparison,
the diffusion method was outperformed only by SwISS and its execution time was
not much higher than in the other experiments.

Table 5: High-dimensional logistic regression discrepancies and average execution time. Method | Mah | IAD | Skew | Training | Sampling
---|---|---|---|---|---
Consensus | 4.73 (0.17) | 0.25 (0.00) | 0.02 (0.00) | - | -
SwISS |  3.36 (0.17) |  0.17 (0.01) |  0.01 (0.00) | - | -
Gaussian | 4.73 (0.17) | 0.25 (0.00) | 0.02 (0.00) | - | -
Semiparametric | 6.94 (0.25) | 0.31 (0.01) | 0.03 (0.00) | - | 733s
Gaussian process | 7.41 (0.17) | 0.32 (0.02) | 0.03 (0.00) | 92s | 12025s
Diffusion | 4.02 (0.17) | 0.21 (0.00) | 0.02 (0.00) | 189s | 10s

##  Appendix C Link between diffusion formulation and parametric Gaussian
approximation

In the proposed merging algorithm, we obtain a sequence of density estimates
for each subposterior that interpolates between a Gaussian prior distribution
$\hat{p}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{1}$ and the final non-
Gaussian diffusion approximation
$\hat{p}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{0}$. Using the
reparameterisation in Section 3.2, the Gaussian priors are $N(\mu_{s},V_{s})$,
where $\mu_{s}$ and $V_{s}$ are the sample mean and covariance respectively of
the samples $\\{\theta^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}\\}$
generated from shard $s$.

The Gaussian prior distribution $\hat{p}_{1}$ for the full posterior is the
product of those of the subposteriors, and has the form $N(\mu,V)$ where:

$\displaystyle\mu=V\sum_{s}V_{s}^{-1}\mu_{s}\quad\mbox{and}\quad
V=\Big{[}\sum_{s}V_{s}^{-1}\Big{]}^{-1}\,.$ (15)

This is exactly the parametric Gaussian approximation of Neiswanger et al.
(2014).

##  Appendix D Failure of SDE sampling for compositions of diffusions

By adding the component score functions $\nabla\log
p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{t}$ together, we obtain the
the score of the product distribution

$\displaystyle\tilde{p}_{t}(\theta_{t})$
$\displaystyle\propto\prod_{s}p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}_{t}(\theta_{t})\,,$
(16)

which at time $t=0$ is exactly equal to the score of the target product
distribution $p^{\text{full}}$. However, this relationship does not hold for
$t>0$ since we cannot interchange the order of adding noise to a distribution
with multiplying densities together:

$\displaystyle p^{(\text{full})}_{t}(\theta_{t})$
$\displaystyle=\int_{\Theta}p_{t|0}(\theta_{t}|\theta_{0})\;p^{\text{full}}(\theta_{0})\,\mathrm{d}\theta_{0}\propto\int_{\Theta}p_{t|0}(\theta_{t}|\theta_{0})\prod_{s}p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(\theta_{0})\,\mathrm{d}\theta_{0}$
(17)
$\displaystyle\not\equiv\tilde{p}_{t}(\theta_{t})\propto\prod_{s}\int_{\Theta}p_{t|0}(\theta_{t}|\theta_{0})\;p^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}(\theta_{0})\,\mathrm{d}\theta_{0}\,.$

This relationship fails even when the target distributions
$p_{0}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}$ are Gaussian. For linear
SDEs the marginal distribution $p_{t|0}$ is $N(m(x_{0},t),S(t)^{2})$, where
$m$ is linear in $x_{0}$, i.e. $m(x,t)=m_{1}(t)x+m_{2}(t)$. So, when the
target $p_{0}$ is Gaussian with parameters $\mu$ and $V$, $p_{t}$ will also be
Gaussian and we can compute its parameters as $m(\mu,t)$ and
$m_{1}(t)^{2}V+S(t)^{2}$, since $m$ is linear and $x_{t}$ can be expressed as
$m(x_{0},t)+S(t)Z$ where $Z$ is a standard normal.

So, assuming
$p_{0}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}=N(\mu_{s},V_{s})$, we can
obtain $p^{\text{full}}_{t}$ by taking a product of Gaussian densities and
then adding noise, i.e.,

$\displaystyle p^{\text{full}}_{t}$
$\displaystyle\propto\phi(m(\mu,t),m_{1}(t)^{2}V+S(t)^{2})\,,\ \mbox{where}$
(18) $\displaystyle\mu$
$\displaystyle=V\sum_{s}V_{s}^{-1}\mu_{s}\,,\quad\mbox{and}\quad
V=\Big{[}\sum_{s}V_{s}^{-1}\Big{]}^{-1}\,.$ (19)

However, $\tilde{p}_{t}$ is a product of the noised densities
$p_{t}^{{\scriptscriptstyle(}s{\scriptscriptstyle)}}$ and will be Gaussian
with parameters:

$\displaystyle\tilde{\mu}=\tilde{V}\sum_{s}[m_{1}^{2}(t)V_{s}+S(t)^{2}]^{-1}m(\mu_{s})\,,\quad\tilde{V}=\Big{[}\sum_{s}[m_{1}(t)^{2}V_{s}+S(t)^{2}]^{-1}\Big{]}^{-1}\,.$
(20)

So, unless the $\mu_{s}$ are zero, the mean of $\tilde{p}_{t}$ will typically
be different to that of $p^{\text{full}}_{t}$. Its variance will also be
incorrect - even if all of the $V_{s}$ are equal, we obtain
$\tilde{V}=m_{1}(t)^{2}V+\frac{1}{S}S(t)^{2}$, so taking the product in this
way effectively shrinks the noise added by the diffusion.

Hence, while the score sum does get closer to the desired score function
$\nabla\log p^{(\text{full})}_{t}$ as time approaches zero, this is not
sufficient to correct errors accumulated earlier in sampling by using an
inaccurate score function estimate. For $t<1$, the density induced by solving
the reverse SDE from 1 to $t$ using the score sum $\nabla\log\tilde{p}_{t}$ is
not the same as $p^{(\text{full})}_{t}$ or even $\tilde{p}_{t}$, since this
would still result in the correct density at time 0. This is because
$\tilde{p}_{t}$ is not obtained by the noising SDE whose coefficients are used
in the reverse SDE.

The problem of mismatch between score function and induced density is also
encountered to some extent in the usual diffusion model sampling procedure,
where numerical integration error causes generated samples to not follow the
‘true’ density induced by the reverse SDE. Song et al. (2021) noticed this and
proposed using a step of ULA after each SDE integration step to correct the
distribution. They called this predictor-corrector sampling since it uses the
MCMC sampler to correct the distribution of the samples proposed by the SDE
solver. The annealed sampling procedure described in Section 3.5 can be seen
as corrector-only sampling in the framework of Song et al. (2021) and is
similar to the method used by Song and Ermon (2019) for sampling in the
discrete-time formulation of the variance exploding SDE.


#  Khinchin-type inequalities via Hadamard’s factorisation

Alex Havrilla, Piotr Nayar and Tomasz Tkocz Carnegie Mellon University;
Pittsburgh, PA 15213, USA. Research supported in part by NSF grant
DMS-1955175. University of Warsaw, Banacha 2, 02-097 Warszawa, Poland.
Research supported in part by the National Science Centre, Poland, grant
2018/31/D/ST1/01355. Carnegie Mellon University; Pittsburgh, PA 15213, USA.
Email<EMAIL_ADDRESS>Research supported in part by NSF grant
DMS-1955175.

(18th February 2020)

###### Abstract

We prove Khinchin-type inequalities with sharp constants for type L random
variables and all even moments. Our main tool is Hadamard’s factorisation
theorem from complex analysis, combined with Newton’s inequalities for
elementary symmetric functions. Besides the case of independent summands, we
also treat ferromagnetic dependencies in a nonnegative external magnetic field
(thanks to Newman’s generalisation of the Lee-Yang theorem). Lastly, we
compare the notions of type L, ultra sub-Gaussianity (introduced by Nayar and
Oleszkiewicz) and strong log-concavity (introduced by Gurvits), with the
latter two being equivalent.

2010 Mathematics Subject Classification. Primary 60E15; Secondary 30D15,
26D15.

Key words. Khinchin inequality, moment comparison, log-concave sequence,
Laguerre-Pólya class

##  1 Introduction and results

Motivated by Ising models of Lee-Yang type, Newman in [13] made the following
definition: a (real-valued) random variable $X$ is of type $\mathscr{L}$, if
(i) for some positive constants $C$ and $C^{\prime}$,
$|\mathbb{E}\exp(zX)|\leq C\exp(C^{\prime}|z|^{2})$ for all $z\in\mathbb{C}$
and (ii) $\mathbb{E}\exp(zX)$ is even with only pure imaginary zeros. Note
that (ii) implies that $X$ is symmetric, that is $-X$ has the same
distribution as $X$. Perhaps the simplest nontrivial example is a Rademacher
random variable (a symmetric random sign) taking the values $\pm 1$ with
probability $\frac{1}{2}$. Newman proved that if $X$ is of type $\mathscr{L}$,
then for an even integer $p\geq 2$,

$\|X\|_{p}\leq\|G\|_{p}\|X\|_{2},$

where $G$ is a standard Gaussian random variable (see [14] or Theorem 5 in
[13]). As usual, $\|Y\|_{p}=(\mathbb{E}|Y|^{p})^{1/p}$ denotes the $p$-norm of
a random variable $Y$.

The goal of this short note is to extend this result to a comparison
inequality for moments of arbitrary even orders. As explained later, the
assumption of evenness in (ii) may be easily dropped (it is just a matter of
cosmetic shifting). Thus we shall say that a random variable $X$ is of type
$\mathscr{L}^{\prime}$ if it satisfies (i) and (ii’) $\mathbb{E}\exp(zX)$ has
only pure imaginary zeros. Our main result is the following theorem.

######  Theorem 1.

Let $X$ be a random variable of type $\mathscr{L}^{\prime}$. Then for every
even integers $2\leq p\leq q$, we have

$\|X\|_{q}\leq\frac{\|G\|_{q}}{\|G\|_{q}}\|X\|_{p},$ (1)

where $G$ is a standard Gaussian random variable.

The usefulness of this result in Khinchin inequalities for sums of independent
random variables comes from the following evident fact.

######  Remark 2.

Finite sums of independent type $\mathscr{L}$ (resp. $\mathscr{L}^{\prime}$)
random variables are of type $\mathscr{L}$ (resp. $\mathscr{L}^{\prime}$).

Thus the above theorem immediately gives classical Khinchin inequalities with
sharp constants for even moments, the result due to Nayar and Oleszkiewicz
from [10]. Furthermore, it automatically extends to Hilbert space coefficients
by isometrical embeddings into $L_{p}$ spaces based on standard Gaussians
(Orlicz-Szarek’s argument, see Remark 3 in [21]).

######  Corollary 3.

Let $(H,\|\cdot\|)$ be a separable (real or complex) Hilbert space. If
$X_{1},\ldots,X_{n}$ are independent type $\mathscr{L}^{\prime}$ random
variables, then for every vectors $v_{1},\ldots,v_{n}$ in $H$, the sum
$X=\sum_{j=1}^{n}X_{j}v_{j}$ satisfies (1) for all positive even integers
$p\leq q$, where we denote $\|X\|_{p}=(\mathbb{E}\|X\|^{p})^{1/p}$.

Thanks to Newman’s generalisation of the Lee-Yang theorem and Griffiths’
method of a “ghost” spin (see [3, 12, 15]), we moreover obtain (1) when $X$ is
a positive linear combination of random variables with ferromagnetic
dependencies $J=(J_{jk})$ in a nonnegative external magnetic field
$h=(h_{j})$.

######  Corollary 4.

Let $\mu_{1},\ldots,\mu_{n}$ be Borel probability measures on $\mathbb{R}$,
each one of type $\mathscr{L}$. Suppose that $(X_{1},\dots,X_{n})$ is a random
vector in $\mathbb{R}^{n}$ whose law $\rho$ on $\mathbb{R}^{n}$ is of the form

$\mathrm{d}\rho(x_{1},\ldots,x_{n})=Z^{-1}\exp\left(\sum_{j=1}^{n}h_{j}x_{j}+\sum_{j,k=1}^{n}J_{jk}x_{j}x_{k}\right)\mathrm{d}\mu_{1}(x_{1})\dots\mathrm{d}\mu_{n}(x_{n})$
(2)

with $h_{j}\geq 0$, $J_{jk}\geq 0$ for all $j,k\leq n$, where $Z$ is the
normalising constant. Then for every nonnegative $a_{1},\ldots,a_{n}$, the sum
$X=\sum_{j=1}^{n}a_{j}X_{j}$ satisfies (1) for all positive even integers
$p\leq q$.

For a certain subclass of the type $\mathscr{L}$ random variables (cf. the
section below devoted to examples), we are also able to derive sharp moment
comparison between the second and $p$-th moment, $p\geq 3$.

######  Theorem 5.

Let $X$ be type $\mathscr{L}$ random variable with characteristic function of
the form $\phi_{X}(t)=e^{-at^{2}/2}\prod_{j=1}^{\infty}(1-b_{j}t^{2})$ with
$a>0$, $b_{j}\geq 0$, $\sum b_{j}\leq a$. Let
$\sigma=\sqrt{\operatorname{Var}(X)}$. Then for every $p\geq 3$,

$\mathbb{E}|\sigma Z_{1}|^{p}\leq\mathbb{E}|X|^{p}\leq\mathbb{E}|\sigma
Z_{0}|^{p},$

where $Z_{0}$ is a standard Gaussian random variable and $Z_{1}$ is a random
variable with density $(2\pi)^{-1/2}x^{2}e^{-x^{2}/2}$.

First we prove Theorem 1 (relying on the product form of the moment generating
function $\mathbb{E}\exp(zX)$ and exploiting Newton’s inequalities for the
elementary symmetric functions) and show how to obtain Corollaries 3 and 4.
Then we provide several examples of type $\mathscr{L}$ random variables and
prove Theorem 5 (using techniques from [2, 11], relying on convexity
properties of certain functions related to $|\cdot|^{p}$, hence the
restriction $p\geq 3$). We conclude with a discussion comparing type
$\mathscr{L}$, ultra sub-Gaussianity (introduced in [10]) and strong log-
concavity (introduced in [4]).

##  2 Even moments

We recall that a nonnegative sequence $(a_{n})_{n\geq 0}$ is called log-
concave if it is supported on a contiguous interval (which may be infinite)
and satisfies $a_{n}^{2}\geq a_{n-1}a_{n+1}$, $n\geq 1$. We remark that log-
concavity implies the inequalities $a_{k}a_{l}\geq a_{k+j}a_{l-j}$ for all
$k\geq l$, $0\leq j\leq l$. We refer for instance to the classical survey [20]
by Stanley.

###  2.1 Proof of Theorem 1

_Type $\mathscr{L}$ vs. $\mathscr{L}^{\prime}$._ We begin with a
straightforward observation.

######  Lemma 6.

If $X$ is a random variable of type $\mathscr{L}^{\prime}$, then there is
$c\in\mathbb{R}$ such that $X+c$ is of type $\mathscr{L}$.

###### Proof.

As in Proposition 2 in [13], by the Hadamard factorization theorem,

$\mathbb{E}\exp(zX)=\exp(bz^{2}/2+cz)\prod_{j}(1+\alpha_{j}z^{2}/2)$ (3)

for some $c\in\mathbb{R}$, $b\geq 0$ and
$0<\alpha_{1}\leq\alpha_{2}\leq\ldots$ with $\sum\alpha_{j}<\infty$ (the set
$\\{\alpha_{j}\\}$ may be empty, finite or infinite countable). Thus $Y=X-c$
is of type $\mathscr{L}$. ∎

_Reduction to type $\mathscr{L}$._ Let $X$ be a random variable of type
$\mathscr{L}^{\prime}$. Thus $Y=X-c$ is of type $\mathscr{L}$. Since $Y$ is a
symmetric random variable,

$\mathbb{E}|X|^{p}=\mathbb{E}|Y+c|^{p}=\mathbb{E}|Y+c\varepsilon|^{p}$ (4)

where $\varepsilon$ is a Rademacher random variable independent of $Y$. Note
that $Y+c\varepsilon$ is of type $\mathscr{L}$ (as a sum of two independent
random variables of type $\mathscr{L}$). Assuming the theorem holds for all
random variables of type $\mathscr{L}$, we thus conclude it for $X$ because
$X$ has the same moments as $Y+c\varepsilon$.

_Proof for type $\mathscr{L}$._ Let $X$ be a random variable of type
$\mathscr{L}$. Then (3) holds with $c=0$. Applying it to $z=\sqrt{2t}$ with
$t\geq 0$, we obtain

$\sum_{n=0}^{\infty}\frac{\mathbb{E}X^{2n}}{(2n)!}2^{n}t^{n}=\mathbb{E}\exp(\sqrt{2t}X)=\exp(bt)\prod_{j}(1+\alpha_{j}t).$

Note that $\mathbb{E}G^{2n}=\frac{(2n)!}{2^{n}n!}=(2n-1)!!$. Moreover,
$\prod_{j}(1+\alpha_{j}t)=\sum_{k=0}^{\infty}\sigma_{k}t^{k}$, where
$\sigma_{k}=\sigma_{k}(\\{\alpha_{j}\\})=\sum_{A\subset\\{\alpha_{j}\\}}\prod_{a\in
A}a$ is the $k$th elementary symmetric function. Therefore we arrive at

$\sum_{n=0}^{\infty}\frac{\mathbb{E}X^{2n}}{n!\mathbb{E}G^{2n}}t^{n}=\sum_{n=0}^{\infty}\left(\sum_{k=0}^{n}\frac{b^{n-k}}{(n-k)!}\sigma_{k}\right)t^{n}.$

For $n\geq 0$, set

$r_{n}=\frac{\mathbb{E}X^{2n}}{\mathbb{E}G^{2n}}.$ (5)

By comparing the coefficients,

$r_{n}=n!\sum_{k=0}^{n}\frac{b^{n-k}}{(n-k)!}\sigma_{k}=\sum_{k=0}^{n}\binom{n}{k}b^{n-k}\sigma_{k}k!.$

Claim. The sequence $(\sigma_{k}k!)_{k\geq 0}$ is log-concave.

###### Proof.

Fix $k\geq 1$. By Newton’s inequalities,

$\frac{\sigma_{k}(\\{\alpha_{j}\\}_{j=1}^{n})^{2}}{\binom{n}{k}^{2}}\geq\frac{\sigma_{k-1}(\\{\alpha_{j}\\}_{j=1}^{n})}{\binom{n}{k-1}}\frac{\sigma_{k+1}(\\{\alpha_{j}\\}_{j=1}^{n})}{\binom{n}{k}}.$

Letting $n\to\infty$ gives the claim. ∎

Walkup’s theorem (see, e.g. [5, 9, 10, 22]) which says that the binomial
convolution preserves log-concavity yields that the sequence $(r_{n})_{n\geq
0}$ is also log-concave. In particular, by monotonicity of slopes, the
sequence $\frac{\log r_{n}-\log r_{0}}{n}$ is nonincreasing. In view of
$r_{0}=1$, this gives that the sequence $(r_{n}^{1/n})_{n\geq 0}$ is
nonincreasing, hence (1) and the proof of Theorem 1 is complete. $\square$

Incidentally, since $(r_{n})_{n\geq 0}$ being log-concave means that $X$ is
ultra sub-Gaussian (see Section 5 and [10]), we have obtained the following
result.

######  Theorem 7.

If a random variable $X$ is of type $\mathscr{L}^{\prime}$, then the sequence
$(r_{n})_{n=0}^{\infty}$ defined in (5) is log-concave. In particular, type
$\mathscr{L}$ random variables are ultra sub-Gaussian.

###### Proof.

It has been shown above that if $X$ is of type $\mathscr{L}$, then $X$ is
ultra sub-Gaussian. Consequently, by Lemma 6 and (4), the same holds if $X$ is
of type $\mathscr{L}^{\prime}$. ∎

###  2.2 Proof of Corollary 3

First suppose that $H$ is real Hilbert space. Fix an orthonormal basis
$e_{1},e_{2},\ldots$ in $H$. Let $g_{1},g_{2},\ldots$ be i.i.d. standard
Gaussian random variables (independent of the $X_{j}$). Consider the map

$H\ni x\mapsto c_{q}\sum_{k}\\!\left\langle x,e_{k}\right\rangle\\!g_{k}\in
L_{q}(\Omega,\mathbb{P}).$

Since sums of independent Gaussians are Gaussian, this map with
$c_{q}=\|G\|_{q}^{-1}$ is an isometrical embedding. Given vectors
$v_{1},\ldots,v_{n}$ in $H$, for appropriate scalars $a_{j}$, we thus have

$\left\|\sum_{j}X_{j}v_{j}\right\|^{q}=\mathbb{E}_{g}\left|\sum_{j}a_{j}X_{j}\right|^{q},\qquad
a_{j}=c_{q}\\!\left\langle v_{j},\sum_{k}g_{k}e_{k}\right\rangle\\!.$ (6)

Using this twice, Theorem 1 applied to $\sum_{j}a_{j}X_{j}$ and Minkowski’s
inequality, we arrive at

$\displaystyle\mathbb{E}_{X}\left\|\sum_{j}X_{j}v_{j}\right\|^{q}$
$\displaystyle=\mathbb{E}_{g}\mathbb{E}_{X}\left|\sum_{j}a_{j}X_{j}\right|^{q}$
$\displaystyle\leq\left(\frac{\|G\|_{q}}{\|G\|_{p}}\right)^{q}\mathbb{E}_{g}\left(\mathbb{E}_{X}\left|\sum_{j}a_{j}X_{j}\right|^{p}\right)^{q/p}$
$\displaystyle\leq\left(\frac{\|G\|_{q}}{\|G\|_{p}}\right)^{q}\left(\mathbb{E}_{X}\left(\mathbb{E}_{g}\left|\sum_{j}a_{j}X_{j}\right|^{q}\right)^{p/q}\right)^{q/p}$
$\displaystyle=\left(\frac{\|G\|_{q}}{\|G\|_{p}}\right)^{q}\left(\mathbb{E}_{X}\left\|\sum_{j}X_{j}v_{j}\right\|^{p}\right)^{q/p}.$

This finishes the proof when $H$ is real. For the complex case, the only
required modification of the above argument is in the definition of the
isometrical embedding which now reads

$H\ni x\mapsto c_{q}\sum_{k}\Big{(}\text{Re}\\!\left\langle
x,e_{k}\right\rangle\\!g_{k}+\text{Im}\\!\left\langle
x,e_{k}\right\rangle\\!g_{k}^{\prime}\Big{)}\in L_{q}(\Omega,\mathbb{P}),$

where $g_{1},g_{1}^{\prime},g_{2},g_{2}^{\prime},\ldots$ are i.i.d. standard
Gaussian random variables. Identity (6) holds with
$a_{j}=c_{q}\Big{(}\text{Re}\\!\left\langle
v_{j},\sum_{k}g_{k}e_{k}\right\rangle\\!+\text{Im}\\!\left\langle
v_{j},\sum_{k}g_{k}^{\prime}e_{k}\right\rangle\\!\Big{)}$ and the rest of the
argument is unchanged. $\square$

###  2.3 Proof of Corollary 4

Let $(X_{1},\ldots,X_{n})$ be a random vector with distribution given by (2)
and let $\varepsilon$ be an independent Rademacher random variable. By
Proposition 6 from [15], the vector
$(Y_{0},Y_{1},\dots,Y_{n})=(\varepsilon,\varepsilon X_{1},\dots,\varepsilon
X_{n})$ has distribution $\rho^{\prime}$ of the form

$\mathrm{d}\rho^{\prime}(x_{0},x_{1},\ldots,x_{n})=Z^{\prime-1}\exp\left(\sum_{j,k=0}^{n}J_{jk}^{\prime}x_{j}x_{k}\right)\mathrm{d}\mu_{0}(x_{0})\mathrm{d}\mu_{1}(x_{1})\dots\mathrm{d}\mu_{n}(x_{n}),$

where $\mu_{0}$ is the distribution of $\varepsilon$, $J_{0,0}^{\prime}=0$,
$J_{0,k}^{\prime}=J_{k,0}^{\prime}=h_{k}/2$, $J_{jk}^{\prime}=J_{jk}$,
$j,k\geq 1$, so of the form (2) with $h\equiv 0$. Therefore, by Theorem 2 from
[14], for every $a_{0},a_{1},\dots,a_{n}\geq 0$, the sum
$S=\sum_{j=0}^{n}a_{j}Y_{j}=a_{0}\varepsilon+\sum_{j=1}^{n}a_{j}\varepsilon
X_{j}$ is of type $\mathscr{L}$ and in particular, $S$ satisfies (1). Hence,
taking $a_{0}=0$ yields that $\sum_{j=1}^{n}a_{j}X_{j}$ also satisfies (1).
$\square$

##  3 Examples of type $\mathscr{L}$ random variables

We list some examples of probability distributions of type $\mathscr{L}$. This
is based on several well-known sufficient conditions for Fourier transforms of
positive functions to have only real zeros. We mention merely in passing that
the development of such conditions, related to the study of the Laguerre-Pólya
class, was initiated mainly with connections to the Riemann hypothesis and
refer to the comprehensive survey [1]. In what follows, $X$ is a symmetric
random variable.

  1. (a)

Let $X$ be integer-valued with $\mathbb{P}\left(X=0\right)=p_{0}$ and
$\mathbb{P}\left(X=-k\right)=\mathbb{P}\left(X=k\right)=p_{k}$, $k=1,\dots,n$
for nonnegative $p_{0},\dots,p_{n}$ with $p_{0}+2\sum_{k=1}^{n}p_{k}=1$.

If $\frac{1}{2}p_{0}\leq p_{1}\leq\dots\leq p_{n}$, then
$\mathbb{E}\cos(zX)=p_{0}+\sum_{k=1}^{n}(2p_{k})\cos(kz)$ has only real zeros,
as it follows from the Eneström-Kakeya theorem (see, e.g. Problem III.204 in
[16]). As a result, $X$ is of type $\mathscr{L}$. In particular, if $X$ is
uniform on $\\{-n,\dots,1,1,\dots,n\\}$ with a possible atom at $0$ satisfying
$\mathbb{P}\left(X=0\right)\leq\frac{1}{n+1}$, then $X$ is of type
$\mathscr{L}$ (the last class was also recently studied in [6] in the context
of Khinchin inequalities).

By the symmetry of $X$, the polynomial $Q(w)=\mathbb{E}w^{X+n}$ is self-
inversive (the sequence of its coefficients is a palindrome, in other words,
$w^{2n}Q(1/w)=Q(w)$). In particular, all its roots are symmetric with respect
to the unit circle, that is if $w_{0}$ is a root of $Q$, then so is $1/w_{0}$.
For instance, if for some $\alpha\geq 1$,

$\frac{1}{2}p_{0}^{\alpha}+\sum_{k=1}^{n-1}p_{k}^{\alpha}\leq\left(\frac{2}{n-2}\right)^{\alpha-1}p_{n}^{\alpha},$

where $n$ is the number of nonzero coefficients of $Q$, then $Q$ has zeros
only on the unit circle, so $X$ is of type $\mathscr{L}$, which follows from
the main result of [18]. In general, the Jury test provides an efficient
algorithm to determine whether $Q$ has zeros only on the unit circle, given
its coefficients, see [8].

  2. (b)

Let $X$ take values in $[-1,1]$ and have a density $f$ (which is even). Each
of the following conditions implies that $X$ is of type $\mathscr{L}$.

    1. (i)

$f$ is nondecreasing on $(0,1)$.

    2. (ii)

$f$ is $C^{2}$ with $f^{\prime}<0$ and $f^{\prime\prime}<0$ on $(0,1)$.

    3. (iii)

$f(t)=h(t)^{\alpha}$, where $\alpha>-1$ and $h$ is an entire even function
which is real-valued on the real line, $h(1)=0$, $h(0)>0$ and $h^{\prime}(iz)$
is in the Laguerre-Pólya class. In particular,
$f(t)=\text{const}\cdot(1-t^{2m})^{\alpha}$ for a nonnegative integer $m$.

Moreover, if $X$ has a density on $\mathbb{R}$ of the following form, then it
is of type $\mathscr{L}$.

    1. (iv)

$f(t)=\text{const}\cdot e^{-t^{2m}}$ for a nonnegative integer $m$.

    2. (v)

$f(t)=(2\pi)^{-1/2}e^{-t^{2}/2}(1-b+bt^{2})$, $0\leq b\leq 1$.

Condition (i) is justified again by the Eneström-Kakeya theorem combined with
a limit argument (see, e.g. Problem III.205 in [16]), (ii) is due to Pólya
(see, e.g. Problem V.173 in [17]), (iii) is due to Ilieff (see [7]), which
gives (iv) by a limit argument (see also Problem V.170 and 171 in [17]), (v)
is justified by a direct computation of the moment generating function which
is $(1+bz^{2})e^{z^{2}/2}$. Moreover, if the density of $X$ is of the form
$f(t)=\text{const}\cdot e^{-|t|^{\alpha}}$ with $\alpha\geq 2$,
$\alpha\notin\\{2,4,\dots\\}$, then its characteristic function has infinitely
many non-real zeros, in particular $X$ is not of type $\mathscr{L}$ (see the
solution of Problem V.171 in [17]).

We refer to [13] for additional examples important in statistical mechanics.

##  4 Type $\mathscr{L}$ random variables with “enough Gaussianity”

The goal in this section is to prove Theorem 5. We need a little bit of
preparation. For a real parameter $b$, consider the following function

$\phi_{b}(t)=e^{-t^{2}/2}(1-bt^{2}),\qquad t\in\mathbb{R},$

which is in $L_{1}$. By the Fourier inversion formula, we have,

$f_{b}(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\phi_{b}(t)e^{-itx}\mathrm{d}t=(1-b+bx^{2})\frac{e^{-x^{2}/2}}{\sqrt{2\pi}},\qquad
x\in\mathbb{R}.$

Clearly, $f_{b}$ is nonnegative on $\mathbb{R}$ if and only if $b\in[0,1]$.
Thus $\phi_{b}$ is the characteristic function of a random variable precisely
for every $b\in[0,1]$ and then its density is $f_{b}$ and variance is the
coefficient $2\phi_{b}(it)[t^{2}]=1+2b$. Throughout, $Z_{b}$ is a random
variable with characteristic function $\phi_{b}$. In particular, $Z_{0}$ is
standard Gaussian and $Z_{1}$ has the bimodal density
$(2\pi)^{-1/2}x^{2}e^{-x^{2}/2}$.

######  Lemma 8.

For every $a>0$ and $b_{1},b_{2},\ldots\geq 0$ with $\sum b_{j}\leq a$, the
function

$e^{-at^{2}/2}\prod_{j}(1-b_{j}t^{2})$

is the characteristic function of a type $\mathscr{L}$ random variable.

###### Proof.

Let $X_{1},X_{2},\ldots$ be i.i.d. copies of $Z_{1}$ and let $Z_{0}$ be an
independent standard Gaussian. Then $\sum\sqrt{b_{j}}X_{j}+\sqrt{a-\sum
b_{j}}Z_{0}$ (which converges in $L_{2}$ if infinitely many $b_{j}$ are
positive) is the desired random variable. ∎

######  Lemma 9.

For $\lambda\in(0,1)$, let $g_{\lambda}$ be the density of
$\sqrt{\lambda}X_{1}+\sqrt{1-\lambda}X_{2}$, where $X_{1},X_{2}$ are
independent copies of $Z_{1}$. Then for every
$0<\lambda_{1}<\lambda_{2}<\frac{1}{2}$, the function
$g_{\lambda_{2}}-g_{\lambda_{1}}$ on $(0,+\infty)$ has exactly two zeros and
the sign pattern $+-+$.

###### Proof.

By a direct computation,

$g_{\lambda}(x)=\Big{(}x^{2}+\lambda(1-\lambda)(3-6x^{2}+x^{4})\Big{)}\frac{e^{-x^{2}/2}}{\sqrt{2\pi}}$

so $g_{\lambda_{2}}-g_{\lambda_{1}}$ has the same sign as
$(\lambda_{2}-\lambda_{1})(1-\lambda_{1}-\lambda_{2})(3-6x^{2}+x^{4})$. ∎

######  Lemma 10.

Let $X_{1},X_{2},\dots$ be i.i.d. copies of $Z_{1}$ and let $Y$ be a symmetric
random variable independent of the $X_{j}$. Then the function

$\Psi(b_{1},\dots,b_{n})=\mathbb{E}|\sqrt{b_{1}}X_{1}+\dots+\sqrt{b_{n}}X_{n}+Y|^{p}$

is Schur-concave on $[0+\infty)^{n}$.

###### Proof.

We use the technique of interlacing densities (see, e.g. [2] or [11]). Let
$h(x)=|x+1|^{p}+|x-1|^{p}$. It suffices to show that for every
$0<\lambda_{1}<\lambda_{2}<\frac{1}{2}$, we have

$\int_{0}^{\infty}h(x)(g_{\lambda_{2}(x)}-g_{\lambda_{1}}(x))\mathrm{d}x\geq
0,$

where $g_{\lambda}$ is as in Lemma 9. For arbitrary $\alpha,\beta$,
$\int(\alpha
x^{2}+\beta)(g_{\lambda_{2}(x)}-g_{\lambda_{1}}(x))\mathrm{d}x=0$, so the
desired inequality is equivalent to

$\int_{0}^{\infty}\tilde{h}(|x|)(g_{\lambda_{2}(x)}-g_{\lambda_{1}}(x))\mathrm{d}x\geq
0,$

with $\tilde{h}(x)=h(x)+\alpha x^{2}+\beta$. Let $x_{1},x_{2}$ be the zeros of
$g_{\lambda_{2}(x)}-g_{\lambda_{1}}(x)$. Choose $\alpha$ and $\beta$ such that
$\tilde{h}$ has zeros at $x_{1}$ and $x_{2}$. Since for $p\geq 3$,
$\tilde{h}(\sqrt{x})$ is convex on $(0,+\infty)$, $\tilde{h}$ on $(0,+\infty)$
has no other zeros and the sign pattern $+-+$. Thus the integrand is pointwise
nonnegative, hence the result. ∎

Theorem 5 is a rather straighforward consequence of the above lemma.

###### Proof of Theorem 5.

By a standard approximation argument of truncating the infinite product if
necessary, we can assume that only finitely many $b_{j}$ are nonzero, say
$b_{1},\dots,b_{n}>0$. By homogeneity, we can also assume that $\sum b_{j}=1$.
Set $a=1+v$ with $v\geq 0$. Then $X$ has the same distribution as
$\sum\sqrt{b_{j}}X_{j}+\sqrt{v}Z_{0}$, where $X_{1},X_{2},\dots$ are i.i.d.
copies of $Z_{1}$ and $Z_{0}$ is an independent standard Gaussian. Note that

$\sigma^{2}=\operatorname{Var}(X)=\sum
b_{j}\operatorname{Var}(X_{j})+v\operatorname{Var}(Z_{0})=\operatorname{Var}(X_{1})+v=3+v.$

By Lemma 10,

$\mathbb{E}|X|^{p}\leq\mathbb{E}\left|\sum_{j=1}^{n}\frac{1}{\sqrt{n}}X_{j}+\sqrt{v}Z_{0}\right|^{p}.$

Moreover, also by Lemma 10, the right hand side is nondecreasing as $n$
increases. Letting $n\to\infty$ and invoking the central limit theorem gives
the upper bound.

For the lower bound, again thanks to Lemma 10,

$\mathbb{E}|X|^{p}\geq\mathbb{E}|Z_{1}+\sqrt{v}Z_{0}|^{p}.$

The random variable $\sigma^{-1/2}(Z_{1}+\sqrt{v}Z_{0})$ has the same
distribution as $(1+2b)^{-1/2}Z_{b}$ with $b=\frac{1}{1+v}$, by comparing the
characteristic functions. By a direct computation (using the density of
$Z_{b}$),

$\mathbb{E}|(1+2b)^{-1/2}Z_{b}|^{p}=\frac{2^{p/2}\Gamma\left(\frac{p+1}{2}\right)}{\sqrt{\pi}}\frac{1+pb}{(1+2b)^{p/2}}$

and we check that the right hand side as a function of $b$ with $p\geq 3$
fixed is decreasing, which finishes the proof. ∎

##  5 Ultra sub-Gaussianity and related notions

Here we only consider real-valued random variables. Recall that after [10], a
random variable $X$ is called ultra sub-Gaussian if it is symmetric, has all
moments finite and the sequence $(a_{n})_{n=0}^{\infty}$ defined by $a_{0}=1$,
$a_{n}=\frac{\mathbb{E}X^{2n}}{\mathbb{E}G^{2n}}$, $n\geq 1$, is log-concave.
Note that then the function

$f_{X}(t)=\mathbb{E}\exp(\sqrt{t}X)=\sum_{n=0}^{\infty}\frac{\sqrt{t}^{2n}\mathbb{E}X^{2n}}{(2n)!}=\sum_{n=0}^{\infty}\frac{\mathbb{E}X^{2n}}{(2n-1)!!}\frac{t^{n}}{2^{n}n!}=\sum_{n=0}^{\infty}\frac{a_{n}}{n!}(t/2)^{n}$

is $C^{\infty}(\mathbb{R})$ and the inequalities $a_{n}^{2}\geq
a_{n-1}a_{n+1}$, $n\geq 1$, are equivalent to

$\big{(}f_{X}^{(n)}(0)\big{)}^{2}\geq f_{X}^{(n-1)}(0)f_{X}^{(n+1)}(0),\qquad
n\geq 1.$ (7)

After [4], an entire function $f\colon\mathbb{C}\to\mathbb{C}$ is called
strongly log-concave if for every $n\geq 0$, $f^{(n)}$ is either zero or $\log
f^{(n)}$ is concave on $\mathbb{R}_{+}$. Thus we shall say that a random
variable is strongly log-concave if it is symmetric,
$z\mapsto\mathbb{E}\exp(zX)$ defines an entire function and $f_{X}^{(n)}$ is
log-concave for every $n\geq 0$, where $f_{X}(t)=\mathbb{E}\exp(\sqrt{t}X)$.
Then, in particular, inequalities (7) hold and thus $X$ is also ultra sub-
Gaussian. Remarkably, the converse is also true by the following result of
Gurvits (see Proposition 4.2(ii) in [4]).

######  Theorem 11 (Gurvits, [4]).

If a nonnegative sequence $(a_{n})_{n=0}^{\infty}$ is log-concave, then the
power series $f(t)=\sum_{n=0}^{\infty}\frac{a_{n}}{n!}t^{n}$ defines a
$C^{\infty}$ function which is log-concave on $\mathbb{R}_{+}$.

Gurvits obtains this as a by-product of his characterisation of _shifts_
preserving log-concavity. He also indicates a very different, short argument
(see [4, 5]): by Shephard’s theorem (see [19]), if a finite sequence
$(a_{n})_{n=0}^{n_{0}}$ is log-concave, then there are simplices $K_{1},K_{2}$
in $\mathbb{R}^{N}$ such that $p_{N}(t)=\sum_{n=0}^{N}\binom{N}{n}a_{n}t^{n}$
is of the form $\text{vol}_{N}(K_{1}+tK_{2})$. By the Brunn-Minkowski
inequality, $p_{N}(t)^{1/N}$ is concave, which gives
$(p_{N}(t)^{\prime})^{2}\geq\frac{N}{N-1}p_{N}(t)\cdot
p_{N}(t)^{\prime\prime}$. Fix $s>0$, let $t=s/N$ with $N\to\infty$ to get
$(p^{\prime}(s))^{2}\geq p^{\prime}(s)p^{\prime\prime}(s)$ with
$p(s)=\sum_{n=0}^{n_{0}}\frac{a_{n}}{n!}s^{n}$, allowing to conclude the
result. We present a third, rather elementary proof.

###### Proof of Theorem 11.

Thanks to log-concavity, $(a_{n})$ is majorised by a geometric sequence, so
the radius of convergence of $\sum\frac{a_{n}}{n!}t^{n}$ is $\infty$, hence
$f$ is $C^{\infty}$ on $\mathbb{R}$.

We have $f^{(k)}(t)=\sum_{n\geq 0}\frac{a_{n+k}}{n!}t^{n}$. The log concavity
of $f$ is equivalent to $(f^{\prime})^{2}\geq f^{\prime\prime}f$, that is

$\left(\sum_{n\geq 0}\frac{a_{n+1}}{n!}t^{n}\right)^{2}-\left(\sum_{n\geq
0}\frac{a_{n}}{n!}t^{n}\right)\left(\sum_{n\geq
0}\frac{a_{n+2}}{n!}t^{n}\right)\geq 0.$

The left hand side is a power series with coefficient in front of $t^{n}$
equal to

$\sum_{k+l=n}\frac{a_{k+1}}{k!}\cdot\frac{a_{l+1}}{l!}-\sum_{k+l=n}\frac{a_{k}}{k!}\cdot\frac{a_{l+2}}{l!}.$

It is enough to show that this expression is non-negative. Multiplying both
sides by $n!$, we see that it is enough to prove the inequality

$\sum{n\choose k}a_{k+1}a_{n-k+1}-\sum{n\choose k}a_{k}a_{n-k+2}\geq 0$

with the usual convention that ${n\choose k}=0$ if $k<0$ or $k>n$. The sums
are over $\mathbb{Z}$. Shifting the summation index in the first sum we get an
equivalent form

$\sum\left({n\choose k-1}-{n\choose k}\right)a_{k}a_{n-k+2}\geq 0.$

For $n=0$ this is just $a_{1}^{2}\geq a_{0}a_{2}$ and for $n=1$, it is
$2a_{1}a_{2}\geq a_{0}a_{3}+a_{1}a_{2}$, which is just $a_{1}a_{2}\geq
a_{0}a_{3}$. By log-concavity, both hold. Thus in what follows we can assume
that $n\geq 2$.

Let $b_{k}=a_{k}a_{n-k+2}$ (we set $a_{k}=0$ for $k<0$). Our goal is to prove

$\sum\left({n\choose k-1}-{n\choose k}\right)b_{k}\geq 0.$

Note that $b_{k}=b_{n-k+2}$. From log concavity we know that the sequence
$b_{k}$ is nondecreasing, $b_{k-1}\leq b_{k}$ for $k\leq\frac{n}{2}+1$. Since
$b_{k}$ is symmetric about $k=\frac{n}{2}+1$, we can write the desired
inequality in an equivalent form

$\left({n\choose\frac{n}{2}}-{n\choose\frac{n}{2}-1}\right)b_{\frac{n}{2}+1}+\sum_{k<\frac{n}{2}+1}\left(2{n\choose
k-1}-{n\choose k}-{n\choose k-2}\right)b_{k}\geq 0$

with the first term present only if $n$ is even. To prove this, we set

$s_{k}=2{n\choose k-1}-{n\choose k}-{n\choose k-2},\qquad k<\frac{n}{2}+1$

and additionally when $n$ is even,

$s_{n/2+1}={n\choose\frac{n}{2}}-{n\choose\frac{n}{2}-1}$

and use the following elementary lemma (Lemma 6 in [10]).

######  Lemma 12 (Nayar-Oleszkiewicz, [10]).

If $b_{0}\leq b_{1}\leq\ldots\leq b_{m}$ and real numbers
$s_{0},s_{1},\ldots,s_{m}$ are such that $\sum_{k=0}^{m}s_{k}=0$ and the
sequence $(\operatorname{sgn}(s_{k}))$ is nondecreasing, then
$\sum_{k=0}^{m}s_{k}b_{k}\geq 0$.

The lemma is applied with $m=\lfloor n/2+1\rfloor$ and the sequence $(s_{k})$
defined above (which plainly satisfies $\sum_{k=0}^{m}s_{k}=0$). It remains to
check that $(sgn(s_{k}))$ is a nondecreasing sequence. For $n\geq 2$ the
inequality $s_{k}\leq 0$ is equivalent to

$-4k^{2}+4kn+8k-n^{2}-3n-2\leq 0.$

The maximum of this parabola is attained for $k=\frac{n}{2}+1$ and thus for
$k\leq\frac{n}{2}+1$ this quadratic function is increasing. This shows that
$(\operatorname{sgn}(s_{k}))$ is non-decreasing, which finishes the proof. ∎

This gives that the classes of strongly log-concave and ultra sub-Gaussian
random variables are the same. As observed earlier in Theorem 7, type
$\mathscr{L}$ random variables are ultra sub-Gaussian. The converse is false:
it is enough to consider $X$ with $3$ atoms:
$\mathbb{P}\left(X=0\right)=1-\beta$, $\mathbb{P}\left(X=\pm
1\right)=\beta/2$. Then $\mathbb{E}\exp(zX)=(1-\beta)+\beta\cosh(z)$. The
equation $(1-\beta)+\beta\cosh(z)=0$ has only purely imaginary zeroes only if
$\beta>1/2$. On the other hand, this random variables is ultra sub-Gaussian
for any $\beta\in(0,1)$. Another examples are random variables with densities
proportional to $e^{-|x|^{\alpha}}$: they are ultra sub-Gaussian for every
$\alpha\geq 2$ (see Corollary 2 in [10]), whereas they are of type
$\mathscr{L}$ only for integral $\alpha$ (see Section 3). This discussion is
summarised in Figure 1.

USGSLCtype $\mathscr{L}$// Figure 1: Implications between ultra sub-Gaussian
(USG), strongly log-concave (SLC) and type $\mathscr{L}$ random variables.

## References

  * [1] Dimitrov, D. K., Rusev, P. K., Zeros of entire Fourier transforms. East J. Approx. 17 (2011), no. 1, 1–110.
  * [2] Eskenazis, A., Nayar, P., Tkocz, T., Sharp comparison of moments and the log-concave moment problem. Adv. Math. 334 (2018), 389–416.
  * [3] Griffiths, R., Correlations in Ising Ferromagnets. II. External Magnetic Fields. J. Math. Phys. 8, 484 (1967).
  * [4] Gurvits, L., On multivariate Newton-like inequalities. Advances in combinatorial mathematics, 61–78, Springer, Berlin, 2009.
  * [5] Gurvits, L., A short proof, based on mixed volumes, of Liggett’s theorem on the convolution of ultra-logconcave sequences. Electron. J. Combin. 16 (2009), no. 1, Note 5, 5 pp.
  * [6] Havrilla, A., Tkocz, T., Sharp Khinchin-type inequalities for symmetric discrete uniform random variables. Preprint (2019), arXiv:1912.13345, to appear in Israel J. Math.
  * [7] Ilieff, L., Über trigonometrische Integrale, welche ganze Funktionen mit nur reellen Nullstellen darstellen. Acta Math. Acad. Sci. Hungar. 6 (1955), 191–194.
  * [8] Jury, E. I., Theory and Application of the z-Transform Method, John Wiley and Sons, New York, NY, 1964.
  * [9] Liggett, T., Ultra logconcave sequences and negative dependence. J. Combin. Theory Ser. A 79 (1997), no. 2, 315–325.
  * [10] Nayar, P., Oleszkiewicz, K., Khinchine type inequalities with optimal constants via ultra log-concavity. Positivity 16 (2012), no. 2, 359–371.
  * [11] Nayar, P., Zwara, S., Sharp variance-entropy comparison for nonnegative gaussian quadratic forms. Preprint (2020), arXiv:2005.11705.
  * [12] Newman, C. M., Zeros of the partition function for generalized Ising systems. Comm. Pure Appl. Math. 27 (1974), 143–159.
  * [13] Newman, C. M., Inequalities for Ising models and field theories which obey the Lee-Yang theorem. Comm. Math. Phys. 41 (1975), 1–9.
  * [14] Newman, C. M., An extension of Khintchine’s inequality. Bull. Amer. Math. Soc. 81 (1975), no. 5, 913–915.
  * [15] Newman, C, M., Moment inequalities for ferromagnetic Gibbs distributions. J. Math. Phys. 16 (1975), no. 9, 1956–1959.
  * [16] Pólya, G., Szegö, G., Problems and theorems in analysis. I. Series, integral calculus, theory of functions. Translated from the German by Dorothee Aeppli. Reprint of the 1978 English translation. Classics in Mathematics. Springer-Verlag, Berlin, 1998.
  * [17] Pólya, G., Szegö, G., Problems and theorems in analysis. II. Theory of functions, zeros, polynomials, determinants, number theory, geometry. Translated from the German by C. E. Billigheimer. Reprint of the 1976 English translation. Classics in Mathematics. Springer-Verlag, Berlin, 1998.
  * [18] Sinclair, C. D., Vaaler, J. D., Self-inversive polynomials with all zeros on the unit circle. Number theory and polynomials, 312–321, London Math. Soc. Lecture Note Ser., 352, Cambridge Univ. Press, Cambridge, 2008.
  * [19] Shephard, G. C., Inequalities between mixed volumes of convex sets. Mathematika 7 (1960), 125–138.
  * [20] Stanley, R., Log-concave and unimodal sequences in algebra, combinatorics, and geometry. Graph theory and its applications: East and West (Jinan, 1986), 500–535, Ann. New York Acad. Sci., 576, New York Acad. Sci., New York, 1989.
  * [21] Szarek, S., On the best constant in the Khintchine inequality. Stud. Math. 58, 197–208 (1976).
  * [22] Walkup, D., Pólya sequences, binomial convolution and the union of random sets. J. Appl. Probability 13 (1976), no. 1, 76–85.