# Population Monte Carlo with Normalizing Flow

Soumyasundar Pal, Antonios Valkanas, and Mark Coates Soumyasundar Pal is with
Huawei Noah’s Arc Lab, Montréal, Québec, Canada. Antonios Valkanas and Mark
Coates are with the Dept. of Electrical and Computer Engineering, McGill
University, Montréal, Québec, Canada. Corresponding author: Soumyasundar Pal.
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Code to reproduce our experiments is available at
https://github.com/soumyapal91/AIS_flow

###### Abstract

Adaptive importance sampling (AIS) methods provide a useful alternative to
Markov Chain Monte Carlo (MCMC) algorithms for performing inference of
intractable distributions. Population Monte Carlo (PMC) algorithms constitute
a family of AIS approaches which adapt the proposal distributions iteratively
to improve the approximation of the target distribution. Recent work in this
area primarily focuses on ameliorating the proposal adaptation procedure for
high-dimensional applications. However, most of the AIS algorithms use simple
proposal distributions for sampling, which might be inadequate in exploring
target distributions with intricate geometries. In this work, we construct
expressive proposal distributions in the AIS framework using normalizing flow,
an appealing approach for modeling complex distributions. We use an iterative
parameter update rule to enhance the approximation of the target distribution.
Numerical experiments show that in high-dimensional settings, the proposed
algorithm offers significantly improved performance compared to the existing
techniques.

###### Index Terms:

Bayesian inference, Importance sampling, Monte Carlo methods, Normalizing
flow, Population Monte Carlo

##  I Introduction

Monte Carlo (MC) methods [1, 2] have numerous applications in statistical
signal processing, communications, and machine learning. For example, in the
Bayesian setting, the statistical inference of an unknown parameter is carried
out by computing the posterior distribution conditioned on the observed data
[3]. However, in many practical problems of interest, the posterior
distribution does not have an analytic form, and approximation via MC methods
becomes necessary.

Markov Chain Monte Carlo (MCMC) [4] algorithms have been a popular approach
for addressing posterior inference problems. They involve constructing a
Markov process whose stationary distribution is the same as the target
posterior. Several efficient MCMC algorithms [5, 6] have been developed which
obtain state-of-the-art results in many high dimensional scenarios. Adaptive
versions of MCMC algorithms [7, 8] are often used for accelerating the mixing
of the Markov chain. Despite the effectiveness of the MCMC techniques, there
are several disadvantages of these algorithms. In higher dimensions, they
often result in very high computational costs and/or slow mixing. In this
setting, an importance sampling (IS) based methodology can provide an
alternative approach.

Importance sampling is a class of Monte Carlo methods, where the target
distribution is approximated by the weighted samples from another proposal
distribution [1]. In high-dimensional applications, the design of an effective
proposal distribution is crucial since a greater mismatch between the target
and the proposal distributions leads to a high variance of the IS estimator.
Several proposal distributions can be used for efficient exploration of state
space of complicated target distributions; this approach is referred to as
multiple IS (MIS) [9].

Although the use of multiple proposals can be beneficial in some cases,
adaptation of the proposal distributions [10] can be more advantageous for
inference of higher dimensional target distributions. In adaptive IS (AIS)
schemes, iterative adaptation of the proposal distributions is performed to
improve the approximation of the target distribution. Generally, an iteration
of an AIS scheme is composed of three different steps [10]: a) sampling from a
set of proposals (sampling), b) computing importance weights for the samples
(weighting), and c) updating (adapting) the parameters of the proposals for
the next iteration.

There are numerous variants of the AIS algorithms, and the differences among
them arise from consideration of how the three aforementioned steps are
implemented. For example, the standard Population Monte Carlo (PMC) algorithm
[11] draws only one sample from each proposal at each iteration, whereas the
more recent schemes [12, 13, 14, 15] use multiple samples. The weighting of
the samples is not uniform in the literature either. [11] uses individual
proposals in the denominator of the importance weight computation step,
whereas [12, 13] use a deterministic mixture (DM-PMC) to compute the weights,
as this leads to lower variance estimates (albeit at the cost of additional
complexity). A detailed comparison of related algorithms in terms of their
similarities and differences is provided in [10].

From a practical point of view, the most crucial step in an AIS algorithm is
proposal adaptation. Earlier works like [11] and its variants [12, 16] use
resampling of the samples to adapt the location parameters of the proposal
distributions, which is inefficient in higher dimensions. More effective
resampling strategies are considered in [12]. Intuitively the goal of the
adaptation step is to move the proposals in a way such that they are placed in
the high-density regions of the target distribution. One way to achieve this
is to use the gradient of the target to adapt the location parameters of the
proposals [14, 13, 17, 18], although other strategies like minimizing KL
divergence [19], employing MCMC steps [20, 21, 22], performing moment matching
[23] or using stochastic optimization [15] can also be used. While these
techniques offer improvement in performance, the application of AIS algorithms
in high dimensions remains challenging.

