Asymptotically exact data augmentation: models, properties and algorithms

TL;DR

This paper introduces AXDA, a unified framework for data augmentation that maintains the original distribution asymptotically, offering statistical benefits and efficient inference algorithms, with theoretical analysis and practical illustrations.

Contribution

It proposes a systematic approach to asymptotically exact data augmentation, unifying existing models and analyzing their statistical properties and computational efficiency.

Findings

AXDA models preserve the original distribution asymptotically.

Theoretical bounds on approximation quality are established.

AXDA enables efficient inference algorithms with proven properties.

Abstract

Data augmentation, by the introduction of auxiliary variables, has become an ubiquitous technique to improve convergence properties, simplify the implementation or reduce the computational time of inference methods such as Markov chain Monte Carlo ones. Nonetheless, introducing appropriate auxiliary variables while preserving the initial target probability distribution and offering a computationally efficient inference cannot be conducted in a systematic way. To deal with such issues, this paper studies a unified framework, coined asymptotically exact data augmentation (AXDA), which encompasses both well-established and more recent approximate augmented models. In a broader perspective, this paper shows that AXDA models can benefit from interesting statistical properties and yield efficient inference algorithms. In non-asymptotic settings, the quality of the proposed approximation is assessed with several theoretical results. The latter are illustrated on standard statistical problems. Supplementary materials including computer code for this paper are available online.

Tables

Table 1. Table 1: Closed-form expressions of m 2 subscript 𝑚 2 m_{2} appearing in ( 20 ) for multivariate generalizations of standard kernels where d 𝑑 d denotes the dimension.

	Gaussian	Cauchy	Laplace	Dirichlet	Uniform	Triangular	Epanechnikov
$m_2$	$\sqrt{d}$	-	$\sqrt{2d}$	-	$\sqrt{d/3}$	$\sqrt{d/6}$	$\sqrt{d/5}$

Table 2. Table 2: Illustration of the bound derived in ( 40 ) for the marginal posterior π ρ subscript 𝜋 𝜌 \pi_{\rho} depicted in Section 5.2 . The (1- α 𝛼 \alpha )-credibility intervals 𝒞 α subscript 𝒞 𝛼 \mathcal{C}_{\alpha} and 𝒞 α ρ superscript subscript 𝒞 𝛼 𝜌 \mathcal{C}_{\alpha}^{\rho} are the highest posterior density regions associated to each density with α = 0.05 𝛼 0.05 \alpha=0.05 .

$\rho$	${𝒞}_{\alpha }$	${𝒞}_{\alpha }^{\rho }$	${\int }_{{𝒞}_{\alpha }^{\rho }}\pi \left({\theta }_1\right)d{\theta }_1$	${ℐ}_{\alpha }^{\rho }$
${10}^{-3}$	[-0.47,1.24]	[-0.47,1.24]	0.95	[0.949,0.951]
${10}^{-2}$	idem	[-0.47,1.24]	0.95	[0.948,0.952]
${10}^{-1}$	idem	[-0.47,1.24]	0.95	[0.88,1]
${10}^0$	idem	[-0.47,1.37]	0.96	[0.34,1]

Table 3. Table 3: Examples of classical kernels K 𝐾 K that can be used to define an appropriate density κ ρ subscript 𝜅 𝜌 \kappa_{\rho} verifying Property 1 in the main paper.

name	support	$K\left(u\right)$
Gaussian	$ℝ$	$\frac{1}{\sqrt{2\pi }}\exp \left(-u^2/2\right)$
Cauchy	$ℝ$	$\frac{1}{\pi \left(1+u^2\right)}$
Laplace	$ℝ$	$\frac{1}{2}\exp \left(-\left|u\right|\right)$
Dirichlet	$ℝ$	$\frac{{\sin }^2\left(u\right)}{\pi u^2}$
Uniform	$[-1,1]$	$\frac{1}{2}{\mathrm{𝟙}}_{\left|u\right|\le 1}$
Triangular	$[-1,1]$	$\left(1-\left|u\right|\right){\mathrm{𝟙}}_{\left|u\right|\le 1}$
Epanechnikov	$[-1,1]$	$\frac{3}{4}\left(1-u^2\right){\mathrm{𝟙}}_{\left|u\right|\le 1}$

Table 4. Table 4: Examples of divergence functions ϕ italic-ϕ \phi that can be used to define an appropriate density κ ρ subscript 𝜅 𝜌 \kappa_{\rho} verifying Property 1 in the main paper.

name	$\Theta$	$\varphi (z,\theta )$
Squared loss	$ℝ$	${\left(z-\theta \right)}^2$
Absolute loss	$ℝ$	$\left|z-\theta \right|$
Logistic loss	$[0,1]$	$z\log \left(\frac{z}{\theta }\right)+\left(1-z\right)\log \left(\frac{1-z}{1-\theta }\right)$
Itakura-Saito divergence	${ℝ}_{+}$	$\frac{z}{\theta }-\log \left(\frac{z}{\theta }\right)-1$
Kullback-Leibler divergence	$[0,1]$	$z\log \left(\frac{z}{\theta }\right)$

Table 5. Table 5: Lipschitz loss functions f 𝑓 f used in standard statistical learning problems. Their domain of definition is denoted 𝒟 f subscript 𝒟 𝑓 \mathcal{D}_{f} and y 𝑦 y stands for an observation. The notation “reg.” stands for regression.

