Augmented Message Passing Stein Variational Gradient Descent

SVGD approximates an intractable unknown target distribution $p(\mathbf{x})$ where $\mathbf{x}=[x_{1},...,x_{D}]^{T}\in\mathcal{X}\subseteq\mathbb{R}^{D}$ with the best $q^{*}(\mathbf{x})$ by minimizing the Kullback-Leibler (KL) divergence (Liu and Wang (2016)). Here, $D$ is the dimension of the target distribution and

SVGD takes a group of particles $\left\{\mathbf{x}\right\}^{m}_{i=0}$ from the initial distribution $q_{0}(\mathbf{x})$, and after a series of smooth transforms, these particles finally converge to the target distribution $p(\mathbf{x})$. Each smooth transformation can be expressed by $\mathbf{T(x)}=\mathbf{x}+\epsilon\Phi(\mathbf{x})$, where $\epsilon$ is the step size and $\Phi:\mathcal{X}\to\mathbb{R}^{D}$ is the transformation direction. Here, $\mathcal{X}$ is the collection of the particles. Let $\mathcal{H}$ denote the set of positive definite kernels $k\left(\cdot,\cdot\right)$ in the reproducing kernel Hilbert space (RKHS) and denote $\mathcal{H}^{D}=\mathcal{H}\times\cdots\times\mathcal{H}$, where $\times$ is the Cartesian product. The steepest descent direction is obtained by minimizing the KL divergence,

where $q_{[T]}$ means particles take the distribution $q$ after $\mathbf{T(x)}$ mapping, $\mathcal{A}_{p}$ is the Stein operator given by $\mathcal{A}_{p}\phi(\mathbf{x})=\nabla_{\mathbf{x}}\log p(\mathbf{x})\phi(\mathbf{x})^{T}+\nabla_{\mathbf{x}}\phi(\mathbf{x})$. SVGD updates the particles $\left\{\mathbf{x}\right\}_{i=1}^{m}$ drawn from the initial distribution $q_{0}(\mathbf{x})$ by

The steepest descent direction is given by

The kernel function can be chosen as the RBF kernel $k(\mathbf{x},\mathbf{y})=\exp(-\|\mathbf{x}-\mathbf{y}\|_{2}^{2}/(2h))$ (Liu and Wang (2016)) or the IMQ kernel $k(\mathbf{x},\mathbf{y})=1/\sqrt{1+\|\mathbf{x}-\mathbf{y}\|_{2}^{2}/(2h)}$ (Gorham and Mackey (2017)). Equation (1) can be divided into two parts: the driving force term $\left[k\left(\boldsymbol{x}^{\prime},\boldsymbol{x}\right)\nabla_{\boldsymbol{x}^{\prime}}\log p\left(\boldsymbol{x}^{\prime}\right)\right]$ and the repulsive force term $\nabla_{\boldsymbol{x}^{\prime}}k\left(\boldsymbol{x}^{\prime},\boldsymbol{x}\right)$. It has been demonstrated that SVGD falls under the high-dimension curse (Zhuo et al. (2018); Liu (2017)). Related research shows that there exists a negative correlation between the problem dimension and the repulsive force of SVGD (Ba et al. (2021)). The influence of repulsive force is mainly related to the dimension of the target distribution.

MP-SVGD. The Message Passing SVGD (MP-SVGD) (Zhuo et al. (2018); Wang et al. (2018)) is proposed to reduce the dimension of the target distribution by identifying the Markov blanket for problems with a known graph structure. For dimension index $d$, its Markov blanket $\Gamma_{d}=\cup\{F:d\in F\}\backslash\{d\}$ contains neighborhood nodes of $d$ such that $p\left(x_{d}\mid\mathbf{x}_{\neg d}\right)=p\left(x_{d}\mid\mathbf{x}_{\Gamma_{d}}\right)$, where $F=\{1,...,D\}$ and $\neg d=\{1,...,D\}\backslash\{d\}$. However, it is necessary to rely on the sparse correlation of variables of the target distribution in order to obtain good results.