A major limitation of the existing AIS algorithms is their use of simple
proposal distributions. Most of these methods use isotropic Gaussian
distributions as proposals for ease of sampling and density evaluation.
Despite the presence of an effective adaptation rule, such a choice might be
sub-optimal for the traversal of high-dimensional target distributions with
multi-modality, complex correlation structures, and highly non-uniform scales
across different dimensions. Some existing methods [13, 17] address this issue
to some extent by updating the scale parameters of the proposals using the
local curvature of the target, along with the locations. However, such
strategies scale poorly to high dimensions since they require the computation
of the inverse of the Hessian matrix of the logarithm of the target
distribution.

Normalizing Flows (NF) [24, 25, 26] represent a class of expressive
probability distributions, which enables efficient and exact sampling as well
as density evaluation. Various flow models have widespread applications in
generative modelling [24, 27], density estimation [28, 25, 26], and modelling
conditional distributions [29, 30, 31]. Several approaches [32, 33, 34, 35,
36] combine NF models with annealed importance sampling and MCMC for efficient
sampling from intractable distributions.

In this paper, we present a novel PMC algorithm, which uses normalizing flow
to parameterize effective proposal distributions. The incorporation of
normalizing flow enhances the sampling capability of the PMC algorithm
remarkably which results in impressive performance in high dimensions with
favorable complexity.

The paper is organized as follows. Sec. II states the inference task that we
address. Sec. III reviews the AIS framework and normalizing flow. Sec. IV
introduces the proposed PMC algorithm using normalizing flow and Sec. V
presents and discusses the results of numerical simulation experiments.
Concluding remarks are summarized in Sec. VI.

##  II Problem Statement

Although Monte Carlo methods can be applied for sampling from any probability
distribution, in practice they are often used to target inference of posterior
distributions in a Bayesian setting. Let us consider the problem of
statistical inference of an unknown parameter $\mathbf{x}\in\mathbb{R}^{d}$ in
the presence of some observed data $\mathbf{z}\in\mathbb{R}^{d_{z}}$. The
prior distribution of $\mathbf{x}$ is denoted by $p(\mathbf{x})$ and the
observed data $\mathbf{z}$ is related to parameter $\mathbf{x}$ through a
probabilistic generative model, specified by the likelihood function
$p(\mathbf{z}|\mathbf{x})$. The posterior distribution of $\mathbf{x}$ is
computed using Bayes’ rule:

$\displaystyle\tilde{\pi}(\mathbf{x}|\mathbf{z})=\frac{p(\mathbf{x})p(\mathbf{z}|\mathbf{x})}{p(\mathbf{z})}\propto\pi(\mathbf{x}|\mathbf{z})\,.$
(1)

Here, $\pi(\mathbf{x}|\mathbf{z}){=}p(\mathbf{x})p(\mathbf{z}|\mathbf{x})$ is
the unnormalized posterior distribution of $\mathbf{x}$ and
$p(\mathbf{z}){=}\int p(\mathbf{x})p(\mathbf{z}|\mathbf{x})d\mathbf{x}>0$ is a
non-negative function of $\mathbf{z}$, which is typically unknown. So, in
general, the posterior distribution $\tilde{\pi}(\mathbf{x}|\mathbf{z})$ is
known only up to a normalizing constant. In many problems, Bayesian decision-
making involves computing the expectation of a function $f$ w.r.t. the
posterior distribution, given as

$\displaystyle\mathbf{E}_{\tilde{\pi}}[f]=\int
f(\mathbf{x})\tilde{\pi}(\mathbf{x}|\mathbf{z})d\mathbf{x}=\frac{1}{p(\mathbf{z})}\int
f(\mathbf{x})\pi(\mathbf{x}|\mathbf{z})d\mathbf{x}\,.$ (2)

Here $f$ is an integrable function with respect to
$\tilde{\pi}(\mathbf{x}|\mathbf{z})$. However, in most cases the integral in
Eq. (2) cannot be computed analytically because $\pi$ and/or $f$ are
complicated and/or $p(\mathbf{z})$ is unknown. In these situations, our goal
is to approximate the posterior distribution by its samples. From here onward,
we denote the target (posterior) distribution as $\pi(\mathbf{x})$ to simplify
notation.

##  III Background

###  III-A Generic AIS Algorithm

Let $\pi(\mathbf{x})$ be the target distribution with a possibly unknown
normalizing constant. We assume that there are $N$ proposal distributions,
used at each of $J$ total iterations of the algorithm. Let
$q_{n}(\cdot;\theta_{n}^{j})$ be the $n$-th proposal at the $j$-th iteration,
where $\theta_{n}^{j}$ typically consists of a location parameter
$\mu_{n}^{j}$ and a scale parameter $\Sigma_{n}^{j}$. In most AIS algorithms,
Gaussian proposals are used for sampling, which implies that $\mu_{n}^{j}$ and
$\Sigma_{n}^{j}$ are the mean vector and the covariance matrix of $q_{n}$
respectively. In each iteration, we draw $K$ samples
$\\{\mathbf{x}_{n,j}^{(k)}\\}_{k{=}1}^{K}{\sim}q_{n}(\cdot;\theta_{n}^{j})$
from each of the proposals. The generic AIS [10] method is summarized in
Algorithm 1.

Algorithm 1 Generic AIS algorithm

