Scalable expectation propagation
for generalized linear models

Niccolò Anceschi1<EMAIL_ADDRESS>
Augusto Fasano2<EMAIL_ADDRESS>
Beatrice Franzolini3<EMAIL_ADDRESS>
Giovanni Rebaudo4<EMAIL_ADDRESS>1Duke University, Durham, USA

2Catholic University of the Sacred Heart, Milan, Italy

3Bocconi University, Milan, Italy

4University of Torino, Turin, Italy

Abstract

Generalized linear models (glms) arguably represent the standard approach for
statistical regression beyond the Gaussian likelihood scenario. When Bayesian
formulations are employed, the general absence of a tractable posterior
distribution has motivated the development of deterministic approximations,
which are generally more scalable than sampling techniques. Among them,
expectation propagation (ep) showed extreme accuracy, usually higher than many
variational Bayes solutions. However, the higher computational cost of ep
posed concerns about its practical feasibility, especially in high-dimensional
settings. We address these concerns by deriving a novel efficient formulation
of ep for glms, whose cost scales linearly in the number of covariates $p$.
This reduces the state-of-the-art $\mathcal{O}(p^{2}n)$ per-iteration
computational cost of the ep routine for glms to
$\mathcal{O}(pn\cdot\min\\{p,n\\})$, with $n$ being the sample size. We also
show that, for binary models and log-linear glms approximate predictive means
can be obtained at no additional cost. To preserve efficient moment matching
for count data, we propose employing a combination of log-normal Laplace
transform approximations, avoiding numerical integration. These novel results
open the possibility of employing ep in settings that were believed to be
practically impossible. Improvements over state-of-the-art approaches are
illustrated both for simulated and real data. The efficient ep implementation
is available at https://github.com/niccoloanceschi/EPglm.

Key Words: Generalized Linear Model, Expectation Propagation, Bayesian
Inference, Extended Multivariate Skew-Normal Distribution, Unified Skew-Normal
Distribution.

Short/Running Title: Scalable EP for GLMs.

##  1 Introduction

Generalized linear models (glms) arguably represent the most popular off-the-
shelf tool for data regression for a wide variety of data, with binary and
count responses representing their main application areas. The lack of a
tractable posterior distribution in their Bayesian formulation has motivated
notable efforts to develop computational methods that can efficiently scale to
high-dimensional settings (large $p$) or large datasets (large $n$), with the
former providing a primary area of interest. See, for instance, the review
paper by [12] about the big $p$ problem faced by the available computational
methods for Bayesian binary regression computations. Among the scalable
deterministic approximate methods, a central role has certainly been played by
variational Bayes (vb) approximations [9, e.g.,], which are usually much
faster than sampling-based techniques. This comes however at the cost of
possible accuracy issues in certain settings, making some of the vb
approximations a sound replacement for the exact posterior only in some
scenarios. This is sub-par compared to having an approximation scheme offering
reliable performances regardless of the specific setting considered.
Empirically, expectation propagation (ep) [37, 39] has demonstrated an
improved quality in posterior approximation for a wide variety of models,
including binary [12], tobit [2], Poisson [53], heteroscedastic, lasso-
penalized, and quantile regression [60], as well as for hierarchical models
more in general [51]. Indeed, when testing its performance with logit
likelihood and various priors on datasets of moderate size, the empirical
approximation error turned out to be so small that [12] conclude that it seems
that ep may be safely used as a complete replacement of sampling-based methods
(Section 5.1.1). However, such improved accuracy comes at a higher
computational cost than vb methods, thus motivating questions about the
applicability of this excellent approximate method in large $p$ settings,
namely $p\gg 1000$ [12]. To address this issue we develop a novel ep
formulation for Bayesian glms with multivariate Gaussian prior which scales
linearly in the number of covariates $p$. This improves over current state-of-
the-art implementations by multiple orders of magnitude and makes ep a
feasible option in glms with dimensions well beyond the above threshold,
paving the way to the employment of ep in settings where available sampling
methods are impractical and other available approximations are less accurate.

In glms, the effect of the $p$-dimensional vector of predictors ${\bf
x}=(x_{1},\ldots,x_{p})^{\intercal}$ on the distribution of the response
variable $y$ is modelled through a link function $g(\cdot)$, which relates the
linear predictor $\eta={\bf x}^{\intercal}{\bm{\beta}}$ to the mean of the
response variable via the relation
${\mathbb{E}}[y\mid{\bm{\beta}}]=g^{-1}(\eta)$, with ${\bm{\beta}}$ a
$p$-dimensional unknown parameter. For binary regression models with Bernoulli
likelihood, classical link functions are represented by the logit or probit
links, while for count data default specifications combine a Poisson or
Negative Binomial likelihood with $\log(\cdot)$ link function. Other models
have also been studied, including, for instance, Gamma regression for non-
negative continuous responses, where a Gamma likelihood is usually combined
with a $\log(\cdot)$ link function. When considering their Bayesian
formulation, the model is completed by specifying a prior distribution for the
parameter ${\bm{\beta}}$, with the zero-mean spherical multivariate Gaussian
distribution constituting a default choice. Unfortunately, a tractable closed-
form posterior distribution for the glms is generally lacking beyond the
Gaussian regression.

State-of-the-art Bayesian approaches usually rely on posterior sampling
methods — see, e.g., [31, 41, 18, 54, 52] and the seminal contribution by [1]
for binary regression models, or [24, 11, 15] for the Poisson-like models — or
more scalable approximate posterior inference — see, e.g., [32, 13, 19, 20]
for binary regression models, or [4, 53] for Poisson regression. Most of the
above-mentioned sampling methods rely on mcmc schemes, often introducing a
data augmentation step, in the spirit of the well-known Gibbs sampler for the
probit model originally developed by [1], who exploited the augmented
variables to obtain conjugate full conditionals. These approaches make in
principle posterior inference feasible for a wide variety of models, but the
possible high correlation between the augmented data and the parameters can
cause convergence issues and thus unreliable results, especially in unbalanced
datasets [33]. Other sampling alternatives, like the no-U-turn Hamiltonian
sampler by [30] or adaptive Metropolis-Hastings techniques [29, 42], are
possible. Although these approaches might mitigate autocorrelation among
samples and improve mixing, they remain computationally impractical for high-
dimensional problems [18, 20]. A notable exception is represented by [18], who
developed an i.i.d. sampler for the posterior distribution in Bayesian probit
models. This leverages on the representation of the $\textsc{sun}_{p,n}$
posterior [3] as a linear combination of a $p$-variate Gaussian and an
$n$-variate truncated Gaussian, making it particularly well suited for high-
dimensional problems. However, the need to sample from truncated multivariate
Gaussians whose dimension equals the number of the observations represents the
computational challenge of such an algorithm: despite recent promising results
on the topic [40, 10, 27, 54], i.i.d. posterior samples are hard to get when
$n$ is above a few hundred.

More scalable deterministic approximations were thus developed to circumvent
the computational issues mentioned above. However, in high-dimensional
settings, some of the most popular approximate methods, like mean-field vb or
Laplace approximation, may turn out to be unreliable in high dimensions [20,
23, 2]. In binary and multinomial Bayesian regression, the above contributions
overcame such drawbacks by developing a partially-factorized variational (pfm-
vb) approximation of the posterior distribution which, for any fixed $n$,
converges to the true posterior distribution as $p$ diverges under reasonable
assumptions. Moreover, pfm-vb is always guaranteed to improve over mean-field
approximations [13, 28]. The computational costs are orders of magnitude
smaller than the ones of sampling methods, and, empirically, the approximate
posterior moments closely match the ones obtained via i.i.d. sampling for
$p/n\geq 2$. Other variational approximations have been developed for the
logit [32] and Poisson [4] to get fast posterior approximations. See also [12]
for a complete review of computational methods for Bayesian binary regression,
covering also Laplace and improved nested Laplace approximation (inla) [45],
as well as other exact methods whose review is beyond the scope of the paper.
The overall accuracy of the vb approaches is still an open area of research
and in many cases, it is very setting-specific. In general, ep provides better
posterior approximations across scenarios for glms [12, 2, 53]. This comes,
however, at a higher computational cost, which makes current implementations
not feasible for settings where $p>1000$. We show that an efficient derivation
of ep is possible, having linear cost in $p$ and reducing the
$\mathcal{O}(p^{2}n)$ per-iteration cost reported in [12] to
$\mathcal{O}(pn\cdot\min\\{p,n\\})$. This improvement makes ep feasible in
settings that are out of reach for state-of-the-art implementations. The
resulting efficient implementation can be fully leveraged to provide accurate
approximations of the marginal likelihood, which is often overlooked in the
literature. Besides such structural enhancement, we achieve further
computational improvement for count data, by proposing a novel approximate
moment matching step for Poisson regression. Similarly to standard practice in
the logit case, this avoids the need for numerical integration, benefiting
both in terms of efficiency and stability.

The rest of the article is organized as follows. Section 2 reviews the ep
algorithm and its employment to compute the approximate posterior moments and
marginal likelihood. Section 3 details our efficient ep routines for various
glms, while Section 4 shows multiple approximations of the hybrid moments
required by the ep algorithm, avoiding numerical quadrature. We also show that
in several relevant cases, the required hybrid moments and/or the predictive
means are available in closed form. Section 5 reports the results of our
proposed algorithm over simulated data, while Section 6 presents the
performance of the proposed method over real data applications. Section 7
contains concluding remarks and further interesting research questions arising
from our findings. The code to reproduce the results is available at
https://github.com/niccoloanceschi/EPglm.

##  2 Expectation propagation

In this section, we briefly review the principles and the main steps of
expectation propagation (ep) [37, 39] for Bayesian inference in a generic
regression model

$y_{i}\mid{\bm{\beta}}\overset{ind}{\sim}p(y_{i}\mid{\bm{\beta}})\quad
i=1,\ldots,n,\quad{\bm{\beta}}\sim p({\bm{\beta}}),$ (1)

where the likelihood contribution $p(y_{i}\mid{\bm{\beta}})$ depends on the
unknown parameter ${\bm{\beta}}$ and a known covariate vector ${\bf
x}_{i}\in{\mathbb{R}}^{p}$, which is thus removed from the conditioning part.
Throughout the paper, we will assume univariate $y_{i}$’s and the well-
established zero-mean spherical multivariate Gaussian prior for the parameter
${\bm{\beta}}$, although computations trivially adapt to more general cases.
We will then take $p({\bm{\beta}})=\phi_{p}({\bm{\beta}};\nu^{2}{\bf I}_{p})$,
where ${\bf I}_{p}$ is the identity matrix of dimension $p$ and $\phi_{p}({\bf
t};{\bf S})$ denotes the density of a $p$-variate Gaussian random variable
with mean $\bm{0}_{p}$ and covariance matrix ${\bf S}$, evaluated at ${\bf
t}$.

###  2.1 Approximation of posterior moments

Calling ${\bf y}=(y_{1},\ldots,y_{n})$, the ep framework approximates the
posterior distribution

$p({\bm{\beta}}\mid{\bf y})\propto
p({\bm{\beta}})\prod_{i=1}^{n}p(y_{i}\mid{\bm{\beta}})$

with a density $q_{{\textsc{ep}}}({\bm{\beta}})$ having the same factorization
form

$q_{{\textsc{ep}}}({\bm{\beta}})\propto
q_{0}({\bm{\beta}})\prod_{i=1}^{n}q_{i}({\bm{\beta}}),$

where the global approximation $q_{{\textsc{ep}}}({\bm{\beta}})$ and the
factors (also referred to as sites) $q_{i}({\bm{\beta}})$, $i=0,\ldots,n$, are
taken from a possibly unnormalized exponential family, as clarified below. The
key operating principle of the ep scheme is to refine iteratively each
approximate site $q_{i}({\bm{\beta}})$ against the exact one
$p(y_{i}\mid{\bm{\beta}})$, exploiting the information coming from the
remaining factors $\\{q_{j}({\bm{\beta}})\\}_{j\neq i}$. This helps in
targeting regions of high posterior density and, thus, constructing a globally
accurate approximation $q_{{\textsc{ep}}}({\bm{\beta}})$.

When the parameter ${\bm{\beta}}$ is continuous, as in glms, a typical choice,
which will also be used hereafter, is to take Gaussian-like sites
$q_{i}({\bm{\beta}})\propto\exp\\{-\frac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}_{i}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf r}_{i}\\}$ for $i=0,\ldots,n$,
with the natural parameters ${\bf r}_{i}$ and ${\bf Q}_{i}$ to be determined
by the ep algorithm optimization rules. As a consequence, the global
approximation $q_{{\textsc{ep}}}({\bm{\beta}})$ is Gaussian, with
$q_{{\textsc{ep}}}({\bm{\beta}})=\phi_{p}({\bm{\beta}}-{\bm{\xi}}_{{\textsc{ep}}},{\bm{\Omega}}_{{\textsc{ep}}})$,
where ${\bm{\xi}}_{{\textsc{ep}}}={\bf Q}_{{\textsc{ep}}}^{-1}{\bf
r}_{{\textsc{ep}}}$ and ${\bm{\Omega}}_{{\textsc{ep}}}={\bf
Q}_{{\textsc{ep}}}^{-1}$ are obtained from the global natural parameters ${\bf
r}_{{\textsc{ep}}}=\sum_{i=0}^{n}{\bf r}_{i}$, ${\bf
Q}_{{\textsc{ep}}}=\sum_{i=0}^{n}{\bf Q}_{i}$. If the prior $p({\bm{\beta}})$
has the same form of $q_{0}({\bm{\beta}})$, i.e. $p({\bm{\beta}})$ is a
multivariate Gaussian density, then we can set
$q_{0}({\bm{\beta}})=p({\bm{\beta}})$ and compute only the natural parameters
${\bf r}_{i}$ and ${\bf Q}_{i}$ for sites $i=1,\ldots,n$. It is worth noting
that the sites $q_{i}({\bm{\beta}})$ provide an approximation of the
likelihood term $p(y_{i}\mid{\bm{\beta}})$ and so $q_{i}({\bm{\beta}})$ is not
a probability density in ${\bm{\beta}}$. Hence, it generally does not hold
that $\int q_{i}({\bm{\beta}})d{\bm{\beta}}=1$. Likewise, the contributions
${\bf Q}_{i}$’s of each site $i=1,\dots,n$ to the global precision matrix
${\bf Q}_{{\textsc{ep}}}$ are not explicitly constrained to be positive
definite. ${\bf Q}_{{\textsc{ep}}}$ must instead be positive definite, since
the global approximation $q_{{\textsc{ep}}}({\bm{\beta}})$, representing the
ep surrogate of the exact posterior distribution, is a valid probability
density in ${\bm{\beta}}$.

In principle, one could think of targeting a global approximation
$q_{{\textsc{ep}}}({\bm{\beta}})$ minimizing the Kullback-Leibler divergence
[35] $\textsc{kl}\big{[}p({\bm{\beta}}\mid{\bf
y})||q_{{\textsc{ep}}}({\bm{\beta}})\big{]}$. However, this would involve the
computation of integrals with respect to the posterior distribution, which is
usually intractable or computationally impractical. ep overcomes this issue by
iteratively updating the natural parameters ${\bf r}_{i}$ and ${\bf Q}_{i}$ of
site $i$ minimizing a different kl divergence, informed by the remaining
$\\{{\bf r}_{j}\\}_{j\neq i}$ and $\\{{\bf Q}_{i}\\}_{j\neq i}$, which instead
is computationally tractable. More specifically, at each iteration ep updates
each site $i=1,\ldots,n$ as follows. Keeping the parameters for all sites
$j\neq i$ fixed, first compute the so-called cavity density

$q_{-i}({\bm{\beta}})\propto\dfrac{q_{{\textsc{ep}}}({\bm{\beta}})}{q_{i}({\bm{\beta}})}\propto\exp\left\\{-\frac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}_{-i}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf r}_{-i}\right\\},$

where ${\bf r}_{-i}=\sum_{j\neq i}{\bf r}_{j}$ and ${\bf Q}_{-i}=\sum_{j\neq
i}{\bf Q}_{j}$. Contrary to the local sites, $q_{-i}({\bm{\beta}})$ is a
proper distribution in ${\bm{\beta}}$. As such, defining $\Psi({\bf r},{\bf
Q})=\int\exp\left\\{-\dfrac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf r}\right\\}d{\bm{\beta}}$ so that
$\log\Psi({\bf r},{\bf Q})=\dfrac{1}{2}{\bf r}^{\intercal}{\bf Q}^{-1}{\bf
r}+\dfrac{p}{2}\log(2\pi)-\dfrac{1}{2}\log|{\bf Q}|$, one gets

$q_{-i}({\bm{\beta}})=\dfrac{1}{\Psi({\bf r}_{-i},{\bf
Q}_{-i})}\exp\left\\{-\frac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}_{-i}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf
r}_{-i}\right\\}=\phi_{p}({\bm{\beta}}-{\bf Q}_{-i}^{-1}{\bf r}_{-i},{\bf
Q}_{-i}^{-1}).$

Similarly, one defines the so-called hybrid (or tilted) distribution

$h_{i}({\bm{\beta}})=\dfrac{1}{Z_{h_{i}}}p(y_{i}\mid{\bm{\beta}})q_{-i}({\bm{\beta}}),$

with $Z_{h_{i}}=\int
p(y_{i}\mid{\bm{\beta}})q_{-i}({\bm{\beta}})d{\bm{\beta}}$ and update the
parameters ${\bf r}_{i}$ and ${\bf Q}_{i}$ of site $i$ to ${\bf
r}_{i}^{{\textsc{new}}}$ and ${\bf Q}_{i}^{{\textsc{new}}}$ so that the
resulting global approximation
$q_{{\textsc{ep}}}^{{\textsc{new}}}({\bm{\beta}})\propto
q_{i}^{{\textsc{new}}}({\bm{\beta}})q_{-i}({\bm{\beta}})$ minimizes
$\textsc{kl}\big{[}h_{i}({\bm{\beta}})||q_{{\textsc{ep}}}^{{\textsc{new}}}({\bm{\beta}})\big{]}$.
When the $q_{i}({\bm{\beta}})$ are taken from an exponential family, the
latter minimization is achieved by matching the moments of the sufficient
statistics of $q_{{\textsc{ep}}}^{{\textsc{new}}}({\bm{\beta}})$ computed with
respect to $h_{i}({\bm{\beta}})$ and
$q_{{\textsc{ep}}}^{{\textsc{new}}}({\bm{\beta}})$ [46, 8]. Calling
${\bm{\mu}}_{h_{i}}={\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]$ and
${\bm{\Sigma}}_{h_{i}}={\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]$,
the ep moment-matching conditions for site $i$ translate into imposing the
equality