Since SVGD forms a dynamic system where particles interact with each other, one can no longer treat the converged particles as independent and identically distributed. Therefore, we need to use the mathematical tool “mixing”, cf. Bradley (2005) for more information.

Mixing. Let $\left\{\mathbf{x}_{m},m\geq 1\right\}$ be a sequence of random variables over some probability space $(\Omega,F,P)$, where $\sigma(\mathbf{x}_{i},i\leq j)$ denotes the $\sigma$-algebra. For any two $\sigma$-fields $\sigma\left(\mathbf{x}_{i},i\leq j\right),\ \sigma\left(\mathbf{x}_{i},i\geq j+k\right)$ with any $k\geq 1$, let

where

is the maximum taken over all finite partitions $\left(A_{i}\right)_{i\in I}$ and $\left(B_{i}\right)_{i\in J}$ of $\Omega$. A family $\left\{\mathbf{x}_{m},m\geq 1\right\}$ of random variables will be said to be absolutely regular (or $\beta-$mixing) if $\lim_{m\to\infty}\beta_{m}=0$, where the coefficients of absolute regularity $\beta_{m}$ are defined in Equation (2) (Banna et al. (2016)). These coefficients quantify the strength of dependence between the $\sigma$-algebra generated by $(\mathbf{x}_{i})_{1\leq i\leq k}$ and the one generated by $(\mathbf{x}_{i})_{i\geq k+m}$ for all $k\in\mathbb{N}^{*}$. $\beta_{m}$ tends to zero while $m$ goes to infinity implies that the $\sigma$-algebra generated by $(\mathbf{x}_{i})_{i\geq k+m}$ is less and less dependent on $\sigma\left(\mathbf{x}_{i},i\leq k\right)$ (Bradley (2005)).

A1 (Fixed points). Although the convergence of SVGD under finite particles is still an open problem, it can be found in many experimental studies that all particles converge to fixed points (Ba et al. (2021)). In this paper, we also assume that for SVGD with finite number of particles, these particles eventually converge to the target distribution and approach fixed points.

A2 ($m$-dependent of particles)

We assume that the fixed points of SVGD satisfy the $m$-dependent assumption. In Ba et al. (2021), only weak correlation between SVGD particles is reported. We perform numerical verification in Appendix B and leave the rigorous proof for our future work. Under Assumptions A1-A2 and the zero-mean of the target distribution, we analyze the variance collapse of SVGD quantitatively in the following part.

We first give the upper bound of the concentrated particles. The main tool used here is the Jensen gap (Gao et al. (2019)).

We leave this proof in Appendix A. For most of the sampling-based inference methods, such as MCMC, VI, et al., samples spread across the whole sample space but with extremely small probability for some samples. However, Proposition 1 shows that SVGD with finite particles are confined to a certain range, although this range may expand as the number of particles increases. With the RBF kernel, this upper bound is related to the trace of the covariance of the target distribution. Under the IMQ kernel $k(\mathbf{x},\mathbf{y})=1/\sqrt{1+\|\mathbf{x}-\mathbf{y}\|_{2}^{2}/(2h)}$ conditions, this range will be further reduced (Gorham and Mackey (2017)).

For independent and identically distributed (i.i.d.) samples, the variance can be estimated using the Bernstein inequality or the Lieb inequality. However, these random matrix results typically require the particles to be i.i.d, which is no longer satisfied by SVGD due to its interacting update. To analyze the variance of particles from SVGD, we assume these particles are $m$-dependent. For a sequence of random variables, the $m$-dependent assumption implies that they satisfy the $\beta$-mixing condition (Bradley (2005)). We obtain the following results based on the Bernstein inequality of dependent random matrices (Banna et al. (2016)).