1: Input: $\pi(\mathbf{x})$, $N,J,K$, $\\{\theta_{n}^{1}\\}_{n=1}^{N}$

2: Output:
$\\{\omega_{n,j}^{(k)},\mathbf{x}_{n,j}^{(k)}\\},n{=}1:N,k{=}1:K,j{=}1:J$

3: for $j{=}1:J$ do

4: Sampling: draw $K$ samples from each of $N$ proposals,
$\\{x_{n,j}^{(k)}\\}_{k{=}1}^{K}{\sim}q_{n}(\cdot;\theta_{n}^{j})$

5: Weighting: compute importance weight $\omega_{n,j}^{(k)}$ of
$\mathbf{x}_{n,j}^{(k)}$

6: Adaptation: update the proposal parameters:
$\\{\theta_{n}^{j}\\}_{n{=}1}^{N}\rightarrow\\{\theta_{n}^{j+1}\\}_{n{=}1}^{N}$

7: end for

###  III-B Normalizing Flow

Normalizing Flow (NF) [24, 25, 26] approaches are capable of constructing
complex distributions by transforming a simple probability density (which is
easy to sample from and easy to evaluate) via a sequence of invertible and
differentiable mappings. These invertible transformations are designed in such
a way that the forward mappings, their Jacobian matrices (and their
determinants), and the inverse mappings can be computed efficiently. The
change of variable formula is used to evaluate the density of the transformed
random variable.

Suppose $\mathbf{x}^{\prime}{\sim}q(\mathbf{x}^{\prime})$ is mapped to
$\mathbf{x}{=}\mathcal{T}(\mathbf{x}^{\prime})$ via an invertible function
$\mathcal{T}^{\phi}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$. Here, $\phi$
denotes the parameters of $\mathcal{T}^{\phi}$. Different NF models use
different neural architectures to parameterize $\mathcal{T}^{\phi}$. The
density of $\mathbf{x}$ is computed in closed form as:

$\displaystyle
q(\mathbf{x})=q(\mathbf{x}^{\prime})\lvert\text{det}(\mathcal{J}_{\mathcal{T}^{\phi}}(\mathbf{x}^{\prime}))\rvert^{-1}\,,$
(3)

where $\mathcal{J}_{\mathcal{T}^{\phi}}(\cdot)\in\mathbb{R}^{d\times d}$ is
the Jacobian of the transformation and $\text{det}(\cdot)$ denotes the
determinant of a matrix. We design $\mathcal{T}^{\phi}(\cdot)$ using the real-
NVP [25] architecture, composed of two coupling layers in an alternating
pattern as follows:

$\displaystyle x_{\lceil\frac{d}{2}\rceil+1:d}$
$\displaystyle{=}x^{\prime}_{\lceil\frac{d}{2}\rceil+1:d}\odot\exp\Big{(}s_{1}^{\phi_{1}}\big{(}x^{\prime}_{1:\lceil\frac{d}{2}\rceil}\big{)}\Big{)}+t_{1}^{\phi_{2}}\big{(}x^{\prime}_{1:\lceil\frac{d}{2}\rceil}\big{)}\,,$
$\displaystyle x_{1:\lceil\frac{d}{2}\rceil}$
$\displaystyle=x^{\prime}_{1:\lceil\frac{d}{2}\rceil}\odot\exp\Big{(}s_{2}^{\phi_{3}}\big{(}x_{\lceil\frac{d}{2}\rceil+1:d}\big{)}\Big{)}+t_{2}^{\phi_{4}}\big{(}x_{\lceil\frac{d}{2}\rceil+1:d})\,.$
(4)

Here, $\odot$ denotes Hadamard product, and
$\big{(}s_{1}^{\phi_{1}}(\cdot),t_{1}^{\phi_{2}}(\cdot)\big{)}$ and
$\big{(}s_{2}^{\phi_{3}}(\cdot),t_{2}^{\phi_{4}}(\cdot)\big{)}$ are the scale
and translation functions of the first and second coupling layers
respectively. Each of these functions is parameterized by a separate _feed-
forward network_ (FFN) and $\phi=\\{\phi_{1:4}\\}$ denotes the set of all
weight matrices and bias vectors of these FFNs.

##  IV PMC with NF

###  IV-A Proposal Step

We design the $n$-th proposal distribution $q_{n}$ by constructing a
normalizing flow $\mathcal{T}_{n}^{\phi_{n}}(\cdot)$ with a multivariate
normal base distribution with mean $\mu_{n}$ and covariance $\Sigma_{n}$. In
other words, sampling from $q_{n}$ is carried out by first sampling
$\mathbf{x}^{\prime}{\sim}\mathcal{N}(\mu_{n},\Sigma_{n})$ and then computing
$\mathbf{x}{=}\mathcal{T}_{n}^{\phi_{n}}(\mathbf{x}^{\prime})$. The proposal
is evaluated as follows:

$\displaystyle
q_{n}(\mathbf{x};\theta_{n})=\mathcal{N}(\mathbf{x}^{\prime};\mu_{n},\Sigma_{n})\lvert\text{det}(\mathcal{J}_{\mathcal{T}_{n}^{\phi_{n}}}(\mathbf{x}^{\prime}))\rvert^{-1}\,,$
(5)