$\begin{cases}\left({\bf Q}_{-i}+{\bf Q}_{i}^{{\textsc{new}}}\right)^{-1}({\bf
r}_{-i}+{\bf r}_{i}^{{\textsc{new}}})={\bm{\mu}}_{h_{i}}\\\ \left({\bf
Q}_{-i}+{\bf
Q}_{i}^{{\textsc{new}}}\right)^{-1}={\bm{\Sigma}}_{h_{i}},\end{cases}$

which implies

$\begin{cases}{\bf r}_{i}^{{\textsc{new}}}=\left({\bf Q}_{-i}+{\bf
Q}_{i}^{{\textsc{new}}}\right){\bm{\mu}}_{h_{i}}-{\bf r}_{-i}\\\ {\bf
Q}_{i}^{{\textsc{new}}}={\bm{\Sigma}}_{h_{i}}^{-1}-{\bf Q}_{-i}.\end{cases}$
(2)

At each iteration, the moment-matching conditions (2) are enforced
sequentially for each site $i=1,\ldots,n$, until convergence is eventually
reached. It is well-known that ep might be prone to converge issues in some
settings [38]. When this happens, damping techniques can be applied to reach
convergence, although usually at a lower rate [51]. After convergence is
reached, one can compute the ep global approximation
$q_{{\textsc{ep}}}({\bm{\beta}})=\phi_{p}({\bm{\beta}}-{\bm{\xi}}_{{\textsc{ep}}},{\bm{\Omega}}_{{\textsc{ep}}})$,
where ${\bm{\xi}}_{{\textsc{ep}}}={\bf Q}_{{\textsc{ep}}}^{-1}{\bf
r}_{{\textsc{ep}}}$ and ${\bm{\Omega}}_{{\textsc{ep}}}={\bf
Q}_{{\textsc{ep}}}^{-1}$ and ${\bf r}_{{\textsc{ep}}}=\sum_{i=0}^{n}{\bf
r}_{i}$, ${\bf Q}_{{\textsc{ep}}}=\sum_{i=0}^{n}{\bf Q}_{i}$.

To completely specify the ep algorithm, initial values of the site parameters
need to be provided. Following common practice employed in the case of
multivariate Gaussian prior $p({\bm{\beta}})=\phi_{p}({\bf m},{\bf C})$ [51],
we consider the case in which $q_{{\textsc{ep}}}({\bm{\beta}})$ is initialized
to the prior distribution, that is we initialize ${\bf r}_{0}={\bf C}^{-1}{\bf
m}$, ${\bf Q}_{0}={\bf C}^{-1}$, while ${\bf r}_{i}=\bm{0}$, ${\bf
Q}_{i}=\bm{0}$ for $i=1,\ldots,n$.

###  2.2 Approximation of the marginal likelihood

ep can be used as well to obtain an estimation of the marginal likelihood
$p({\bf
y})=\int\prod_{i=1}^{n}p(y_{i}\mid{\bm{\beta}})p({\bm{\beta}})d{\bm{\beta}}$
at no additional cost [51, 37]. To this end, it becomes necessary to specify a
proportionality constant $Z_{i}$ for each site approximation
$q_{i}({\bm{\beta}})$, so that

$q_{i}({\bm{\beta}})=\dfrac{1}{Z_{i}}\exp\left\\{-\dfrac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}_{i}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf r}_{i}\right\\}.$

Such $Z_{i}$’s were not strictly required to compute the moments of the global
posterior approximation $q_{{\textsc{ep}}}({\bm{\beta}})$. It turns out that
the former task is straightforward when the normalizing constant of the hybrid
distribution $Z_{h_{i}}$ can be computed in closed form (as is the case for
the probit model). Since $q_{0}({\bm{\beta}})=p({\bm{\beta}})$, it holds $\log
Z_{0}=\log\Psi({\bf r}_{0},{\bf Q}_{0})$. Recall that each
$q_{i}({\bm{\beta}})$, for $i=1,\ldots,n$, is not a probability density in
${\bm{\beta}}$, but it is rather an approximation of the likelihood term
$p(y_{i}\mid{\bm{\beta}})$. Adhering to such rigorous representation in terms
of global normalization constant $Z_{{\textsc{ep}}}$ and local proportionality
constants $Z_{i}$, the overall ep approximation can be rewritten as
$q_{{\textsc{ep}}}^{{\textsc{new}}}({\bm{\beta}})=\dfrac{1}{Z_{\textsc{ep}}^{{\textsc{new}}}}q_{i}^{{\textsc{new}}}({\bm{\beta}})q_{-i}({\bm{\beta}})$,
where, analogously to $Z_{h_{i}}$, one has

$\begin{split}Z_{\textsc{ep}}^{{\textsc{new}}}&=\int
q_{-i}({\bm{\beta}})q_{i}^{{\textsc{new}}}({\bm{\beta}})d{\bm{\beta}}=\dfrac{1}{\Psi({\bf
r}_{-i},{\bf Q}_{-i})\,Z_{i}^{{\textsc{new}}}}\int
e^{-\frac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}_{\textsc{ep}}^{{\textsc{new}}}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf
r}_{\textsc{ep}}^{{\textsc{new}}}}d{\bm{\beta}}\\\ &=\dfrac{\Psi({\bf
r}_{\textsc{ep}}^{{\textsc{new}}},{\bf
Q}_{\textsc{ep}}^{{\textsc{new}}})}{\Psi({\bf r}_{-i},{\bf
Q}_{-i})\,Z_{i}^{{\textsc{new}}}}.\end{split}$

Such reformulation is not immediately useful per se, since the the values of
these $Z_{i}^{{\textsc{new}}}$’s are not determined by the ep moment-matching
conditions. Accordingly, an additional constraint must be fixed to uniquely
identify them, and, in turn, to produce an approximation of the marginal
likelihood. This can be done by extending the moment matching procedure to the
so-called zeroth order moments, imposing that
$Z_{\textsc{ep}}^{{\textsc{new}}}=Z_{h_{i}}$. Intuitively, this can be
interpreted as a way of enforcing $q_{i}^{{\textsc{new}}}({\bm{\beta}})$ to be
the best approximation of $p(y_{i}\mid{\bm{\beta}})$, such that $\int
q_{-i}({\bm{\beta}})q_{i}^{{\textsc{new}}}({\bm{\beta}})d{\bm{\beta}}=\int
q_{-i}({\bm{\beta}})p(y_{i}\mid{\bm{\beta}})d{\bm{\beta}}$. Accordingly, it is
trivial to see that

$\log Z_{i}^{{\textsc{new}}}=\log\Psi({\bf
r}_{{\textsc{ep}}}^{{\textsc{new}}},{\bf
Q}_{{\textsc{ep}}}^{{\textsc{new}}})-\log\Psi({\bf r}_{-i},{\bf Q}_{-i})-\log
Z_{h_{i}}.$ (3)

Upon convergence of the main ep routine, the ep approximation of the marginal
likelihood can be computed by replacing each $p(y_{i}\mid{\bm{\beta}})$ with
the corresponding approximate site the factor $q_{i}({\bm{\beta}})$, so that

$m_{{\textsc{ep}}}({\bf
y})=\int\prod_{i=1}^{n}q_{i}({\bm{\beta}})q_{0}({\bm{\beta}})d{\bm{\beta}}=\dfrac{\int\exp\left\\{{\textstyle-\dfrac{1}{2}{\bm{\beta}}^{\intercal}{\bf
Q}_{{\textsc{ep}}}{\bm{\beta}}+{\bm{\beta}}^{\intercal}{\bf
r}_{{\textsc{ep}}}}\right\\}d{\bm{\beta}}}{\Psi({\bf r}_{0},{\bf
Q}_{0})\prod_{i=1}^{n}Z_{i}}=\dfrac{\Psi({\bf r}_{{\textsc{ep}}},{\bf
Q}_{{\textsc{ep}}})}{\Psi({\bf r}_{0},{\bf Q}_{0})\prod_{i=1}^{n}Z_{i}}.$

Equivalently, its logarithm reads

$\log m_{{\textsc{ep}}}({\bf y})=\log\Psi({\bf r}_{{\textsc{ep}}},{\bf
Q}_{{\textsc{ep}}})-\log\Psi({\bf r}_{0},{\bf Q}_{0})-\sum_{i=1}^{n}\log
Z_{i}.$ (4)

##  3 Efficient updates in GLMs

###  3.1 Posterior moments and marginal likelihood with linear cost in $n$

The ep routine presented above applies regardless of the details of the
specific $p(y_{i}\mid{\bm{\beta}})$ considered. However, the implementation of
such a scheme for glms greatly benefits from a general structure common to
most of them. This is because the associated likelihood contributions depend
on the parameters ${\bm{\beta}}$ only thorough the inner product ${\bf
x}_{i}^{\intercal}{\bm{\beta}}$ so that
$p(y_{i}\mid{\bm{\beta}})=\ell_{i}({\bf x}_{i}^{\intercal}{\bm{\beta}})$. Here
the subscript $i$ conveys the dependence of $\ell_{i}$ on the response $y_{i}$
and possibly on additional nuisance parameters $\upsilon_{i}$ of the glm,
although this is irrelevant to the structure we want to elicit in the ep
approximate (i.e., $\ell_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})=p(y_{i}\mid{\bf
x}_{i}^{\intercal}{\bm{\beta}},\upsilon_{i})$). In the remainder, possible
additional parameters $\upsilon_{i}$ are considered fixed constants and thus
omitted from the conditioning.

The crucial resulting feature is that the corresponding approximate Gaussian
site depends on ${\bm{\beta}}$ only via ${\bf x}_{i}^{\intercal}{\bm{\beta}}$
as well [49, 48]. This can be proved by studying the moment-generating
function (mgf) of the hybrid distribution, as detailed in Appendix S.1.
Defining
${\bm{\Omega}}_{i}={\operatorname{var}}_{q_{-i}({\bm{\beta}})}[{\bm{\beta}}]={\bf
Q}_{-i}^{-1}$ and
${\bm{\xi}}_{i}={\mathbb{E}}_{q_{-i}({\bm{\beta}})}[{\bm{\beta}}]={\bf
Q}_{-i}^{-1}{\bf r}_{-i}$, this entails a rank one update of the global ep
parameters since the first two moments of the hybrid distribution take the
following form.

######  Lemma 1.

For any model of the form (1) with $p(y_{i}\mid{\bm{\beta}})=\ell_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})$, the moments of the hybrid distribution
$h_{i}({\bm{\beta}})$ needed for the ep updates (2) take the form

$\begin{split}{\bm{\mu}}_{h_{i}}&={\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf
x}_{i}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}\big{)}\\\
{\bm{\Sigma}}_{h_{i}}&={\bm{\Omega}}_{i}+{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}\big{)}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\;.\end{split}$

This reduces the dimensionality of the model-dependent moment matching since
the latter entails solely the moments induced by the hybrid
$h_{i}({\bm{\beta}})$ on the linear subspace identified by the projection on
${\bf x}_{i}$ [60]. Crucially, this means that the update of the $i$-th site
via conditions (2) does not require the direct computation of
${\bm{\Sigma}}_{h_{i}}^{-1}$. Indeed, the natural parameters associated with
the resulting updated site $q_{i}({\bm{\beta}})$ are parametrized by scalar
quantities, as clarified in the following proposition.

######  Proposition 1.

For any model of the form (1) with $p(y_{i}\mid{\bm{\beta}})=\ell_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})$, the updated parameters of the $i^{th}$ ep
site, $i=1,\ldots,n$, have the form ${\bf
r}_{i}^{{\textsc{new}}}=m_{i}^{{\textsc{new}}}{\bf x}_{i}$ and ${\bf
Q}_{i}^{{\textsc{new}}}=k_{i}^{{\textsc{new}}}{\bf x}_{i}{\bf
x}_{i}^{\intercal}$, with $m_{i}^{{\textsc{new}}}$ and
$k_{i}^{{\textsc{new}}}$ scalar quantities defined as

$\begin{split}k_{i}^{{\textsc{new}}}&=\operatorname{var}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]^{-1}-({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}\\\\[2.0pt]
m_{i}^{{\textsc{new}}}&=\operatorname{var}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]^{-1}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
r}_{-i}\;.\end{split}$

Hence, one can implement ep without storing and updating the $p\times p$
matrices ${\bf Q}_{i}$, but rather storing only scalar quantities $k_{i}$ and
$m_{i}$ within

$q_{i}({\bm{\beta}})\propto\exp\Big{\\{}-\frac{1}{2}k_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})^{2}+m_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})\Big{\\}}\;,\qquad\operatorname{for}\;i=1,\ldots,n.$

As mentioned before, the global approximation is typically initialized to
match the prior distribution, which boils down to setting all $k_{i}$’s and
$m_{i}$’s to zero. At each iteration, $k_{i}$ and $m_{i}$ are then updated
according to Proposition 1. The associated updates involve the potentially
large matrix inversion ${\bm{\Omega}}_{i}={\bf Q}_{-i}^{-1}$. Nonetheless,
computational efficiency can be preserved thanks to Woodbury’s identity, which
gives

${\bm{\Omega}}_{i}={\bf Q}_{-i}^{-1}=\left({\bf Q}_{{\textsc{ep}}}-k_{i}{\bf
x}_{i}{\bf
x}_{i}^{\intercal}\right)^{-1}={\bm{\Omega}}_{{\textsc{ep}}}+\dfrac{k_{i}}{1-k_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}}\left({\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}\right)\left({\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}\right)^{\intercal},$

and avoids any explicit $\mathcal{O}(p^{3})$ matrix inversion. Likewise, the
updated global covariance
${\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}=({\bf Q}_{-i}+{\bf
Q}_{i}^{{\textsc{new}}})^{-1}={\bm{\Sigma}}_{h_{i}}$ can be updated without
matrix inversion as in Lemma 1.

The same structure can be exploited in the computation of the approximate
marginal likelihood, presented in Section 2, while matching the same reduction
in computational costs. To this end, the proportionality constant $Z_{i}$ of
each site needs to be updated at each iteration together with the natural
parameters ${\bf r}_{i}$ and ${\bf Q}_{i}$ according to the rule (3). In the
case of glms, some simplifications occur also in this respect. The computation
of the normalizing constant for the hybrid distribution benefits as well of
the dimensionality reduction, since now
$Z_{h_{i}}=\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i},\,{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}d\eta_{i}$. Moreover, the
structure of ${\bf r}_{\textsc{ep}}$ and ${\bf Q}_{\textsc{ep}}$ can be
further leveraged, leading to the update reported in the proposition below.

######  Proposition 2.

For any model of the form (1) with $p(y_{i}\mid{\bm{\beta}})=\ell_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})$, the updated value of the log-proportionality
constant for the $i$-th ep site, $\log Z_{i}$, $i=1,\ldots,n$, equals

$\begin{split}\log
Z_{i}^{{\textsc{new}}}=\dfrac{1}{2}\bigg{[}&\dfrac{1}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}}\left(2m_{i}^{{\textsc{new}}}{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}+(m_{i}^{{\textsc{new}}})^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}-k_{i}^{{\textsc{new}}}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{2}\right)\\\
&-\log(1+k_{i}^{{\textsc{new}}}{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})\bigg{]}-\log Z_{h_{i}}.\end{split}$

After the algorithm has converged, one can compute the final ep approximate
marginal likelihood thanks to (4), which under a spherical Gaussian prior
$p({\bm{\beta}})=\phi_{p}({\bm{\beta}},\nu^{2}{\bf I}_{p})$ simplifies to

$\log m_{{\textsc{ep}}}({\bf y})=\dfrac{1}{2}\big{[}\,{\bf
r}_{{\textsc{ep}}}^{\intercal}{\bm{\xi}}_{{\textsc{ep}}}-\log|{\bf
Q}_{{\textsc{ep}}}|-2p\log\nu\,\big{]}-\sum_{i=1}^{n}\log Z_{i}.$

It is worth noting that the $\mathcal{O}(p^{3})$ direct computation of
$\log|{\bf Q}_{{\textsc{ep}}}|$ can be avoided. After initializing $\log|{\bf
Q}_{{\textsc{ep}}}|$ to $\log|{\bf Q}_{0}|=-2p\log\nu$, the log-determinant
can be updated coherently with the moment-matching step time of each site $i$.
Indeed, one has

$\begin{split}\log|{\bf Q}_{{\textsc{ep}}}^{{\textsc{new}}}|&=\log\big{|}{\bf
Q}_{{\textsc{ep}}}\big{[}{\bf
I}_{p}+(k_{i}^{{\textsc{new}}}-k_{i}){\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}{\bf x}_{i}^{\intercal}\big{]}\big{|}\\\ &=\log\big{|}{\bf
Q}_{{\textsc{ep}}}\big{|}+\log\big{[}1+(k_{i}^{{\textsc{new}}}-k_{i}){\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf x}_{i}\big{]}.\end{split}$
(5)

Initialization:

${\bm{\Omega}}_{{\textsc{ep}}}=\nu^{2}{\bf I}_{p};\,\ {\bf
r}_{{\textsc{ep}}}=\bm{0};\,\ \log|{\bf Q}_{{\textsc{ep}}}|=-2p\log\nu;$

$k_{i}=0$, $m_{i}=0$, and $\log Z_{i}=0$ for $i=1,\ldots,n$.

Optimization:

for  _$\,t\,$ from $\,1\,$ until convergence _ do

for  _$\,i\,$ from $\,1\,$ to $\,n\,$_ do

Cavity distribution

[1.5pt]
${\bm{\Omega}}_{i}={\bm{\Omega}}_{{\textsc{ep}}}+k_{i}/\left(1-k_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}\right)\left({\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}\right)\left({\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}\right)^{\intercal}$

[2pt] ${\bf r}_{-i}={\bf r}_{{\textsc{ep}}}-m_{i}{\bf x}_{i}$

[-4pt] Hybrid distribution

[1.5pt] $Z_{h_{i}}=\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i},\,{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}d\eta_{i}$

[2pt] $\mu_{i}={\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]$

[2pt] $\rho_{i}^{2}=\operatorname{var}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]$

[-3pt] $i$-th site approximation

[1.5pt] $k_{i}=1\,/\,\rho_{i}^{2}-1\,/\,{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}$