Proposition 2 shows that the main factors that affect the upper bound of the variance error include the inter-particle correlation, the number of particles, the dimension of the target distribution, and the trace of its covariance matrix. According to Ba et al. (2021), it can be considered that $tr(\Sigma_{m})\leq tr(\Sigma)$ should be true for SVGD, therefore the constant $c$ in Proposition 1 can be replaced by $\frac{2}{h}e^{\frac{-c_{0}tr(\Sigma)}{h}}$.

The update direction $\Delta$ of SVGD is given by

The log derivative term in the update rule $\mathbb{E}_{\boldsymbol{x}^{\prime}\sim q}\left[k(\boldsymbol{x}^{\prime},\boldsymbol{x})\nabla_{\boldsymbol{x}^{\prime}}\log p(\boldsymbol{x}^{\prime})\right]$ corresponds to the driving force that guides particles toward the high-likelihood region. The (Stein) score function $\nabla_{\boldsymbol{x}^{\prime}}\log p\left(\boldsymbol{x}^{\prime}\right)$ is a vector field describing the target distribution. $\nabla_{\boldsymbol{x}^{\prime}}k\left(\boldsymbol{x}^{\prime},\boldsymbol{x}\right)$ provides a repulsive force to prevent particles from aggregating. However, as the dimension of the target distribution increases, this repulsive force gradually decreases (Zhuo et al. (2018)). This causes SVGD to fall under the curse of high dimensions. How to effectively reduce the dimension has become the guiding ideology to make SVGD overcome the curse of dimensionality.

Our concern is the problem with the continuous graphical model, i.e., the target distribution that satisfies the following form $p(\mathbf{x})\propto\exp\left[\sum_{s\in\mathcal{S}}\psi\left(\mathbf{x}_{s}\right)\right]$ where $\mathcal{S}$ is the family of index sets $s\subseteq\{1,\ldots,D\}$ that specifies the Markov structure. For any index $i$, its Markov blanket is represented by $\mathcal{N}_{i}:=\cup\{s:s\subset\mathcal{S},i\in s\}\backslash\{i\}$. According to Wang et al. (2018), one can transform the global kernel function into a local kernel function and this local kernel function is just related to $\mathcal{C}_{i}:=\mathcal{N}_{i}\cup\{i\}$ for any $i$. Under this transformation, the dimension is reduced from $D$ to $|\mathcal{C}_{i}|$, where $|\mathcal{C}_{i}|$ is the size of the set $\mathcal{C}_{i}$. Then, SVGD becomes

where $s(x_{i}^{\ell})=\nabla_{x_{i}^{\ell}}\log p(x_{i}^{\ell})$. Such method is known as the message passing SVGD (MP-SVGD) (Zhuo et al. (2018); Liu (2017)). However, MP-SVGD needs to know the graph structure of the target distribution in advance to determine the Markov blanket and the sparsity of this graph is needed in order to achieve dimension reduction.

Inspired by MP-SVGD, we propose the augmented MP-SVGD which is suitable for more complex graph structures. We keep the previous symbols but redefine them for clarity. We assume $p(\mathbf{x})$ can be factorized as $p(\mathbf{x})\propto\prod_{F\in\mathcal{F}}\psi_{F}\left(\mathbf{x}_{F}\right)$ where $F\subset\{1,\ldots,D\}$ is the index set. Partition $\Gamma_{d}\subset\{1,\ldots,D\}\backslash\{d\}$ and $\mathbf{S_{d}}\subset\{1,\ldots,D\}\backslash\{d\}$ such that $\Gamma_{d}\cap\mathbf{S_{d}}=\emptyset$ and $\Gamma_{d}\cup\mathbf{S_{d}}=\{1,\ldots,D\}\backslash\{d\}$. Let $\mathbf{x}_{\mathbf{S_{d}}}=[x_{i},\ ...\ ,\ x_{j}],\ i,...,j\in\mathbf{S_{d}}$. Similarly, $\mathbf{x}_{\Gamma_{d}}=[x_{i},\ ...\ ,\ x_{j}],\ i,...,j\in\Gamma_{d}$.