The parameters of the $n$-th proposal density are defined as
$\theta_{n}{=}\\{\mu_{n},\Sigma_{n},\phi_{n}\\}$ and
$\theta{=}\\{\theta_{n}\\}_{n{=}1}^{N}$ denotes all proposal parameters.
Utilization of normalizing flow allows $q_{n}$ to be non-Gaussian and
expressive.

###  IV-B Weighting Step

Let $\theta_{n}^{j}$ be the value of $\theta_{n}$ at the $j$-th iteration. We
use the deterministic mixture (DM) weighting [16] scheme, as it reduces the
variance of the unnormalized AIS estimator, i.e., we compute:

$\displaystyle\omega_{n,j}^{(k)}\propto\frac{\pi(\mathbf{x}_{n,j}^{(k)})}{\frac{1}{N}\displaystyle{\sum_{l=1}^{N}}q_{l}(\mathbf{x}_{n,j}^{(k)};\theta_{l}^{j})}\,.$
(6)

###  IV-C Adaptation Step

We adapt the proposal parameters $\theta$ by performing a mini-batch gradient
descent step towards minimizing the KL divergence between the mixture proposal
distribution
$q(\mathbf{x};\theta)=\frac{1}{N}\sum_{n{=}1}^{N}q_{n}(\mathbf{x};\theta_{n})$
and the target density $\pi(\mathbf{x})$ w.r.t. $\theta$. Since this KL
divergence cannot be evaluated in closed form, at the $j$-th iteration, we
compute an unbiased Monte Carlo estimate of
$\ell(\theta^{j})=\text{KL}\big{(}q(\mathbf{x};\theta^{j})\lvert\rvert\pi(\mathbf{x})\big{)}$
as follows:

$\displaystyle\widehat{\ell(\theta^{j})}$
$\displaystyle=\frac{-1}{NK}\sum_{n=1}^{N}\sum_{k=1}^{K}\log\omega_{n,j}^{(k)}\,.$
(7)

The update rule for $\theta$ is:

$\displaystyle\theta^{j+1}=\theta^{j}-\epsilon^{j}\nabla_{\theta}\widehat{\ell(\theta)}\rvert_{\theta=\theta^{j}}\,,$
(8)

where $\epsilon^{j}$ denotes the learning rate at the $j$-th iteration. In
practice, the stochastic gradient used in Eq. (8) is computed using automatic
differentiation.

##  V Numerical Experiments

We conduct simulation experiments for two scenarios. The goal of the first
experiment is to draw samples from a high-dimensional Gaussian mixture target
distribution, whereas the second experiment considers a high-dimensional
Bayesian logistic regression problem. The proposed Normalizing Flow (NF)-PMC
algorithm is compared with PMC [11], Global Resampling (GR)-PMC [12], Local
Resampling (LR)-PMC [12], Scaled Langevin (SL)-PMC [13], Optimized (O)-PMC
[17], GRadient-based Adaptive Multiple Importance Sampler [18], Hamiltonian
Adaptive Importance Sampler (HAIS) [21], and Variational Adaptive Population
Importance Sampler (VAPIS) [15] algorithms.

GR-PMC and LR-PMC algorithms are two variants of the DM-PMC [16] technique.
Langevin dynamics [37] is used for designing the adaptation rule of SL-PMC,
whereas O-PMC and HAIS algorithms utilize Newton’s optimization and
Hamiltonian Monte Carlo [38, 5] respectively. GRAMIS adapts the location
parameters of the proposals using a physics-inspired repulsion term, which
facilitates coordinated exploration of the target density. Location parameters
of the deterministic mixture proposal distribution are adapted by minimizing
the per-sample variance of the IS estimator of the normalizing constant via
stochastic gradient variational inference [39] for the VAPIS algorithm.

For each algorithm, except PMC, we use $J{=}50$ iterations, in each of which
$K{=}10$ samples are randomly drawn from each of the $N{=}100$ proposal
distributions. Since the PMC algorithm uses $K{=}1$ sample from each proposal
at each iteration, we use $J{=}50\times 10{=}500$ iterations for PMC, so that
all of the algorithms generate the same number of samples in total. For each
setting, we conduct 100 Monte Carlo trials.

Deterministic mixture (DM) [16] weights are computed for all algorithms except
PMC, which uses the standard weighting. As in the original papers, we follow
the same backtracking strategy for tuning the step size for SL-PMC and O-PMC.
The burn-in period in HAIS is set to 10, and we use 50 leapfrog steps. The
RMSprop111http://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf
optimizer with a decaying learning rate schedule is employed for VAPIS and NF-
PMC because of its faster convergence compared to the SGD method. The scale
and translation functions in each coupling layer are parameterized by FFNs
with 2 hidden layers with 8 neurons and $tanh$ activation for the
implementation of NF-PMC. Additionally, we treat $\Sigma_{n}$ as
hyperparameters and share the NF parameters $\\{\phi_{n}\\}_{n{=}1}^{N}$
across all proposals to reduce parameter complexity.