[2pt] $m_{i}=\mu_{i}\,/\,\rho_{i}^{2}-{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}\,/\,{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}$

[3pt] $\log Z_{i}=\dfrac{1}{2}\bigg{[}\dfrac{2m_{i}{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}+m_{i}^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}-k_{i}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{2}}{1+k_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}-\log(1+k_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})\bigg{]}$

[3pt] $-\log Z_{h_{i}}$

[-4pt] Global approximation

${\bf r}_{{\textsc{ep}}}={\bf r}_{-i}+m_{i}{\bf x}_{i}$

[2pt]
${\bm{\Omega}}_{{\textsc{ep}}}={\bm{\Omega}}_{i}+\rho_{i}^{2}({\bm{\Omega}}_{i}{\bf
x}_{i})({\bm{\Omega}}_{i}{\bf x}_{i})^{\intercal}$

[2pt] $\log|{\bf Q}_{{\textsc{ep}}}|=\log|{\bf
Q}_{{\textsc{ep}}}|+\log[1+\delta_{k_{i}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf x}_{i}]$

Final computation of quantities of interest:

${\bm{\xi}}_{{\textsc{ep}}}={\bm{\Omega}}_{{\textsc{ep}}}{\bf
r}_{{\textsc{ep}}}$

[2pt] $\log m_{{\textsc{ep}}}({\bf y})=\tfrac{1}{2}\big{(}{\bf
r}_{{\textsc{ep}}}^{\intercal}{\bm{\xi}}_{{\textsc{ep}}}-\log|{\bf
Q}_{{\textsc{ep}}}|-2p\log\nu\big{)}-\sum_{i=1}^{n}\log Z_{i}$

Output: $({\bm{\xi}}_{{\textsc{ep}}},{\bm{\Omega}}_{{\textsc{ep}}},\log
m_{{\textsc{ep}}}({\bf y}))$

Algorithm 1 ep glm \- $\mathcal{O}(p^{2}n)$ cost per iteration

Combining all the structural reformulations presented in this section one gets
Algorithm 1. Except for the computation of the approximate marginal
likelihoods, Algorithm 1 extends to the general glms framework the ep
algorithm presented in [12] and implemented in the R package EPGLM for probit
and logit regressions. However, the exploitation of simplifications and
iterative updates following from Propositions 1–2 and equation (5) leads to
substantial computational advantages (see Section 6). Note that if the
interest is only on posterior moments approximation, the $Z_{i}$ are not
needed, and neither $|{\bf Q}_{{\textsc{ep}}}|$ is, so in such a case the
algorithm can be implemented avoiding the corresponding lines. Noticeably,
Algorithm 1 has per-iteration algebraic cost of $\mathcal{O}(p^{2}n)$, since
the update of each site $i=1,\ldots,n$ requires the computation of quadratic
forms involving $p\times p$ matrices and the computation of the product
between $p\times p$ matrices and $p$-dimensional vectors. This has to be
complemented with a further $\mathcal{O}(n\,\mathcal{C}_{\textsc{mm}})$ cost
per iteration, where $\mathcal{C}_{\textsc{mm}}$ represent the cost of each
model-dependent moment matching computations. The implementation from
Algorithm 1 is efficient as long as $p$ is of the order of a few hundred, and
in particular when $p<n$. However, although explicit $p\times p$ matrix
inversions are avoided, this computational cost might be computationally
inefficient or impractical in high-dimensional settings. To overcome this
issue, in Section 3.2 we derive in full detail a general implementation of ep
for Bayesian glms (1) having per-iteration cost $\mathcal{O}(pn^{2})$. The key
to achieving such improved scaling lies in performing a simultaneous update of
sufficient functionals of the global covariance ${\bf Q}_{\textsc{ep}}$ and
all local ones ${\bf Q}_{i}$ at every site update, as first noted in [2] for a
specific class of models generalizing the probit regression.

###  3.2 Posterior moments and marginal likelihood with linear cost in $p$

A close inspection of Algorithm 1 shows that the $p\times p$ dimensional
matrices ${\bm{\Omega}}_{{\textsc{ep}}}$ and ${\bm{\Omega}}_{i}$ needed for
the optimization of the $i$-th site enter into the computations always by
multiplication with the vector ${\bf x}_{i}$. This includes both the algebraic
component of updates and the computation of the reduced hybrid moments. As a
consequence, the whole ep routine can be written directly in terms of the
$p$-dimensional vectors ${\bf w}_{i}={\bm{\Omega}}_{i}{\bf x}_{i}$ and ${\bf
v}_{i}={\bm{\Omega}}_{{\textsc{ep}}}{\bf x}_{i}$, $i=1,\ldots,n$. In turn,
such quantity can be efficiently computed and updated in a vectorized form, as
shown in Proposition 3. Notice that, at each iteration, updating the value of
$k_{i}$ to $k_{i}^{{\textsc{new}}}$ for any $i$ induces a new value of ${\bf
Q}_{i}$ and, as a consequence, of ${\bm{\Omega}}_{{\textsc{ep}}}$. This means
that each time a site $i$ is updated, one has to recompute the corresponding
new values of all the ${\bf v}_{j}$’s, $j=1,\ldots,n$, as reported in
Proposition 3, whose proof can be found in the appendix.

######  Proposition 3.

For any model of the form (1) with $p(y_{i}\mid{\bm{\beta}})=\ell_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})$, and calling ${\bf
w}_{i}={\bm{\Omega}}_{i}{\bf x}_{i}$ and ${\bf
v}_{i}={\bm{\Omega}}_{{\textsc{ep}}}{\bf x}_{i}$, $i=1,\ldots,n$, it holds

${\bf w}_{i}=d_{i}{\bf v}_{i},$

where $d_{i}=(1-k_{i}{\bf x}_{i}^{\intercal}{\bf v}_{i})^{-1}$ and $k_{i}$
refers to the value before site $i$ is updated. After site $i$ is updated, it
holds that, for $j=1,\ldots,n$,

${\bf v}_{j}^{{\textsc{new}}}={\bf v}_{j}-c_{i}({\bf x}_{i}^{\intercal}{\bf
v}_{j}){\bf v}_{i},$

where we introduced
$c_{i}=(k_{i}^{{\textsc{new}}}-k_{i})/(1+(k_{i}^{{\textsc{new}}}-k_{i}){\bf
x}_{i}^{\intercal}{\bf v}_{i})$. Defining the $p\times n$ matrix ${\bf
V}=[{\bf v}_{1},{\bf v}_{2},\dots,{\bf
v}_{n}]={\bm{\Omega}}_{{\textsc{ep}}}{\bf X}^{\intercal}$, with ${\bf X}=({\bf
x}_{1},\ldots,{\bf x}_{n})^{\intercal}$, this update can be written in matrix
form as

${\bf V}^{{\textsc{new}}}={\bf V}-c_{i}{\bf v}_{i}{\bf x}_{i}^{\intercal}{\bf
V}.$ (6)

Update (6) is the most expensive computation performed at each site update,
having cost $\mathcal{O}(pn)$. As a consequence, this alternative ep
implementation has per-iteration algebraic-component cost
$\mathcal{O}(pn^{2})$, contrary to the $\mathcal{O}(p^{2}n)$ of Algorithm 1.
This has great practical relevance in large-$p$ settings, as it leads to
massive computational gains over routinely implemented functions, see Section
6. In particular, the linear cost in $p$ makes ep feasible in settings with
$p\gg 1000$, moving forward the frontier of scenarios in which ep can be
effectively used and thus relaxing one of the critical problems mentioned in
Section 7.3 of [12].
The initialization of the newly-introduced quantities is trivial when
${\bm{\Omega}}_{{\textsc{ep}}}$ is initialized to the prior covariance matrix
${\bf Q}_{0}^{-1}$, as standard practice. Considering the spherical Gaussian
prior case, where ${\bf Q}_{0}^{-1}=\nu^{2}{\bf I}_{p}$, this translates into
simply initializing ${\bf V}$ to $\nu^{2}{\bf X}^{\intercal}$. Furthermore, it
is sufficient to specify the initialization and updates only for ${\bf V}$,
since each ${\bf w}_{i}$ is proportional to the $i$-th column ${\bf v}_{i}$ of
${\bf V}$.

Contrarily to Algorithm 1, the global covariance matrix
${\bm{\Omega}}_{{\textsc{ep}}}$ is not readily available once the procedure
has converged. Nonetheless, it can still obtained efficiently as

${\bm{\Omega}}_{{\textsc{ep}}}=\nu^{2}{\bf I}_{p}-\nu^{2}{\bf V}{\bf K}{\bf
X},$ (7)

leading to a post-processing cost of $\mathcal{O}(p^{2}n)$. Indeed, defined
${\bf K}={\text{diag}}(k_{1},\ldots,k_{n})$, one has
${\bm{\Omega}}_{{\textsc{ep}}}^{-1}={\bf Q}_{{\textsc{ep}}}={\bf
Q}_{0}+\sum_{i=1}^{n}k_{i}{\bf x}_{i}{\bf x}_{i}^{\intercal}=\nu^{-2}{\bf
I}_{p}+{\bf X}^{\intercal}{\bf K}{\bf X}$. Calling ${\bm{\Lambda}}=({\bf
I}_{n}+\nu^{2}{\bf K}{\bf X}{\bf X}^{\intercal})^{-1}$, Woodbury’s identity
applied to the former gives ${\bm{\Omega}}_{{\textsc{ep}}}={\bf
Q}_{{\textsc{ep}}}^{-1}=\nu^{2}{\bf I}_{p}-\nu^{4}{\bf
X}^{\intercal}{\bm{\Lambda}}{\bf K}{\bf X}$. This in turn allows to elicit

${\bf V}={\bm{\Omega}}_{{\textsc{ep}}}{\bf X}^{\intercal}=\nu^{2}{\bf
X}^{\intercal}\big{[}{\bf I}_{n}-\nu^{2}{\bm{\Lambda}}{\bf K}{\bf X}{\bf
X}^{\intercal}\big{]}=\nu^{2}{\bf
X}^{\intercal}{\bm{\Lambda}}\big{[}{\bm{\Lambda}}^{-1}-\nu^{2}{\bf K}{\bf
X}{\bf X}^{\intercal}\big{]}=\nu^{2}{\bf X}^{\intercal}{\bm{\Lambda}}.$

If the interest is only in marginal posterior means and variances, the post-
processing cost can be further reduced to $\mathcal{O}(pn)$ by computing only
the diagonal in (7). The whole routine is summarized in Algorithm 2, which is
then well-suited for large-$p$ settings, and in particular cases where $p>n$.

Initialization:

${\bf r}_{{\textsc{ep}}}=\bm{0};\,\ {\bf V}=\left[{\bf v}_{1},\ldots,{\bf
v}_{n}\right]=\nu^{2}{\bf X}^{\intercal};\,\ \log|{\bf
Q}_{{\textsc{ep}}}|=-2p\log\nu;$

$k_{i}=0$, $m_{i}=0$, and $\log Z_{i}=0$ for $i=1,\ldots,n$.

Optimization:

for  _$\,t\,$ from $\,1\,$ until convergence _ do

for  _$\,i\,$ from $\,1\,$ to $\,n\,$_ do

Cavity distribution

[1.5pt] ${\bf w}_{i}=(1-k_{i}{\bf x}_{i}^{\intercal}{\bf v}_{i})^{-1}{\bf
v}_{i}$

[2pt] ${\bf r}_{-i}={\bf r}_{{\textsc{ep}}}-m_{i}{\bf x}_{i}$

[-4pt] Hybrid distribution

[1.5pt] $Z_{h_{i}}=\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i},\,{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}d\eta_{i}$

[2pt] $\mu_{i}={\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]$

[2pt] $\rho_{i}^{2}=\operatorname{var}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]$

[-3pt] $i$-th site approximation

[1.5pt] $\delta_{k_{i}}=\big{(}1\,/\,\rho_{i}^{2}-1\,/\,{\bf
x}_{i}^{\intercal}{\bf w}_{i}\big{)}-k_{i}$

[2pt] $k_{i}=\delta_{k_{i}}+k_{i}$

[2pt] $m_{i}=\mu_{i}\,/\,\rho_{i}^{2}-{\bf r}_{-i}^{\intercal}{\bf
w}_{i}\,/\,{\bf x}_{i}^{\intercal}{\bf w}_{i}$

[3pt] $\log Z_{i}=\dfrac{1}{2}\bigg{[}\dfrac{2m_{i}{\bf
r}_{-i}^{\intercal}{\bf w}_{i}+m_{i}^{2}{\bf x}_{i}^{\intercal}{\bf
w}_{i}-k_{i}({\bf r}_{-i}^{\intercal}{\bf w}_{i})^{2}}{1+k_{i}{\bf
x}_{i}^{\intercal}{\bf w}_{i}}-\log(1+k_{i}{\bf x}_{i}^{\intercal}{\bf
w}_{i})\bigg{]}-\log Z_{h_{i}}$

[3pt] Global approximation

${\bf r}_{{\textsc{ep}}}={\bf r}_{-i}+m_{i}{\bf x}_{i}$

[2pt] ${\bf V}={\bf V}-{\bf
v}_{i}\left[\delta_{k_{i}}/\left(1+\delta_{k_{i}}{\bf x}_{i}^{\intercal}{\bf
v}_{i}\right)\right]{\bf x}_{i}^{\intercal}{\bf V}$

[2pt] $\log|{\bf Q}_{{\textsc{ep}}}|=\log|{\bf
Q}_{{\textsc{ep}}}|+\log[1+\delta_{k_{i}}{\bf x}_{i}^{\intercal}{\bf v}_{i}]$

Final computation of quantities of interest:

${\bm{\Omega}}_{{\textsc{ep}}}=\nu^{2}{\bf I}_{p}-\nu^{2}{\bf V}{\bf K}{\bf
X}$, with ${\bf K}={\text{diag}}(k_{1},\ldots,k_{n})$

[2pt] ${\bm{\xi}}_{{\textsc{ep}}}={\bf r}_{{\textsc{ep}}}-\nu^{2}{\bf V}({\bf
K}{\bf X}{\bf r}_{{\textsc{ep}}})$

[2pt] $\log m_{{\textsc{ep}}}({\bf y})=\tfrac{1}{2}\big{(}{\bf
r}_{{\textsc{ep}}}^{\intercal}{\bm{\xi}}_{{\textsc{ep}}}-\log|{\bf
Q}_{{\textsc{ep}}}|-2p\log\nu\big{)}-\sum_{i=1}^{n}\log Z_{i}$

Output: $({\bm{\xi}}_{{\textsc{ep}}},{\bm{\Omega}}_{{\textsc{ep}}},\log
m_{{\textsc{ep}}}({\bf y}))$

Algorithm 2 Efficient ep glm for large $p$ \- $\mathcal{O}(pn^{2})$ cost per
iteration

##  4 Hybrids moments in GLMs: some relevant cases

The efficient structure for the ep updates described above still left
unaddressed the computation of the univariate hybrid moments, associated with
the specific model of interest. In the current section, we provide the details
of such computations for several notable glms, including probit, logit,
Poisson, and gamma regression. Even for such fundamental constructions, the
required functionals are often unavailable in closed form. Probit models
constitute an eminent exception, in that the hybrid distributions belong to a
tractable family of distributions and the exact ep updates are readily
available at no additional computational cost [2, 21, 22].
For any model lacking these desirable properties, $Z_{h_{i}}$,
${\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf x}_{i}^{\intercal}{\bm{\beta}}]$ and
$\operatorname{var}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]$ could be computed via one-dimensional
numerical integration, benefiting from the dimensionality reduction induced by
the inner product with ${\bf x}_{i}$. Choosing a large number $\mathcal{T}$ of
subintervals for the domain of the desired integral typically provides an
accurate approximation of the exact value. However, this would result in a
$\mathcal{O}(\mathcal{T})$ cost $\mathcal{C}_{\textsc{mm}}$ per each moment-
matching step, besides being prone to potential numerical issues. To preserve
the efficient nature of the ep scheme, we here provide model-specific closed-
form approximations of the required moments, cutting down
$\mathcal{C}_{\textsc{mm}}$ to $\mathcal{O}(1)$.

In practice, this boils down to supplying an accurate approximation
$\widetilde{Z}_{h_{i}}$ of the normalizing constants

$Z_{h_{i}}=\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}d\eta_{i}$

and working out its derivative [12]. Here $\lambda_{i}={\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}={\bf w}_{i}^{\intercal}{\bf
r}_{-i}$ and $\rho_{i}^{2}={\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}={\bf w}_{i}^{\intercal}{\bf x}_{i}$, depending on which version of the
ep updates is being considered. Indeed, the required moments can be further
obtained from the derivatives of $\log Z_{h_{i}}$ [47] via

$\begin{split}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]&=\lambda_{i}+\rho_{i}^{2}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}\\\ \operatorname{var}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]&=\rho_{i}^{2}+\rho_{i}^{2}\left(2\frac{\partial}{\partial\rho_{i}^{2}}\log
Z_{h_{i}}-\Big{(}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}\Big{)}^{2}\right)\rho_{i}^{2}\;,\end{split}$ (8)

as proven in Appendix S.2. In turn, the above expressions can be directly
combined with the updates from Propostion 1. Defining
$\Theta_{i}=\frac{\partial}{\partial\lambda_{i}}\log Z_{h_{i}}$ and
$\Delta_{i}=2\frac{\partial}{\partial\rho_{i}^{2}}\log
Z_{h_{i}}-\Big{(}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}\Big{)}^{2}$, simply algebraic calculations lead to

$k_{i}^{{\textsc{new}}}=-\Delta_{i}\big{/}\big{(}1+\rho_{i}^{2}\Delta_{i}\big{)},\qquad\qquad
m_{i}^{{\textsc{new}}}=k_{i}^{{\textsc{new}}}\big{(}\lambda_{i}-\Delta_{i}^{-1}\Theta_{i}\big{)}\;.$

Accordingly, any differentiable approximation $\widetilde{Z}_{h_{i}}$ for
$Z_{h_{i}}$ naturally induces approximate updates for the natural parameters
of the Gaussian ep sites. The quality of the resulting approximations will be
highly specific to the model considered, the nature of $\widetilde{Z}_{h_{i}}$
and the order of magnitude of the parameters $\lambda_{i}$ and $\rho_{i}^{2}$.

###  4.1 Probit regression