Consider the probability graph model illustrated by Figure 1, $p(\mathbf{x})$ can be represented as $p(\mathbf{x})=p(x_{d}|\mathbf{x}_{\Gamma_{d}},\mathbf{x}_{\mathbf{S_{d}}})p(\mathbf{x}_{\Gamma_{d}}|\mathbf{x}_{\mathbf{S_{d}}})p(\mathbf{x}_{\mathbf{S_{d}}})$. Our method relies on the key observation of Zhuo et al. (2018) that

To minimize $\mathrm{KL}(q\|p)$, we adopt a two-stage procedure so that $\text{prKL}(q,p,d)$ and $\text{seKL}(q,p,d)$ are optimized alternatively. At the first stage, $\text{prKL}(q,p,d)$ is further decomposed into

We can fix $\mathbf{x}_{\mathbf{S_{d}}}$ to apply the local kernel function to the second part of Equation (4) to minimize $\text{prKL}(q,p,d)$,

This optimization procedure is given by Proposition 3 and we leave the proof in the appendix.

At the second stage, $\mathbf{x}_{\Gamma_{d}}$ and $\mathbf{x}_{\mathbf{S_{d}}}$ are fixed while only $\mathbf{x}_{d}$ is updated. We can further decompose $\mathrm{seKL}(q,p,d)$ into three parts via the convexity of the KL divergence,

where $C=\mathrm{KL}\left[\frac{1}{q(\mathbf{x}_{\neg d})}\|\frac{1}{p(\mathbf{x}_{\neg d})}\right]$ is a positive constant due to the fact that $\mathbf{x}_{\neg d}=\mathbf{x}_{\Gamma_{d}}\cup\mathbf{x}_{\mathbf{S_{d}}}$ is fixed. Therefore,

The optimization procedure is described by Proposition 4 and we also leave the proof in the appendix.

For $d\in\left\{1,2,\dots,D\right\}$, $x_{d}$ is updated through the above two-stage procedure, which reduces the dimension of the original problem from $D$ to the size of $min(|\mathbf{x}_{\mathbf{S_{d}}}|,\ |\mathbf{x}_{\Gamma_{d}}|)$. In this way, we are able to solve the problem of inferring more complex target distributions that traditional MP-SVGD failed. This comparison is illustrated in Experiment 2. Moreover, we do not need to know the real probability graph structure of the target distribution in advance.

The key start of our AUMP-SVGD is to choose $\Gamma_{d}$ or $\mathbf{S_{d}}$. This problem can be formulated in the following form, for $m$ particles with each $x\in\mathbb{R}^{D}$, denote the ensemble matrix of these particles by $X$, i.e., $X\in\mathbb{R}^{m\times D}$. Let $r=|\mathbf{x}_{\Gamma_{d}}|$ and $\mathcal{M}(\mathrm{r},m\times D)$ be the set of sub-matrices $\Gamma_{x}\in\mathbb{R}^{m\times r}$ of $X$. Determining the set $\Gamma_{d}$ corresponds to select the sub-matrix $X_{r}$ from the set $\mathcal{M}(\mathrm{r},m\times D)$ to minimize $\mathbf{E}\|\Sigma_{m}-\Sigma_{x}\|$, where $\Sigma_{m}$ is the empirical covariance matrix of particles and $\Sigma_{x}$ is the empirical covariance matrix of sub-ensembles. The upper bound of $\mathbf{E}\|\Sigma_{m}-\Sigma_{x}\|$ is already given by Proposition 2 and it is related to $tr(\Sigma_{x})$. Therefore, we choose the sub-matrix with the smallest $tr(\Sigma_{x})$ to ensure that $\mathbb{E}\|\Sigma_{m}-\Sigma_{x}\|$ is minimized. This corresponds to select $X_{r}\in\mathbb{R}^{m\times r}$ from $X$ to get minimal $\operatorname{tr}\left(X_{r}X_{r}^{T}\right)$. Since $\operatorname{tr}\left(X_{r}X_{r}^{T}\right)=\operatorname{tr}\left(X_{r}^{T}X_{r}\right)$ , we just need to reorder each column of $X$ by the 2-norm and choose $r$ columns with smallest 2-norm. The computational complexity for this part is $\mathcal{O}\left(Dm^{2}+D\log D\right)$. Finally, we give the complete form of our AUMP-SVGD by Algorithm 1.