name	problem	${𝒟}_f$	$f(y;t)$
hinge	SVM	$\{-1,1\}\times ℝ$	$\max (0,1-yt)$
Huber	robust reg.	$ℝ\times ℝ$	$\{\begin{array}{c}{\left(y-t\right)}^2/\left(2\delta \right)\text{if}\left|y-t\right|\le \delta \\ \left|y-t\right|-\delta /2\text{otherwise, where}\delta >0\end{array}$
logistic	logistic reg.	$\{-1,1\}\times ℝ$	$\log \left(1+\exp \left(-yt\right)\right)$
pinball	quantile reg.	$ℝ\times ℝ$	$\tau \max (0,t-y)+\left(1-\tau \right)\max (0,y-t)$, $\tau \in (0,1)$

Equations

$$\pi (\mathbf{\theta })\propto \exp \left(-f(\mathbf{\theta })\right),\text{or}\pi (\mathbf{y}|\mathbf{\theta })\propto \exp \left(-f(y;\mathbf{\theta })\right),$$

$${\int }_{\mathcal{Z}}\pi (\mathbf{\theta },\mathbf{z})\mathrm{d}\mathbf{z}=\pi (\mathbf{\theta }).$$

$${\pi }_{\rho }(\mathbf{\theta })={\int }_{\mathcal{Z}}{\pi }_{\rho }(\mathbf{\theta },\mathbf{z})\mathrm{d}\mathbf{z},$$

$${\pi }_{\rho }(\mathbf{\theta },\mathbf{z})=\pi (\mathbf{z}){\kappa }_{\rho }(\mathbf{z},\mathbf{\theta }),$$

$$d_{\psi }(\mathbf{z},\mathbf{\theta })=\psi (\mathbf{z})-\psi (\mathbf{\theta })-\nabla \psi (\mathbf{\theta }{)}^T(\mathbf{z}-\mathbf{\theta }).$$

$$y_i|\mathbf{\theta }\sim p(y_i|g^{-1}(h({\mathbf{x}}_i,\mathbf{\theta })),{\sigma }^2),\forall i\in [n],$$

$$\mathbf{\theta }\sim \mathcal{N}(\mathbf{\theta }|0_d,{\nu }^2{\mathbf{I}}_d),$$

$$h({\mathbf{x}}_i,\mathbf{\theta })=\sum _{j=1}^k{\theta }_{j}B({\mathbf{x}}_i,{\mathbf{k}}_j),$$

$$y_i|z_i\sim p(y_i|g^{-1}(z_i),{\sigma }^2),\forall i\in [n],$$

$$z_i|\mathbf{\theta }\sim \mathcal{N}(z_i|h({\mathbf{x}}_i,\mathbf{\theta }),{\rho }^2),\forall i\in [n],$$

$$\mathbf{\theta }\sim \mathcal{N}(\mathbf{\theta }|0_d,{\nu }^2{\mathbf{I}}_d).$$

$$\pi (\mathbf{\theta }|\mathbf{y},\mathbf{X})\propto \nu (\mathbf{\theta })\prod _{b=1}^B\prod _{i\in \text{node }b}\exp \left(-f_i(y_i;h({\mathbf{x}}_i,\mathbf{\theta }))\right).$$

$${\pi }_{\rho }(\mathbf{\theta },\mathbf{z}|\mathbf{y},\mathbf{X})\propto \nu (\mathbf{\theta })\prod _{b=1}^B\prod _{i\in \text{node }b}\exp \left(-f_i(y_i;z_i)\right){\kappa }_{\rho }(z_i,h({\mathbf{x}}_i,\mathbf{\theta })).$$

$$\pi (\mathbf{y}|\mathbf{\theta })\propto \prod _{i=1}^n\pi (y_i|\mathbf{\theta }),$$

$${\pi }_{\rho }(\mathbf{y},{\mathbf{z}}_{1\text{mathchar}58n}|\mathbf{\theta })\propto \prod _{i=1}^n\pi (y_i|{\mathbf{z}}_i){\kappa }_{\rho }({\mathbf{z}}_i,\mathbf{\theta }).$$

$${\pi }_{\rho }(\mathbf{y}|\mathbf{\theta })\triangleq {\int }_{\Theta }{\pi }_{\rho }(\mathbf{y},\mathbf{z}|\mathbf{\theta })\mathrm{d}\mathbf{z}={\int }_{\Theta }\pi (\mathbf{z}|\mathbf{\theta }){\kappa }_{\rho }(\mathbf{z},\mathbf{y})\mathrm{d}\mathbf{z}.$$

$${\mathbb{E}}_{{\pi }_{\rho }}(\theta )$$

$${\mathrm{var}}_{{\pi }_{\rho }}(\theta )$$

$$W_p(\pi ,{\pi }_{\rho })={\left(\underset{\mu }{\min }\left\{{\int }_{{\mathbb{R}}^d}{\int }_{{\mathbb{R}}^d}{\Vert \mathbf{\theta }-\mathbf{z}\Vert }_2^p\mathrm{d}\mu (\mathbf{z},\mathbf{\theta });\mu \in \Gamma ({\pi }_{\rho },\pi )\right\}\right)}^{1/p},$$

$$W_p(\pi ,{\pi }_{\rho })\le \rho m_p.$$

$${\kappa }_{\rho }(\mathbf{z},\mathbf{\theta })\propto \exp \left(-\frac{d_{\psi }(\mathbf{z},\mathbf{\theta })}{\rho }\right).$$

$${\pi }_{\rho }(\mathbf{\theta })-\pi (\mathbf{\theta })=\mathcal{O}(\sqrt{\rho }).$$

$${\pi }_{\rho }(\mathbf{\theta })-\pi (\mathbf{\theta })=\frac{\rho }{2}\mathrm{Trace}\left({\mathbf{H}}_{\pi }(\mathbf{\theta }){\mathbf{H}}_{d_{\psi }}(\mathbf{\theta }{)}^{-1}\right)+o(\rho ).$$

$${\kappa }_{\rho }(\mathbf{z},\mathbf{\theta })=\mathcal{N}(\mathbf{z}|\mathbf{\theta },{\rho }^2{\mathbf{I}}_d).$$

$${\Vert {\pi }_{\rho }-\pi \Vert }_{\mathrm{TV}}\le 1-{\Delta }_d(\rho ),$$

$${\Delta }_d(\rho )=\frac{D_{-d}(L_f\rho )}{D_{-d}(-L_f\rho )}.$$

$$D_{-d}(z)=\frac{\exp (-z^2/4)}{\Gamma (d)}{\int }_0^{+\infty }{e}^{-xz-x^2/2}{x}^{d-1}\mathrm{d}x.$$

$${\Vert {\pi }_{\rho }-\pi \Vert }_{\mathrm{TV}}\le \rho L_f\frac{2\sqrt{2}\Gamma \left(\frac{d+1}{2}\right)}{\Gamma \left(\frac{d}{2}\right)}+o(\rho ),$$

$${\Vert {\pi }_{\rho }-\pi \Vert }_{\mathrm{TV}}\le 1-\frac{1}{(1+2{\rho }^2{M}_f{)}^{d/2}}\left(1-\frac{{\rho }^4{M}_f{\mathsf{M}}_f}{1+2{\rho }^2{M}_f}\right).$$

$${\Vert {\pi }_{\rho }-\pi \Vert }_{\mathrm{TV}}\le {\rho }^{2}d{M}_f+o({\rho }^2).$$

$$\pi (\mathbf{\theta })\propto \prod _{j=1}^J{\pi }_j(\mathbf{\theta })\propto \exp \left(-\sum _{j=1}^J{f}_j(\mathbf{\theta })\right),$$

$${\pi }_{\rho }(\mathbf{\theta },{\mathbf{z}}_{1\text{mathchar}58J})\propto \prod _{j=1}^J{\pi }_j({\mathbf{z}}_j){\kappa }_{{\rho }_j}({\mathbf{z}}_j;\mathbf{\theta })\propto \exp \left(-\sum _{j=1}^J{f}_j({\mathbf{z}}_j)+\frac{1}{2{\rho }_j^2}{\Vert {\mathbf{z}}_j-\mathbf{\theta }\Vert }_2^2\right).$$

Full text

Keywords: Approximation, auxiliary variables, divide-and-conquer, Bayesian inference, robustness.

1 Introduction

Starting at least from the 1960s with the seminal paper of Hartley (1958) on the expectation-maximization (EM) algorithm, introducing auxiliary variables has been a widely adopted strategy to derive iterative algorithms able to deal with possibly complicated inference problems. Indeed, either by coming from statistical physics (Swendsen and Wang 1987) or by the broad statistical community (Dempster et al. 1977), auxiliary (also called latent) variables have been used to improve (Duane et al. 1987; Edwards and Sokal 1988; Marnissi et al. 2018) and/or simplify (Tanner and Wong 1987; Doucet et al. 2002) inference methods, such as maximum likelihood (ML) estimation or simulation-based ones. Insightful reviews of these methods were conducted by Besag and Green (1993); van Dyk and Meng (2001); Tanner and Wong (2010). Among many others, slice sampling and half-quadratic (HQ) methods are archetypal instances of such auxiliary variable-based methods. These methods, by introducing auxiliary variables, appear to be an interesting alternative when sampling cannot be performed directly from a target distribution $\pi$. Nonetheless, the superiority of simulation-based algorithms based on data augmentation (DA) over classical Markov chain Monte Carlo (MCMC) methods without DA is not obvious as pointed out by Polson (1996); Damien et al. (1999). DA methods have been found to be slower than single-site update approaches in some cases (Hurn 1997) and some improvements have been derived to cope with these problems such as partial decoupling (Higdon 1998) or the introduction of a working parameter (Meng and van Dyk 1997). Moreover, DA techniques are often used on a case-by-case basis (Geman and Reynolds 1992; Albert and Chib 1993; Geman and Yang 1995; Polson et al. 2013) and could not be applied in general scenarios due to the absence of exact DA schemes yielding an efficient inference and low computation costs.

Similarly to approximate Bayesian computation (ABC) methods to circumvent intractable likelihoods (Beaumont et al. 2002; Sisson et al. 2018b), these limitations can be tackled by considering approximate DA schemes that become exact asymptotically. For instance, inspired from the variable splitting technique used in the alternating direction method of multipliers (ADMM) (Boyd et al. 2011), Vono et al. (2019) and Rendell et al. (2020) recently and independently proposed a novel and broad Bayesian inference framework that can circumvent limitations of exact DA approaches. By introducing a collection of instrumental (also called splitting) variables, the aforementioned authors considered the inference from an approximate probability distribution which can be simpler, more efficient and distributed over multiple computational workers (e.g., machines or kernels).

This paper aims at deeply investigating a broad framework coined asymptotically exact data augmentation (AXDA) which encompasses previously proposed special instances such as approximate models used in Vono et al. (2019); Rendell et al. (2020), among others. More precisely, Section 2 details how such models can be built in a quasi-systematic and simple way which is highly appreciable compared to the case-by-case search of computationally efficient DA schemes. In Section 3, we revisit some already-proposed special instances of AXDA models in order to show the potential benefits of AXDA on specific examples and to exhibit interesting properties which can be generally inherited by AXDA approaches. In Section 4, we assess quantitatively the bias of AXDA models with non-asymptotic theoretical results by considering Wasserstein and total variation distances. Then, Section 5 illustrates the previous theoretical results and the benefits of the proposed methodology on several statistical problems. In order to facilitate the use of AXDA, we eventually point out that the supplementary material involves a dedicated section (Section 5) presenting how such models can be instantiated to perform efficient inference through classical simulation-based, variational Bayes (VB), optimization or expectation-maximization (EM) methods. The proofs are also given in the supplementary material, see Section 1.

2 Asymptotically exact data augmentation

This section introduces AXDA schemes that aim to circumvent exact DA main issue: the art (van Dyk and Meng 2001) of finding the exact DA associated to a statistical model and its inference limitations. For sake of simplicity, with little abuse, we shall use the same notations for a probability distribution and its associated probability density function (pdf).

2.1 Motivations

In this paper, we are interested in performing the inference of a variable of interest $\boldsymbol{\theta}\in\Theta\subseteq\mathbb{R}^{d}$, where $\Theta$ is a closed convex set and $\mathrm{dim}(\Theta)=d$, by relying on a probability distribution with density $\pi$ writing

[TABLE]

where the potential $f$ taking values in $\mathbb{R}$ is such that $\pi$ defines a proper, bounded and continuous probability distribution. For sake of generality, note that $\pi$ in (1) shall describe various quantities. First, with a little abuse of notations, $\pi(\boldsymbol{\theta})$ may simply refer to a pdf associated to the random variable $\boldsymbol{\theta}$, e.g., its prior distribution $\pi(\boldsymbol{\theta})$ or its posterior distribution $\pi(\boldsymbol{\theta})\triangleq\pi(\boldsymbol{\theta}|\mathbf{y})$ when referring to a set of observations denoted by $\mathbf{y}$. Depending on the problem, we also allow $\pi$ to stand for a likelihood function $\pi(\mathbf{y}|\boldsymbol{\theta})$. We will work under this convention and write explicitly the form of $\pi$ when required. For sake of simplicity and clarity, only the case corresponding to $\pi(\boldsymbol{\theta})$ will be detailed in this section. The application of the proposed methodology to $\pi(\mathbf{y}|\boldsymbol{\theta})$ is very similar and can be retrieved by a straightforward derivation.

We consider situations where direct inference from (1) is difficult because intractable or computationally prohibitive. To overcome these issues, an option is to rely on exact DA which introduces some auxiliary variables stacked into a vector $\mathbf{z}\in\mathcal{Z}\subseteq\mathbb{R}^{k}$ and defines a new density, simpler to handle, such that

[TABLE]

Much research has been devoted to these models in order to simplify an inference task or to improve the convergence properties of direct inference approaches (e.g., slice sampling and HQ methods introduced in Section 1). Nonetheless, these approaches have several limitations. Indeed, finding a convenient form for the augmented density in order to satisfy (2) while leading to efficient algorithms generally requires some knowledge and can even be impossible in some cases (Geman and Yang 1995). For instance, the mixture representation of a binomial likelihood function based on the Polya-Gamma distribution has been used to derive a promising Gibbs sampler for logistic regression problems (Polson et al. 2013). Nonetheless, even if this algorithm has been proved to be uniformly ergodic by Choi and Hobert (2013), the corresponding ergodicity constant depends exponentially on the number of observations $n$ and on the dimension of the regression coefficients vector $d$.

To tackle these limitations, we propose to relax the constraint (2) and consider an approximate DA model. This will permit the choice of an augmented density with more flexibility, fix the issues associated to the initial model and make inference more efficient in some cases. To this purpose, Section 2.2 presents the so-called AXDA framework which embeds approximate DA models controlled by a positive scalar parameter $\rho$. These models become asymptotically exact when $\rho$ tends towards 0. Of course, some assumptions will be required on the approximate augmented density to guarantee a good approximation. The quality of this approximation will be assessed in Section 4 with non-asymptotic theoretical results.

2.2 Model

Instead of searching for an exact data augmentation scheme (2), some auxiliary variables $\mathbf{z}$ can be introduced in order to define an approximate but asymptotically exact probability distribution. One possibility is to introduce an augmented distribution depending on a parameter $\rho>0$ and such that the associated marginal density defined by

[TABLE]

satisfies the following property.

Property 1**.**

For all $\boldsymbol{\theta}\in\Theta$, $\lim_{\rho\rightarrow 0}\pi_{\rho}(\boldsymbol{\theta})=\pi(\boldsymbol{\theta})$ .

By applying Scheffé’s lemma (Scheffé 1947), this property yields the convergence in total variation, that is $\left\|\pi_{\rho}-\pi\right\|_{\mathrm{TV}}\rightarrow 0$ as $\rho\rightarrow 0$. A natural question is: how to choose the augmented density in (3) such that Property 1 is met? In this paper, we assume that $\mathcal{Z}=\Theta$ and investigate AXDA schemes associated to an initial density (1) and defined by the approximate augmented density

[TABLE]

where $\kappa_{\rho}$ is such that (4) defines a proper density.

Remark 1**.**

When $\pi$ stands for a product of $J$ densities, that is $\pi\propto\prod_{j=1}^{J}\pi_{j}$, the proposed approximate model can naturally be generalized to $\pi_{\rho}(\boldsymbol{\theta},\mathbf{z}_{1\mathrel{\mathop{\mathchar 58\relax}}J})\propto\prod_{j=1}^{J}\pi_{j}(\mathbf{z}_{j})\kappa_{\rho}(\mathbf{z}_{j},\boldsymbol{\theta})$. Such a generalization will for instance be considered in Sections 3.1 and 3.2.

The introduction of the proposed model (4) is aimed at avoiding a case-by-case search of an appropriate augmented approach. Although there might exist other marginal densities $\pi_{\rho}$ satisfying Property 1, we restrict our analysis to models where $\kappa_{\rho}(\cdot,\boldsymbol{\theta})$ weakly converges towards the Dirac measure at $\boldsymbol{\theta}$ as $\rho\rightarrow 0$ (Aguirregabiria et al. 2002). This is a sufficient condition to satisfy Property 1. In the sequel, we will call AXDA any approach based on (4) and satisfying these properties.

A natural choice for $\kappa_{\rho}$ is to consider a standard kernel $K$ (Wand and Jones 1995). Based on the latter, we define for all $\mathbf{z},\boldsymbol{\theta}\in\Theta$, $\kappa_{\rho}(\mathbf{z},\boldsymbol{\theta})\propto_{\mathbf{z}}\rho^{-d}K(\rho^{-1}(\boldsymbol{\theta}-\mathbf{z}))$ (Dang and Ehrhardt 2012). Beyond standard kernels but motivated by the same idea of measuring the discrepancy between $\mathbf{z}$ and $\boldsymbol{\theta}$, one can also build on divergence functions widely used in the optimization literature to define a potentially asymmetric density $\kappa_{\rho}$ such that for all $\mathbf{z},\boldsymbol{\theta}\in\Theta$, $\kappa_{\rho}(\mathbf{z},\boldsymbol{\theta})\propto_{\mathbf{z}}\exp(-\rho^{-1}\phi(\mathbf{z},\boldsymbol{\theta}))$ where $\phi$ is a strictly convex function w.r.t. $\mathbf{z}$ admitting a unique minimizer $\mathbf{z}^{*}=\boldsymbol{\theta}$ (Ben-Tal et al. 2001; Krichene et al. 2015; Fellows et al. 2019). Specific instances of such potentials are Bregman divergences such as the logistic loss and the Kullback-Leibler divergence, see Definition 1.

Definition 1** (Bregman divergence).**

Let $\psi$ a continuously-differentiable and strictly convex function defined on a closed convex set. The Bregman divergence associated to $\psi$ is defined by

[TABLE]

Additional details associated to standard kernels and Bregman divergences are given in Section 2 in the supplementary material.

3 Benefits of AXDA by revisiting existing models

Before providing theoretical guarantees for AXDA models, this section proposes to review some important state-of-the-art works from the AXDA perspective described in Section 2. We do not pretend to give new insights about these approaches. We rather use them to illustrate potential benefits that can be gained by resorting to the proposed framework. For sake of clarity, these benefits are directly highlighted in the title of the following sections before being discussed in the latter.

3.1 Tractable posterior inference

This first section illustrates how an AXDA approach can alleviate the intractability of an initial posterior distribution $\pi$ and significantly aid in the computations.

To this purpose, we consider the case where the posterior distribution $\pi$ is intractable. Such a model for instance appears when $\pi$ involves a constraint on some set (Liechty et al. 2009), admits a non-standard potential function such as the total variation norm (Chambolle et al. 2010; Pereyra 2016; Vono et al. 2019) or yields complicated conditional posterior distributions (Holmes and Mallick 2003). To simplify the inference, the aforementioned authors have considered special instances of AXDA by relying on an additional level involving latent variables $\mathbf{z}$, leading a hierarchical Bayesian model. In these cases, AXDA has been invoked in order to move a difficulty to the conditional posterior of $\mathbf{z}$ where it can be dealt with more easily by using standard inference algorithms, see Section 5 in the supplementary material for more details. The following example, derived from Holmes and Mallick (2003), illustrates this idea.

Example 1**.**

Let $\mathbf{y}\in\mathbb{R}^{n}$ be a set of observations and $\mathbf{X}=(\mathbf{x}_{1},\ldots,\mathbf{x}_{n})^{T}\in\mathbb{R}^{n\times d}$ a design matrix filled with covariates. We consider a generalized non-linear model which writes

[TABLE]

where $p$ belongs to the exponential family and has mean $g^{-1}(h(\mathbf{x}_{i},\boldsymbol{\theta}))$ and variance $\sigma^{2}$ where $g$ is a link function. As in classical regression problems, we are interested in infering the regression coefficients $\boldsymbol{\theta}\in\mathbb{R}^{d}$. In the sequel, we set the non-parametric model $h$ to be

[TABLE]

where $B(\mathbf{x}_{i},\mathbf{k}_{j})$ is a non-linear function of $\mathbf{x}_{i}$ (e.g., regression splines) and $\mathbf{k}_{j}$ is the knot location of the $j$-th basis. The difficulty here is the non-linearity of $h$ which, combined with the non-Gaussian likelihood, rules out the use of efficient simulation schemes to sample from the posterior $\pi(\boldsymbol{\theta}|\mathbf{y})$. In order to mitigate this issue, Holmes and Mallick (2003) proposed to rely on an additional level which boils down to consider the approximate model (4). More specifically, the aforementioned authors treated the non-linear predictor $h$ as a Gaussian random latent variable which leads to the approximate model

[TABLE]

Here, AXDA has been applied only to the likelihood function with $\kappa_{\rho}$ chosen as the univariate normal distribution (10) leading to a smoothed likelihood function. The main advantage of relying on such a model is that the posterior conditional distribution $\pi_{\rho}(\boldsymbol{\theta}|\mathbf{z},\mathbf{X})$, with $\mathbf{z}=[z_{1},\ldots,z_{n}]^{T}$, is now a multivariate normal distribution. In addition, by moving the difficulty induced by $h$ to the conditional posterior of $z_{i}$, we are now dealing with a generalized linear model where standard techniques can be applied (Albert and Chib 1993; Polson et al. 2013).

Beyond the widely-used Gaussian choice for $\kappa_{\rho}$ (Holmes and Mallick 2003; Liechty et al. 2009; Barbos et al. 2017; Vono et al. 2019), more general AXDA approaches can be built by taking inspiration from these works. To this purpose, we recommand to choose $\kappa_{\rho}$ w.r.t. the prior and likelihood at stake. For instance, when a Poisson likelihood function and a complex prior distribution on its intensity $\boldsymbol{\theta}$ are considered, one option for $\phi$ (see Section 2.2) would be an Itakura-Saito divergence since it preserves the positivity constraint on $\boldsymbol{\theta}$ and yields the well-known Gamma-Poisson model (Canny 2004).

3.2 Distributed inference

When data are stored on multiple machines and/or one is interested in respecting their privacy, this section illustrates how AXDA can be resorted to perform distributed computations.

Let consider observed data $\{y_{i},\mathbf{x}_{i}\}_{i=1}^{n}$, where $\mathbf{x}_{i}$ stands for the covariates associated to observation $y_{i}$, which are distributed among $B$ nodes within a cluster. By adopting a prior $\nu(\boldsymbol{\theta})$ and by assuming that the likelihood can be factorized w.r.t. the $B$ nodes, the posterior distribution of the variable of interest $\boldsymbol{\theta}$ writes

[TABLE]

Such models classically appear in statistical machine learning when generalized linear models (GLMs) (Dobson and Barnett 2008) are considered. In these cases, $h(\mathbf{x}_{i},\boldsymbol{\theta})=\mathbf{x}_{i}^{T}\boldsymbol{\theta}$. Due to the distributed environment, sampling efficiently from (12) is challenging and a lot of “divide-and-conquer” approaches have been proposed in the past few years to cope with this issue (Wang and Dunson 2013; Scott et al. 2016). These methods launch independent Markov chains on each node $b$ and then combine the outputs of these local chains to obtain an approximation of the posterior of interest (12). Nonetheless, the averaging schemes used to combine the local chains might lead to poor approximations when $\pi$ is high-dimensional and non-Gaussian. Instead, considering a special instance of AXDA circumvents the previously mentioned drawbacks by introducing local auxiliary variables on each node such that

[TABLE]

The posterior distribution of the auxiliary variables conditionally to $\boldsymbol{\theta}$ only depends on the data available at a given node. Based on this nice property, the joint posterior can be sampled efficiently with a Gibbs sampler, see Rendell et al. (2020) for a comprehensive review. We emphasize that the benefits described in this section for Monte Carlo sampling also hold when one wants to use other types of algorithms (e.g., expectation-maximization or variational Bayes), see Section 5 in the supplementary material.

3.3 Robust inference

By noting that classical robust hierarchical models fall into the proposed framework, this section shows that AXDA is also a relevant strategy to perform robust inference by coping with model misspecification by modeling additional sources of uncertainty.

Considering a well-chosen demarginalization procedure is known to yield robustness properties in some cases (Robert and Casella 2004). Some approaches took advantage of this idea in order to build robust hierarchical Bayesian models w.r.t. possible outliers in the data. For instance, such models can be built by allowing each observation to be randomly drawn from a local statistical model, as described in the recent review of Wang and Blei (2018). This “localization” idea is illustrated in Figure 1.

Many of these models can be viewed as particular instances of AXDA. Indeed, assume that $n$ data points $y_{i}$ are independently and identically distributed (i.i.d.) defining the likelihood function

[TABLE]

where $\boldsymbol{\theta}\in\Theta$ is a common parameter. Applying AXDA as described in Section 2 by introducing $n$ $d$-dimensional auxiliary variables stacked into the vector $\mathbf{z}_{1\mathrel{\mathop{\mathchar 58\relax}}n}$ leads to the augmented likelihood

[TABLE]

The statistical model defined by (15) implies a hierarchical Bayesian model similar to the localized one depicted on Figure 1(b) and corresponds in general to an approximation of the initial one, see Example 2.

Example 2**.**

Robust logistic regression – Assume that for all $i\in[n]$, $\pi(y_{i}|\boldsymbol{\theta})=\mathcal{B}\left(\sigma\left(\mathbf{x}_{i}^{T}\boldsymbol{\theta}\right)\right)$, where $\mathcal{B}$ stands for the Bernoulli distribution, $\sigma$ for the sigmoid function, $\mathbf{x}=\left[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\right]$ for the transpose of the design matrix and $\boldsymbol{\theta}$ for the regression coefficients vector to infer. Then as proposed by Wang and Blei (2018), one can robustify the inference by assuming that each observation $y_{i}$ is drawn from a local and independent model $\mathcal{B}\left(\sigma\left(\mathbf{x}_{i}^{T}\mathbf{z}_{i}\right)\right)$ associated to an auxiliary parameter $\mathbf{z}_{i}\sim\mathcal{N}(\boldsymbol{\theta},\rho^{2}\mathbf{I}_{d})$. In this case, $\kappa_{\rho}(\mathbf{z},\boldsymbol{\theta})\propto\prod_{i=1}^{n}\mathcal{N}(\mathbf{z}_{i}\ |\ \boldsymbol{\theta},\rho^{2}\mathbf{I}_{d})$.

Beyond the convenient Gaussian prior $\kappa_{\rho}$, its choice can be motivated by robust loss functions. In the statistical machine learning literature, the absolute or Huber losses are of common use (She and Owen 2011). In Bayesian linear inverse problems considered in the signal processing community, it is classical to approximate a complicated forward physical model in order to yield tractable computations. If the latter can be written as $\mathbf{y}=h(\boldsymbol{\theta})+\boldsymbol{\epsilon}$, with $\boldsymbol{\epsilon}\sim\pi(\boldsymbol{\epsilon})$, then introducing a latent variable $\mathbf{z}\sim\kappa_{\rho}(\mathbf{z},h(\boldsymbol{\theta}))$ such that $\mathbf{y}=\mathbf{z}+\boldsymbol{\epsilon}$ allows to take into consideration the model approximation. In those cases, one can set $\kappa_{\rho}$ to be the distribution of the modeling error which could be adjusted thanks to some expertise.

3.4 Inheriting sophisticated inference schemes from ABC

Finally, this section shows that AXDA models, by sharing strong connections with ABC, might inherit sophisticated algorithms to sample from (4).

ABC stands for a family of methods that permit to cope with intractable likelihoods by sampling from the latter instead of evaluating them. In a nutshell, if one’s goal is to infer a parameter $\boldsymbol{\theta}$ based on a posterior of interest, the simplest ABC rejection sampler is as follows. At iteration $t$, draw a candidate $\boldsymbol{\theta}^{(t)}$ from the prior, generate pseudo-observations $\mathbf{z}$ from the likelihood given this candidate and accept $\boldsymbol{\theta}^{(t)}$ if $\mathbf{z}=\mathbf{y}$ where $\mathbf{y}$ is the observations vector. Many more sophisticated ABC samplers have been derived. We refer the interested reader to the recent review by Sisson et al. (2018a) for more information about ABC methods.

Among a huge literature on ABC (also called likelihood-free) methods, noisy ABC approaches proposed and motivated by Fearnhead and Prangle (2012) and Wilkinson (2013) are strongly related to AXDA. Indeed, only comparing the underlying models, AXDA with observation splitting is equivalent to noisy ABC. To see this, let $\pi(\mathbf{y}|\boldsymbol{\theta})$ stand for an intractable likelihood. Noisy ABC replaces the exact inference based on $\pi$ by considering the pseudo-likelihood with density

[TABLE]

This density has exactly the same formulation as the one defined in (4) except that noisy ABC splits the observations $\mathbf{y}$ instead of the parameter of interest $\boldsymbol{\theta}$. Capitalizing on this equivalence property, also pointed out by Rendell et al. (2020), one can derive efficient algorithms for AXDA from the ABC framework. For instance, Rendell et al. (2020) recently built on the works of Beaumont et al. (2002); Del Moral et al. (2012) in the ABC context to propose a bias correction approach and a sequential Monte Carlo (SMC) algorithm avoiding the tuning of the tolerance parameter $\rho$. Obviously, many other inspirations from ABC can be considered, such as the parallel tempering approach of Baragatti et al. (2013) among others, to make the inference from an AXDA model more flexible and efficient.

4 Theoretical guarantees

By building on existing approaches, Section 3 showed that AXDA can be used in quite general and different settings depending on ones motivations. In order to further promote the use of such approximate augmented models, this section goes beyond the empirical bias analysis performed by previous works and provides quantitative bounds on the error between the initial and the approximate model. More precisely, for a fixed tolerance parameter $\rho>0$, non-asymptotic results on the error associated to densities, potentials and credibility regions are derived. We will assume all along this section that $\Theta=\mathbb{R}^{d}$. The proofs of the results of this section can be found in Section 1 of the supplementary material.

4.1 Results for standard kernels

In this section, we consider the case $\kappa_{\rho}(\mathbf{z},\boldsymbol{\theta})\propto\rho^{-d}K(\rho^{-1}(\boldsymbol{\theta}-\mathbf{z}))$ where $K$ is a kernel, see Section 2.2. Under this model, $\pi_{\rho}$ stands for the convolution of $\pi$ and $\kappa_{\rho}$ and the following results hold.

Proposition 1**.**

Let $\pi\in L^{1}$. The marginal with density $\pi_{\rho}$ in (3) has the following properties.

1. i)