Probit regression represents a prominent example of models for which the
hybrids belong to a known family of distributions, further allowing for
closed-form ep updates. This can be seen as a relevant special case of a
broader class of models considered in [2], where the ep algorithm further
simplifies, avoiding the need to work with the possibly computationally
challenging sun hybrid distributions, as shown below. In particular, the form
of the ep updates can be obtained by leveraging results on the efficiently
tractable multivariate extended skew-normal (sn) random variables (see [57]).
In probit models, the binary responses $y_{i}\in\\{0,1\\}$ are modeled as
Bernoulli random variables with success probability $\Phi({\bf
x}_{i}^{\intercal}{\bm{\beta}})$, for each $i=1,\ldots,n$. Accordingly, each
likelihood contribution takes the form

$p(y_{i}\mid{\bm{\beta}})=\big{(}\Phi({\bf
x}_{i}^{\intercal}{\bm{\beta}})\big{)}^{y_{i}}\big{(}1-\Phi({\bf
x}_{i}^{\intercal}{\bm{\beta}})\big{)}^{1-y_{i}}=\Phi\big{(}(2y_{i}-1){\bf
x}_{i}^{\intercal}{\bm{\beta}}\big{)}\;.$ (9)

This allowed [18] to first show that the posterior distribution under Gaussian
priors belongs to the sun family. The same rationale can be extended to the
analysis of the ep hybrid distributions $h_{i}({\bm{\beta}})$, resulting in
the following proposition.

######  Proposition 4.

For the probit likelihood (9), the hybrid density $h_{i}({\bm{\beta}})$ is the
density of a multivariate extended skew-normal distribution
$\textsc{sn}_{p}({\bm{\xi}}_{i},{\bm{\Omega}}_{i},{\bm{\alpha}}_{i},\tau_{i})$
(see [57]), with

$\begin{split}{\bm{\omega}}_{i}&=\left[\text{diag}\left({\bm{\Omega}}_{i}\right)\right]^{1/2}\hskip
50.0pt\bar{{\bm{\Omega}}}_{i}={\bm{\omega}}_{i}^{-1}{\bm{\Omega}}_{i}{\bm{\omega}}_{i}^{-1}\\\
{\bm{\alpha}}_{i}&=(2y_{i}-1){\bm{\omega}}_{i}{\bf x}_{i}\hskip
55.0pt\tau_{i}=(2y_{i}-1)(1+{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1/2}{\bf x}_{i}^{\intercal}{\bm{\xi}}_{i}\;.\end{split}$

and ${\bm{\xi}}_{i}={\bf Q}_{-i}^{-1}{\bf r}_{-i}$ and ${\bm{\Omega}}_{i}={\bf
Q}_{-i}^{-1}$, as before. Its normalizing constant is
$Z_{h_{i}}=\Phi(\tau_{i})$, while the updates for one dimensional parameters
in the $i$-th approximate site take the form

$\begin{split}k_{i}^{{\textsc{new}}}&=-\zeta_{2}(\tau_{i})/\big{(}1+{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}+\zeta_{2}(\tau_{i}){\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}\\\\[5.0pt]
m_{i}^{{\textsc{new}}}&=\zeta_{1}(\tau_{i})s_{i}+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
r}_{-i}+k_{i}^{{\textsc{new}}}\zeta_{1}(\tau_{i})s_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\;,\end{split}$

where $s_{i}=(2y_{i}-1)(1+{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1/2}$, $\zeta_{1}(x)=\phi(x)/\Phi(x)$ and
$\zeta_{2}(x)=-\zeta_{1}(x)^{2}-x\zeta_{1}(x)$.

All the above holds even under the reformulation from Section 3.2. This
amounts simply to replacing ${\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}$ and ${\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}$ with ${\bf
w}_{i}^{\intercal}{\bf x}_{i}$ and ${\bf w}_{i}^{\intercal}{\bf r}_{-i}$,
respectively.

####  4.1.1 Closed-form predictive probabilities

The analytic tractability of the probit carries over to predictions for new
observations. Indeed, the Gaussian form of the global ep approximation also
has a very convenient side effect that the resulting approximate predictive
probabilities for new observations admit a simple closed form. Considering a
new observation $y_{{\textsc{new}}}$ having covariate vector ${\bf
x}_{{\textsc{new}}}$, by analogy with the exact predictive probability
$\Pr{}[y_{{\textsc{new}}}=1\mid{\bf
y}]={\mathbb{E}}[y_{{\textsc{new}}}\mid{\bf
y}]={\mathbb{E}}_{p({\bm{\beta}}\mid{\bf y})}\big{[}\Phi\big{(}{\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\beta}}\big{)}\big{]}$, one gets the ep
approximation counterpart

$\Pr{}_{\textsc{ep}}[y_{{\textsc{new}}}=1\mid{\bf
y}]={\mathbb{E}}_{q_{{\textsc{ep}}}({\bm{\beta}})}\big{[}\Phi\big{(}{\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\beta}}\big{)}\big{]}=\Phi\big{(}\big{(}1+u\big{)}^{-1/2}{\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\xi}}_{\textsc{ep}}\big{)},\vspace*{-0.2cm}$
(10)

where $u={\bf x}_{{\textsc{new}}}^{\intercal}{\bm{\Omega}}_{\textsc{ep}}{\bf
x}_{{\textsc{new}}}$ and the last equality in (10) follows by Lemma 7.1 in
[57]. The evaluation of (10) does not pose any computational difficulties and
has a cost that is always dominated by the ep algorithm, as the only
computationally relevant part is the computation of the quadratic form $u$.
This can be computed directly when using Algorithm 1 (usually, when $p<n$), as
${\bm{\Omega}}_{\textsc{ep}}$ is directly returned by the algorithm, and $u$
can be computed at cost $\mathcal{O}(p^{2})$. On the other hand, when using
Algorithm 2 (usually, when $p>n$), this direct computation can be avoided. For
simplicity, assume once more a spherical Gaussian prior
$p({\bm{\beta}})=\phi_{p}({\bm{\beta}},\nu^{2}{\bf I}_{p})$, so that
$q_{0}({\bm{\beta}})$ is fixed and with natural parameter ${\bf r}_{0}=\bm{0}$
and ${\bf Q}_{0}=\nu^{-2}{\bf I}_{p}$. Then, exploiting (7), it can be seen
that $u=\nu^{2}\left[{\bf x}_{{\textsc{new}}}^{\intercal}{\bf
x}_{{\textsc{new}}}-\big{(}{\bf V}^{\intercal}{\bf
x}_{{\textsc{new}}}\right)^{\intercal}{\bf K}\left({\bf X}{\bf
x}_{{\textsc{new}}}\right)\big{]}$, which can be computed at cost
$\mathcal{O}(pn)$. As a consequence of this argument, the closed-form
approximation of the exact predictive probability (10) can be efficiently
computed at cost $\mathcal{O}(p\cdot\min\\{p,n\\})$ from the ep parameters.

###  4.2 Logistic regression

Similarly to the probit case, logistic regression describes binary responses
$y_{i}\in\\{0,1\\}$ as Bernoulli random variables, for each $i=1,\ldots,n$,
although now the success probabilities $(1+e^{-{\bf
x}_{i}^{\intercal}{\bm{\beta}}})^{-1}$ are mediated by the logit link. In such
a case, the likelihood terms become

$p(y_{i}\mid{\bm{\beta}})=\Big{(}\frac{1}{1+e^{-{\bf
x}_{i}^{\intercal}{\bm{\beta}}}}\Big{)}^{y_{i}}\Big{(}\frac{1}{1+e^{{\bf
x}_{i}^{\intercal}{\bm{\beta}}}}\Big{)}^{1-y_{i}}=\frac{1}{1+e^{-(2y_{i}-1){\bf
x}_{i}^{\intercal}{\bm{\beta}}}}\;.$

Contrary to the probit case, the hybrid moments for the logit exact sites do
not admit closed-form expressions. Following [12], we propose a two-fold
application of the approximation introduced by [36]

$\begin{split}\int\frac{1}{1+e^{-\eta}}\phi\big{(}\eta-\lambda,\,\rho^{2}\big{)}d\eta\approx\int\Phi\left(\eta\sqrt{\pi/8}\right)\phi\big{(}\eta-\lambda,\,\rho^{2}\big{)}d\eta&=\Phi\left(\frac{\lambda}{\sqrt{8\pi+\rho^{2}}}\right)\\\
&\approx\frac{1}{1+e^{-\lambda/\sqrt{1+\rho^{2}\pi/8}}}\end{split}$

Accordingly, the hybrid normalizing constants can be approximated as

$Z_{h_{i}}\approx\widetilde{Z}_{h_{i}}\coloneq\frac{1}{1+e^{-(2y_{i}-1)\lambda_{i}/\sqrt{1+\rho_{i}^{2}\pi/8}}}\;.$

Exploiting the fact that $\frac{\partial}{\partial
s}\big{(}(1+e^{-s})^{-1}\big{)}=(1+e^{-s})^{-1}\big{(}1-(1+e^{-s})^{-1}\big{)}$,
one gets

$\begin{split}\frac{\partial}{\partial\lambda_{i}}\log\widetilde{Z}_{h_{i}}&=\frac{(2y_{i}-1)}{\sqrt{1+\rho_{i}^{2}\pi/8}}\,(1-\widetilde{Z}_{h_{i}})\\\
\frac{\partial}{\partial\rho_{i}^{2}}\log\widetilde{Z}_{h_{i}}&=-\frac{1}{2}\frac{\pi}{8}\frac{\lambda_{i}}{\sqrt{1+\rho_{i}^{2}\pi/8}}\,\frac{\partial}{\partial\lambda_{i}}\log\widetilde{Z}_{h_{i}}\;.\end{split}$

An approximate implementation of the ep moment matching is thus obtained
replacing $Z_{h_{i}}$, $\frac{\partial}{\partial\lambda_{i}}\log Z_{h_{i}}$
and $\frac{\partial}{\partial\rho_{i}^{2}}\log Z_{h_{i}}$ respectively with
$\widetilde{Z}_{h_{i}}$,
$\frac{\partial}{\partial\lambda_{i}}\log\widetilde{Z}_{h_{i}}$ and
$\frac{\partial}{\partial\rho_{i}^{2}}\log\widetilde{Z}_{h_{i}}$.

###  4.3 Poisson regression

We now move beyond the modeling of binary responses, focusing instead on count
data $y_{i}\in\mathbb{N}$. In particular, we focus on Poisson regression under
the logarithmic link. This corresponds to a likelihood of the form

$p(y_{i}\mid{\bm{\beta}})=\frac{\big{(}e^{{\bf
x}_{i}^{\intercal}{\bm{\beta}}}\big{)}^{y_{i}}e^{-(e^{{\bf
x}_{i}^{\intercal}{\bm{\beta}}})}}{y_{i}!}\;,$

in turn, associated with hybrid normalizing constants

$Z_{h_{i}}=\frac{1}{y_{i}!}\,e^{y_{i}\lambda_{i}+\frac{1}{2}\rho_{i}^{2}y_{i}^{2}}\int
e^{-e^{\eta_{i}}}\phi\big{(}\eta_{i}-(\lambda_{i}+\rho_{i}^{2}y_{i}),\,\rho_{i}^{2}\big{)}d\eta_{i}\;.$
(11)

The integral within the expression above does not admit an exact solution.
Indeed, it is easy to see that it coincides with the Laplace transform
$\mathcal{L}(s)$ of a log-normal random variable, evaluated at a specific
value of its argument $s$. More details on the connection with the log-normal
Laplace transform are provided in Appendix S.2. Several dedicated
approximations have been developed over the years. As a default for the ep
updates, we propose exploiting the contribution by [56], which gives

$\widetilde{Z}_{h_{i}}=\frac{1}{y_{i}!}\,e^{y_{i}\lambda_{i}+\frac{1}{2}\rho_{i}^{2}y_{i}^{2}}\,\frac{e^{-\frac{1}{\rho_{i}^{2}}\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}\left(1+\frac{1}{2}\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}\right)}}{\sqrt{1+\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}}}\;.$

Here $\mathcal{W}(x)$ the so-called Lambert W function, which is defined only
implicitly as the solution of $\mathcal{W}(x)e^{\mathcal{W}(x)}=x$, but that
can be evaluated efficiently in most programming languages. As proved in
Appendix S.2, the derivatives of the logarithm of $\widetilde{Z}_{h_{i}}$ read

$\begin{split}\frac{\partial}{\partial\lambda_{i}}\log\widetilde{Z}_{h_{i}}&=y_{i}-\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}\left(\frac{1}{\rho_{i}^{2}}+\frac{1}{2}\Big{(}1+\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}\Big{)}^{-2}\right)\\\
\frac{\partial}{\partial\rho_{i}^{2}}\log\widetilde{Z}_{h_{i}}&=\frac{y_{i}^{2}}{2}+\frac{1}{\rho_{i}^{2}}\frac{\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}}{\rho_{i}^{2}}\Big{(}1+\frac{1}{2}\mathcal{W}\big{(}\rho_{i}^{2}e^{\lambda_{i}+\rho_{i}^{2}y_{i}}\big{)}\Big{)}\\\
&\quad+\frac{1+\rho_{i}^{2}y_{i}}{\rho_{i}^{2}}\Big{(}\frac{\partial}{\partial\lambda_{i}}\log\widetilde{Z}_{h_{i}}-y_{i}\Big{)}\;.\end{split}$

Our empirical insights suggest that the above provides an accurate
approximation as long as $y_{i}>0$ or $\lambda_{i}$ is large enough. Indeed,
[56] essentially apply Laplace’s method to the Laplace transform
$\mathcal{L}(s_{i})$ of the log-normal distribution. Such an approximation is
exact only in the limit $s_{i}\rightarrow\infty$. In the context of ep Poisson
regression, $\mathcal{L}(s_{i})$ needs to be evaluated at
$s_{i}=\exp{(\lambda_{i}+\rho_{i}^{2}y_{i})}$ for the $i^{th}$ site.
Observations with $y_{i}=0$ and small $\lambda_{i}$ clearly fall outside the
region of large $s_{i}$. To preserve accuracy in all settings, we propose
employing the alternative approximation scheme by [59, 43] for sites with
$y_{i}=0$ and $\lambda_{i}$ smaller than some heuristic threshold. The details
on the resulting approximate moment matching are provided in Appendix S.2.

####  4.3.1 Closed-form predictive means

The combination of the logarithmic link with the Gaussian ep posterior
approximation also brings a very convenient closed-form expression for the
predictive mean of a new observation $y_{{\textsc{new}}}$ having covariate
vector ${\bf x}_{{\textsc{new}}}$, that is
${\mathbb{E}}[y_{{\textsc{new}}}\mid{\bf y}]$. Note that this provides the
analog of the closed-form predictive probabilities derived for the probit
model in Section 4.1.1 as, in that case, $\Pr{}[y_{{\textsc{new}}}=1\mid{\bf
y}]={\mathbb{E}}[y_{{\textsc{new}}}\mid{\bf y}]$. When a logarithmic link is
employed, as in the default choice for the Poisson model (although the same
argument holds true for any glms with a logarithmic link), one has

${\mathbb{E}}[y_{{\textsc{new}}}\mid{\bf
y}]={\mathbb{E}}_{p({\bm{\beta}}\mid{\bf
y})}\left[{\mathbb{E}}_{p(y_{{\textsc{new}}}\mid{\bm{\beta}})}[y_{{\textsc{new}}}]\right]={\mathbb{E}}_{p({\bm{\beta}}\mid{\bf
y})}\left[e^{{\bf x}_{{\textsc{new}}}^{\intercal}{\bm{\beta}}}\right].$ (12)