Multivariate Gaussian. The first example is a $D$-dimensional multivariate Gaussian $p(\mathbf{x})=\mathcal{N}\left(0,I_{D}\right)$. For each method, 100 particles are initialized from $q_{0}(\mathbf{x})=\mathcal{N}\left(10,I_{D}\right)$.

Spaceship Mixture. The target in the second experiment is a $D$-dimensional mixture of two correlated Gaussian distributions $p(\mathbf{x})=0.5\mathcal{N}\left(x;\ \mu_{1},\ \Sigma_{1}\right)+0.5\mathcal{N}\left(x;\ \mu_{2},\ \Sigma_{2}\right)$. The mean $\mu_{1}$, $\mu_{2}$ of each Gaussian have components equal to 1 in the first two coordinates and 0 otherwise. The covariance matrix admits a correlated block diagonal structure. The mixture hence manifests as a “spaceship” density margin in the first two dimensions (see Figure 2).

It can be seen from Figure 2 that for the high-dimensional inference, particles from SVGD aggregate, which leads to a high-dimensional curse (Zhuo et al. (2018); Liu and Wang (2016)). However, AUMP-SVGD can estimate the true probability distribution well in these high-dimensional situations. We calculate the energy distance and the mean-square error (MSE) $\mathbb{E}\|\Sigma_{m}-\Sigma\|_{2}$ between the samples from the inference algorithm and the real samples. The energy distance is given by $D^{2}(F,G)=2\mathbb{E}\|X-Y\|-\mathbb{E}\left\|X-X^{\prime}\right\|-\mathbb{E}\left\|Y-Y^{\prime}\right\|$, where $F$ and $G$ are the cumulative distribution function (CDF) of $X$ and $Y$, respectively. $X^{\prime}$ and $Y^{\prime}$ denote an independent and identically distributed (i.i.d.) copy of $X$ and $Y$ (Rizzo and Székely (2016)). 10 experiments are performed and the averaged results are given in Figure 3.

Figure 3 demonstrates the gradual expansion of the error difference of SVGD as the dimension increases. For the sparse problem, when the graph structure is already known, MP-SVGD achieves comparable outcomes with S-SVGD and PSVGD-2. In the above example, PSVGD achieves its best results when the problem dimension is reduced to 2. The AUMP-SVGD with a set size of $\Gamma_{d}$ of 1 or 3 outperforms other methods. In Experiment 1, we demonstrate that AUMP-SVGD yields a variance estimate equivalent to that of MP-SVGD under the simplified graph structure. In Experiment 2, as the correlation between different dimensions of the target becomes strong or the density of the graph increases, our algorithm performs better than other SVGD variants. Furthermore, our approach exhibits superior variance estimation compared with SVGD, MP-SVGD, S-SVGD, and PSVGD. In practice, the covariance matrix of the target distribution may not be sparse, making it challenging to capture this structure. As a result, the effectiveness of MP-SVGD is significantly limited. However, AUMP-SVGD can still attain stable and superior results.

Non-sparse experiment We set the dimension of the target distribution to 50 and systematically transfer the sparse covariance matrix to a non-sparse one by augmenting the correlations between different dimensions along the main diagonal. The results are presented in Figure  4.

As the density of the probability map increases, the discrepancy between MP-SVGD and the actual target distribution gradually magnifies. Eventually, MP-SVGD succumbs to the curse of dimensionality. This phenomenon arises due to the fact that the dispersion force of particles in MP-SVGD primarily depends on the size of the Markov blanket. With increasing density, MP-SVGD encounters the same challenges associated with high-dimensional scenarios as observed in SVGD. However, as illustrated in Figure 4, regardless of the density of the target distribution’s graph structure, AUMP-SVGD remains unaffected by variance collapse. This is attributed to the repulsion force of particles in AUMP-SVGD being influenced by the artificially chosen Markov blanket size, underscoring the superiority of our algorithm over MP-SVGD.