Let $\pi$ stand for a pdf associated to the random variable $\boldsymbol{\theta}$ and $\mathbb{E}_{\kappa_{\rho}}(X)=0$. Then, the expectation and variance under $\pi_{\rho}$ are given by

[TABLE] 2. ii)

$\mathrm{supp}(\pi_{\rho})\subseteq S$* where $S$ is the closure of $\{\mathbf{x}+\mathbf{z};\mathbf{x}\in\mathrm{supp}(\pi),\mathbf{z}\in\mathrm{supp}(\kappa_{\rho})\}$. The notation $\mathrm{supp(h)}=\{\mathbf{x}\in\mathcal{X}\ |\ h(\mathbf{x})\neq 0\}$ refers to the support of a function $h\mathrel{\mathop{\mathchar 58\relax}}\mathcal{X}\rightarrow\mathbb{R}$.* 3. iii)

If both $\pi$ and $\kappa_{\rho}$ are log-concave, then $\pi_{\rho}$ is log-concave. 4. iv)

If $\kappa_{\rho}\in\mathcal{C}^{\infty}(\mathbb{R}^{d})$ and $|\partial^{k}\kappa_{\rho}|$ is bounded for all $k\geq 0$, then $\pi_{\rho}$ is infinitely differentiable w.r.t. $\boldsymbol{\theta}$.