Thus, by replacing the true posterior distribution with its ep approximation
in (12), the ep approximate predictive mean for observation
$y_{{\textsc{new}}}$ is
${\mathbb{E}}_{{\textsc{ep}}}[y_{{\textsc{new}}}\mid{\bf
y}]={\mathbb{E}}_{q_{{\textsc{ep}}}({\bm{\beta}})}\left[e^{{\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\beta}}}\right].$ Now, calling $\eta={\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\beta}}$, since
$q_{{\textsc{ep}}}({\bm{\beta}})=\phi_{p}({\bm{\beta}}-\xi_{{\textsc{ep}}},{\bm{\Omega}}_{{\textsc{ep}}})$,
one has that, according to the ep approximation of the posterior distribution
of ${\bm{\beta}}$, $\eta\sim N(m_{{\textsc{new}}},s_{{\textsc{new}}}^{2})$,
with $m_{{\textsc{new}}}={\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\xi}}_{{\textsc{ep}}}$ and
$s_{{\textsc{new}}}^{2}={\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{{\textsc{new}}}$. Consequently, $e^{\eta}$ has log-normal distribution
with parameters $m_{{\textsc{new}}}$ and $s_{{\textsc{new}}}^{2}$ and we
immediately obtain

${\mathbb{E}}_{{\textsc{ep}}}[y_{{\textsc{new}}}\mid{\bf
y}]={\mathbb{E}}_{\eta\sim
N(m_{{\textsc{new}}},s_{{\textsc{new}}}^{2})}\left[e^{\eta}\right]=\exp\\{m_{{\textsc{new}}}+s_{{\textsc{new}}}^{2}/2\\},$

where the quadratic form ${\bf
x}_{{\textsc{new}}}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{{\textsc{new}}}$ can be efficiently computed at cost
$\mathcal{O}(p\cdot\min{p,n})$, as explained in Section 4.1.1, without
affecting the overall cost of the ep procedure.

###  4.4 Gamma regression

Gamma regression is widely used glm to model continuous, non-negative, and
positively skewed data. As the name points out, it postulates gamma
distributed responses
$y_{i}\mid{\bm{\beta}}\overset{ind}{\sim}\mathcal{G}a(\upsilon_{i},\varrho({\bf
x}_{i}^{\intercal}{\bm{\beta}})/\upsilon_{i})$. Here we focus on the
logarithmic link $\varrho({\bf x}_{i}^{\intercal}{\bm{\beta}})=e^{{\bf
x}_{i}^{\intercal}{\bm{\beta}}}$ since it is the standard choice in the gamma
likelihood glm, while we will consider the additional shape parameter
$\upsilon_{i}>0$ to be fixed. This leads to a likelihood of the form

$p(y_{i}\mid{\bm{\beta}})=\frac{1}{\Gamma(\upsilon_{i})}\Big{(}\frac{\upsilon_{i}}{e^{{\bf
x}_{i}^{\intercal}{\bm{\beta}}}}\Big{)}^{\upsilon_{i}}y_{i}^{\upsilon_{i}-1}e^{-\upsilon_{i}y_{i}e^{-{\bf
x}_{i}^{\intercal}{\bm{\beta}}}}\;,$

reflected in hybrid normalizing constants

$Z_{h_{i}}=\frac{\upsilon_{i}^{\upsilon_{i}}}{\Gamma(\upsilon_{i})}\,y_{i}^{\upsilon_{i}-1}\,e^{-\upsilon_{i}\lambda_{i}+\frac{1}{2}\rho_{i}^{2}\upsilon_{i}^{2}}\int
e^{-\upsilon_{i}y_{i}e^{\eta_{i}}}\phi\big{(}\eta_{i}-(\rho_{i}^{2}\upsilon_{i}-\lambda_{i}),\rho_{i}^{2}\big{)}d\eta_{i}\;.$

Similarly to Poisson regression, the required integral takes the form of the
log-normal Laplace transform. Accordingly, the same heuristics exploited
before can be adapted to approximate moment matching in gamma regression as
well.

##  5 Simulation Studies

###  5.1 Probit regression

In this section, we demonstrate the advantages of our proposed efficient ep
through a simulation study. We apply the probit regression model (9), with
$\nu^{2}=25$, to multiple simulated datasets, with $n=500$ and
$p=125,\,250,\,500,\,1000$ and $2000$. The $p-1$ covariates, excluding the
intercept, are generated independently from a standard normal distribution and
are then standardized to have zero mean and $0.5$ standard deviation,
following [26] as in [12]. Our investigation focuses on the performance of ep
when utilizing the efficient implementations outlined in Algorithm 1 and
Algorithm 2, choosing automatically the optimal one based on the
dimensionality of the input data. These implementations are used to compute
posterior moments and predictive probabilities (for $\tilde{n}=100$ test
units) when $p<n$ and $p\geq n$, respectively. We refer to this implementation
as ep-eff in the following discussion.

We compare the ep-eff approach with the pfm-vb approximation introduced in
[20], which has higher approximation accuracy than the mean-field variational
approach presented in [13]. This comparison is made in terms of both running
time and the quality of the approximation. We measure the latter in terms of
the goodness of the approximation of the $p$ posterior means and standard
deviations and the $\tilde{n}$ predictive probabilities. For each of these
three quantities, we compute the median absolute difference (together with the
first and third quartile) with respect to those computed via $10000$ mcmc
samples obtained using the rstan [50] implementation of the No-U-Turn sampler
(nuts), which performs Hamiltonian Monte Carlo (hmc) sampling. These samples
were generated by $5$ chains of $2000$ samples each, after a burn-in of
$1000$. For more information on the hmc sampler in probit regression, please
refer to [12].

As presented in Table 1, the nuts method is considerably slower, even though
running five shorter chains in parallel. This is why we considered $p$ only up
to $2000$ in order to be able to get posterior samples from the nuts method,
which is used as a benchmark. To illustrate the computational advantages in
comparison to standard ep implementations, we also compare the running time
needed to obtain the ep approximations using the R function EPprobit from the
EPGLM package, which implements the ep derivations reported in [12]. As shown
in Table 1, ep-eff dramatically reduces the computational effort to
approximate the posterior moments and predictive probabilities, compared to
the standard EPprobit in high dimensions. This results in a decrease in
running time by more than three orders of magnitude when $p=2000$, with
computational gains increasing as expected with higher values of $p$. The ep-
eff running times, although generally much lower than the ones of EPprobit,
are still higher than the ones of pfm-vb for $p>n$. For $p<n$ instead, the
lower number of iterations needed to reach convergence by ep make ep-eff
faster than pfm-vb.

Table 1: Running time, in seconds, to compute posterior moments and predictive probabilities for $\tilde{n}=100$ test observations with the ep approximation as in Algorithms 1 and 2 (ep-eff), with the ep approximation computed via the R function EPprobit (EPprobit), with the pfm-vb approximation (pfm-vb) and with the rstan implementation of the No-U-Turn hmc sampler (nuts) for probit regression with $n=500$ and $\nu^{2}=25$.  | p
---|---
| Method | 125 | 250 | 500 | 1000 | 2000
Running time | ep-eff |  0.55 |  1.13 |  2.06 |  3.88 |  8.15
(seconds) | EPprobit |  6.06 |  38.36 |  232.92 |  1867.35 |  21206.49
| pfm-vb |  2.88 |  1.90 |  3.25 |  1.21 |  1.81
| nuts |  152.13 |  491.08 |  1971.50 |  7745.48 |  20415.65

Nevertheless, when examining the quality of the approximation of the two
posterior moments and predictive probabilities (Figure 1), ep-eff consistently
provides accurate approximations across different dimensions of $p$, while
pfm-vb gets similar accuracy only for $p\gtrsim 2n$. Combined with the
improved accuracy of ep in all scenarios, this highlights the importance of
the proposed efficient implementations for ep, making it computationally
feasible in almost every practical scenario including the challenging high-
dimensional and/or large sample-size settings where previously available
implementations are either impractical or inaccurate.

Finally, evaluation of the ep approximate predictive probabilities is a direct
by-product of the computation of optimal posterior approximation, thanks to
the closed-form expression derived in Section 4.1.1. Accordingly, the
associate additional computation time is minimal in the case of ep. In
contrast, the equivalent step in pfm-vb requires a (slightly) higher cost, due
to the lack of a closed-form expression of the predictive probabilities, which
are computed from Equation (11) in [20] sampling from $n$ independent
univariate truncated normals. Although this can be done efficiently and
requires no more than one second in the considered scenarios, the availability
of a closed form for ep makes this step computationally faster. For the sake
of completeness, we present the running times for the various methods in
calculating only the posterior moments in Table S.1 of the Appendix. The
comparison mirrors that of the running times for predictive probabilities in
Table 1. All computations are carried out in R, using plain R code, on an
eight-core i9 a Lenovo ThinkStation P330 with 16Gb RAM (Windows 11 Pro).

Figure 1: For varying $p$, median absolute difference between the $p$
posterior means, the standard deviations and the $\tilde{n}=100$ predictive
probabilities resulting from $10000$ hmc (nuts) samples and the ones arising
from ep-eff and pfm-vb for probit regression with $n=500$ and $\nu^{2}=25$.
Shaded areas denote the first and third quartiles.

###  5.2 Poisson regression

We now show the performance of the proposed ep efficient algorithm for
Bayesian Poisson regression under Gaussian prior. This constitutes a
widespread popular Bayesian glm, for which, however, many open challenges
remain open due to the lack of a tractable posterior distribution and possible
mixing issues in the available sampling methods (see also [15]). We generated
synthetic datasets with the same dimensions used in Section 5.1 so that
sampling methods could also be run and used as a benchmark. Specifically, we
generated five datasets fixing the number of observations $n=500$ and held-out
units $\tilde{n}=100$, while we varied $p$ over the values
$125,250,500,1000,2000$. The generation process for such datasets is as
follows. The $p-1$ covariate values for each observation were preliminarily
generated independently from a standard Gaussian. Then, for each variable, we
centered the generated covariate values by subtracting the variable empirical
mean and rescaled them by dividing by $2\sqrt{p-1}\widehat{sd}_{j}$, where
$\widehat{sd}_{j}$ is the estimated standard deviation for variable $j$. This
is done to have a constant variance of the linear predictor generating the
data across the multiple scenarios, since the parameters $\beta_{j}$,
$j=2,\ldots,p$ are uniformly sampled from the interval $[-5,5]$, while
$\beta_{1}$ was fixed to $5$. The model was completed by taking a dimension-
adjusted prior on the parameters $\nu^{2}=250/p$.

In line with Section 5.2, we used $10000$ hmc samples obtained from nuts,
implementing the Bayesian Poisson model in rstan. Also in this case, the
samples were obtained from $5$ parallel chains of $2000$ samples each, after a
burn-in of $1000$. We thus compared the accuracy of the resulting parameters’
posterior means and standard deviations and of the predictive means for the
$\tilde{n}=100$ held out units, i.e. ${\mathbb{E}}[{\bf
y}_{{\textsc{new}}}\mid{\bf y}]$, resulting from the efficient ep
implementation reported in Algorithms 1 and 2 (ep-eff) and from the
Metropolis-Hastings (mh) scheme developed by [15], implemented in the R
package bpr. This was recently introduced to get posterior samples for Poisson
regression models with possibly different choices for the prior distribution,
including the Gaussian one. bpr is particularly efficient for smaller $p$
dimensions, which is indeed the focus of [15]. Conversely, it may become more
computationally demanding in higher dimensions due to the intrinsic
limitations of mh schemes in high dimensions, as we can also see in Table 2
and Figure 2.

Table 2: Running time, in seconds, to compute posterior moments and predictive means of the linear predictor for $\tilde{n}=100$ test observations with the ep approximation as in Algorithms 1 and 2 (ep-eff), with the Metropolis-Hastings sampler from the R package bpr (bpr) and with the rstan implementation of the No-U-Turn hmc sampler (nuts) for Poisson regression with $n=500$ and $\nu^{2}=250/p$.  | p
---|---
| Method | 125 | 250 | 500 | 1000 | 2000
Running time | ep-eff |  0.94 |  2.48 |  3.14 |  5.56 |  12.09
(seconds) | bpr |  164.62 |  727.21 |  4021.31 |  25370.04 |  184153.96
| nuts |  62.08 |  138.64 |  580.82 |  1082.71 |  1672.56
Figure 2: For varying $p$, median absolute difference between the $p$
posterior means, the standard deviations and the $\tilde{n}=100$ predictive
probabilities resulting from $10000$ hmc (nuts) samples and the ones arising
from ep-eff and from $10000$ samples from bpr [15] for Poisson regression with
$n=500$ and $\nu^{2}=250/p$. Shaded areas denote the first and third
quartiles.

As shown in Table 2, the nuts method is considerably slower, despite the
parallel execution of five chains. As expected, ep-eff significantly minimizes
the computational effort required to approximate the posterior moments and
predictive means, outperforming the standard sampling methods available. The
computational advantages increase as the problem’s dimensionality grows (i.e.,
higher values of $p$), resulting in a decrease in running time by more than
two orders of magnitude when $p=2000$. However, when assessing the quality of
the approximation of the two posterior moments and predictive means (Figure
2), ep-eff consistently provides accurate approximations across different
dimensions of $p$, proving to be practically indistinguishable from nuts. In
conjunction with the improved accuracy of ep in all scenarios, this
underscores the importance of the proposed efficient implementations for ep.
It renders inference and prediction via glms computationally feasible in
almost every practical scenario, including the challenging high-dimensional
and large sample-size settings where previously available implementations are
impractical. Finally, note that for ep-eff we rely on a basic and easily
reproducible R implementation. Optimized computational routines relying on C++
coding can be further considered as it is done in nuts and bpr.

##  6 Illustrations

###  6.1 Probit regression

As discussed in [12], state-of-the-art ep implementations for Bayesian binary
regression provide accurate results but are feasible in small-to-moderate $p$
settings, while they tend to become rapidly impractical in large $p$ contexts,
such as $p>1000$. As a first example to illustrate the improvements provided
by the ep algorithm for high-dimensional glms (Algorithm 2), we consider a
large $p$ study to quantify the accuracy and computational costs, expressed as
running times, in a real-world application with binary responses. More
precisely, we demonstrate the performance of the ep approximation, derived in
Section 4.1, in a medical application that was previously used in [20] for
comparisons among some state-of-the-art algorithms for approximating posterior
and predictive distributions under probit regression.

Following [20], we model the presence ($y_{i}=1$) or absence ($y_{i}=0$) of
Alzheimer’s disease in its early stages as a function of demographic data,
genotype, and assay results via an interpretable probit regression. We further
increase flexibility by considering pairwise interactions in the design
matrix, ${\bf X}$, thus obtaining $p=9036$ predictors, including the
intercept, for $333$ subjects. The original dataset is available in the R
library AppliedPredictiveModeling [14]. The predictors have been standardized
to have mean $0$ and standard deviation $0.5$, before being included in the
probit models, following a common recommendation in probit regression [26, 12,
20]. As for the simulation study, we perform Bayesian inference by employing
independent Gaussian priors with mean $0$ and variance $\nu^{2}=25$ for each
$\beta_{j}$ $(j=1,\ldots,9036)$ [26]. These priors are updated with the probit
likelihood of $n=300$ units, after excluding $33$ subjects to study the
behavior of the predictive probabilities in such large $p$ settings, along
with the performance of the overall approximation of the posterior.

Figure 3: Results for the Alzheimer data under probit regression with
$\nu^{2}=25$. Boxplot of the differences of the
${\mathbb{E}}[\beta_{1}\mid{\bf y}],\ldots,{\mathbb{E}}[\beta_{p}\mid{\bf
y}]$, $\log\big{(}\sqrt{\text{var}[\beta_{1}\mid{\bf
y}]}\,\big{)},\ldots,\log\big{(}\sqrt{\text{var}[\beta_{p}\mid{\bf
y}]}\,\big{)}$, and $\Pr[y_{n+1}=1\mid{\bf y}],\ldots,\Pr[y_{n+100}=1\mid{\bf
y}]$ computed with ep and the pfm-vb solutions, using the inferences obtained
via $10000$ i.i.d. samples from the exact sun as benchmark.

We compare the accuracy of ep, implemented according to Algorithm 2, and the
pfm-vb approximation [20] considered in Section 5.1 for recovering functionals
of interest of the posterior distribution as well as the predictive
probabilities for the $33$ held-out units. The dimensions of the problem —
i.e., large $p$ and moderate $n$ — make the i.i.d. sampler developed by [18]
feasible, and thus a preferable solution over mcmc schemes. Indeed, since the
nuts hmc sampler discussed in Section 5.1 takes more than six hours to run in
such an example, comparisons are conducted using as a benchmark the results
obtained from $10^{4}$ i.i.d. samples from [18]. As illustrated in Section 5.1
and by theoretical results in [20], we are focusing on the best scenario for
the pfm-vb, i.e., $p\gg n$ since we have already shown in simulation study
that ep outperform pfm-vb in term of accuracy in $n\geq p$ settings. In such a
high-dimensional setting, with $p\gg 1000$, however, routinely-implemented
functions, such as EPprobit in the EPglm R package, are not practically
usable, as they require more than six hours to run (see also [20]). This also
highlights the improvements of the proposed efficient ep implementation
derived in Algorithm 2, which instead requires only a few seconds, see Table
3.

The associate results are illustrated in Figure 3, where we show the boxplots
of the differences between the posterior means ${\mathbb{E}}[\beta_{1}\mid{\bf
y}],\ldots,{\mathbb{E}}[\beta_{p}\mid{\bf y}]$, log standard deviations
$\log\big{(}\sqrt{\text{var}[\beta_{1}\mid{\bf
y}]}\,\big{)},\ldots,\allowbreak\log\big{(}\sqrt{\text{var}[\beta_{p}\mid{\bf
y}]}\,\big{)}$, and predictive probabilities $\Pr[y_{n+1}=1\mid{\bf
y}],\ldots,\Pr[y_{n+100}=1\mid{\bf y}]$ obtained via the i.i.d. sampler and
the ones resulting from the two approximate methods. As expected, both
approximate solutions are extremely accurate for all the quantities of
interest, which are almost indistinguishable from the ones obtained with the
i.i.d. sampler. ep shows even better accuracy for the posterior means and
predictive probabilities.

Finally, Table 3 provides insights into the computational time of the
different algorithms. The approximated methods enable the computation of the
desired quantities in a fraction of the time required by the i.i.d. sampler,
which takes approximately $25$ minutes. Although slightly slower than pfm-vb,
which takes only approximately $3$ seconds to produce the results, the
proposed implementation of ep developed in Algorithm 2 takes about $12$
seconds to run, paving the way to the use of ep in high-dimensional Bayesian
glms and showing massive improvements in running times over the available ep
routine for the probit model (EPprobit), which could not run within six hours.

Table 3: Running time, in seconds, to compute posterior moments and predictive probabilities for $\tilde{n}=33$ test observations with the Monte Carlo sampler (mc), with the rstan implementation of the No-U-Turn hmc sampler (nuts), with the ep approximation as in Algorithms 1 and 2 (ep-eff), with the ep approximation computed via the R function EPprobit (EPprobit), and with the pfm-vb approximation (pfm-vb) for probit regression on the Alzheimer dataset ($n=300$) with $\nu^{2}=25$.  | Method |  |
---|---|---|---
| nuts | EPprobit | mc | ep-eff | pfm-vb
Running time (in seconds) | > 21600 | > 21600 | 1541.76 | 12.31 | 3.17

###  6.2 Poisson regression

To further illustrate the improvements provided by the ep algorithm for high-
dimensional glms (Algorithm 2), we undertake a $p\gg n$ study to quantify the
accuracy and computational costs, measured as running times, in a real-world
application with integer responses. Specifically, we showcase the performance
of the ep approximation for glms, derived in Algorithm 2 (with the efficient
hybrid moment approximation computed in Section 4), in a sports analytics
application that was previously used in [15] for comparisons among state-of-
the-art sampling algorithms under Poisson log-linear regression. We analyze
data of match scores from the Italian Serie A 2020–2021 season, which are
publicly available at http://www.football-data.co.uk. The number of goals is
considered as the variable of interest ($y_{i}\in\mathbb{N}$), while the team,
several betting odds (for different betting types and bookmakers), and an
indicator of whether the team is playing at home are included as covariates.
We enhance the flexibility in the design matrix ${\bf X}$, by considering
pairwise interactions between the teams and the other covariates in the source
dataset, thus obtaining $p=1762$ predictors, including the intercept, for
$760$ observations.

As for the simulation studies with count responses, we conduct Bayesian
inference by employing independent Gaussian priors with mean $0$ and variance
$\nu^{2}=250/p$ for each $\beta_{j}$ $(j=1,\ldots,1762)$. These priors are
updated with the Poisson likelihood (under logarithm link) of $n=700$ units,
after excluding $60$ observations to study the behavior of the predictive
means and posterior moments in such $p\gg n$ setting. In such a high-
dimensional scenario, both bpr and rstan incurred insurmountable mixing
issues. To obtain adequate mcmc inference, it was necessary to modify the
metric within the momentum term of hmc [7]. In particular, we run nuts via the
shell interface CmdStan of stan [25], initializing the mass matrix to the
approximate posterior variance learned by ep. Even so, nuts is considerably
slower than ep, despite running five shorter chains in parallel as done in the
simulation studies. To produce the results reported in Figure 4, the running
time of nuts is $\bf{46154.41}$ secs compared to the $\bf{19.56}$ secs of our
epp-eff implementation.

Figure 4: Results for the football data under Poisson regression with
$\nu^{2}=250/1762$. Scatterplot of posterior expectations and standard
deviations of the parameters and predictive means, i.e.,
${\mathbb{E}}[y_{n+1}\mid{\bf y}],\ldots,{\mathbb{E}}[y_{n+60}\mid{\bf y}]$,
for $60$ held-out units computed via nuts (x-axis) and via the efficient ep
implementation reported in Algorithm 2 (y-axis).

Results are depicted in Figure 4, where we present the scatter plots of the
posterior means ${\mathbb{E}}[\beta_{1}\mid{\bf
y}],\ldots,{\mathbb{E}}[\beta_{p}\mid{\bf y}]$, log standard deviations
$\log\big{(}\sqrt{\text{var}[\beta_{1}\mid{\bf
y}]}\,\big{)},\ldots,\,\allowbreak\log\big{(}\sqrt{\text{var}[\beta_{p}\mid{\bf
y}]}\,\big{)}$, and predictive means ${\mathbb{E}}[y_{n+1}\mid{\bf
y}],\ldots,{\mathbb{E}}[y_{n+60}\mid{\bf y}]$ obtained via the nuts sampler
and the ones resultant from the efficient ep approximation derived in
Algorithm 2 (with the efficient hybrid moment approximation for Poisson
computed in Section 4). Despite being computationally much faster, the ep
solution is extremely accurate for all the quantities of interest, which are
almost identical to the ones obtained with the nuts sampler in most of the
cases.

##  7 Discussion

This contribution develops a general and scalable implementation of ep for the
Bayesian glms under spherical Gaussian prior, making such approximation
feasible in scenarios that were not practically possible before [12]. Such an
algorithm provides an approximation of the posterior moments and the marginal
likelihood, which may be of practical use for model selection [34]. Moreover,
being the global approximation Gaussian, this also facilitates the computation
of ep predictive mean ${\mathbb{E}}[y_{{\textsc{new}}}\mid{\bf y}]$, which is
available in closed-form at no additional cost for the probit model (where it
coincides with the predictive probability $\Pr[y_{{\textsc{new}}}=1\mid{\bf
y}]$), and any glm with a logarithmic link (that is the canonical choice for
Poisson and Gamma likelihoods). We further improved computational scalability
by considering tailored approximations for the moment-matching step, while
preserving accuracy. As shown on simulated and real data in Sections 5 and 6,
this decreased cost leads to a massive reduction of the running time compared
to currently available implementations of ep, paving the way to applying such
a method in high-dimensional settings. Considering the quality of the
approximation empirically observed in the present work, in [12], and in [2],
this represents a major contribution to Bayesian data regression with
responses beyond the continuous case. In accordance with the robust evidence
in the literature, ep seems, in general, the preferable choice among
approximate methods. Finally, these empirical results can also motivate future
theoretical research about the goodness of the approximation provided by ep,
for which some preliminary findings on the asymptotic behavior have been
derived in [16] and [17].

## References

  * [1] James H. Albert and Siddhartha Chib “Bayesian analysis of binary and polychotomous response data” In _J. Am. Stat. Assoc._ 88 Taylor & Francis, 1993, pp. 669–679
  * [2] Niccolò Anceschi, Augusto Fasano, Daniele Durante and Giacomo Zanella “Bayesian conjugacy in probit, tobit, multinomial probit and extensions: a review and new results” In _J. Am. Stat. Assoc._ 118 Taylor & Francis, 2023, pp. 1451–1469
  * [3] Reinaldo B. Arellano-Valle and Adelchi Azzalini “On the unification of families of skew-normal distributions” In _Scand. J. Stat._ 33 Wiley, 2006, pp. 561–574
  * [4] Simon R Arridge, Kazufumi Ito, Bangti Jin and Chen Zhang “Variational Gaussian approximation for Poisson data” In _Inverse Probl._ 34 IOP Publishing, 2018, pp. 1–29
  * [5] Søren Asmussen, Jens Ledet Jensen and Leonardo Rojas-Nandayapa “On the Laplace transform of the lognormal distribution” In _Methodol. Comput. Appl. Probab._ 18 Springer ScienceBusiness Media LLC, 2014, pp. 441–458
  * [6] Adelchi Azzalini and Antonella Capitanio “The Skew-Normal and Related Families” UK: Cambridge Univ. Press, 2014
  * [7] Michael Betancourt and Leo C. Stein “The geometry of Hamiltonian Monte Carlo” In _Preprint at arXiv: 1112.4118_ , 2011
  * [8] Christopher M. Bishop “Pattern Recognition and Machine Learning” Springer, 2006
  * [9] David M Blei, Alp Kucukelbir and Jon D. McAuliffe “Variational inference: A review for statisticians” In _J. Am. Stat. Assoc._ 112, 2017, pp. 859–877
  * [10] Zdravko I Botev “The normal law under linear restrictions: simulation and estimation via minimax tilting” In _J. R. Stat. Soc. Series B Stat. Methodol._ 79 Oxford Univ. Press, 2017, pp. 125–148
  * [11] Charalampos Chanialidis, Ludger Evers, Tereza Neocleous and Agostino Nobile “Efficient Bayesian inference for COM-Poisson regression models” In _Stat. Comput._ 28 Springer, 2018, pp. 595–608
  * [12] Nicolas Chopin and James Ridgway “Leave Pima Indians alone: binary regression as a benchmark for Bayesian computation” In _Stat. Sci._ 32, 2017, pp. 64–87
  * [13] Guido Consonni and Jean-Michel Marin “Mean-field variational approximate Bayesian inference for latent variable models” In _Comput. Stat. Data Anal._ 52, 2007, pp. 790–798
  * [14] Rebecca Craig-Schapiro, Max Kuhn, Chengjie Xiong, Eve H. Pickering, Jingxia Liu, Thomas P. Misko, Richard J. Perrin, Kelly R. Bales, Holly Soares, Anne M. Fagan and David M. Holtzman “Multiplexed immunoassay panel identifies novel CSF biomarkers for Alzheimer’s disease diagnosis and prognosis” In _PloS one_ 6 Public Library of Science San Francisco, USA, 2011, pp. 1–15
  * [15] Laura D’Angelo and Antonio Canale “Efficient posterior sampling for Bayesian Poisson regression” In _J. Comput. Graph. Stat._ 32 Taylor & Francis, 2023, pp. 917–926
  * [16] Guillaume Dehaene and Simon Barthelmé “Bounding errors of expectation-propagation” In _Adv. Neural Inf. Process. Syst._ 28, 2015
  * [17] Guillaume Dehaene and Simon Barthelmé “Expectation propagation in the large data limit” In _J. R. Stat. Soc. Series B Stat. Methodol._ 80, 2018, pp. 199–217
  * [18] D. Durante “Conjugate Bayes for probit regression via unified skew-normal distributions” In _Biometrika_ 106, 2019, pp. 765–779
  * [19] Daniele Durante and Tommaso Rigon “Conditionally conjugate mean-field variational Bayes for logistic models” In _Stat. Sci._ 34 Institute of Mathematical Statistics, 2019, pp. 472–485
  * [20] A. Fasano, D. Durante and G. Zanella “Scalable and accurate variational Bayes for high-dimensional binary regression models” In _Biometrika_ 109, 2022, pp. 901–919
  * [21] Augusto Fasano, Niccolò Anceschi, Beatrice Franzolini and Giovanni Rebaudo “Efficient computation of predictive probabilities in probit models via expectation propagation” In _Book of Short Papers - CLADAG 2023_ , 2023, pp. 449–452
  * [22] Augusto Fasano, Niccolò Anceschi, Beatrice Franzolini and Giovanni Rebaudo “Efficient expectation propagation for posterior approximation in high-dimensional probit models” In _Book of Short Papers - SIS 2023_ , 2023, pp. 1133–1138
  * [23] Augusto Fasano and Daniele Durante “A class of conjugate priors for multinomial probit models which includes the multivariate normal one” In _J. Mach. Learn. Res._ 23, 2022, pp. 1–16
  * [24] Sylvia Frühwirth-Schnatter, Rudolf Frühwirth, Leonhard Held and Håvard Rue “Improved auxiliary mixture sampling for hierarchical models of non-Gaussian data” In _Stat. Comput._ 19 Springer, 2009, pp. 479–492
  * [25] Jonah Gabry, Rok Češnovar and Andrew Johnson “cmdstanr: R Interface to ’CmdStan”’, 2023
  * [26] Andrew Gelman, Aleks Jakulin, Maria Grazia Pittau and Yu-Sung Su “A weakly informative default prior distribution for logistic and other regression models” In _Ann. Appl. Stat._ 2 Institute of Mathematical Statistics, 2008, pp. 1360–1383
  * [27] Marc G. Genton, David E. Keyes and George Turkiyyah “Hierarchical decompositions for the computation of high-dimensional multivariate normal probabilities” In _J. Comput. Graph. Stat._ 27 Taylor & Francis, 2018, pp. 268–277
  * [28] Mark Girolami and Simon Rogers “Variational Bayesian multinomial probit regression with Gaussian process priors” In _Neural Comput._ 18 MIT Press, 2006, pp. 1790–1817
  * [29] Heikki Haario, Eero Saksman and Johanna Tamminen “An adaptive Metropolis algorithm” In _Bernoulli_ JSTOR, 2001, pp. 223–242
  * [30] Matthew D. Hoffman and Andrew Gelman “The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.” In _J. Mach. Learn. Res._ 15, 2014, pp. 1593–1623
  * [31] Chris C. Holmes and Leonhard Held “Bayesian auxiliary variable models for binary and multinomial regression” In _Bayesian Anal._ 1, 2006, pp. 145–168
  * [32] Tommi S Jaakkola and Michael I Jordan “Bayesian parameter estimation via variational methods” In _Stat. Comput._ 10 Springer, 2000, pp. 25–37
  * [33] James E. Johndrow, Aaron Smith, Natesh Pillai and David B. Dunson “MCMC for imbalanced categorical data” In _J. Am. Stat. Assoc._ 114 Taylor & Francis, 2019, pp. 1394–1403
  * [34] Robert E. Kass and Adrian E. Raftery “Bayes factors” In _J. Am. Stat. Assoc._ 90 Taylor & Francis, 1995, pp. 773–795
  * [35] Solomon Kullback and Richard A Leibler “On information and sufficiency” In _Ann. Stat._ 22, 1951, pp. 79–86
  * [36] David J.. MacKay “The evidence framework applied to classification networks” In _Neural Comput._ 4, 1992, pp. 720–736
  * [37] Thomas P. Minka “Expectation propagation for approximate Bayesian inference” In _Proceedings of the Seventeenth Conference on Uncertainty in Artificial Intelligence_ San Francisco, USA: Morgan Kaufmann Publishers Inc., 2001, pp. 362–369
  * [38] Thomas P. Minka “Power EP”, 2004
  * [39] Thomas Peter Minka “A Family of Algorithms for Approximate Bayesian Inference”, 2001
  * [40] Ari Pakman and Liam Paninski “Exact hamiltonian monte carlo for truncated multivariate gaussians” In _J. Comput. Graph. Stat._ 23 Taylor & Francis, 2014, pp. 518–542
  * [41] Nicholas G. Polson, James .G Scott and Jesse Windle “Bayesian inference for logistic models using Pólya–Gamma latent variables” In _J. Am. Stat. Assoc._ 108 Taylor & Francis, 2013, pp. 1339–1349
  * [42] Gareth O. Roberts and Jeffrey S. Rosenthal “Examples of adaptive MCMC” In _J. Comput. Graph. Stat._ 18 Taylor & Francis, 2009, pp. 349–367
  * [43] A.. Rossberg “Laplace transforms of probability distributions and their inversions are easy on logarithmic scales” In _J. Appl. Probab._ 45, 2008, pp. 531–541
  * [44] A.. Rossberg, T. Amemiya and K. Itoh “Accurate and fast approximations of moment-generating functions and their inversion for log-normal and similar distributions” Citeseer, 2008
  * [45] Håvard Rue, Sara Martino and Nicolas Chopin “Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations” In _J. R. Stat. Soc. Series B Stat. Methodol._ 71 Oxford University Press, 2009, pp. 319–392
  * [46] Matthias Seeger “Bayesian Gaussian Process Models: PAC-Bayesian Generalisation Error Bounds and Sparse Approximations”, 2003
  * [47] Matthias Seeger “Expectation propagation for exponential families” Univ. of California at Berkeley, 2005
  * [48] Matthias Seeger “Bayesian inference and optimal design for the sparse linear model” In _J. Mach. Learn. Res._ 9, 2008, pp. 759–813
  * [49] Matthias Seeger and Matthias Gerwinn “Bayesian inference for sparse generalized linear models” In _Mach. Learn: ECML_ , 2007, pp. 298–309
  * [50] Stan Development Team “RStan: the R interface to Stan” R package version 2.32.3, 2023
  * [51] Aki Vehtari, Andrew Gelman, Tuomas Sivula, Pasi Jylänki, Dustin Tran, Swupnil Sahai, Paul Blomstedt, John P. Cunningham, David Schiminovich and Christian P. Robert “Expectation propagation as a way of life: a framework for Bayesian inference on partitioned data” In _J. Mach. Learn. Res._ 21, 2020, pp. 577–629
  * [52] Gregor Zens, Sylvia Frühwirth-Schnatter and Helga Wagner “Ultimate Pólya gamma samplers – efficient MCMC for possibly imbalanced binary and categorical data” In _J. Am. Stat. Assoc._ in press Taylor & Francis, 2023
  * [53] Chen Zhang, Simon Arridge and Bangti Jin “Expectation propagation for Poisson data” In _Inverse Probl._ 35 IOP Publishing, 2019, pp. 1–27
  * [54] Zhenyu Zhang, Akihiko Nishimura, Paul Bastide, Xiang Ji, Rebecca P. Payne, Philip Goulder, Philippe Lemey and Marc A. Suchard “Large-scale inference of correlation among mixed-type biological traits with phylogenetic multivariate probit models” In _Ann. Appl. Stat._ 15, 2021, pp. 230–251
  * [55] Jackson Zhou, John T. Ormerod and Clara Grazian “Fast expectation propagation for heteroscedastic, lasso-penalized, and quantile regression” In _J. Mach. Learn. Res._ 24, 2023, pp. 1–39

## References

  * [56] Søren Asmussen, Jens Ledet Jensen and Leonardo Rojas-Nandayapa “On the Laplace transform of the lognormal distribution” In _Methodol. Comput. Appl. Probab._ 18 Springer ScienceBusiness Media LLC, 2014, pp. 441–458
  * [57] Adelchi Azzalini and Antonella Capitanio “The Skew-Normal and Related Families” UK: Cambridge Univ. Press, 2014
  * [58] Kaare Brandt Petersen and Michael Syskind Pedersen “The Matrix Cookbook”, 2008, pp. 1–72
  * [59] A.. Rossberg, T. Amemiya and K. Itoh “Accurate and fast approximations of moment-generating functions and their inversion for log-normal and similar distributions” Citeseer, 2008
  * [60] Jackson Zhou, John T. Ormerod and Clara Grazian “Fast expectation propagation for heteroscedastic, lasso-penalized, and quantile regression” In _J. Mach. Learn. Res._ 24, 2023, pp. 1–39

Supplementary Materials for
Scalable expectation propagation
for generalized linear models

Niccolò Anceschi1<EMAIL_ADDRESS>
Augusto Fasano2<EMAIL_ADDRESS>
Beatrice Franzolini3<EMAIL_ADDRESS>
Giovanni Rebaudo4<EMAIL_ADDRESS>1Duke University, Durham, USA

2Catholic University of the Sacred Heart, Milan, Italy

3Bocconi University, Milan, Italy

4University of Torino, Turin, Italy

##  S.1 Derivation of the results in Section 3

#### Proof of Lemma 1 and Proposition 1

Adapting the rationale from [60], the (mgf) $\mathcal{M}_{h_{i}}({\bf t})$ of
the hybrid distribution $h_{i}({\bm{\beta}})$ reads

$\begin{split}\mathcal{M}_{h_{i}}({\bf
t})&={\mathbb{E}}_{h_{i}({\bm{\beta}})}\big{[}e^{{\bm{\beta}}^{\intercal}{\bf
t}}\big{]}=\int e^{{\bm{\beta}}^{\intercal}{\bf
t}}\textstyle{\frac{1}{Z_{h_{i}}}}\ell_{i}({\bf
x}_{i}^{\intercal}{\bm{\beta}})\phi_{p}({\bm{\beta}}-{\bm{\xi}}_{i},{\bm{\Omega}}_{i})d{\bm{\beta}}\\\
&=\textstyle{\frac{1}{Z_{h_{i}}}}e^{{\bm{\xi}}_{i}^{\intercal}{\bf
t}+{\textstyle\frac{1}{2}}{\bf t}^{\intercal}{\bm{\Omega}}_{i}{\bf
t}}\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf t}),{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}d\eta_{i}\;,\end{split}$