The next example is a benchmark that is often used to test inference methods in high dimensions (Detommaso et al. (2018); Chen and Ghattas (2020); Liu et al. (2022)). We consider a stochastic process $u:[0,1]\rightarrow\mathbb{R}$ governed by

where $t\in(0,1]$, the forcing term $x=\left(x_{t}\right)_{t\geq 0}$ follows a Brownian motion so that $x\sim\mathcal{N}(0,B)$ with $B\left(t,t^{\prime}\right)=\min\left(t,t^{\prime}\right)$. The noisy data $\mathbf{z}=\left(z_{t_{1}},\ldots,z_{t_{50}}\right)^{T}\in\mathbb{R}^{50}$ at 50 equi-spaced time points with $t_{i}=0.02i$, where $z_{t_{i}}=u_{t_{i}}+\epsilon$ for $\epsilon\sim\mathcal{N}\left(0,\sigma^{2}\right)$ with $\sigma=0.1$. The objective is to use $z$ to infer the forcing term $x$ and thus the state of the solution $u$. The results are given in Figure 5 where the shadow interval depicts the mean plus/minus the standard deviation.

As depicted by Figure 5, it is evident that in the case of 50 dimensions, both SVGD and S-SVGD exhibit certain deviations from the ground truth while S-SVGD exhbits excessively large variances in numerous tests. Consequently, SVGD and S-SVGD prove inadequate in effectively addressing the conditional diffusion model. Conversely, our AUMP-SVGD demonstrates satisfactory performance with the size of $\Gamma_{d}$ of 5 and 10.

We investigate a Bayesian logistic regression model from Liu and Wang (2016) applied to the Covertype dataset from Asuncion and Newman (2007). We use 70% data for training and 30% for testing. We compare AUMP-SVGD with SVGD, S-SVGD, PSVGD, MP-SVGD, and Hamiltonian Monte Carlo(HMC) and the number of generated samples ranges from 100 to 500. Each experiment uses ten different random seeds and the error of each value does not exceed 0.01. We verify the impact of different sampling methods on the prediction accuracy and the results are given by Table 1.

Table 1 demonstrates that our AUMP-SVGD has a prediction accuracy rate similar to that of PSVGD, and as the number of particles increases, our algorithm is more accurate than other algorithms.

According to the fixed points assumption (A1),

Then,

In Equation (5), $\nabla{\mathbf{x}}\log p\left(\mathbf{x}\right)$ and $k\left(\mathbf{x}_{i},\mathbf{x}\right)$ are independent and,

Associated with the Jensen gap (Gao et al. [2019]), we get

which in turn gives

since $\left|\mathbf{x}_{i}-\mathbb{E}\mathbf{x}\right|_{2}^{2}=\mathrm{tr}\left((\mathbf{x}_{i}-\mathbb{E}\mathbf{x})(\mathbf{x}_{i}-\mathbb{E}\mathbf{x})^{T}\right)\leq c_{0}\mathrm{tr}(\Sigma)$, where $c_{0}$ is a positive number.

According to Proposition 1:

Apply the expectation version of the Bernstein inequality (Banna et al. [2016]) for the sum of mean zero random matrices $\mathbf{x}_{i}\mathbf{x}_{i}^{T}-\Sigma$ and we obtain,

and $M$ is any number chosen such that,

To bound $M$ is simple:

which completes the proof.

Note that

Now we derive the optimality condition for Equation (3). Note that,

Following the proof of Theorem 3.1 by Liu and Wang [2016], we have

According to Liu and Wang [2016], we can show that the optimal solution is given by $\frac{\phi_{\Gamma_{d}}^{*}}{\|\phi_{\Gamma_{d}}^{*}\|_{\mathcal{H}_{\Gamma_{d}}}}$, where