TABLE I: Mean and standard deviation of MSE in the mean estimation of the Gaussian mixture along with avg. runtime, the lowest MSE in each row is shown in bold. Alg. | PMC | GR-PMC | LR-PMC | SL-PMC | O-PMC | GRAMIS | HAIS | VAPIS | NF-PMC
---|---|---|---|---|---|---|---|---|---
MSE |
$\sigma{=}1$ |  36.24$\pm$5.21 |  32.49$\pm$4.32 |  38.65$\pm$4.62 |  18.21$\pm$7.91 |  17.08$\pm$7.75 |  13.46$\pm$7.88 |  20.71$\pm$2.80 |  29.26$\pm$2.88 |  10.49$\pm$2.00*
$\sigma{=}2$ |  61.65$\pm$14.17 |  46.46$\pm$5.35 |  39.72$\pm$6.82 |  13.96$\pm$8.34 |  59.31$\pm$10.57 |  40.71$\pm$12.28 |  22.60$\pm$3.05 |  30.63$\pm$3.42 |  10.38$\pm$1.94*
$\sigma{=}3$ |  54.97$\pm$21.49 |  67.45$\pm$7.52 |  41.14$\pm$3.44 |  28.86$\pm$8.19 |  28.86$\pm$8.19 |  28.10$\pm$8.20 |  26.80$\pm$3.46 |  34.11$\pm$3.69 |  10.89$\pm$1.86*
Average Runtime (seconds/iteration) |
| 0.31 | 0.59 | 0.56 | 1.92 | 2.33 | 2.37 | 11.22 | 0.70 | 0.77
TABLE II: Mean and standard deviation of relative MSE in estimating the weight vector $\mathbf{x}$ and test log-likelihood in the Bayesian logistic regression along with avg. runtime, the lowest relative MSE and the highest test log-likelihood in each row are shown in bold. Alg. | PMC | GR-PMC | LR-PMC | SL-PMC | O-PMC | GRAMIS | HAIS | VAPIS | NF-PMC
---|---|---|---|---|---|---|---|---|---
Relative MSE
$\sigma{=}1$ |  0.72$\pm$0.18 |  0.66$\pm$0.15 |  0.95$\pm$0.12 |  0.51$\pm$0.12 |  0.48$\pm$0.13 |  0.60$\pm$0.07 |  0.49$\pm$0.20 |  0.96$\pm$0.07 |  0.42$\pm$0.13*
$\sigma{=}2$ |  1.35$\pm$0.24 |  0.82$\pm$0.19 |  1.35$\pm$0.19 |  0.51$\pm$0.15 |  0.46$\pm$0.14 |  0.46$\pm$0.11 |  0.49$\pm$0.20 |  0.95$\pm$0.08 |  0.42$\pm$0.12*
$\sigma{=}3$ |  2.40$\pm$0.53 |  1.32$\pm$0.24 |  1.26$\pm$0.17 |  0.59$\pm$0.17 |  0.60$\pm$0.15 |  0.47$\pm$0.11 |  0.51$\pm$0.20 |  0.99$\pm$0.08 |  0.41$\pm$0.11*
Test Log-Likelihood
$\sigma{=}1$ |  -1.66$\pm$0.35 |  -1.37$\pm$0.26 |  -2.17$\pm$0.47 |  -1.00$\pm$0.21 |  -0.92$\pm$0.17 |  -0.81$\pm$0.15 |  -1.02$\pm$0.18 |  -1.96$\pm$0.47 |  -0.75$\pm$0.16*
$\sigma{=}2$ |  -2.39$\pm$0.54 |  -1.90$\pm$0.41 |  -3.17$\pm$0.79 |  -1.20$\pm$0.25 |  -1.09$\pm$0.23 |  -0.88$\pm$0.20 |  -1.12$\pm$0.20 |  -1.97$\pm$0.46 |  -0.76$\pm$0.15*
$\sigma{=}3$ |  -3.00$\pm$0.67 |  -2.59$\pm$0.53 |  -2.75$\pm$0.69 |  -1.44$\pm$0.32 |  -1.46$\pm$0.33 |  -1.03$\pm$0.24 |  -1.21$\pm$0.24 |  -2.12$\pm$0.48 |  -0.74$\pm$0.15*
Average Runtime (seconds/iteration) |
| 0.23 | 0.61 | 0.58 | 1.12 | 1.52 | 1.56 | 3.23 | 0.69 | 0.73

###  V-A Gaussian Mixture Target Distribution