for any ${\bf t}\in{\mathbb{R}}^{p}$, with ${\bm{\Omega}}_{i}={\bf
Q}_{-i}^{-1}$ and ${\bm{\xi}}_{i}={\bf Q}_{-i}^{-1}{\bf r}_{-i}$. The last
identity follows from $\eta_{i}={\bm{\xi}}_{i}^{\intercal}{\bm{\theta}}$ and
the linear subspace properties of Gaussian measure. Accordingly, the moments
of the hybrid can be obtained as

$\begin{split}{\bm{\mu}}_{h_{i}}&={\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]=\nabla_{\bf
t}\mathcal{M}_{h_{i}}({\bf 0})\\\
{\bm{\Sigma}}_{h_{i}}&=\operatorname{var}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]=\nabla_{\bf
t}^{2}\mathcal{M}_{h_{i}}({\bf 0})-\nabla_{\bf t}\mathcal{M}_{h_{i}}({\bf
0})\nabla_{\bf t}^{\intercal}\mathcal{M}_{h_{i}}({\bf 0})\;.\end{split}$

For ease of notation, let us define $\mathcal{F}_{h_{i}}({\bf
t})=\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf t}),{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}d\eta_{i}$ and
$\mathcal{G}_{h_{i}}({\bf t})=e^{{\bm{\xi}}_{i}^{\intercal}{\bf
t}+{\textstyle\frac{1}{2}}{\bf t}^{\intercal}{\bm{\Omega}}_{i}{\bf t}}$, so
that