Proposition 1 permits to draw several conclusions about the inference based on $\pi_{\rho}$. Firstly, the infinite differentiability of $\pi_{\rho}$ (Property iv)) implies that it stands for a smooth approximation of $\pi$, see Figure 5 in Section 5.2. Secondly, Property i) of Proposition 1 is reassuring regarding the inference task. Indeed, if $\pi$ stands for a prior distribution, then considering the approximation $\pi_{\rho}$ simply corresponds to a more diffuse prior knowledge around the same expected value, see Section 5.2. Thus, more weight will be given to the likelihood if a posterior distribution is derived with this prior. On the other hand, if $\pi$ stands for a likelihood, then considering the approximation $\pi_{\rho}$ yields the opposite behavior: the likelihood becomes less informative w.r.t. the prior. This idea is directly related to robust hierarchical Bayesian models discussed in Section 3.3.

We now provide quantitative bounds on the approximation implied by considering the marginal $\pi_{\rho}$ instead of $\pi$. For $p\geq 1$, we define the $p$-Wasserstein distance between $\pi$ and $\pi_{\rho}$ by

[TABLE]

where $\Gamma(\pi_{\rho},\pi)$ is the set of probability distributions $\mu(\boldsymbol{\theta},\mathbf{z})$ with marginals $\pi_{\rho}$ and $\pi$ w.r.t. $\boldsymbol{\theta}$ and $\mathbf{z}$, respectively. Under mild assumptions on the kernel $K$, Proposition 2 gives a simple and practical upper bound on (19).

Proposition 2**.**

Assume that $\pi_{\rho}$ in (3) stands for a pdf associated to the variable $\boldsymbol{\theta}$. Let $p\geq 1$ such that $m_{p}\triangleq\left(\displaystyle\int_{\mathbb{R}^{d}}\left\|\mathbf{u}\right\|_{2}^{p}K(\mathbf{u})\mathrm{d}\mathbf{u}\right)^{1/p}<\infty$. Then, we have

[TABLE]

Note that (20) holds without assuming additional assumptions on the initial density $\pi$ such as infinite differentiability. If the latter is assumed w.r.t. the parameter of interest $\boldsymbol{\theta}$, then one can estimate the bias $\pi-\pi_{\rho}$ with a Taylor expansion of $\pi$ similarly to bias analysis in ABC, see Sisson et al. (2018b). Table 1 gives closed-form expressions of $m_{2}$ for the multivariate generalizations of standard kernels. One can see that the constant $m_{2}$ has the same dependence w.r.t. the dimension $d$ for the considered standard kernels $K$. Hence, in high-dimensional scenarios, the approximation quality will be more affected by an inappropriate value for the tolerance parameter $\rho$ rather than by the choice of $K$. In Section 5, we illustrate Proposition 2 with numerical experiments.

4.2 Pointwise bias for Bregman divergences

In complement to Section 4.1 where $\kappa_{\rho}$ was built using kernels, we now analyze the bias induced by considering $\pi_{\rho}$ when $\kappa_{\rho}$ is derived from a Bregman divergence $d_{\psi}$ (see Definition 1), that is

[TABLE]

Under regularity assumptions on both $\pi$ and $\kappa_{\rho}$, Proposition 3 shows the dependence of the pointwise bias $\pi_{\rho}-\pi$ w.r.t. to the tolerance parameter $\rho$ when the latter is sufficiently small.

Proposition 3**.**

Assume that $\pi$ is analytic and twice differentiable on $\mathbb{R}^{d}$ and so does $d_{\psi}$ w.r.t. its first argument. Let $\boldsymbol{\theta}\in\mathbb{R}^{d}$ such that both $\mathbf{H}_{\pi}(\boldsymbol{\theta})$ and $\mathbf{H}_{d_{\psi}}(\boldsymbol{\theta})^{-1}$ exist and are continuous, where $\mathbf{H}_{\pi}(\boldsymbol{\theta})$ is the Hessian matrix of $\pi$ and \mathbf{H}_{d_{\psi}}(\boldsymbol{\theta})\triangleq\frac{\partial^{2}d_{\psi}(\mathbf{z},\boldsymbol{\theta})}{\partial\mathbf{z}^{2}}\Bigr{|}_{\begin{subarray}{c}\mathbf{z}=\boldsymbol{\theta}\end{subarray}} is the Hessian matrix associated to $d_{\psi}(\cdot,\boldsymbol{\theta})$. Then, if

• •

$\left\|\mathbf{H}_{\pi}\right\|\leq C<\infty$**

• •

$\left\|\mathbf{H}_{d_{\psi}}\right\|\geq c>0$,

it follows that

[TABLE]

In addition, if we have $\displaystyle\int_{\mathbb{R}^{d}}\mathbf{u}\kappa_{\rho}(\boldsymbol{\theta}-\sqrt{\rho}\mathbf{u},\boldsymbol{\theta})\mathrm{d}\mathbf{u}=\mathbf{0}_{d}$, then

[TABLE]

Note that when $\psi(\mathbf{z})=\left\|\mathbf{z}\right\|_{2}^{2}/2$, $\kappa_{\rho}$ stands for a Gaussian smoothing kernel, see Section 4.1. In that case, we have the sanity check that the dependence w.r.t. $\rho$ of the bias between $\pi$ and $\pi_{\rho}$ in (23) is the same as the one derived by Sisson et al. (2018b) when interpreting $\kappa_{\rho}$ as a kernel.

4.3 A detailed non-asymptotic analysis for Gaussian smoothing

The previous sections gave quantitative approximation results for a large class of densities $\kappa_{\rho}$ built either via a kernel or a Bregman divergence. In this section, we provide complementary results by restricting our analysis on the case

[TABLE]

This particular yet convenient assumption will allow to complement and sharpen results of Section 4.1 by deriving quantitative bounds which take into account the regularity properties of $f$. Furthermore, these bounds can be extended to a sum of potential functions $f=\sum_{i}f_{i}$ and used to assess the bias associated to both log-densities and credibility regions. This analysis is also motivated by the fact that the Gaussian smoothing case has been widely advocated in the literature since it generally leads to simple inference steps (Holmes and Mallick 2003; Giovannelli 2008; Liechty et al. 2009; Dümbgen and Rufibach 2009), and can be related to both the ADMM in optimization (Boyd et al. 2011; Vono et al. 2019) and the approximation involved in proximal MCMC methods (Pereyra 2016; Durmus et al. 2018; Salim et al. 2019). Unfortunately, a straigthforward generalization of the proof techniques used in the sequel does not give informative upper bounds for smoothing associated to other Bregman divergences.

4.3.1 Assumptions

To derive non-asymptotic bounds between quantities related to $\pi_{\rho}$ defined in (3) and $\pi$ in (1), some complementary assumptions on $f=-\log\pi$ will be required. They are detailed hereafter. For simplicity and with a little abuse of notations, we also denote here by $f(\boldsymbol{\theta})$ the potential associated to (1) when $\pi(\mathbf{y}|\boldsymbol{\theta})$ stands for a likelihood.

1. ($A_{1}$)

$f$ is $L_{f}$-Lipschitz w.r.t. $\left\|\cdot\right\|_{2}$, that is $\exists\ L_{f}\geq 0$ such that for all $\boldsymbol{\theta},\boldsymbol{\eta}\in\mathbb{R}^{d}$, $|f(\boldsymbol{\theta})-f(\boldsymbol{\eta})|\leq L_{f}\left\|\boldsymbol{\theta}-\boldsymbol{\eta}\right\|_{2}$. When $\pi$ is a likelihood, it is further assumed that $L_{f}$ is independent of $\mathbf{y}$. 2. ($A_{2}$)

$f$ is continuously differentiable and has an $M_{f}$-Lipschitz continuous gradient w.r.t. $\left\|\cdot\right\|_{2}$, that is $\exists\ M_{f}\geq 0$ such that for all $\boldsymbol{\theta},\boldsymbol{\eta}\in\mathbb{R}^{d}$, $\left\|\nabla f(\boldsymbol{\theta})-\nabla f(\boldsymbol{\eta})\right\|_{2}\leq M_{f}\left\|\boldsymbol{\theta}-\boldsymbol{\eta}\right\|_{2}$. 3. ($A_{3}$)

$f$ is convex, that is for every $\alpha\in[0,1]$, $\boldsymbol{\theta},\boldsymbol{\eta}\in\mathbb{R}^{d}$, $f(\alpha\boldsymbol{\theta}+(1-\alpha)\boldsymbol{\eta})\leq\alpha f(\boldsymbol{\theta})+(1-\alpha)f(\boldsymbol{\eta})$. 4. ($A_{4}$)

$\mathsf{M}_{f}=\displaystyle\int_{\mathbb{R}^{d}}\left\|\nabla f(\boldsymbol{\theta})\right\|_{2}^{2}\pi(\boldsymbol{\theta})\mathrm{d}\boldsymbol{\theta}<\infty$.

Assumptions $(A_{1})$, $(A_{2})$ and $(A_{3})$ on the potential $f$ stand for standard regularity assumptions in the optimization literature and cover a large class of functions $f$ (Beck and Teboulle 2009; Bolte et al. 2014). In the broad statistical community, $(A_{1})$ has been used by Durmus et al. (2018) to derive non-asymptotic bounds on the total variation distance between probability distributions while $(A_{2})$ stands for a sufficient condition to have a strong solution to the overdamped Langevin stochastic differential equation (Durmus and Moulines 2017).

Under the previous assumptions (not used all at once), non-asymptotic upper bounds on the total variation distance between $\pi_{\rho}$ and $\pi$ are derived in Section 4.3.2. Then, Sections 4.3.3 and 4.3.4 take advantage of this bound to state theoretical properties on the potential functions and credibility regions.

4.3.2 Non-asymptotic bounds on the total variation distance

In this section, we make additional regularity assumptions on the potential $f$ in order to show quantitative results depending explicitly on regularity constants associated to $f$. Two different cases will be considered, namely Lipschitz potentials, and differentiable, gradient-Lipschitz and convex ones.

Lipschitz potential – When the potential function $f$ is assumed to be Lipschitz continuous but not necessarily continuously differentiable, the following result holds.

Theorem 1**.**

Let a potential function $f$ satisfy $(A_{1})$. Then,

[TABLE]

where

[TABLE]

The function $D_{-d}$ is a parabolic cylinder function defined for all $d>0$ and $z\in\mathbb{R}$ by

[TABLE]

As expected from Property 1, note that this bound tends towards zero when $\rho\rightarrow 0$. Additionally, this bound depends on few quantities that can be computed, bounded or approximated in real applications: the dimension of the problem $d$, the Lipschitz constant $L_{f}$ associated to the regularized potential $f$ and the tolerance parameter $\rho$. In the limiting case $\rho\rightarrow 0$, the following equivalent function for the upper bound derived in (25) holds.

Corollary 1**.**

In the limiting case $\rho\rightarrow 0$, we have:

[TABLE]

where for all $z>0$ as $\Gamma(z)=\displaystyle\int_{0}^{+\infty}x^{z-1}e^{-x}\mathrm{d}x$.

Under some regularity conditions (here Lipschitz continuity) on the potential function $f$, Proposition 1 states that $\left\|\pi_{\rho}-\pi\right\|_{\mathrm{TV}}$ grows at most linearly w.r.t. the parameter $\rho$ and w.r.t. $L_{f}$ when $\rho$ is sufficiently small. Moreover, using Stirling-like approximations when $d$ is large in the equivalence relation (28) may give a mild dependence on the dimensionality of the problem in $\mathcal{O}\left(L_{f}d^{1/2}\right)$. Potential functions verifying the hypothesis of Theorem 1 are common in machine learning and signal/image processing problems, see Section 3 in the online supplementary material. As an archetypal example, the sparsity promoting potential function defined for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$ by $f(\boldsymbol{\theta})=\tau\left\|\boldsymbol{\theta}\right\|_{1}$ with $\tau>0$ is Lipschitz continuous with Lipschitz constant $L_{f}=\tau\sqrt{d}$ and satisfies Theorem 1 and Proposition 1. In this case, the dependence of (28) is linear w.r.t. $d$ when $d$ is large and $\rho$ is small. Note also that continuously differentiable functions on a compact set are Lipschitz continuous.

Convex and gradient-Lipschitz potential – We now show a complementary result by assuming $f$ to be convex and continuously differentiable with a Lipschitz-continuous gradient.

Theorem 2**.**

Let a potential function $f$ satisfy $(A_{2})$, $(A_{3})$ and $(A_{4})$. Then, when $\pi$ stands for a pdf associated to $\boldsymbol{\theta}$, we have:

[TABLE]

In the limiting case $\rho\rightarrow 0$, the upper bound in (29) has a simpler expression as shown hereafter.

Corollary 2**.**

In the limiting case $\rho\rightarrow 0$, we have:

[TABLE]

Note that the dependences w.r.t. both $\rho$ and $d$ in Corollary 1 and 2 are similar to the ones found by Nesterov and Spokoiny (2017) for optimization purposes.

Figure 2 gives the behavior of the upper bounds in (25) and (29) w.r.t. the dimensionality $d$ of the problem ranging from $1$ to $10^{6}$ and as a function of $\rho$ in log-log scale. The linear (resp. quadratic) relation between this upper bound and $\rho$ shown in (28) (resp. (30)) is clearly observed for small values of $\rho$. Nonetheless, these upper bounds are not a silver bullet. Indeed, as expected, for a fixed value of the parameter $\rho$, the approximation error increases as the dimension $d$ grows. Thus, these bounds suffer from the curse of dimensionality and become non-informative in high-dimension if $\rho$ is not sufficiently small.

Theorem 1 is easily extended to the case where the initial density $\pi$ is expressed as a product of several terms. If $\pi$ stands for the pdf associated to the variable $\boldsymbol{\theta}$, this boils down to considering

[TABLE]

where for all $j\in[J]$, $f_{j}\mathrel{\mathop{\mathchar 58\relax}}\mathbb{R}^{d}\rightarrow\mathbb{R}$, and a natural generalization of AXDA when applied to each $\pi_{j}$, which writes

[TABLE]

Under this product form, we have the following corollary.

Corollary 3**.**

For all $j\in[J]$, let $f_{j}$ satisfy $(A_{1})$. Then,

[TABLE]

where $\Delta_{d}^{(j)}(\rho_{j})=D_{-d}(L_{f_{j}}\rho_{j})/D_{-d}(-L_{f_{j}}\rho_{j})$.

Unfortunately, Theorem 2 cannot be extended to the multiple splitting scenario. We are nevertheless confident that quantitative bounds can be found with different proof techniques but this task goes beyond the scope of this paper.

4.3.3 Uniform bounds on potentials

From an optimization point of view, it is quite common to consider potential functions associated to densities. For such applications, we give hereafter a quantitative uniform bound on the difference between the potential functions associated to $\pi$ and $\pi_{\rho}$. Similarly to the definition of the potential function $f$ in (1), we define the potential function $f_{\rho}$ associated to the approximate marginal $\pi_{\rho}$ in (3), for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$, by

[TABLE]

By considering a Gaussian smoothing kernel $\kappa_{\rho}$, the potential $f_{\rho}$ becomes

[TABLE]

Note that $f_{\rho}(\boldsymbol{\theta})$ appears as a regularized version of $f(\boldsymbol{\theta})$.

Proposition 4**.**

Let $f$ satisfy $(A_{1})$. Then, for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$,

[TABLE]

with

[TABLE]

and

[TABLE]

It is easily observed that these bounds are informative in the limiting case $\rho\rightarrow 0$ since they both tend towards zero.

4.3.4 Uniform bounds on credibility regions

When $\pi$ stands for the density associated to a posterior distribution, one advantage of Bayesian analysis is its ability to derive the underlying probability distribution of the variable of interest $\boldsymbol{\theta}$ and thereby to provide credibility information under this distribution. This uncertainty information is particularly relevant and essential for real-world applications. Since the marginal $\pi_{\rho}$ stands for an approximation of the original target distribution $\pi$, it is important to control the credibility regions under $\pi_{\rho}$ w.r.t. those drawn under $\pi$. The control in total variation distance given by Theorem 1 is already a good indication. However, it is possible to quantify more precisely the difference between the credible regions (Robert 2001) with confidence level ($1-\alpha$) under $\pi_{\rho}$ and $\pi$, as stated below.

Proposition 5**.**

Let $\pi$ be a posterior distribution associated to $\boldsymbol{\theta}$ and $f$ such that $(A_{1})$ is verified. Let $\mathcal{C}_{\alpha}^{\rho}$ an arbitrary $(1-\alpha)$-credibility region under $\pi_{\rho}$, that is $\mathbb{P}_{\pi_{\rho}\text{}}\left(\boldsymbol{\theta}\in\mathcal{C}_{\alpha}^{\rho}\right)=1-\alpha$ with $\alpha\in(0,1)$. Then,

[TABLE]

where $N_{\rho}$ is defined in (39).

Proposition 5 states that the coverage of $\pi$ under $\mathcal{C}_{\alpha}^{\rho}$ can be determined for a fixed value of $\rho$. Thus, it is even possible to obtain a theoretical comprehensive description of $\mathcal{C}_{\alpha}^{\rho}$ w.r.t. the initial target density $\pi$ before conducting an AXDA-based inference. The bounds in (40) permit to choose a parameter $\rho$ in order to ensure a prescribed coverage property. The behavior of these bounds w.r.t. $\rho$ is the same as in Section 4.3.2, i.e., linear behavior w.r.t. $\rho$ when this parameter is sufficiently small.

5 Numerical illustrations

This section illustrates the quantitative results shown in Sections 4.1 and 4.3 on three different examples which classically appear in statistical signal processing and machine learning. As shown in Table 1, the bias induced by considering $\pi_{\rho}$ is mostly driven by the value of the tolerance parameter $\rho$ rather than by the choice of $\kappa_{\rho}$. Hence, for simplicity, most of the numerical illustrations hereafter consider the case where $\kappa_{\rho}$ is a Gaussian smoothing kernel. Additional illustrations can be found in the online supplementary material.

5.1 Multivariate Gaussian example

We start by performing a sanity check with the simple case where $\pi$ stands for a multivariate Gaussian density that is

[TABLE]

where $\mathbf{\Sigma}$ is assumed to be positive definite. If $\kappa_{\rho}(\cdot,\boldsymbol{\theta})$ is taken to be Gaussian density with mean $\boldsymbol{\theta}$ and covariance matrix $\rho^{2}\mathbf{I}_{d}$, then one can show that

[TABLE]

In particular, let consider the univariate setting, that is $\Theta=\mathbb{R}$, $\mathbf{\Sigma}=\sigma^{2}$. In this case, the variance under $\pi_{\rho}$ is $\sigma^{2}+\rho^{2}$ and simply corresponds to the variance under $\pi$ inflated by a factor $\rho^{2}$. Therefore, the approximation will be reasonable if $\rho^{2}/\sigma^{2}$ is sufficiently small, see Figure 3. In this Figure, we also show the approximation induced by considering a uniform kernel instead of a Gaussian one. The smoothing via the uniform kernel performs slightly better than Gaussian smoothing due to its lower variance. In both cases, the approximation is reasonable for small $\rho$ although $\pi_{\rho}$, built with a uniform kernel, no longer belongs to the Gaussian family.

In order to illustrate the proposed upper bounds on both 2-Wasserstein and total variation distances, we consider a covariance matrix $\mathbf{\Sigma}$ which stands for a squared exponential matrix commonly used in applications involving Gaussian processes (Higdon 2007) and which writes

[TABLE]

where $a=1.5$, $s_{i,i\in[d]}$ are regularly spaced scalars on $[-3,3]$ and $\delta_{ij}=1$ if $i=j$ and zero otherwise.

Figure 4 shows the behavior of the quantitative bounds derived in Proposition 2 and Theorem 2 for $d\in\{10,100\}$. The Gaussian case allows to compute exactly $W_{2}(\pi,\pi_{\rho})$ by noting that $W_{2}^{2}(\pi,\pi_{\rho})=\mathrm{Trace}(\mathbf{\Sigma}+\rho^{2}\mathbf{I}_{d}-2\rho\mathbf{\Sigma}^{1/2})$. On the other hand, $\left\|\pi-\pi_{\rho}\right\|_{\mathrm{TV}}$ has been estimated by using a Monte Carlo approximation. One can note that the general upper bound on the 2-Wasserstein distance is quite conservative for small $\rho$ since it does not catch the behavior in $\mathcal{O}(\rho^{2})$ when $\rho$ is small. This is essentially due to the fact that this bound only assumes a finite moment property and does not require any regularity assumptions on $\pi$ such as differentiability or strong convexity of its potential. On the contrary, the bound on the total variation distance, derived under stronger assumptions, manages to achieve a rate of the order $\mathcal{O}(\rho^{2})$ for small $\rho$.

5.2 Sparse linear regression

We study here a generalized version of the least absolute shrinkage and selection operator (lasso) regression problem analyzed by Park and Casella (2008). We assume a standard linear regression problem where centered observations $\mathbf{y}\in\mathbb{R}^{n}$ are related to the unknown parameters $\boldsymbol{\theta}\in\mathbb{R}^{d}$ via the model $\mathbf{y}=\mathbf{X}\boldsymbol{\theta}+\boldsymbol{\varepsilon}$, where $\mathbf{X}\in\mathbb{R}^{n\times d}$ stands for a known standardized design matrix and $\boldsymbol{\varepsilon}\sim\mathcal{N}(\mathbf{0}_{n},\sigma^{2}\mathbf{I}_{n})$. By considering a generalized Laplacian prior distribution for $\boldsymbol{\theta}$, the target posterior distribution has density for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$,

[TABLE]

where $g(\mathbf{B}\boldsymbol{\theta})=\tau\left\|\mathbf{B}\boldsymbol{\theta}\right\|_{1}$ with $\tau>0$ and $\mathbf{B}\in\mathbb{R}^{k\times d}$ an arbitrary matrix acting on $\boldsymbol{\theta}$. The choice of such a prior may promote a form of sparsity (lasso). For instance, this matrix $\mathbf{B}$ might stand for a $p$-th order difference operator (Bredies et al. 2010) which is highly used in signal and image processing problems. As an archetypal example, the case $p=1$ leads to the well-known total variation regularization function (Chambolle et al. 2010) used to recover piecewise constant signals.

Note that because of the presence of the matrix $\mathbf{B}$, finding an exact data augmentation leading to an efficient sampling scheme is not possible for the general case $\mathbf{B}\neq\mathbf{I}_{d}$. Instead, an AXDA model makes the posterior sampling task possible. Indeed, by regularizing the prior with a Gaussian term, the joint density $\pi_{\rho}$ writes

[TABLE]

By resorting to a Gibbs algorithm to sample from (45), one can now use a simple data augmentation scheme (Park and Casella 2008) to sample from the $\mathbf{z}$-conditional. On the other hand, sampling from the $\boldsymbol{\theta}$-conditional, which is a multivariate Gaussian distribution, can be undertaken efficiently with state-of-the-art approaches (Papandreou and Yuille 2010; Barbos et al. 2017; Marnissi et al. 2018).

In this specific case, the potential $g_{\rho}$ associated to the smoothed prior distribution (see (35)) has a closed-form expression given for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$, by

[TABLE]

with $a(\boldsymbol{\theta})=\sqrt{\pi\rho^{2}/2}\exp\left(-(\mathbf{b}_{i}^{T}\boldsymbol{\theta})^{2}/(2\rho^{2})\right)$, $b(\boldsymbol{\theta})=\sqrt{\rho^{2}/2}(\tau-\mathbf{b}_{i}^{T}\boldsymbol{\theta}/\rho^{2})$, $c(\boldsymbol{\theta})=\sqrt{\rho^{2}/2}(\tau+\mathbf{b}_{i}^{T}\boldsymbol{\theta}/\rho^{2})$ and $\mathbf{b}_{i}\in\mathbb{R}^{d}$ standing for the $i$-th row of $\mathbf{B}$. Note that in more general cases where $g_{\rho}$ has no closed form, one can estimate it by a Monte Carlo approximation.

Figure 5 shows the behavior of the regularized potential $g_{\rho}$ defined in (47) for several values of the parameter $\rho$ along with the associated smoothed prior and posterior distributions. For simplicity and pedagogical reasons, the univariate case corresponding to $\boldsymbol{\theta}=\theta_{1}\in\mathbb{R}$ and $\mathbf{B}=1$ has been considered. The regularization parameter $\tau$ has been set to $\tau=1$. The contours of the shaded area correspond to $g+L_{\rho}$ and $g+U_{\rho}$. The potential $g_{\rho}$ is a smooth approximation of the potential $g$ associated to the initial prior as expected, see Property iv) in Proposition 1. Note that the inequalities derived in (36) are verified. Although this approximation seems similar to the Moreau-Yosida regularization of a non-smooth potential function (Combettes and Pesquet 2011), the rationale behind this approximation is different. Indeed, the Moreau-Yosida envelope stands for a particular instance of the infimal convolution between two convex functions (an initial potential and a Gaussian one). On the other hand, $g_{\rho}$ is the potential associated to a smoothed density obtained by convolution with a Gaussian kernel. In addition, the third row of Figure 5 shows the form of the posterior of $\theta_{1}$ defined in (45) for $y=1$, $x=2$ and $\sigma=1$ and derived from the smoothed prior distributions shown in Figure 5. For sufficiently small values of $\rho$, the marginal $\pi_{\rho}$ stands for a quite accurate approximation of the original target $\pi$.