We consider the task of approximating a mixture of $P$ Gaussian densities in
$\mathbb{R}^{d}$, given as
$\pi(\mathbf{x}){=}\sum_{p=1}^{P}\alpha_{p}\mathcal{N}(\mathbf{x},\mathbf{m}_{p},\mathbf{Q}_{p})$,
where $\alpha_{p},\mathbf{m}_{p}$, and $\mathbf{Q}_{p}$ denote the proportion,
mean, and covariance of the $p$-th mixture component. We set $d{=}200$ and
$P{=}5$. In each trial, a different target density is considered by sampling
$(\alpha_{1},\cdots,\alpha_{p}){\sim}Dir(10.\mathbf{1}_{P})$,
$\mathbf{m}_{p}\sim Unif([-10,10]^{d})$, and setting
$\mathbf{Q}_{p}{=}\widetilde{\mathbf{Q}}_{p}+2\mathbb{I}_{d}$, where
$\widetilde{\mathbf{Q}}_{p}\sim\mathcal{IW}(\mathbb{I}_{d},d)$. Here
$Dir(10.\mathbf{1}_{P})$, $Unif(\mathcal{S})$, and
$\mathcal{IW}(\mathcal{B},d_{f})$ denote the Dirichlet distribution of order
$P$ with parameters $(10,\cdots,10)$, uniform distribution over the set
$\mathcal{S}$, and inverse Wishart distributions with scale matrix
$\mathcal{B}$ and degree of freedom $d_{f}$ respectively. The goal is to
estimate the mean of this multi-modal target distribution.

In each trial, the proposal means
$\\{\mu_{n}^{1}\in\mathbb{R}^{d}\\}_{n{=}1:N}$ are initialized by sampling
from $Unif([-10,10]^{d})$ for all algorithms except the proposed NF-PMC
approach. We conduct three different experiments by initializing
$\Sigma_{n}^{1}{=}\sigma^{2}\mathbf{I}_{d\times d}$ with
$\sigma\in\\{1,2,3\\}$ for $1\leqslant n\leqslant N$. We use the same
$\\{\mu_{n}^{1},\Sigma_{n}^{1}\\}$ for initializing the base distributions of
the real-NVP and Xavier initialization [40] for coupling layer parameters. The
initial learning rate is set to 0.25 and 0.005 for the VAPIS and the NF-PMC
algorithms respectively. HAIS uses a fixed step size of 0.005.

Table I shows the average MSE in the estimation of the mean of $\mathbf{x}$.
The result from the best-performing approach in each experiment is marked with
an asterisk if it is significantly better at the 0.5% level compared to all
other baselines based on the Wilcoxon signed-rank test. We observe that PMC,
GR-PMC, and LR-PMC algorithms struggle in this high-dimensional, multi-modal
setting because of the inefficiency of resampling-based proposal adaptation.
Although SL-PMC, O-PMC, GRAMIS, HAIS, and VAPIS methods achieve improved
performance, they are significantly outperformed by the proposed NF-PMC
algorithm in all cases. Moreover, NF-PMC exhibits the lowest variability of
MSE for different choices of $\sigma$ among all algorithms.

###  V-B Bayesian Logistic Regression

We consider a logistic regression problem where the binary response
$y_{i}\in\\{0,1\\}$ is related to the predictor
$\mathbf{z}_{i}\in\mathbb{R}^{d}$ via the following likelihood function:

$\displaystyle
p(y_{i}|\mathbf{z}_{i},\mathbf{x})=q_{i}^{y_{i}}(1-q_{i})^{(1-y_{i})},\text{
}q_{i}=(1+e^{-\mathbf{z}_{i}^{T}\mathbf{x}})^{-1}\,.$ (9)

We assume that the prior distribution of the weight vector
$\mathbf{x}\in\mathbb{R}^{d}$ is
$\mathcal{N}(\mathbf{0},\zeta^{2}\mathbf{I})$. Conditioned on $N_{train}$
observations $(\mathbf{z}_{i},y_{i})$, our objective is to approximate the
posterior distribution of $\mathbf{x}$, which can be evaluated up to a
constant as:

$\displaystyle\tilde{\pi}(\mathbf{x}|\\{\mathbf{z}_{i},y_{i}\\}_{i=1}^{N_{train}})\propto
p(\mathbf{x})\prod_{i=1}^{N_{train}}p(y_{i}|\mathbf{z}_{i},\mathbf{x})\,.$
(10)

In each trial, the weight vector $\mathbf{x}$ is sampled from its prior
distribution. The artificial predictors
$\mathbf{z}_{i}{=}[z_{i,1},z_{i,2},...z_{i,d}]^{T}$ are simulated as follows:
$z_{i,1}{=}1$ and $z_{i,\ell}{\sim}\mathcal{N}(0,\delta^{2})$. We generate
$y_{i}$ according to the model in Eq. (9). $N_{train}{=}500$ artificial data
points are used for training and we construct a held-out test set using
another $N_{test}{=}500$ points. We set $d{=}200,\zeta{=}10$ and
$\delta{=}0.1$. For all algorithms, we follow the same initialization strategy
as in Section V-A. The initial learning rate is set to 0.5 and 0.05 for the
VAPIS and the NF-PMC algorithms respectively. We use a fixed step size of 0.05
for HAIS.

In addition to reporting the relative MSE in estimating the weight vector
$\mathbf{x}$ as a performance metric, we estimate the average predictive log-
likelihood for the held-out test data using the Monte Carlo approximation of
the posterior distribution of $\mathbf{x}$. From Table II, we observe that in
each setting, the proposed PF-PMC algorithm obtains the lowest average
relative MSE for the estimation of the weight vector and the highest average
test log-likelihood among all competing techniques.

##  VI Conclusion