$\begin{split}\nabla_{\bf t}\mathcal{M}_{h_{i}}({\bf
t})&=\textstyle{\frac{1}{Z_{h_{i}}}}\mathcal{F}_{h_{i}}({\bf t})\nabla_{\bf
t}\mathcal{G}_{h_{i}}({\bf
t})+\textstyle{\frac{1}{Z_{h_{i}}}}\mathcal{G}_{h_{i}}({\bf t})\nabla_{\bf
t}\mathcal{F}_{h_{i}}({\bf t})\\\ \nabla_{\bf t}^{2}\mathcal{M}_{h_{i}}({\bf
t})&=\textstyle{\frac{1}{Z_{h_{i}}}}\mathcal{F}_{h_{i}}({\bf t})\nabla_{\bf
t}^{2}\mathcal{G}_{h_{i}}({\bf
t})+\textstyle{\frac{1}{Z_{h_{i}}}}\mathcal{G}_{h_{i}}({\bf t})\nabla_{\bf
t}^{2}\mathcal{F}_{h_{i}}({\bf t})\\\
&\quad+\textstyle{\frac{1}{Z_{h_{i}}}}\nabla_{\bf t}\mathcal{F}_{h_{i}}({\bf
t})\nabla_{\bf t}^{\intercal}\mathcal{G}_{h_{i}}({\bf
t})+\textstyle{\frac{1}{Z_{h_{i}}}}\nabla_{\bf t}\mathcal{G}_{h_{i}}({\bf
t})\nabla_{\bf t}^{\intercal}\mathcal{F}_{h_{i}}({\bf t})\;.\end{split}$

Cleary, it holds that $\mathcal{F}_{h_{i}}({\bf 0})=Z_{h_{i}}$ and
$\mathcal{G}_{h_{i}}({\bf 0})=1$, while it can be seen that $\nabla_{\bf
t}\mathcal{G}_{h_{i}}({\bf t})=({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf
t})\mathcal{G}_{h_{i}}({\bf t})$ and $\nabla_{\bf
t}^{2}\mathcal{G}_{h_{i}}({\bf t})={\bm{\Omega}}_{i}\mathcal{G}_{h_{i}}({\bf
t})+({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf
t})({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf
t})^{\intercal}\mathcal{G}_{h_{i}}({\bf t})$. Accordingly, one has

$\nabla_{\bf t}\mathcal{G}_{h_{i}}({\bf
0})={\bm{\xi}}_{i}\qquad\qquad\nabla_{\bf t}^{2}\mathcal{G}_{h_{i}}({\bf
0})={\bm{\Omega}}_{i}+{\bm{\xi}}_{i}{\bm{\xi}}_{i}^{\intercal}\;.$

Conversely,

$\begin{split}\nabla_{\bf t}\mathcal{F}_{h_{i}}({\bf
t})&=\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf t}),{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}\,\\\
&\qquad\Big{(}{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}-{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
t}\big{)}\Big{)}d\eta_{i}\\\ \nabla_{\bf t}^{2}\mathcal{F}_{h_{i}}({\bf
t})&=-\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf t}),{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}\,{\bm{\Omega}}_{i}{\bf
x}_{i}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\,d\eta_{i}\\\
&\quad\,+\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}({\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf t}),{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\big{)}\,\\\
&\qquad\quad\Big{(}{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}-{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
t}\big{)}\Big{)}\\\ &\qquad\quad\Big{(}{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}\big{(}\eta_{i}-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}-{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
t}\big{)}\Big{)}^{\intercal}d\eta_{i}\;,\end{split}$

so that

$\begin{split}\nabla_{\bf t}\mathcal{F}_{h_{i}}({\bf
0})&=Z_{h_{i}}{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}\big{)}\\\ \nabla_{\bf
t}^{2}\mathcal{F}_{h_{i}}({\bf 0})&=-Z_{h_{i}}{\bm{\Omega}}_{i}{\bf
x}_{i}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}+Z_{h_{i}}{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}\;\cdot\\\
&\quad\;\cdot\Big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}\big{[}({\bf
x}_{i}^{\intercal}{\bm{\beta}})^{2}\big{]}-2{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]{\bf x}_{i}^{\intercal}{\bm{\xi}}_{i}+({\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i})^{2}\Big{)}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\;.\end{split}$

Putting it all together, one gets

$\begin{split}\nabla_{\bf t}&\mathcal{M}_{h_{i}}({\bf
0})={\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}\big{)}\\\ \nabla_{\bf
t}^{2}&\mathcal{M}_{h_{i}}({\bf
0})={\bm{\xi}}_{i}{\bm{\xi}}_{i}^{\intercal}+{\bm{\Omega}}_{i}-{\bm{\Omega}}_{i}{\bf
x}_{i}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\\\ &+{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\Big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}\big{[}({\bf
x}_{i}^{\intercal}{\bm{\beta}})^{2}\big{]}-2{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]{\bf x}_{i}^{\intercal}{\bm{\xi}}_{i}+({\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i})^{2}\Big{)}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\\\ &+\big{(}{\bm{\Omega}}_{i}{\bf
x}_{i}{\bm{\xi}}_{i}^{\intercal}+{\bm{\xi}}_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\big{)}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}\big{)}\end{split}$

and

$\begin{split}\nabla_{\bf t}^{2}&\mathcal{M}_{h_{i}}({\bf 0})-\nabla_{\bf
t}\mathcal{M}_{h_{i}}({\bf 0})\nabla_{\bf
t}^{\intercal}\mathcal{M}_{h_{i}}({\bf
0})={\bm{\Omega}}_{i}-{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\\\ &+{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\Big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}\big{[}({\bf
x}_{i}^{\intercal}{\bm{\beta}})^{2}\big{]}-{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]^{2}\Big{)}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\;.\end{split}$

or equivalently that

$\begin{split}{\bm{\mu}}_{h_{i}}&={\bm{\xi}}_{i}+{\bm{\Omega}}_{i}{\bf
x}_{i}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf
x}_{i}^{\intercal}{\bm{\xi}}_{i}\big{)}\\\
{\bm{\Sigma}}_{h_{i}}&={\bm{\Omega}}_{i}+{\bm{\Omega}}_{i}{\bf x}_{i}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}\big{(}{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}\big{)}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}\;.\end{split}$

Exploiting Woodbury’s identity, the natural parameters for the updated global
ep approximant read

$\begin{split}{\bf
Q}_{\textsc{ep}}^{{\textsc{new}}}&=\big{(}{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]\big{)}^{-1}\\\
&={\bm{\Omega}}_{i}^{-1}-{\bf x}_{i}\Big{(}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})\big{(}{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}\big{)}^{-1}({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})\Big{)}^{-1}{\bf x}_{i}^{\intercal}\\\ &={\bm{\Omega}}_{i}^{-1}-{\bf
x}_{i}{\bf x}_{i}^{\intercal}\big{(}({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}-{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]^{-1}\big{)}\end{split}$

and

$\begin{split}{\bf
r}_{\textsc{ep}}^{{\textsc{new}}}&=\big{(}{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]\big{)}^{-1}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]\\\
&={\bm{\Omega}}_{i}^{-1}{\bm{\xi}}_{i}+{\bf x}_{i}{\bf
x}_{i}^{\intercal}\big{(}{\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]^{-1}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bf
x}_{i}^{\intercal}{\bm{\beta}}]-({\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}{\bf x}_{i}^{\intercal}{\bm{\xi}}_{i}\big{)}\;.\end{split}$

Together with the definition of the cavity parameters as
${\bm{\Omega}}_{i}={\bf Q}_{-i}^{-1}$ and ${\bm{\xi}}_{i}={\bf
Q}_{-i}^{-1}{\bf r}_{-i}$, this gives the one-dimensional updates for the
$i$-th site, reported in Proposition 1.

#### Proof of Proposition 2

Specifying equation (3),

$\displaystyle\log Z_{i}^{{\textsc{new}}}$ $\displaystyle=\dfrac{1}{2}{\bf
r}_{{\textsc{ep}}}^{{\textsc{new}}\
\intercal}{\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}{\bf
r}_{{\textsc{ep}}}^{{\textsc{new}}}+\dfrac{p}{2}\log(2\pi)-\dfrac{1}{2}\log|{\bf
Q}_{{\textsc{ep}}}^{{\textsc{new}}}|$ (S.11)
$\displaystyle\quad\quad-\dfrac{1}{2}{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
r}_{-i}-\dfrac{p}{2}\log(2\pi)+\dfrac{1}{2}\log|{\bf Q}_{-i}|-\log Z_{h_{i}}$
$\displaystyle=\dfrac{1}{2}\left({\bf r}_{{\textsc{ep}}}^{{\textsc{new}}\
\intercal}{\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}{\bf
r}_{{\textsc{ep}}}^{{\textsc{new}}}-{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
r}_{-i}\right)-\dfrac{1}{2}\left(\log|{\bf
Q}_{{\textsc{ep}}}^{{\textsc{new}}}|-\log|{\bf Q}_{-i}|\right)-\log
Z_{h_{i}}.$

Now, for the first parenthesis, exploiting Woodbury’s identity, one gets

$\begin{split}&={\bf r}_{{\textsc{ep}}}^{{\textsc{new}}\
\intercal}{\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}{\bf
r}_{{\textsc{ep}}}^{{\textsc{new}}}-{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}\\\ &=({\bf
r}_{-i}+m_{i}^{{\textsc{new}}}{\bf x}_{i})^{\intercal}({\bf
Q}_{-i}+k_{i}^{{\textsc{new}}}{\bf x}_{i}{\bf x}_{i}^{\intercal})^{-1}({\bf
r}_{-i}+m_{i}^{{\textsc{new}}}{\bf x}_{i})-{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}\\\ &=({\bf
r}_{-i}+m_{i}^{{\textsc{new}}}{\bf
x}_{i})^{\intercal}\left[{\bm{\Omega}}_{i}-\dfrac{k_{i}^{{\textsc{new}}}}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}{\bm{\Omega}}_{i}{\bf
x}_{i}{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}\right]({\bf
r}_{-i}+m_{i}^{{\textsc{new}}}{\bf x}_{i})\\\
&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}\\\ &={\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}+2m_{i}^{{\textsc{new}}}{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}+(m_{i}^{{\textsc{new}}})^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\\\
&\quad\quad-\dfrac{k_{i}^{{\textsc{new}}}}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}\Big{[}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{2}+2m_{i}^{{\textsc{new}}}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})({\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})\\\
&\qquad\qquad\qquad\qquad\qquad\qquad+(m_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{2}\Big{]}-{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf r}_{-i}\\\
&=2m_{i}^{{\textsc{new}}}{\bf r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}\left[1-\dfrac{k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}}\right]-\dfrac{k_{i}^{{\textsc{new}}}}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{2}\\\
&\quad\quad+(m_{i}^{{\textsc{new}}})^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}\left[1-\dfrac{k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}\right]\\\
&=\dfrac{1}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}}\left[2m_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
r}_{-i}+(m_{i}^{{\textsc{new}}})^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}-k_{i}^{{\textsc{new}}}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{2}\right]\end{split}$

As for the second parenthesis, exploiting standard properties of determinants
of matrices [58] it holds

$\begin{split}\log|{\bf Q}_{{\textsc{ep}}}^{{\textsc{new}}}|-\log|{\bf
Q}_{-i}|&=\log|{\bf Q}_{-i}+k_{i}^{{\textsc{new}}}{\bf x}_{i}{\bf
x}_{i}^{\intercal}|-\log|{\bf Q}_{-i}|\\\ &=\log|{\bf Q}_{-i}[{\bf
I}_{p}+k_{i}^{{\textsc{new}}}{\bm{\Omega}}_{i}{\bf x}_{i}{\bf
x}_{i}^{\intercal}]|\\\ &=\log|{\bf Q}_{-i}|+\log|{\bf
I}_{p}+k_{i}^{{\textsc{new}}}{\bm{\Omega}}_{i}{\bf x}_{i}{\bf
x}_{i}^{\intercal}|-\log|{\bf Q}_{-i}|\\\ &=\log(1+k^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}).\end{split}$

Substituting the above derivations in (S.11), one finally obtains

$\begin{split}\log
Z_{i}^{{\textsc{new}}}=\dfrac{1}{2}\bigg{[}&\dfrac{1}{1+k_{i}^{{\textsc{new}}}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}}\left(2m_{i}^{{\textsc{new}}}{\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i}+(m_{i}^{{\textsc{new}}})^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}-k_{i}^{{\textsc{new}}}({\bf
r}_{-i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{2}\right)\\\
&-\log(1+k_{i}^{{\textsc{new}}}{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})\bigg{]}-\log Z_{h_{i}},\end{split}$

as claimed. ∎

#### Proof of Proposition 3

By Woodbury’s identity, one gets

$\begin{split}{\bf w}_{i}&={\bf Q}_{-i}^{-1}{\bf x}_{i}=({\bf Q}-{\bf
Q}_{i})^{-1}{\bf x}_{i}={\bf Q}^{-1}{\bf x}_{i}+(1-k_{i}{\bf
x}_{i}^{\intercal}{\bf Q}^{-1}{\bf x}_{i})^{-1}k_{i}({\bf Q}^{-1}{\bf
x}_{i})({\bf Q}^{-1}{\bf x}_{i})^{\intercal}{\bf x}_{i}\\\ &={\bf
v}_{i}+k_{i}(1-k_{i}{\bf x}_{i}^{\intercal}{\bf v}_{i})^{-1}{\bf v}_{i}{\bf
v}_{i}^{\intercal}{\bf x}_{i}=\big{[}1+(1-k_{i}{\bf x}_{i}^{\intercal}{\bf
v}_{i})^{-1}(k_{i}{\bf x}_{i}^{\intercal}{\bf v}_{i})\big{]}{\bf
v}_{i}=d_{i}{\bf v}_{i},\end{split}$

where $d_{i}=(1-k_{i}{\bf x}_{i}^{\intercal}{\bf v}_{i})^{-1}$.

When site $i$ is updated, ${\bf Q}_{i}$ changes and consequently also
${\bm{\Omega}}_{{\textsc{ep}}}$ does, assuming a new value
${\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}=({\bf Q}_{-i}+{\bf
Q}_{i}^{{\textsc{new}}})^{-1}$. Thus also the ${\bf v}_{j}$’s, $j=1,\ldots,n$,
should be updated. Exploiting the form of
${\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}$, together with Woodbury’s
identity, one gets

$\begin{split}{\bf
v}_{j}^{{\textsc{new}}}&={\bm{\Omega}}_{{\textsc{ep}}}^{{\textsc{new}}}{\bf
x}_{j}=({\bf Q}_{{\textsc{ep}}}-{\bf Q}_{i}+{\bf
Q}_{i}^{{\textsc{new}}})^{-1}{\bf x}_{j}=[{\bf
Q}_{{\textsc{ep}}}+(k_{i}^{{\textsc{new}}}-k_{i}){\bf x}_{i}{\bf
x}_{i}^{\intercal}]^{-1}{\bf x}_{j}\\\
&=[{\bm{\Omega}}_{{\textsc{ep}}}-(k_{i}^{{\textsc{new}}}-k_{i})[1+(k_{i}^{{\textsc{new}}}-k_{i}){\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{i}]^{-1}{\bm{\Omega}}_{{\textsc{ep}}}{\bf x}_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{{\textsc{ep}}}]{\bf x}_{j}\\\
&={\bm{\Omega}}_{{\textsc{ep}}}{\bf
x}_{j}-[(k_{i}^{{\textsc{new}}}-k_{i})^{-1}+{\bf x}_{i}^{\intercal}{\bf
v}_{i}]^{-1}{\bf v}_{i}{\bf x}_{i}^{\intercal}{\bf v}_{j}={\bf
v}_{j}-c_{i}({\bf x}_{i}^{\intercal}{\bf v}_{j}){\bf v}_{i},\end{split}$

where
$c_{i}=(k_{i}^{{\textsc{new}}}-k_{i})/(1+(k_{i}^{{\textsc{new}}}-k_{i}){\bf
x}_{i}^{\intercal}{\bf v}_{i})$.∎

##  S.2 Derivation of the results in Section 4

#### Proof of equation (8)

Recalling that here $\eta_{i}={\bf x}_{i}^{\intercal}{\bm{\beta}}$ and
exploring differentiation under integral sign, one has

$\begin{split}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}&=\frac{1}{Z_{h_{i}}}\frac{\partial
Z_{h_{i}}}{\partial\lambda_{i}}=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\frac{\partial}{\partial\lambda_{i}}\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}d\eta_{i}\\\
&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}\frac{\partial}{\partial\lambda_{i}}\Big{(}-\frac{(\eta_{i}-\lambda_{i})^{2}}{2\rho_{i}^{2}}\Big{)}d\eta_{i}\\\
&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}\frac{(\eta_{i}-\lambda_{i})}{\rho_{i}^{2}}d\eta_{i}=\frac{1}{\rho_{i}^{2}}\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}]-\lambda_{i}\big{)}\;.\end{split}$

A simple rearrangement gives
${\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}]=\lambda_{i}+\rho_{i}^{2}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}$. Similarly, it holds that

$\begin{split}\frac{\partial}{\partial\rho_{i}^{2}}\log
Z_{h_{i}}&=\frac{1}{Z_{h_{i}}}\frac{\partial
Z_{h_{i}}}{\partial\rho_{i}^{2}}=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\frac{\partial}{\partial\rho_{i}^{2}}\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}d\eta_{i}\\\
&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}\Big{(}-\frac{1}{2\rho_{i}^{2}}-\frac{\partial}{\partial\rho_{i}^{2}}\frac{(\eta_{i}-\lambda_{i})^{2}}{2\rho_{i}^{2}}\Big{)}d\eta_{i}\\\
&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}\frac{1}{2\rho_{i}^{2}}\Big{(}\frac{(\eta_{i}-\lambda_{i})^{2}}{2\rho_{i}^{2}}-1\Big{)}d\eta_{i}\\\
&=\frac{1}{2\rho_{i}^{2}}\Big{(}\frac{1}{\rho_{i}^{2}}\,{\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}^{2}]-2\,\frac{1}{\rho_{i}^{2}}\,\lambda_{i}\,{\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}]+\frac{1}{\rho_{i}^{2}}\,\lambda_{i}^{2}\Big{)}-\frac{1}{2\rho_{i}^{2}}\\\
&=\frac{1}{2\rho_{i}^{2}}\Big{(}\frac{1}{\rho_{i}^{2}}\,\operatorname{var}_{h_{i}({\bm{\beta}})}[\eta_{i}]+\frac{1}{\rho_{i}^{2}}\,\big{(}{\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}]-\lambda_{i}\big{)}^{2}\Big{)}-\frac{1}{2\rho_{i}^{2}}\\\
&=\frac{1}{2\rho_{i}^{2}}\Big{(}\frac{1}{\rho_{i}^{2}}\,\operatorname{var}_{h_{i}({\bm{\beta}})}[\eta_{i}]+\rho_{i}^{2}\big{(}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}\big{)}^{2}\Big{)}-\frac{1}{2\rho_{i}^{2}}\;.\end{split}$