Table 2 illustrates the bounds derived in (40) for $\rho\in\left\{1,10^{-1},10^{-2},10^{-3}\right\}$. For each case, the values of the bounds are summarized in the interval

[TABLE]

and the real coverage $\int_{\mathcal{C}_{\alpha}^{\rho}}\pi(\theta_{1})\mathrm{d}\theta_{1}$ is also reported. The (1-$\alpha$)-credibility intervals $\mathcal{C}_{\alpha}$ and $\mathcal{C}_{\alpha}^{\rho}$ have been chosen to be the highest posterior density regions associated to each density with $\alpha=0.05$. Note that the theoretical coverage interval $\mathcal{I}_{\alpha}^{\rho}$ becomes informative only if $\rho$ is sufficiently small which is not surprising since the assumptions on the potential of $\pi_{\rho}$ are weak. Indeed, the form of the density (e.g. symmetry or unimodality) is not taken into account in the derived bounds. Regarding the empirical value of the coverage $\int_{\mathcal{C}_{\alpha}^{\rho}}\pi(\theta_{1})\mathrm{d}\theta_{1}$, we emphasize that the marginal $\pi_{\rho}$ stands for a conservative approximation of $\pi$ in this example. Indeed, in each case, the (1-$\alpha$)-credibility interval under $\pi_{\rho}$ denoted $\mathcal{C}_{\alpha}^{\rho}$ covers at least $100(1-\alpha)$% of the probability mass under $\pi$.

5.3 Illustration on an image inpainting problem

We illustrate here the correctness of the proposed approach on a multidimensional and non-Gaussian example which classically appears in image processing. To this purpose, we consider the observation of a damaged and noisy image $\mathbf{y}\in\mathbb{R}^{n}$ (represented as a vector by lexicographic ordering) related to the unknown original image $\boldsymbol{\theta}\in\mathbb{R}^{d}$ by the linear model

[TABLE]

where $n<d$, $\mathbf{H}\in\mathbb{R}^{n\times d}$ stands for a decimation binary matrix. The dimension $d$ being typically large (e.g., $10^{3}\leq d\leq 10^{9}$), these problems require scalable inference algorithms. Since the matrix $\mathbf{H}$ is not invertible, the linear inverse problem (49) is ill-posed. To cope with this issue, we assign the total variation prior distribution to the unknown parameter $\boldsymbol{\theta}$, leading to the posterior distribution

[TABLE]

where $\tau>0$ is a regularization parameter, $\mathbf{D}\boldsymbol{\theta}=(\mathbf{D}_{1}\boldsymbol{\theta},\mathbf{D}_{2}\boldsymbol{\theta})\in\mathbb{R}^{2\times d}$ is the two-dimensional discrete gradient associated to the image $\boldsymbol{\theta}$ and the notation $\mathbf{M}_{i}$ stands for the $i$-th column of the matrix $\mathbf{M}$, see Chambolle et al. (2010) for more details about the total variation regularization. The presence of the operator $\mathbf{D}$ and the non-differentiability of the total variation norm rule out the use of common data augmentation schemes and simulation-based algorithms (e.g., Hamiltonian and Langevin Monte Carlo methods). Possible surrogates are proximal MCMC methods (Pereyra 2016; Durmus et al. 2018) which replace the non-differentiable posterior distribution by a smooth approximation based on the proximity operator (Combettes and Pesquet 2011) of the total variation norm. However, the latter does not admit a closed-form expression and iterative routines are commonly used to approximate the latter (Chambolle 2004) leading to higher computational costs.

To mitigate these issues, we propose to rely on a particular instance of AXDA by smoothing the total variation prior with a Gaussian term, leading to the approximate joint posterior density

[TABLE]

where $\mathbf{Z}=(\mathbf{z}_{1},\mathbf{z}_{2})\in\mathbb{R}^{2\times d}$. By relying on (51), the inference is now simplified and can be conducted with a Gibbs sampler, see Section 4 in the supplementary material. Since $\mathrm{ker}(\mathbf{H})\cap\mathrm{ker}(\mathbf{D})=\{\mathbf{0}_{d}\}$, the conditional posterior distribution of $\boldsymbol{\theta}$ is a non-degenerate multivariate Gaussian distribution. Samples from the latter can be obtained efficiently with the two-dimensional discrete Fourier transform by exploiting the periodic boundary conditions for $\boldsymbol{\theta}$ (Wang et al. 2008; Marnissi et al. 2018). On the other hand, samples from $\pi_{\rho}(\mathbf{Z}|\boldsymbol{\theta})$ can be drawn efficiently using exact data augmentation, see Kyung et al. (2010). All the inference details are given in Section 4 of the supplementary material.

We illustrate the proposed approximate model $\pi_{\rho}$ by considering the Shepp-Logan phantom magnectic resonance image of size $100\times 100$ ($d=10^{4}$), see Figure 6. We artifically damaged and added noise to this image to build a noisy observation $\mathbf{y}$ consisting of 90% randomly selected pixels of the initial image. The standard deviation of the Gaussian noise and the regularization parameter have been set to $\sigma=7\times 10^{-2}$ (corresponding to a SNR of 58dB) and $\tau=5$, respectively. The tolerance parameter has been set to $\rho=0.1$.

In order to assess the bias of the proposed approach, we implemented the Moreau-Yosida Metropolis-adjusted Langevin algorithm (MYMALA) of Durmus et al. (2018), specifically designed to sample exactly from high-dimensional and non-smooth posterior distributions. For all the MCMC algorithms, the initialization has been set to $\boldsymbol{\theta}^{[0]}=\mathbf{0}_{d}$. We generated $10^{5}$ samples and kept the last $5\times 10^{4}$ ones.

Figure 6 shows the minimum mean square estimate (MMSE) under $\pi_{\rho}$ along with the original image. One can denote that the MMSE under $\pi_{\rho}$ is visually similar to the original image and hence coherent with the reconstruction task. The relative residual error between the former and the MMSE under $\pi$ is of order 2%. The main differences are located on the boundaries of the image, as depicted in the figure on the left which shows the absolute difference between the pixels of the two posterior means $\mathbb{E}_{\pi}(\boldsymbol{\theta})$ and $\mathbb{E}_{\pi_{\rho}}(\boldsymbol{\theta})$.

To emphasize the correctness of the proposed approach beyond the comparision between pointwise estimates, we also paid attention to the comparison between posterior credibility sets induced by both $\pi$ and $\pi_{\rho}$. To this purpose, we considered the highest posterior density region given by

[TABLE]

where $\gamma_{\alpha}\in\mathbb{R}$ is such that $\int_{\mathcal{C}_{\alpha}^{\star}}\pi(\boldsymbol{\theta}|\mathbf{y})\mathrm{d}\boldsymbol{\theta}=1-\alpha$ and $f$ is the potential function associated to $\pi(\boldsymbol{\theta}|\mathbf{y})$.

Figure 7 shows the different values of the scalar summary $\gamma_{\alpha}$ estimated using $\pi$ and the scalar $\gamma_{\alpha}^{\rho}$ estimated using $\pi_{\rho}$ for $\alpha\in[0.01,0.99]$. Note that the approximation error associated to $\gamma_{\alpha}$ is of order 2.6% whatever the value of $\alpha$, which supports the use of $\pi_{\rho}$ to conduct Bayesian uncertainty analysis in this problem. After the burn-in period, the efficiency of the Gibbs algorithm used to sample from $\pi_{\rho}$ has been measured by comparing the effective sample size (ESS) associated to the slowest component of $\boldsymbol{\theta}$ to the one obtained with MYMALA. We found that the two ESS were roughly similar but the cost per iteration of the Gibbs sampler (0.079 sec/iteration) is almost two times lower than that of MYMALA (0.144 sec/iteration)111Both algorithms have been implemented in Matlab with the same level of efficiency.. In addition, the number of iterations required to reach high-probability regions is much less important for the Gibbs sampler than for MYMALA, showing the interest of AXDA, see Figure 7.

6 Conclusion

This paper presented a unifying framework for asymptotically exact data augmentation (AXDA) schemes. AXDA introduces approximate densities in order to simplify the inference. By building on existing works which considered special instances of AXDA, we illustrate potential benefits that can be inherited by the proposed framework such as distributed computations, robustness or sophisticated inference schemes from the ABC literature. On top of these qualitative properties, we derived a set of theoretical guarantees on the bias involved in the proposed methodology. The latter encompass a large class of AXDA models and a detailed non-asymptotic analysis has been done for Gaussian smoothing. These results have been illustrated on several cases that can arise in statistical learning or signal processing showing the broad scope of application of the proposed approach. In practice, we emphasize that AXDA models can remarkably improve the inference task in big data and high-dimensional settings. In summary, at the price of an approximation which comes with theoretical guarantees, AXDA approaches appear to be a general, systematic and efficient way to conduct simple inference in a wide variety of large-scale problems. They provide accurate estimates with relevant confidence intervals that are crucial in many applications, in particular when no ground truth is available.

Acknowledgements

Part of this work has been supported by the ANR-3IA Artificial and Natural Intelligence Toulouse Institute (ANITI) under grant agreement ANITI ANR-19-PI3A-0004.

SUPPLEMENTARY MATERIALS

Appendices:

The supplementary material includes the proofs of Theorems 1 and 2; Propositions 1, 2, 3, 4 and 5; and Corollary 1 and 3. It also includes additional details about standard kernels and Bregman divergences, inference details associated to the image inpainting example in Section 5.3 and derivations of classical inference algorithms to target AXDA models. (supplementary_material.pdf, pdf file)

Package for AXDA:

The computer code associated to the illustrations and experiments described in this paper is also available online. More precisely, the package “AXDA” contains a Python jupyter notebook to reproduce all the tables and figures of the paper and a Matlab-code associated to the image inpainting example along with a README file. (axda.zip, zipped file)

SUPPLEMENTARY MATERIALS

Appendix 1 Proofs

1.1 Proof of Proposition 1

Property i) follows from the fact that $\pi_{\rho}$ stands for a convolution integral between $\pi$ and $\kappa_{\rho}$, i.e. $\pi_{\rho}=\pi\ast\kappa_{\rho}$. Therefore, the expectation and variance under $\pi_{\rho}$ are the sum of the expectations and variances of two independent random variables under $\pi$ and $\kappa_{\rho}$ respectively. Property ii) follows directly from (Folland 1999, Proposition 8.6). Property iii) follows from the fact that log-concavity is preserved by convolution of distributions (Dharmadhikari and Joag-Dev 1988, Theorem 2.18). Finally, Property iv) follows from the dominated convergence theorem since $\pi\in L^{1}$, $\kappa_{\rho}\in\mathcal{C}^{\infty}(\mathbb{R}^{d})$ and for all $k\geq 0$, $|\partial^{k}\kappa_{\rho}|\leq C_{k}$ (Folland 1999, Proposition 8.10).

1.2 Proof of Proposition 2

The proof can be found in Ambrosio et al. (2008, Lemma 7.1.10). Since it is quite short, we recall it hereafter for completeness. We have

[TABLE]

1.3 Proof of Proposition 3

Let $\boldsymbol{\theta}\in\mathbb{R}^{d}$. Since $\pi$ has been assumed to be analytic and twice differentiable with $\mathbf{H}_{\pi}$ being continuous, there exists $\tilde{\boldsymbol{\theta}}$ lying between $\boldsymbol{\theta}$ and $\boldsymbol{\theta}-\sqrt{\rho}\mathbf{u}$ such that

[TABLE]

where $\mathbf{H}_{\pi}$ stands for the Hessian matrix of $\pi$.

It follows that

[TABLE]

We now show that (63) $=O(\sqrt{\rho})$ and (64) $=\mathcal{O}(\rho)$. To this purpose, we use the analyticity and two times differentiability of $d_{\psi}$ w.r.t. to its first argument and the continuity of $\mathbf{H}_{d_{\psi}}$. By definition of the Bregman divergence (see Definition 1 in the main paper), $d_{\psi}(\boldsymbol{\theta},\boldsymbol{\theta})=0$ and \nabla_{\mathbf{z}}d_{\psi}(\mathbf{z},\boldsymbol{\theta})\Bigr{|}_{\begin{subarray}{c}\mathbf{z}=\boldsymbol{\theta}\end{subarray}}=\mathbf{0}_{d} so that, for all $\mathbf{u}\in\mathbb{R}^{d}$,