Instead of the common practice of drawing samples from isotropic Gaussian
distributions in most AIS algorithms, we employ a proposal mechanism using
normalizing flow and adapt its parameters by minimizing a discrepancy measure
between the proposal and target distributions. Experimental results
demonstrate that the superior capability of the normalizing flow proposals in
approximating complex, high-dimensional target distributions results in
significantly improved performance. The computational burden of the proposed
NF-PMC algorithm remains comparable to that of resampling-based PMC variants
(e.g. GR-PMC and LR-PMC) for generating an equal number of samples. In
contrast, several other alternatives require considerably higher execution
time and scale poorly to high dimensions. An interesting future research
direction is to explore minimization of the $\chi^{2}$ divergence instead of
the KL divergence. In some settings, this has been shown to have superior
sampling performance and has desirable theoretical properties [15, 32]. In our
experiments, minimizing the KL divergence was faster, provided more stable
convergence, and achieved better sampling performance. Another potential
avenue is an investigation of the applicability of other normalizing flow
architectures in the AIS framework. For example, observation-driven proposals
using conditional normalizing flow [29, 30, 31] might be beneficial in
Bayesian inference problems, whereas auto-regressive flow [28, 27, 26] might
be better suited for approximating high-dimensional target distributions with
complex dependency structures.