Trivial rearrangements give
$\operatorname{var}_{h_{i}({\bm{\beta}})}[\eta_{i}]=\rho_{i}^{2}+\rho_{i}^{2}\Big{(}2\frac{\partial}{\partial\rho_{i}^{2}}\log
Z_{h_{i}}-\big{(}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}\big{)}^{2}\Big{)}\rho_{i}^{2}$. We note that the expression for the
variance could have been equivalently derived by studying the second
derivative of $\log Z_{h_{i}}$ with respect to $\lambda_{i}^{2}$

$\frac{\partial^{2}}{\partial\lambda_{i}^{2}}\log
Z_{h_{i}}=\frac{\partial}{\partial\lambda_{i}}\Big{(}\frac{\partial}{\partial\lambda_{i}}\log
Z_{h_{i}}\Big{)}=\frac{1}{Z_{h_{i}}}\frac{\partial^{2}Z_{h_{i}}}{\partial\lambda_{i}^{2}}-\Big{(}\frac{1}{Z_{h_{i}}}\frac{\partial
Z_{h_{i}}}{\partial\lambda_{i}}\Big{)}^{2}\;.$

Working out the first term

$\begin{split}\frac{1}{Z_{h_{i}}}\frac{\partial^{2}Z_{h_{i}}}{\partial\lambda_{i}^{2}}&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\frac{\partial^{2}}{\partial\lambda_{i}^{2}}\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}d\eta_{i}\\\
&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\frac{\partial}{\partial\lambda_{i}}\Big{(}\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}\frac{(\eta_{i}-\lambda_{i})}{\rho_{i}^{2}}\Big{)}d\eta_{i}\\\
&=\frac{1}{Z_{h_{i}}}\int\ell_{i}(\eta_{i})\phi\big{(}\eta_{i}-\lambda_{i},\,\rho_{i}^{2}\big{)}\frac{1}{\rho_{i}^{2}}\Big{(}\frac{(\eta_{i}-\lambda_{i})^{2}}{\rho_{i}^{2}}-1\Big{)}d\eta_{i}\\\
&=\frac{1}{\rho_{i}^{2}}\Big{(}\frac{1}{\rho_{i}^{2}}\,{\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}^{2}]-2\,\frac{1}{\rho_{i}^{2}}\,\lambda_{i}\,{\mathbb{E}}_{h_{i}({\bm{\beta}})}[\eta_{i}]+\frac{1}{\rho_{i}^{2}}\,\lambda_{i}^{2}\Big{)}-\frac{1}{\rho_{i}^{2}}\;,\end{split}$

one has that $\frac{\partial}{\partial\rho_{i}^{2}}\log
Z_{h_{i}}=\frac{1}{2}\frac{\partial^{2}}{\partial\lambda_{i}^{2}}\log
Z_{h_{i}}$. Nonetheless, we note that this last equality holds true in the
exact case, but it typically breaks down for the approximations considered in
Section 4. In such cases, dedicated analysis could be performed to assess
empirically which of the two alternatives provides a better approximation of
the exact hybrid moments.

#### Proof of Proposition 4

By definition, the hybrid distribution $h_{i}({\bm{\beta}})$ is proportional
to

$p(y_{i}\mid{\bm{\beta}})q_{-i}({\bm{\beta}})=\Phi((2y_{i}-1){\bf
x}_{i}^{\intercal}{\bm{\beta}})\phi_{p}({\bm{\beta}}-{\bf Q}_{-i}^{-1}{\bf
r}_{-i},{\bf Q}_{-i}^{-1}).$

Recall that ${\bm{\xi}}_{i}={\bf Q}_{-i}^{-1}{\bf r}_{-i}$ and
${\bm{\Omega}}_{i}={\bf Q}_{-i}^{-1}$, and define

$\begin{split}{\bm{\omega}}_{i}&=\left[\text{diag}\left({\bm{\Omega}}_{i}\right)\right]^{1/2}\hskip
50.0pt\bar{{\bm{\Omega}}}_{i}={\bm{\omega}}_{i}^{-1}{\bm{\Omega}}_{i}{\bm{\omega}}_{i}^{-1}\\\
{\bm{\alpha}}_{i}&=(2y_{i}-1){\bm{\omega}}_{i}{\bf x}_{i}\hskip
55.0pt\tau_{i}=(2y_{i}-1)(1+{\bf x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1/2}{\bf x}_{i}^{\intercal}{\bm{\xi}}_{i}\;.\end{split}$

After some trivial algebraic manipulations, one gets that

$p(y_{i}\mid{\bm{\beta}})q_{-i}({\bm{\beta}})=\Phi\left(\tau_{i}(1+{\bm{\alpha}}_{i}^{\intercal}\bar{{\bm{\Omega}}}_{i}{\bm{\alpha}}_{i})^{1/2}+{\bm{\alpha}}_{i}^{\intercal}{\bm{\omega}}_{i}^{-1}({\bm{\beta}}-{\bm{\xi}}_{i})\right)\phi_{p}\left({\bm{\beta}}-{\bm{\xi}}_{i},{\bm{\Omega}}_{i}\right),$

which is the kernel of a multivariate extended skew-normal distribution
$\textsc{sn}_{p}({\bm{\xi}}_{i},{\bm{\Omega}}_{i},{\bm{\alpha}}_{i},\tau_{i})$
(see Equation (5.61) in [57]), with normalizing constant
$Z_{h_{i}}=\Phi(\tau_{i})$. Furthermore, exploiting formulae (5.71) and (5.72)
in [57], the first two moments of $h_{i}({\bm{\beta}})$ are readily available:

$\begin{split}{\bm{\mu}}_{h_{i}}&={\mathbb{E}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]={\bm{\xi}}_{i}+\zeta_{1}(\tau_{i})s_{i}{\bm{\Omega}}_{i}{\bf
x}_{i},\\\
{\bm{\Sigma}}_{h_{i}}&={\operatorname{var}}_{h_{i}({\bm{\beta}})}[{\bm{\beta}}]={\bm{\Omega}}_{i}+\zeta_{2}(\tau_{i})s_{i}^{2}({\bm{\Omega}}_{i}{\bf
x}_{i})({\bm{\Omega}}_{i}{\bf x}_{i})^{\intercal},\end{split}$

This allows for an alternative derivation of the ep updates for the $i$-th
site one-dimensional parameters, without having to resort to the mgf of the
hybrid distribution. Appealing once more to Woodbury’s identity, one has

$\begin{split}{\bf Q}_{i}^{{\textsc{new}}}&={\bm{\Sigma}}_{h_{i}}^{-1}-{\bf
Q}_{-i}\\\
&={\bm{\Omega}}_{i}^{-1}-\zeta_{2}(\tau_{i})s_{i}^{2}(1+\zeta_{2}(\tau_{i})s_{i}^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bm{\Omega}}_{i}^{-1}{\bm{\Omega}}_{i}{\bf
x}_{i})^{-1}{\bm{\Omega}}_{i}^{-1}{\bm{\Omega}}_{i}{\bf x}_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bm{\Omega}}_{i}^{-1}-{\bf Q}_{-i}\\\
&=-\zeta_{2}(\tau_{i})s_{i}^{2}(1+\zeta_{2}(\tau_{i})s_{i}^{2}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf x}_{i}{\bf
x}_{i}^{\intercal}=-(\zeta_{2}(\tau_{i})^{-1}s_{i}^{-2}+{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i})^{-1}{\bf x}_{i}{\bf
x}_{i}^{\intercal}\\\ &=-\dfrac{\zeta_{2}(\tau_{i})}{1+{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}+\zeta_{2}(\tau_{i}){\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}}{\bf x}_{i}{\bf
x}_{i}^{\intercal}=k_{i}^{{\textsc{new}}}{\bf x}_{i}{\bf
x}_{i}^{\intercal},\end{split}$

with $k_{i}^{{\textsc{new}}}=-\zeta_{2}(\tau_{i})/\left(1+{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}+\zeta_{2}(\tau_{i}){\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\right)$. Then, exploiting the
previous result, one gets

$\begin{split}{\bf
r}_{i}^{{\textsc{new}}}&={\bm{\Sigma}}_{h_{i}}^{-1}{\bm{\mu}}_{h_{i}}-{\bf
r}_{-i}\\\ &={\bf Q}_{-i}{\bm{\mu}}_{h_{i}}+{\bf
Q}_{i}^{{\textsc{new}}}{\bm{\mu}}_{h_{i}}-{\bf r}_{-i}={\bf Q}_{-i}{\bf
Q}_{-i}^{-1}{\bf r}_{-i}+\zeta_{1}(\tau_{i})s_{i}{\bf
Q}_{-i}{\bm{\Omega}}_{i}{\bf x}_{i}+{\bf
Q}_{i}^{{\textsc{new}}}{\bm{\mu}}_{h_{i}}-{\bf r}_{-i}\\\
&=\zeta_{1}(\tau_{i})s_{i}{\bf x}_{i}+{\bf
Q}_{i}^{{\textsc{new}}}{\bm{\mu}}_{h_{i}}=\zeta_{1}(\tau_{i})s_{i}{\bf
x}_{i}+k_{i}^{{\textsc{new}}}{\bf x}_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf
r}_{-i}+k_{i}^{{\textsc{new}}}\zeta_{1}(\tau_{i})s_{i}{\bf x}_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}\\\
&=[\zeta_{1}(\tau_{i})s_{i}+k_{i}^{{\textsc{new}}}({\bm{\Omega}}_{i}{\bf
x}_{i})^{\intercal}{\bf
r}_{-i}+k_{i}^{{\textsc{new}}}\zeta_{1}(\tau_{i})s_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}]{\bf
x}_{i}=m_{i}^{{\textsc{new}}}{\bf x}_{i},\end{split}$

where
$m_{i}^{{\textsc{new}}}=\zeta_{1}(\tau_{i})s_{i}+k_{i}^{{\textsc{new}}}({\bm{\Omega}}_{i}{\bf
x}_{i})^{\intercal}{\bf
r}_{-i}+k_{i}^{{\textsc{new}}}\zeta_{1}(\tau_{i})s_{i}{\bf
x}_{i}^{\intercal}{\bm{\Omega}}_{i}{\bf x}_{i}$. ∎

#### Approximation of the log-normal Laplace transform from [56]

In this section, we work out the details of the approximation of the Laplace
transform for log-normal random variable
$\mathcal{L}og\mathcal{N}orm(0,\rho^{2})$, developed in [56]. The latter is
given by the expectation

$\mathcal{L}(s)={\mathbb{E}}\big{[}e^{-se^{X_{0}}}\big{]}=\int_{-\infty}^{\infty}\frac{1}{\sqrt{2\pi\rho^{2}}}\exp{\Big{(}-se^{x_{0}}-\frac{x_{0}^{2}}{2\rho^{2}}\Big{)}}dx_{0}$

where $X_{0}\sim\mathcal{N}(0,\rho^{2})$, which is defined only for $s>0$.
Extensions to expectations of the form
${\mathbb{E}}\big{[}e^{-se^{X}}\big{]}$, for
$X\sim\mathcal{N}(\varphi,\rho^{2})$, can be then trivially obtained as
$\mathcal{L}(se^{\varphi})$. The authors essentially suggest applying the
standard Laplace method to the integral corresponding to the required
expectation. In general, this gives

$\int_{a}^{b}\frac{1}{\sqrt{2\pi}}e^{-\gamma f(x_{0})}dx_{0}=\frac{e^{-\gamma
f(\hat{x}_{0})}}{\sqrt{\gamma
f^{\prime\prime}(\hat{x}_{0})}}\big{(}1+\mathcal{O}(\gamma^{-1})\big{)}\qquad\text{as}\qquad\gamma\rightarrow\infty\;,$

provided that $f(x_{0})$ has a unique global minimum at $\hat{x}_{0}\in(a,b)$.
Nonetheless, this approximation can be very accurate not only in asymptotic
settings but also for relatively small values of $\gamma$ [56]. Although in
our context we lack an explicit parameter $\gamma$, an analogous role can be
assigned to the second derivative of $f(\cdot)$ at its minimum. In the case of
$\mathcal{L}(s)$, $\gamma f(x_{0})$ gets replaced by
$g(x_{0};s,\rho^{2})=se^{x_{0}}+\frac{x_{0}^{2}}{2\rho^{2}}$, whose
derivatives are trivially obtained as
$g^{\prime}(x_{0};s,\rho^{2})=se^{x_{0}}+\frac{x_{0}}{\rho^{2}}$ and
$g^{\prime\prime}(x_{0};s,\rho^{2})=se^{x_{0}}+\frac{1}{\rho^{2}}$. Note that
$g(x_{0};s,\rho^{2})$ is convex in $x_{0}$ if $s>0$. Hence, its global minimum
$\hat{x}_{0}$ is obtained by setting to zero the first derivative, which gives
$-\hat{x}_{0}e^{-\hat{x}_{0}}=s\rho^{2}$. Recalling that the Lambert W
function is defined as the solution of $\mathcal{W}(x)e^{\mathcal{W}(x)}=x$,
one gets that $\hat{x}_{0}=-\mathcal{W}(s\rho^{2})$. Accordingly

$\begin{split}g(\hat{x}_{0};s,\rho^{2})&=\frac{1}{\rho^{2}}\mathcal{W}(s\rho^{2})+\frac{1}{2\rho^{2}}\mathcal{W}(s\rho^{2})^{2}\\\
g^{\prime\prime}(\hat{x}_{0};s,\rho^{2})&=\frac{1+\mathcal{W}(s\rho^{2})}{\rho^{2}}\;.\end{split}$

Putting it all together, we Laplace method applied to $\mathcal{L}(s)$ gives
the approximation

$\widetilde{\mathcal{L}}(s)=\frac{1}{\sqrt{1+\mathcal{W}(s\rho^{2})}}\exp{\Big{(}-\frac{1}{\rho^{2}}\mathcal{W}(s\rho^{2})-\frac{1}{2\rho^{2}}\mathcal{W}(s\rho^{2})^{2}\Big{)}}\;.$

We further highlight that $\frac{\partial}{\partial
x}\mathcal{W}(x)=\frac{\mathcal{W}(x)}{1+\mathcal{W}(x)}\frac{1}{x}$, which is
going to be useful in the next two Sections for the computations of the
derivatives of $\widetilde{\mathcal{L}}(s)$.

#### Adaptation to approximate moments matching in Poisson regression

Exploiting the results presented in the previous appendix, it is easy to see
that the intractable integral in the normalizing constant of ep hybrids for
Poisson regression in (11) coincides with $\mathcal{L}(s_{i})$ for
$s_{i}=\exp{(\lambda_{i}+\rho_{i}^{2}y_{i})}$. Accordingly, we approximate it
with $\widetilde{\mathcal{L}}(s_{i})$. As a consequence, one has that

$\begin{split}\frac{\partial}{\partial\lambda_{i}}\log\widetilde{Z}_{h_{i}}&=\frac{\partial}{\partial\lambda_{i}}\Big{(}-\log(y_{i}!)+\lambda_{i}y_{i}+\frac{1}{2}\rho_{i}^{2}y_{i}^{2}+\log\widetilde{\mathcal{L}}(s_{i})\Big{)}\\\
&=\frac{\partial}{\partial\lambda_{i}}\Big{(}\lambda_{i}y_{i}+\frac{1}{2}\rho_{i}^{2}y_{i}^{2}-\dfrac{1}{\rho_{i}^{2}}\mathcal{W}(s_{i}\rho_{i}^{2})-\dfrac{1}{2\rho_{i}^{2}}\mathcal{W}(s_{i}\rho_{i}^{2})^{2}-\dfrac{1}{2}\log\big{(}1+\mathcal{W}(s_{i}\rho_{i}^{2})\big{)}\Big{)}\\\
&=y_{i}-\Big{(}\dfrac{1}{\rho_{i}^{2}}\big{(}1+\mathcal{W}(s_{i}\rho_{i}^{2})\big{)}+\dfrac{1}{2}\dfrac{1}{1+\mathcal{W}(s_{i}\rho_{i}^{2})}\Big{)}\frac{\partial}{\partial\lambda_{i}}\mathcal{W}(s_{i}\rho_{i}^{2})\\\
&=y_{i}-\Big{(}\dfrac{1}{\rho_{i}^{2}}\big{(}1+\mathcal{W}(s_{i}\rho_{i}^{2})\big{)}+\dfrac{1}{2}\dfrac{1}{1+\mathcal{W}(s_{i}\rho_{i}^{2})}\Big{)}\frac{\mathcal{W}(s_{i}\rho_{i}^{2})}{1+\mathcal{W}(s_{i}\rho_{i}^{2})}\frac{1}{s_{i}\rho_{i}^{2}}\frac{\partial}{\partial\lambda_{i}}(s_{i}\rho_{i}^{2})\\\
&=y_{i}-\mathcal{W}(s_{i}\rho_{i}^{2})\Big{(}\dfrac{1}{\rho_{i}^{2}}+\frac{1}{2}\big{(}1+\mathcal{W}(s_{i}\rho_{i}^{2})\big{)}^{-2}\Big{)}\;,\end{split}$

and similarly

$\begin{split}\frac{\partial}{\partial\rho_{i}^{2}}\log\widetilde{Z}_{h_{i}}&=\frac{\partial}{\partial\rho_{i}^{2}}\Big{(}\frac{1}{2}\rho_{i}^{2}y_{i}^{2}-\dfrac{1}{\rho_{i}^{2}}\mathcal{W}(s_{i}\rho_{i}^{2})-\dfrac{1}{2\rho_{i}^{2}}\mathcal{W}(s_{i}\rho_{i}^{2})^{2}-\dfrac{1}{2}\log\big{(}1+\mathcal{W}(s_{i}\rho_{i}^{2})\big{)}\Big{)}\\\