[TABLE]

where $\boldsymbol{\theta}^{\prime}$ lies between $\boldsymbol{\theta}$ and $\boldsymbol{\theta}-\sqrt{\rho}\mathbf{u}$.

We first prove (64) $=\mathcal{O}(\rho)$. Using (65), we can re-write (64) as

[TABLE]

Since $\lim_{\rho\rightarrow 0}\boldsymbol{\theta}^{\prime}=\boldsymbol{\theta}$ and $\lim_{\rho\rightarrow 0}\tilde{\boldsymbol{\theta}}=\boldsymbol{\theta}$, we will use the dominated convergence theorem using that

[TABLE]

In addition, since $d_{\psi}$ is strictly convex w.r.t. its first argument, $\mathbf{H}_{d_{\psi}}$ is a symmetric and positive-definite matrix. By using that $\left\|\mathbf{H}_{\pi}\right\|\leq C<\infty$ and $\left\|\mathbf{H}_{d_{\psi}}\right\|\geq c>0$, we have:

[TABLE]

which is integrable on $\mathbb{R}^{d}$. From the dominated convergence theorem, it follows that

[TABLE]

Similarly,

[TABLE]

Hence,

[TABLE]

We now prove (63) $=\mathcal{O}(\sqrt{\rho})$. Using (65), it follows that

[TABLE]

Again, since $\left\|H_{d_{\psi}}\right\|$ has been assumed to be lower bounded, it follows from the dominated convergence theorem that

[TABLE]

1.4 Proof of Theorem 1

We are interested in controlling w.r.t. $\rho$ the quantity $\left\|\pi_{\rho}-\pi\right\|_{\mathrm{TV}}$. To this purpose, let assume that $f$ satisfies $(A_{1})$ in the main paper. In the following, we will assume for simplicity reasons that $\pi$ stands for a pdf associated to the random variable $\boldsymbol{\theta}$. The case when $\pi$ is a likelihood is treated right after. Under this convention, it follows

[TABLE]

where $C_{\pi}$ and $C_{\pi_{\rho}\text{}}$ are the normalizing constants associated to $\pi$ and $\pi_{\rho}$, respectively, and

[TABLE]

Note that

[TABLE]

Since $f$ is assumed to be $L_{f}$-Lipschitz, we have

[TABLE]

We make the change of variables $\mathbf{u}=\mathbf{z}-\boldsymbol{\theta}$, which leads to

[TABLE]

Then, with another change of variables $t=\left\|\mathbf{u}\right\|_{2}$, it follows

[TABLE]

This integral admits a closed-form expression (Gradshteyn and Ryzhik 2015, Formula 3.462 1.) by introducing the special parabolic cylinder function $D_{-d}$ defined for all $d>0$ and $z\in\mathbb{R}$ by

[TABLE]

Then,

[TABLE]

where

[TABLE]

Then, with (80) and (85), we also have

[TABLE]

We now use the triangle inequality in (77) which leads to

[TABLE]

The first term in this upper bound writes

[TABLE]

This allows us to bound (88), that is

[TABLE]

Using one more time the $L_{f}$-Lipschitz assumption on $f$, we have for all $\boldsymbol{\theta},\mathbf{z}$,

[TABLE]

With the same changes of variables as above, it follows

[TABLE]

where

[TABLE]

Then we have $1-\dfrac{1}{A(\rho)}\mathcal{K}(\boldsymbol{\theta})\leq 1-\dfrac{B(\rho)}{A(\rho)}$ which combined with (90) yields

[TABLE]

Note: When $\pi=\pi(\mathbf{y}|\boldsymbol{\theta})$ is a likelihood, (77) becomes

[TABLE]

Since $L_{f}$ is assumed to be independent of $\mathbf{y}$, the same type of proof can be followed in this case and yields the same quantitative bound.

1.5 Proof of Corollary 1

The parabolic cylinder function when $d>0$ has the following expression (Gradshteyn and Ryzhik 2015, Formula 9.241 2.)

[TABLE]

In the limiting case when $z\rightarrow 0$, a first order Taylor expansion of $e^{-xz}$ gives

[TABLE]

recording that $\int_{0}^{+\infty}e^{-x^{2}/2}x^{d}\mathrm{d}x=\Gamma((d+1)/2)2^{d/2-1/2}$ (Gradshteyn and Ryzhik 2015, Formula 3.383 11.). Using (98) for $z=\pm\rho L_{f}$ yields

[TABLE]

1.6 Dependence of (99) with respect to the dimension

The gamma function $\Gamma$ can be expressed for all $z>0$ as $\Gamma(z)=\int_{0}^{+\infty}x^{z-1}e^{-x}\mathrm{d}x$. When $z$ is large, Stirling-like approximations give the following equivalent for $\Gamma(z+1/2)$ and $\Gamma(z)$:

[TABLE]

So that when $d$ is large

[TABLE]

1.7 Proof of Theorem 2

We now prove another bound on the TV distance when $f$ satisfies $(A_{2}),(A_{3})$ and $(A_{4})$ in the main paper. The beginning of the proof follows the same lines as in the proof of Theorem 1 above when $\pi$ stands for a pdf associated to $\boldsymbol{\theta}$. Hence, we have from (77) that

[TABLE]

We now use the convexity of $f$ to write for all $\boldsymbol{\theta}\in\mathbb{R}^{d},\mathbf{z}\in\mathbb{R}^{d}$,

[TABLE]

By using (104) and (79), it follows that

[TABLE]

By using again the convexity of $f$, we also have for all $\boldsymbol{\theta}\in\mathbb{R}^{d},\mathbf{z}\in\mathbb{R}^{d}$,

[TABLE]

Then, (106) leads to

[TABLE]

We now use $(A_{2})$ in the main paper which leads to

[TABLE]

We now apply the triangle inequality in (103) which yields

[TABLE]

The absolute value in the first term of (109) can be removed by noting that

[TABLE]

Then (109) becomes

[TABLE]

We now use the fact that $-\exp(-u)\leq u-1$ for all $u\geq 0$ which yields

[TABLE]

The result in Corollary 3 in the main paper comes from a straightforward Taylor expansion of (112).

1.8 Proof of Corollary 3

Equation (79) becomes

[TABLE]

Bounding each term in (113) and following the proof of Theorem 1 detailed above completes the proof.

1.9 Proof of Proposition 4

By using (85) and (93) we have for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$,

[TABLE]

So that

[TABLE]

The result of Proposition 4 follows from the definition of $A(\rho)$ and $B(\rho)$.

1.10 Proof of Proposition 5

By using (85) and (93) it follows, for all $\boldsymbol{\theta}\in\mathbb{R}^{d}$,

[TABLE]

Using (79) yields

[TABLE]

Using (86) and (94) gives

[TABLE]

where the constant $N_{\rho}$ has been defined in (39) in the main paper.

Let $\mathcal{C}_{\alpha}^{\rho}$ an arbitrary $(1-\alpha)$-credibility region under $\pi_{\rho}$. By integrating (119) on $\mathcal{C}_{\alpha}^{\rho}$,

[TABLE]

Since $\mathcal{C}^{\rho}_{\alpha}\subseteq\mathbb{R}^{d}$ and $\int_{\mathbb{R}^{d}}\pi(\boldsymbol{\theta})\mathrm{d}\boldsymbol{\theta}=1$, the upper bound in (120) can be replaced by

$\min\left\{1,\dfrac{N_{\rho}}{D_{-d}(L_{f}\rho)}(1-\alpha)\right\}$.

1.11 Lipschitz loss functions - Dependence w.r.t. the number of observations

Combining Corollary 4 in the main paper with (99), we have:

[TABLE]

Appendix 2 On the choice of $\kappa_{\rho}$

Section 2.2 of the main document provides two distinct ways to choose $\kappa_{\rho}$. We can first consider a kernel $K$ that is a positive function such that $\displaystyle\int_{\mathbb{R}^{d}}K(\mathbf{u})\mathrm{d}\mathbf{u}=1$ and $K(-\mathbf{u})=K(\mathbf{u})$, for all $\mathbf{u}\in\mathbb{R}^{d}$. Based on the latter, we define for all $\mathbf{z},\boldsymbol{\theta}\in\Theta$, $\kappa_{\rho}(\mathbf{z},\boldsymbol{\theta})\propto_{\mathbf{z}}\rho^{-d}K(\rho^{-1}(\boldsymbol{\theta}-\mathbf{z}))$ (Dang and Ehrhardt 2012). Table 3 lists some classical examples of symmetric kernels $K(\cdot)$ which are not necessarily compactly supported. For sake of simplicity, we only define univariate versions of them but they can obviously be generalized in higher dimension. Figure 8 illustrates these kernels. Some of them have for instance been used in ABC approaches (Sisson et al. 2018b).

Another natural choice consists in resorting to a divergence function $\phi$ such that $\kappa_{\rho}(\mathbf{z},\boldsymbol{\theta})\propto_{\mathbf{z}}\exp(-\rho^{-1}\phi(\mathbf{z},\boldsymbol{\theta}))$. Particular instances of these functions are the family of the Bregman divergences, which are ubiquitous tools in signal processing, machine learning and optimization (see Definition 1 of the main document). Table 4 recalls classical examples of Bregman divergences that can be used to define $\phi$. Again, only univariate examples of such potentials are provided but they can be easily extended to the multivariate case.

Note that, reciprocally, tight connections have been already drawn between Bregman divergences and regular exponential family distributions. Indeed, when the function $\kappa_{\rho}$ belongs to the exponential family, its associated potential function defined by $-\log\kappa_{\rho}$ can be expressed as a Bregman divergence up to an additional term (Banerjee et al. 2005, Theorem 4)

[TABLE]

Appendix 3 Illustration for Lipschitz loss functions used in statistical learning

Some of the results of Section 4.3 in the main paper assume that the potential function $f$ associated to $\pi$ is Lipschitz. Interestingly, such Lipschitz functions are used in standard statistical learning problems to evaluate the discrepancy between observations and model outputs (van de Geer 2016).

Table 5 lists some of them along with their definition and associated statistical problems. Note that the absolute loss stands for a particular instance of the pinball loss with $\tau=0.5$. Figure 9 illustrates the form of these losses and associated regularized potentials $f_{\rho}$ with $\rho=1$ obtained via a Monte Carlo approximation.

Without loss of generality, these problems consider a likelihood function that can be written as in (32) in the main paper with

[TABLE]

where for $j\in[n]$, $\mathbf{x}_{j}$ is the feature vector associated with observation $y_{j}$; $f$ is one of the loss functions in Table 5 and $\boldsymbol{\theta}\in\mathbb{R}^{d}$ is the parameter to infer. Since all the loss functions listed in Table 5 are Lipschitz continuous w.r.t. their second argument $t$ with Lipschitz constant equal to 1, the potential $f_{j}$ in (125) is also Lipschitz with constant $L_{f_{j}}=\left\|\mathbf{x}_{j}\right\|_{2}$. Motivated by the robustness properties inherited by AXDA, see Section 3.3 in the main paper, we consider the smoothing of the likelihood contribution associated to each observation $f_{j}$ with a Gaussian kernel. The results of Corollary 3 in the main paper can then be applied to $\pi$ defined in (32) in the main paper.

In practice, to illustrate the behavior of the upper bound in Corollary 3 w.r.t. the number of observations, we fixed the dimension $d$ and considered several values of $n$ ranging from 1 to $10^{4}$. For each $n$, we randomly generated sets of features $\left\{\mathbf{x}_{j}\right\}_{j\in[n]}$ and we normalized the columns of the matrix $\mathbf{X}^{T}=\left[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\right]^{T}$ such that each entry is a random number between [math] and $1$. The latter operation is classical in machine learning and is also called feature scaling.

Figure 10 shows the behavior of the upper bound in Corollary 3 for two values of the dimension $d=10$ and $d=10^{3}$. As expected, the bound becomes less informative for a fixed value of $\rho$ as the number of likelihood approximations increases with the size of the dataset $n$. Nonetheless, the effect of $n$ on the bound is not highly prohibitive. In both cases $d=10$ and $d=10^{3}$, $\rho$ and $n$ appear to be complementary variables: increasing the value of the latter and decreasing the value of the former by the same factor roughly gives the same bound value. Actually, one can show that the dependence of the bound when $\rho$ is small is of the order $\mathcal{O}(n\rho)$ for a fixed dimension $d$, see the supplementary material. Obviously, one can limit this dependence on $n$ by splitting blocks of observations in minibatches instead of splitting each observation. This splitting strategy has for instance been considered by Rendell et al. (2020).

Appendix 4 Inference details for the image inpainting example

In this section, we detail the steps of the Gibbs sampler used to sample from the posterior distribution $\pi_{\rho}(\boldsymbol{\theta},\mathbf{z}|\mathbf{y})$ in Section 5.4 in the main paper.