## References

  * [1] C. P. Robert and G. Casella,  Monte Carlo Statistical Methods,  Springer, 2004.
  * [2] J. S. Liu,  Monte Carlo Strategies in Scientific Computing,  Springer, 2004.
  * [3] J. V. Candy,  Bayesian signal processing: classical, modern, and particle filtering methods, vol. 54,  John Wiley & Sons, 2016.
  * [4] W. R. Gilks, S. Richardson, and D. Spiegelhalter,  Markov Chain Monte Carlo in Practice,  Taylor & Francis, 1995.
  * [5] R. M Neal,  “MCMC using Hamiltonian dynamics,”  in Handbook of Markov Chain Monte Carlo, S. Brooks, A. Gelman, G. L. Jones, and X. Meng, Eds., chapter 5, pp. 113–162. Chapman and Hall/CRC, Boca Raton, USA, 2011.
  * [6] M. Girolami and B. Calderhead,  “Riemann manifold Langevin and Hamiltonian Monte Carlo methods,”  J. Royal Statist. Society, vol. 73, pp. 123 – 214, Mar. 2011.
  * [7] G. O. Roberts and J. S. Rosenthal,  “Examples of adaptive MCMC,”  J. Comput. and Graphical Statist., vol. 18, no. 2, pp. 349–367, 2009.
  * [8] N. Mahendran, Z. Wang, F. Hamze, and N. D. Freitas,  “Adaptive MCMC with Bayesian optimization,”  in Proc. Int. Conf. Artificial Intell. and Statist., La Palma, Canary Islands, Apr. 2012, vol. 22, pp. 751–760.
  * [9] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo,  “Generalized multiple importance sampling,”  Statist. Science, vol. 34, Nov. 2015.
  * [10] M. F. Bugallo, V. Elvira, L. Martino, D. Luengo, J. Miguez, and P. M. Djuric,  “Adaptive importance sampling: the past, the present, and the future,”  IEEE Signal Process. Magazine, vol. 34, no. 4, pp. 60–79, Jul. 2017.
  * [11] O. Cappé, A. Guillin, J. M. Marin, and C. P Robert,  “Population Monte Carlo,”  J. Computat. and Graphical Statist., vol. 13, no. 4, pp. 907–929, 2004.
  * [12] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo,  “Population Monte Carlo schemes with reduced path degeneracy,”  in IEEE Int. Workshop on Computat. Adv. in Multi-Sensor Adaptive Process., Curacao, Dutch Antilles, Dec. 2017, pp. 1–5.
  * [13] V. Elvira and É. Chouzenoux,  “Langevin-based strategy for efficient proposal adaptation in population Monte Carlo,”  in Proc. IEEE Int. Conf. Acoust., Speech and Signal Process., Brighton, UK, May 2019, pp. 5077–5081.
  * [14] V. Elvira, L. Martino, D. Luengo, and J. Corander,  “A gradient adaptive population importance sampler,”  in Proc. IEEE Int. Conf. Acoust., Speech and Signal Process., Brisbane, Australia, Apr. 2015, pp. 4075–4079.
  * [15] Y. El-Laham, P. M. Djurić, and M. F. Bugallo,  “A variational adaptive population importance sampler,”  in IEEE Int. Conf. on Acoust., Speech and Signal Process., Brighton, UK, May 2019, pp. 5052–5056.
  * [16] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo,  “Improving population monte carlo: alternative weighting and resampling schemes,”  Signal Process., vol. 131, pp. 77–91, 2017.
  * [17] V. Elvira and É. Chouzenoux,  “Optimized population Monte Carlo,”  IEEE Trans. Signal Process., vol. 70, pp. 2489–2501, May 2022.
  * [18] V. Elvira, É. Chouzenoux, Ö. D. Akyildiz, and L. Martino,  “Gradient-based adaptive importance samplers,”  J. Franklin Inst., vol. 360, no. 13, pp. 9490–9514, 2023.
  * [19] O. Cappé, R. Douc, A. Guillin, J. M. Marin, and C. P. Robert,  “Adaptive importance sampling in general mixture classes,”  Statist. Comput., vol. 18, pp. 447–459, 2008.
  * [20] L. Martino, V. Elvira, D. Luengo, and J. Corander,  “Layered adaptive importance sampling,”  Statist. Comput., vol. 27, no. 3, pp. 599–623, 2017.
  * [21] A. Mousavi, R. Monsefi, and V. Elvira,  “Hamiltonian adaptive importance sampling,”  IEEE Signal Process. Letters, vol. 28, pp. 713–717, Mar. 2021.
  * [22] F. Llorente, E. Curbelo, L. Martino, V. Elvira, and D. Delgado,  “MCMC‐driven importance samplers,”  Appl. Math. Model., vol. 111, pp. 310–331, 2022.
  * [23] J. M. Cornuet, J. M. Marin, A. Mira, and C. P. Robert,  “Adaptive multiple importance sampling,”  Scandinavian J. Statist., vol. 39, no. 4, pp. 798–812, Dec. 2012.
  * [24] D. J. Rezende and S. Mohamed,  “Variational inference with normalizing flows,”  in Pro. Int. Conf. Machine Learning, Jul. 2015, p. 1530–1538.
  * [25] L. Dinh, J. Sohl-Dickstein, and S. Bengio,  “Density estimation using real NVP,”  in Proc. Int. Conf. Learning Representations, Toulon, France, Apr. 2017.
  * [26] G. Papamakarios, T. Pavlakou, and I. Murray,  “Masked autoregressive flow for density estimation,”  in Proc. Adv. Neural Info. Process. Syst., Long Beach, CA, USA, Dec. 2017, pp. 2338–2347.
  * [27] D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling,  “Improved variational inference with inverse autoregressive flow,”  in Proc. Adv. Neural Info. Process. Syst., Barcelona, Spain, Dec. 2016, pp. 4736–4744.
  * [28] M. Germain, K.l Gregor, I. Murray, and H. Larochelle,  “MADE: Masked autoencoder for distribution estimation,”  in Proc. Int. Conf. Machine Learning, Lille, France, Jul. 2015, pp. 881–889.
  * [29] B. L. Trippe and R. E. Turner,  “Conditional density estimation with Bayesian normalising flows,”  in Proc. Workshop on Bayesian Deep Learning, Adv. Neural Info. Process. Syst., Long Beach, CA, USA, Dec. 2017.
  * [30] C. Winkler, D. Worrall, E. Hoogeboom, and M. Welling,  “Learning likelihoods with conditional normalizing flows,”  ArXiv e-prints: arXiv 1912.00042, Nov. 2019.
  * [31] Y. Lu and B. Huang,  “Structured output learning with conditional generative flows,”  in Proc. AAAI Conf. Artificial Intell., New York, NY, USA, Feb. 2020, pp. 5005–5012.
  * [32] L. I. Midgley, V. Stimper, G. N. C. Simm, B. Schölkopf, and J. M. Hernández-Lobato,  “Flow annealed importance sampling bootstrap,”  in Proc. Int. Conf. Learning Representations, Kigali, Rwanda, May 2023.
  * [33] M. D. Parno and Y. M. Marzouk,  “Transport map accelerated Markov Chain Monte Carlo,”  SIAM/ASA J. Uncertainty Quantification, vol. 6, no. 2, pp. 645–682, 2018.
  * [34] M. Arbel, A. G. D. G. Matthews, and A. Doucet,  “Annealed flow transport Monte Carlo,”  in Proc. Int. Conf. Machine Learning, Virtual, Jul. 2021.
  * [35] M. Gabrié, G. M. Rotskoff, and E. Vanden-Eijnden,  “Adaptive Monte Carlo augmented with normalizing flows,”  Proc. National Academy of Sciences, vol. 119, no. 10, pp. e2109420119, 2022.
  * [36] R. D. P. Grumitt, B. Dai, and U. Seljak,  “Deterministic Langevin Monte Carlo with normalizing flows for Bayesian inference,”  in Proc. Adv. Neural Info. Process. Syst., New Orleans, LA, USA, Dec. 2022.
  * [37] G. O. Roberts and O. Stramer,  “Langevin diffusions and Metropolis-Hastings algorithms,”  Methodology and Comput. Appl. Probab., vol. 4, no. 4, pp. 337–357, Dec 2002.
  * [38] S. Duane, A. D. Kennedy, B. J. Pendleton, and D. Roweth,  “Hybrid Monte Carlo,”  Phys. Lett. B, vol. 195, no. 2, pp. 216 – 222, 1987.
  * [39] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley,  “Stochastic variational inference,”  J. Machine Learning Research, vol. 14, pp. 1303–1347, 2013.
  * [40] X. Glorot and Y. Bengio,  “Understanding the difficulty of training deep feedforward neural networks,”  in Proc. Int. Conf. Artificial Intell. and Statist., Sardinia, Italy, May 2010, pp. 249–256.