4.1 Sampling the auxiliary vector

The conditional distribution associated to the auxiliary variable $\mathbf{Z}=(\mathbf{z}_{1},\mathbf{z}_{2})$ writes

[TABLE]

This conditional distribution can be sampled exactly by using data augmentation. Indeed, one can re-write the distribution involving the non-differentiable potential $\left\|\cdot\right\|_{2}$ as a mixture of normal and gamma distributions (Kyung et al. 2010, Section 3.1). Hence, sampling from (126) can be performed with the following two steps

[TABLE]

Note that all these sampling steps can be performed efficiently by “vectorizing” them.

4.2 Sampling the parameter of interest

The conditional distribution associated to the image to recover $\boldsymbol{\theta}$ writes

[TABLE]

The distribution (127) is a non-degenerate Gaussian distribution $\mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{\theta}},\mathbf{\Sigma}_{\boldsymbol{\theta}})$ where

[TABLE]

Sampling from this multivariate distribution can be done efficiently in $O(d\log d)$ floating point operations by resorting to the two-dimensional discrete Fourier transform. Indeed, under periodic boundary conditions for $\boldsymbol{\theta}$, the matrix $\mathbf{D}^{T}\mathbf{D}$ is a block circulant matrix and hence diagonalizable in the Fourier domain. On the other hand, $\mathbf{H}^{T}\mathbf{H}$ stands for a diagonal matrix with some zeros on the diagonal corresponding to the missing pixels. Since these two matrices cannot be diagonalized in the same domain, we use the auxiliary variable method of Marnissi et al. (2018) to decouple them. Let $\eta\left\|\mathbf{H}^{T}\mathbf{H}\right\|_{\mathrm{S}}<\sigma^{2}$ where $\left\|\mathbf{M}\right\|_{\mathrm{S}}$ is the spectral norm of the matrix $\mathbf{M}$. Then, we have the following two-step sampling scheme

[TABLE]

where

[TABLE]

Appendix 5 Inference algorithms based on AXDA

Motivated by the issues detailed in Section 2, the good expected properties reviewed in Section 3 as well as the theoretical results shown in Section 4 in the main paper, this section shows that AXDA may allow to derive more efficient and distributed inference algorithms ranging from simulation to optimization-based methods. To do so, the potential benefits of AXDA in comparison with direct inference from $\pi$ are presented and discussed. MCMC and VB methods based on AXDA models are detailed to explore the distribution of the parameters to infer. Optimization-based approaches such as ADMM and the EM-algorithm are also derived if maximum likelihood (ML) or maximum a posteriori (MAP) estimates based on the approximate density $\pi_{\rho}$ are desired.

From now on, we assume that $\Theta=\mathbb{R}^{d}$ and we consider a target density with the general form

[TABLE]

Based on this target density, the augmented density $\pi_{\rho}$ is assumed to take the form

[TABLE]

This writing permits to highlight the benefits of using the augmented density $\pi_{\rho}$ instead of $\pi$ for each of the different inference approaches detailed in the sequel.

5.1 Monte Carlo sampling from AXDA

A standard way to sample from the joint density $\pi_{\rho}$ is to consider a Gibbs sampler as in Algorithm 1. This algorithm can make the inference tractable, simpler and/or faster by targetting $\pi_{\rho}$ instead of $\pi$. First, by splitting the initial potential $\sum_{j}f_{j}$, $\pi_{\rho}$ admits simpler and lower-dimensional conditional posterior distributions, each of them possibly defined by a single potential $f_{j}$. Within a Gibbs algorithm, these conditional posteriors yield simpler sampling steps, which may embed efficient dedicated state-of-the-art sampling methods. Second, given the current iterate $\boldsymbol{\theta}^{[t]}$, sampling each auxiliary variable $\mathbf{z}_{j}^{[t]}$ can be performed in an independent and parallel manner for a faster inference. This is of particular interest in big data settings where datasets are stored on multiple kernels or machines (Rendell et al. 2020). In addition, Vono et al. (2019) experimentally showed that considering AXDA-based models can even improve the convergence properties of classical MCMC methods such as Langevin Monte Carlo by embedding them. A detailed description of additional benefits of AXDA simulation-based methods and their illustration on image processing and machine learning problems can be found in Rendell et al. (2020) and Vono et al. (2019).

Example 3**.**

We consider in this example the penalized logistic regression problem. We assume that $n$ binary responses $\mathbf{y}\in\{-1,1\}^{n}$ are observed and correspond to conditionally independent Bernoulli random variables with probability of success $\sigma(\mathbf{x}_{j}^{T}\boldsymbol{\theta})$. The function $\sigma$ is the sigmoid function, $\mathbf{x}_{j}\in\mathbb{R}^{d}$ stands for the feature vector associated to observation $y_{j}$ and $\boldsymbol{\theta}\in\mathbb{R}^{d}$ are the unknown regression coefficients to infer. We consider a zero-mean Gaussian prior distribution on $\boldsymbol{\theta}$ with precision $2\tau$, that is $g(\boldsymbol{\theta})=\tau\left\|\boldsymbol{\theta}\right\|_{2}^{2}$. The target $\pi$ then stands for the posterior distribution of the unknown regression coefficients $\boldsymbol{\theta}$

[TABLE]

By denoting $f_{(n+1)}=g$ and for all $j\in[n]$ with $J=n$, $f_{j}(u)=\log\left[1+\exp\left(y_{j}u\right)\right]$, the posterior distribution in (132) has the form (130) with $J=n+1$. Following the work of Polson et al. (2013), one can derive a promising DA scheme from $\pi$ based on the Polya-Gamma distribution. Hence, a Gibbs sampler can be used to sample from each conditional distribution as detailed by Polson et al. (2013). However, this Gibbs sampler scales poorly in high-dimensional settings as pointed out by Durmus and Moulines (2016). The AXDA alternative is Algorithm 1 with a quadratic potential and resulting conditional distributions

[TABLE]

Thanks to this splitting scheme, the inference is simpler, might be distributed, and sampling from these conditional distributions can be done exactly and efficiently. Indeed, since (133) is univariate and log-concave, one can use adaptive rejection sampling (Gilks and Wild 1992) while sampling the variable of interest $\boldsymbol{\theta}$ from (134) boils down to high-dimensional Gaussian sampling and efficient methods can be applied.

5.2 Variational Bayes inference from AXDA

AXDA can also be a major asset when conducting variational Bayes (VB) inference, providing important benefits such as simplicity and parallelization. VB methods (Bishop and Tipping 2000; Opper and Saad 2001) circumvent the direct sampling from a target density such as $\pi_{\rho}$ by defining an approximation of the latter denoted $\tilde{\pi}_{\rho}$. The best approximation is found by minimizing the Kullback-Leibler (KL) divergence between $\tilde{\pi}_{\rho}$ and $\pi_{\rho}$ restricted to a set of tractable candidates $\tilde{\pi}_{\rho}$. Depending on this set, a lot of VB approximation methods exist in the literature, see Bishop (2006) and Pereyra et al. (2016) for reviews. In this section, we will consider the widely-used mean-field approximation method where the approximate density $\tilde{\pi_{\rho}\text{}}$ is chosen among the set of conditionally independent (w.r.t. $\rho$) densities, that is $\tilde{\pi}_{\rho}(\boldsymbol{\theta},\mathbf{z}_{1\mathrel{\mathop{\mathchar 58\relax}}J})=\tilde{\pi}_{\rho}(\boldsymbol{\theta})\prod_{j=1}^{J}\tilde{\pi}_{\rho}(\mathbf{z}_{j})$. Under this constraint, the optimal choice of the VB approximation is given by

[TABLE]

The above VB-marginals require to compute expectations under each marginal distribution which are often functions of moments under each marginal. Similarly to what has been encountered for Gibbs sampling in the previous section, deriving a VB approach based on $\pi_{\rho}$, instead of $\pi$, yields important benefits for parallel and possibly easier computations. Indeed, the VB-marginal in (136) shows again that each potential $f_{j}$ contributes independently given $\boldsymbol{\theta}$. The updates of expectations under (136) are thereby simplified since (i) the VB-marginals (136) are simpler than those obtained from a mean-field approximation of $\pi$ and (ii) the moments under the latter can be computed in parallel or distributed.

Example 3 (continued)****.

Jaakkola and Jordan (2000) considered a local VB algorithm for the penalized logistic regression problem. Instead of using a local VB approach and finding bounds on each individual function $f_{j}$, the use of $\pi_{\rho}$ instead of $\pi$ permits to consider directly a global VB approach such as the mean-field approximation. In addition, similarly to Algorithm 1, the updates of (136) and the associated expectations can be computed in parallel and efficiently by using state-of-the-art existing methods. For instance, the expectations $\mathbb{E}_{\tilde{\pi}_{\rho}(\mathbf{z}_{j})}$ can be approximated efficiently using rejection sampling.

5.3 Optimizing AXDA meets quadratic penalty methods

Computing the MAP or ML estimate under the AXDA model (131) boils down to solve the optimization problem

[TABLE]

If $\kappa_{\rho}$ stands for a Gaussian kernel, the problem (137) can be viewed as a quadratically penalized formulation of the initial problem $\min_{\boldsymbol{\theta}}\sum_{j=1}^{J}f_{j}(\boldsymbol{\theta})$, see Nocedal and Wright (2006, Section 17.1). As expected, the solution of (137) stands for an approximate solution w.r.t. the initial optimization problem. The associated algorithm is depicted in Algorithm 2. Regarding this algorithm, one can clearly see the benefit of using a variable splitting approach as in AXDA: the initial potential is split into $J$ individual potentials with no operator acting on $\boldsymbol{\theta}$. Therefore, the corresponding minimization problems are simpler (e.g., associated proximity operators may become available) and can be handled in parallel.

We eventually point out that the benefits of Algorithm 2 highlighted previously are also shared with the ADMM (Boyd et al. 2011). Instead of solving the approximate optimization problem (137) which encodes the splitting operation with a quadratic regularization term, the latter builds on Lagrangian duality in order to provide an exact solution to the initial minimization problem $\min f_{j}$.

Example 3 (continued)****.

Computing directly the MAP estimate with classical forward-backward algorithms (e.g., the fast iterative shrinkage-thresholding algorithm (FISTA) (Beck and Teboulle 2009)) associated to $\pi$ is challenging if the observations are distributed among multiple nodes. In addition, proximity operators associated to $f_{j}$ are generally not available in closed-form because of the operators $\mathbf{A}_{j}$. Algorithm 2 permits to tackle these issues by splitting the initial objective function. In particular, the minimization w.r.t. $\mathbf{z}_{j}$ for $j\in[J]$ corresponds to an unidimensional $l_{2}$ regularized logistic regression problem that can be dealt with gradient-based methods with few iterations. The minimization w.r.t. $\boldsymbol{\theta}$ boils down to the solving of a linear system where efficient solvers can be applied. Note that such a splitting scheme avoids the use of stochastic gradient methods.

5.4 Expectation-maximization for AXDA

An EM algorithm under the augmented density $\pi_{\rho}(\boldsymbol{\theta},\mathbf{z})$ will target the MAP or ML estimator, see Algorithm 3. If the expectations in the E-step cannot be evaluated, one can use a Monte Carlo approximation to approximate them (Wei and Tanner 1990). The benefits of using the augmented density $\pi_{\rho}$ instead of $\pi$ are threefolds. Firstly, as pointed out in Section 2 in the main paper, exact DA schemes based on $\pi$ cannot be derived in general cases and corresponding EM algorithms cannot be implemented. Instead, considering $\pi_{\rho}$ gives a quite systematic way of introducing latent variables in the original statistical model. Secondly, the expectations involved in the E-step of Algorithm 3 can be simpler to derive than the expectation under $\pi$. Indeed, the latter involves the whole potential $\sum_{j}f_{j}$ while the former involves regularized parts $f_{j}-\log\kappa_{\rho}$ of this potential separately. Finally, conditionally on $\boldsymbol{\theta}^{[t]}$, the random variables $\mathbf{z}_{j}$ are independent. Thus, each expectation involved in the E-step can be computed in parallel.

Example 3 (continued)****.

Again, following the work of Polson et al. (2013), if the potential $g$ of the prior distribution is quadratic or corresponds to a sparsity-promoting $\ell_{p}$-penalization ($0<p\leq 1$), a simple EM-algorithm can be derived as detailed by Scott and Sun (2013). However, although this EM algorithm can be generalized to an online version, it does not scale to distributed and high-dimensional problems. On the other hand, the E-step of Algorithm 3 can be processed in parallel by computing the $J$ expectations on individual nodes: thanks to the AXDA approach, the algorithm is therefore suited to distributed and high-dimensional scenarios.