§ INTRODUCTION

Sequential Monte Carlo () methods are a class of sampling algorithms that
combine importance sampling and resampling. They have been primarily used as
“particle filters” to solve optimal filtering problems; see, for example,
Cappé et al., 2007 and Doucet and Johansen, 2011 for recent reviews. They are
also used in a static setting where a target distribution is of interest, for
example, for the purpose of Bayesian modeling. This was proposed by
Del Moral et al., 2006 and developed by Peters, 2005 and
Del Moral et al., 2006. This framework involves the construction of a
sequence of artificial distributions on spaces of increasing dimensions which
admit the distributions of interest as particular marginals.

algorithms are perceived as being difficult to implement while general
tools were not available until the development by Johansen, 2009, which
provided a general framework for implementing algorithms. algorithms
admit natural and scalable parallelization. However, there are only parallel
implementations of algorithms for many problem specific applications,
usually associated with specific related researches. Lee et al., 2010
studied the parallelization of algorithms on s with some
generality. There are few general tools to implement algorithms on
parallel hardware though multicore s are very common today and computing
on specialized hardware such as s are more and more popular.

The purpose of the current work is to provide a general framework for
implementing algorithms on both sequential and parallel hardware. There
are two main goals of the presented framework. The first is reusability. It
will be demonstrated that the same implementation source code can be used to
build a serialized sampler, or using different programming models (for
example, and ) to build parallelized samplers for multicore
s. They can be scaled for clusters using with few modifications.
And with a little effort they can also be used to build parallelized samplers
on specialized massive parallel hardware such as s using . The
second is extensibility. It is possible to write a backend for to use
new parallel programming models while reusing existing implementations. It is
also possible to enhance the library to improve performance for specific
applications. Almost all components of the library can be reimplemented by
users and thus if the default implementation is not suitable for a specific
application, they can be replaced while being integrated with other components

§ SEQUENTIAL MONTE CARLO

§.§ Sequential importance sampling and resampling

Importance sampling is a technique which allows the calculation of the
expectation of a function $\varphi$ with respect to a distribution $\pi$ using
samples from some other distribution $\eta$ with respect to which $\pi$ is
absolutely continuous, based on the identity,

\begin{equation}
\Exp_{\pi}[\varphi(X)]
= \int\varphi(x)\pi(x)\intd x
= \int\frac{\varphi(x)\pi(x)}{\eta(x)}\eta(x)\intd x
= \Exp_{\eta}\Bigl[\frac{\varphi(X)\pi(X)}{\eta(X)}\Bigr]
\end{equation}

And thus, let $\{X^{(i)}\}_{i=1}^N$ be samples from $\eta$, then
$\Exp_{\pi}[\varphi(X)]$ can be approximated by

\begin{equation}
\hat\varphi_1 =
\frac{1}{N}\sum_{i=1}^N\frac{\varphi(X^{(i)})\pi(X^{(i)})}{\eta(X^{(i)})}
\end{equation}

In practice $\pi$ and $\eta$ are often only known up to some normalizing
constants, which can be estimated using the same samples. Let $w^{(i)} =
\pi(X^{(i)})/\eta(X^{(i)})$, then we have

\begin{equation}
\hat\varphi_2 =
\frac{\sum_{i=1}^Nw^{(i)}\varphi(X^{(i)})}{\sum_{i=1}^Nw^{(i)}}
\end{equation}

\begin{equation}
\hat\varphi_3 = \sum_{i=1}^NW^{(i)}\varphi(X^{(i)})
\end{equation}

where $W^{(i)}\propto w^{(i)}$ and are normalized such that
$\sum_{i=1}^NW^{(i)} = 1$.

Sequential importance sampling () generalizes the importance sampling
technique for a sequence of distributions $\{\pi_t\}_{t\ge0}$ defined on
spaces $\{\prod_{k=0}^tE_k\}_{t\ge0}$. At time $t = 0$, sample
$\{X_0^{(i)}\}_{i=1}^N$ from $\eta_0$ and compute the weights $W_0^{(i)}
\propto \pi_0(X_0^{(i)})/\eta_0(X_0^{(i)})$. At time $t\ge1$, each sample
$X_{0:t-1}^{(i)}$, usually termed particles in the literature, is
extended to $X_{0:t}^{(i)}$ by a proposal distribution
$q_t(\cdot|X_{0:t-1}^{(i)})$. And the weights are recalculated by $W_t^{(i)}
\propto \pi_t(X_{0:t}^{(i)})/\eta_t(X_{0:t}^{(i)})$ where

\begin{equation}
\eta_t(X_{0:t}^{(i)}) =
\eta_{t-1}(X_{0:t-1}^{(i)})q_t(X_{0:t}^{(i)}|X_{0:t-1}^{(i)})
\end{equation}

and thus

\begin{align}
W_t^{(i)} \propto \frac{\pi_t(X_{0:t}^{(i)})}{\eta_t(X_{0:t}^{(i)})}
&= \frac{\pi_t(X_{0:t}^{(i)})\pi_{t-1}(X_{0:t-1}^{(i)})}
\pi_{t-1}(X_{0:t-1}^{(i)})} \notag\\
&= \frac{\pi_t(X_{0:t}^{(i)})}
\label{eq:si}
\end{align}

and importance sampling estimate of $\Exp_{\pi_t}[\varphi_t(X_{0:t})]$ can be
obtained using $\{W_t^{(i)},X_{0:t}^{(i)}\}_{i=1}^N$.

However this approach fails as $t$ becomes large. The weights tend to become
concentrated on a few particles as the discrepancy between $\eta_t$ and
$\pi_t$ becomes larger. Resampling techniques are applied such that, a new
particle system $\{\bar{W}_t^{(i)},\bar{X}_{0:t}^{(i)}\}_{i=1}^M$ is obtained
with the property,

\begin{equation}
\Exp\Bigl[\sum_{i=1}^M\bar{W}_t^{(i)}\varphi_t(\bar{X}_{0:t}^{(i)})\Bigr] =
\Exp\Bigl[\sum_{i=1}^NW_t^{(i)}\varphi_t(X_{0:t}^{(i)})\Bigr]
\label{eq:resample}
\end{equation}

In practice, the resampling algorithm is usually chosen such that $M = N$ and
$\bar{W}^{(i)} = 1/N$ for $i=1,\dots,N$. Resampling can be performed at each time
$t$ or adaptively based on some criteria of the discrepancy. One popular
quantity used to monitor the discrepancy is effective sample size
(), introduced by Liu and Chen, 1998, defined as

\begin{equation}
\ess_t = \frac{1}{\sum_{i=1}^N (W_t^{(i)})^2}
\end{equation}

where $\{W_t^{(i)}\}_{i=1}^N$ are the normalized weights. And resampling can
be performed when $\ess\le \alpha N$ where $\alpha\in[0,1]$.

The common practice of resampling is to replicate particles with large weights
and discard those with small weights. In other words, instead of generating a
random sample $\{\bar{X}_{0:t}^{(i)}\}_{i=1}^N$ directly, a random sample of
integers $\{R^{(i)}\}_{i=1}^N$ is generated, such that $R^{(i)} \ge 0$ for $i
= 1,\dots,N$ and $\sum_{i=1}^N R^{(i)} = N$. And each particle value
$X_{0:t}^{(i)}$ is replicated for $R^{(i)}$ times in the new particle system.
The distribution of $\{R^{(i)}\}_{i=1}^N$ shall fulfill the requirement of
Equation <ref>. One such distribution is a multinomial
distribution of size $N$ and weights $(W_t^{(i)},\dots,W_t^{(N)})$. See
Douc et al., 2005 for some commonly used resampling algorithms.

§.§ samplers

samplers allow us to obtain, iteratively, collections of weighted samples
from a sequence of distributions $\{\pi_t\}_{t\ge0}$ over essentially any
random variables on some spaces $\{E_t\}_{t\ge0}$, by constructing a sequence
of auxiliary distributions $\{\tilde\pi_t\}_{t\ge0}$ on spaces of increasing
dimensions, $\tilde\pi_t(x_{0:t})=\pi_t (x_t) \prod_{s=0}^{t-1}
L_s(x_{s+1},x_s)$, where the sequence of Markov kernels $\{L_s\}_{s=0}^{t-1}$,
termed backward kernels, is formally arbitrary but critically influences the
estimator variance. See Del Moral et al., 2006 for further details and
guidance on the selection of these kernels.

Standard sequential importance sampling and resampling algorithms can then be
applied to the sequence of synthetic distributions, $\{\tilde\pi_t\}_{t\ge0}$.
At time $t - 1$, assume that a set of weighted particles
$\{W_{t-1}^{(i)},X_{0:t-1}^{(i)}\}_{i=1}^N$ approximating $\tilde\pi_{t-1}$ is
available, then at time $t$, the path of each particle is extended with a
Markov kernel say, $K_t(x_{t-1}, x_t)$ and the set of particles
$\{X_{0:t}^{(i)}\}_{i=1}^N$ reach the distribution $\eta_t(X_{0:t}^{(i)}) =
\eta_0(X_0^{(i)})\prod_{k=1}^tK_t(X_{t-1}^{(i)}, X_t^{(i)})$, where $\eta_0$
is the initial distribution of the particles. To correct the discrepancy
between $\eta_t$ and $\tilde\pi_t$, Equation <ref> is applied and in
this case,

\begin{equation}
W_t^{(i)} \propto \frac{\tilde\pi_t(X_{0:t}^{(i)})}{\eta_t(X_{0:t}^{(i)})}
= \frac{\pi_t(X_t^{(i)})\prod_{s=0}^{t-1}L_s(X_{s+1}^{(i)}, X_s^{(i)})}
\propto \tilde{w}_t(X_{t-1}^{(i)}, X_t^{(i)})W_{t-1}^{(i)}
\end{equation}

where $\tilde{w}_t$, termed the incremental weights, are calculated as,

\begin{equation}
\tilde{w}_t(X_{t-1}^{(i)},X_t^{(i)}) =
\frac{\pi_t(X_t^{(i)})L_{t-1}(X_t^{(i)}, X_{t-1}^{(i)})}
{\pi_{t-1}(X_{t-1}^{(i)})K_t(X_{t-1}^{(i)}, X_t^{(i)})}
\end{equation}

If $\pi_t$ is only known up to a normalizing constant, say $\pi_t(x_t) =
\gamma_t(x_t)/Z_t$, then we can use the unnormalized incremental

\begin{equation}
w_t(X_{t-1}^{(i)},X_t^{(i)}) =
\frac{\gamma_t(X_t^{(i)})L_{t-1}(X_t^{(i)}, X_{t-1}^{(i)})}
{\gamma_{t-1}(X_{t-1}^{(i)})K_t(X_{t-1}^{(i)}, X_t^{(i)})}
\end{equation}

for importance sampling. Further, with the previously normalized
weights $\{W_{t-1}^{(i)}\}_{i=1}^N$, we can estimate the ratio of normalizing
constant $Z_t/Z_{t-1}$ by

\begin{equation}
\frac{\hat{Z}_t}{Z_{t-1}} =
\sum_{i=1}^N W_{t-1}^{(i)}w_t(X_{t-1}^{(i)},X_t^{(i)})
\end{equation}

Sequentially, the normalizing constant between initial distribution $\pi_0$
and some target $\pi_T$, $T\ge1$ can be estimated. See Del Moral et al., 2006
for details on calculating the incremental weights. In practice, when $K_t$ is
invariant to $\pi_t$, and an approximated suboptimal backward kernel

\begin{equation}
L_{t-1}(x_t, x_{t-1}) = \frac{\pi(x_{t-1})K_t(x_{t-1}, x_t)}{\pi_t(x_t)}
\end{equation}

is used, the unnormalized incremental weights will be

\begin{equation}
w_t(X_{t-1}^{(i)},X_t^{(i)}) =
\frac{\gamma_t(X_{t-1}^{(i)})}{\gamma_{t-1}(X_{t-1}^{(i)})}.
\label{eq:inc_weight_mcmc}
\end{equation}

§.§ Other sequential Monte Carlo algorithms

Some other commonly used sequential Monte Carlo algorithms can be viewed as
special cases of algorithms introduced above. The annealed importance sampling
(; Neal, 2001) can be viewed as samplers without resampling.

Particle filters as seen in the physics and signal processing literature, can
also be interpreted as the sequential importance sampling and resampling
algorithms. See Doucet and Johansen, 2011 for a review of this topic. A simple
particle filter example is used in Section <ref>
to demonstrate basic features of the library.

§ USING THE VSMC LIBRARY

§.§ Overview

The library makes use of 's template generic programming to
implement general algorithms. This library is formed by a few major
modules listed below. Some features not included in these modules are
introduced later in context.

Core The highest level of abstraction of samplers. Users
interact with classes defined within this module to create and manipulate
general samplers. Classes in this module include Sampler,
Particle and others. These classes use user defined callback
functions or callable objects, such as functors, to perform problem
specific operations, such as updating particle values and weights. This
module is documented in Section <ref>.

Symmetric Multiprocessing () This is the form of computing most
people use everyday, including multiprocessor workstations, multicore
desktops and laptops. Classes within this module make it possible to write
generic operations which manipulate a single particle that can be applied
either sequentially or in parallel through various parallel programming
models. A method defined through classes of this module can be used by
Sampler as callback objects. This module is documented in
Section <ref>.

Message Passing Interface is the de facto standard
for parallel programming on distributed memory architectures. This module
enables users to adapt implementations of algorithms written for the module such that the same sampler can be parallelized using . In
addition, when used with the module, it allows easy implementation of
hybrid parallelization such as /. In Section <ref>, an example is shown how to extend
existing parallelized samplers using this module.

OpenCL This module is similar to the two above except it eases the
parallelization through , such as for the purpose of General
Purpose Programming (). is a framework for writing
programs that can be execute across heterogeneous platforms. programs can run on either s or s. It is beyond the scope of
this paper to give a proper introduction to , and their use
in . However, later we will demonstrate the relative performance of
this programming model on both s and s in
Section <ref>.

§.§ Obtaining and installation

is a header only library. There is practically no installation step. The
library can be downloaded from <http://zhouyan.github.io/vSMC/vSMC.zip>.
After downloading and unpacking, one can start using by ensuring that
the compiler can find the headers inside the include directory. To
permanently install the library in a system directory, one can simply copy the
contents of the include directory to an appropriate place.

Alternatively, one can use the (2.8 or later) configuration script and
obtain the source by . On a Unix-like system (such as Mac OS X, BSD,
Linux and others),

git clone git://github.com/zhouyan/vSMC.git
cd vSMC
git submodule init
git submodule update
mkdir build
cd build
cmake .. -DCMAKE_BUILD_TYPE=Release
make install

For Unix-like systems, there is also a shell script build.sh that
builds all examples in this paper and produces the results and benchmarks as
in Section <ref>. See documentations in that script for details of
how to change settings for the users' platform.

§.§ Documentation

To build the reference manual, one need , version 1.8.3 or later.
Continuing the last step (still in the build directory), invoking
make docs will create a doc directory inside build, which
contains the references. Alternatively the reference manual can also be
found on <http://zhouyan.github.io/vSMC/doc/html/index.html>. It is beyond
the scope of this paper to document every feature of the library. In
many places we will refer to this reference manual for further information.

§.§ Third-party dependencies

uses counter-based for random number generating. For an
sampler with $N$ particles, constructs $N$ (statistically)
independent streams. It is possible to use millions of such streams
without a huge memory footprint or other performance penalties. Since each
particle has its own independent stream, it frees users from many
thread-safety and statistical independence considerations. It is highly
recommended that the users install this library. Within , these streams are wrapped under engines, and can be replaced by other
compatible engines seamlessly. Users only need to be familiar with
classes defined in  <random> or their equivalents to
use these streams. See the documentation of the corresponding libraries
for details, as well as examples in Section <ref>.

The other third-party dependency is the library. Version 1.49 or later
is required. However, this is actually optional provided that proper features are available in the standard library, for example using with
. The headers of concern are <functional> and
<random>. To instruct to use the standard library headers instead
of falling back to the library, one needs to define the configuration
macros before including any headers. For example,

clang++ -std=c++11 -stdlib=libc++         -DVSMC_HAS_CXX11LIB_FUNCTIONAL=1      -DVSMC_HAS_CXX11LIB_RANDOM=1          -o prog prog.cpp

tells the library to use  <functional> and <random>.
The availability of these headers are also checked by the configuration

§.§ Compiler support

has been tested with recent versions of , , and
. can optionally use some features to improve performance
and usability. In particular, as mentioned before, can use standard libraries instead of the library. At the time of writing,
with has the most comprehensive support of with respect
to standard compliant and feature completion. 4.8 , 2012 and 2013 also have very good support. Note that, provided the library is available, all language and library features are optional.
can be used with any conforming compilers.

§ THE VSMC LIBRARY

§.§ Core module

The core module abstracts general samplers. samplers can be viewed
as formed by a few concepts regardless of the specific problems. The following
is a list of the most commonly seen components of samplers and their
corresponding abstractions.

* A collection of all particle state values, namely
$\{X_t^{(i)}\}_{i=1}^N$. In , users need to define a class, say
T, to abstract this concept. We refer to this as the value
collection. We will slightly abuse the generic name T in this
paper. Whenever a template parameter is mentioned with the name T,
it always refers to such a value collection type unless stated otherwise.

* A collection of all particle state values and their associated
weights. This is abstracted by a Particle<T> object. We refer to
this as the particle collection. A Particle<T> object has
three primary sub-objects. One is the above type T value collection
object. Another is an object that abstracts weights
$\{W_t^{(i)}\}_{i=1}^N$. By default this is a WeightSet object. The
last is a collection of streams, one for each particle. By default
this is an RngSet object.

* Operations that perform tasks common to all samplers to these
particles, such as resampling. These are the member functions of

* A sampler that updates the particles (state values and weights) using
user defined callbacks. This is a Sampler<T> object.

* Monitors that record the importance sampling estimates of certain
functions defined for the values when the sampler iterates. These are
Monitor<T> objects. A monitor for the estimates of $E[h(X_t)]$
computes $h(X_t^{(i)})$ for each $i = 1,\dots,N$. The function value
$h(X_t)$ is allowed to be a vector.

Note that within the core module, all operations are applied to
Particle<T> objects, that is $\{W_t^{(i)},X_t^{(i)}\}_{i=1}^N$, instead
of a single particle. Later we will see how to write operations that can be
applied to individual particles and can be parallelized easily.

§.§.§ Program structures

A program usually consists of the following tasks.

* Define a value collection type T.

* Constructing a Sampler<T> object.

* Configure the behavior of initialization and updating by adding
callable objects to the sampler object.

* Optionally add monitors.

* Initialize and iterate the sampler.

* Retrieve the outputs, estimates and other informations.

In this section we document how to implement each of these tasks. Within the
library, all public classes, namespaces and free functions, are declared
in the namespace vsmc.

§.§.§ The value collection

The template parameter T is a user defined type that abstracts the
value collection. does not restrict how the values shall be actually
stored. They can be stored in memory, spread among nodes of a cluster, in memory or even in a database. However this kind of flexibility comes with a
small price. The value collection does need to fulfill two requirements. We
will see later that for most common usage, provides readily usable
implementations, on top of which users can create problem specific classes.

First, the value collection class T has to provide a constructor of the

T (SizeType N)

where SizeType is some integer type. Since allows one to allocate
the states in any way suitable, one needs to provide this constructor which
Particle<T> can use to allocate them.

Second, the class has to provide a member function named copy that
copies each particle according to replication numbers given by a resampling
algorithm. For the same reason as above has no way to know how it can
extract and copy a single particle when it is doing the resampling. The
signature of this member function may look like

template <typename SizeType>
void copy (std::size_t N, const SizeType *copy_from);

The pointer copy_from points to an array that has $N$ elements,
where $N$ is the number of particles. After calling this member function, the
value of particle i shall be copied from the particle
j = copy_from[i]. In other words, particle i is a child of
particle copy_from[i], or copy_from[i] is the parent of particle
i. If a particle j shall remain itself, then
copy_from[j] == j. How the values are actually copied is user defined.
Note that, the member function can take other forms, as usual when using template generic programming.
The actual type of the pointer copy_from, SizeType, is
Particle<T>::size_type, which depends on the type T. For
example, define the member function as the following is also allowed,

void copy (int N, const std::size_t *copy_from);

provided that Particle<T>::size_type is indeed std::size_t.
However, writing it as a member function template releases the users from
finding the actual type of pointer copy_from and the sometimes
troubling forward declaration issues. Will not elaborate such more technical
issues further.

§.§.§ Constructing a sampler

Once the value collection class T is defined. One can start
constructing samplers. For example, the following line creates a sampler
with $N$ particles

Sampler<T> sampler(N);

The number of particles is the only mandatory argument of the constructor.
There are two optional parameters. The complete signature of the constructor

explicit Sampler<T>::Sampler (Sampler<T>::size_type N,
ResampleScheme scheme = Stratified, double threshold = 0.5);

The scheme parameter is self-explanatory. provides a few built-in
resampling schemes; see the reference manual for a list of them. User defined
resampling algorithms can also be used. See the reference manual for details.
The threshold is the threshold of $\ess/N$ below which a resampling
will be performed. It is obvious that if $\text{\code{threshold}}\ge1$ then
resampling will be performed at each iteration. If
$\text{\code{threshold}}\le0$ then resampling will never be performed. Both
parameters can be changed later. However the size of the sampler can never be
changed after the construction.

§.§.§ Initialization and updating

All the callable objects that initialize, move and weight the particle
collection can be added to a sampler through a set of member functions. All
these objects operate on the Particle<T> object. Because allows one to manipulate the particle collection as a whole, in principle many
kinds of parallelization are possible.

To set an initialization method, one need to implement a function with the
following signature,

std::size_t init_func (Particle<T>    particle, void *param);

or a class with operator() properly defined, such as

struct init_class
std::size_t operator() (Particle<T>    particle, void *param); ;

They can be added to the sampler through

respectively.  std::function or its equivalent
boost::function can also be used. For example,

std::function<std::size_t (Particle<T>    , void *)> init_obj(init_func);

The addition of updating methods is more flexible. There are two kinds of
updating methods. One is simply called move in , and is performed
before the possible resampling at each iteration. These moves usually perform
the updating of the weights among other tasks. The other is called
mcmc, and is performed after the possible resampling. They are often
type moves. Multiple move's or mcmc's are also allowed. In
fact a sampler consists of a queue of move's and a queue of
mcmc's. The move's in the queue can be changed through

Sampler<T>    move (const Sampler<T>::move_type    new_move, bool append);

If append == true then new_move is appended to the existing
(possibly empty) queue. Otherwise, the existing queue is cleared and
new_move is added. The member function returns a reference to the
updated sampler. For example, the following move,

std::size_t move_func (std::size_t iter, Particle<T>    particle);

can be added to a sampler by

sampler.move(move_func, false);

This will clear the (possibly empty) existing queue of move's and set a
new one. To add multiple moves into the queue,

sampler.move(move1, true).move(move2, true).move(move3, true);

Objects of class type with operator() properly defined can also be
used, similarly to the initialization method. The queue of mcmc's can
be used similarly. See the reference manual for other methods that can be used
to manipulate these two queues.

In principle, one can combine all moves into a single move. However, sometimes
it is more natural to think of a queue of moves. For instance, if a
multi-block Metropolis random walk consists of kernels $K_1$ and $K_2$, then
one can implement each of them as functions, say mcmc_k1 and
mcmc_k2, and add them to the sampler sequentially,

sampler.mcmc(mcmc_k1, true).mcmc(mcmc_k2, true);

Then at each iteration, they will be applied to the particle collection
sequentially in the order in which they are added.

§.§.§ Running the algorithm, monitoring and outputs

Having set all the operations, one can initialize and iterate the sampler by

sampler.initialize((void *)param);

The param argument to initialize is optional, with NULL
as its default. This parameter is passed to the user defined init_func.
The iter_num argument to iterate is also optional; the default
is $1$.

Before initializing the sampler or after a certain time point, one can add
monitors to the sampler. The concept is similar to 's monitor
statement, except it does not monitor the individual values but rather the
importance sampling estimates. Consider approximating $\Exp[h(X_t)]$, where
$h(X_t) = (h_1(X_t),\dots,h_m(X_t))$ is an $m$-vector function. The importance
sampling estimate can be obtained by $AW$ where $A$ is an $N$ by $m$ matrix
where $A(i,j) = h_j(X_t^{(i)})$ and $W = (W_t^{(i)},\dots,W_t^{(N)})^T$ is the
$N$-vector of normalized weights. To compute this importance sampling
estimate, one need to define the following evaluation function (or a class
with operator() properly defined),

void monitor_eval (std::size_t iter, std::size_t m,
const Particle<T>    particle, double *res)

and add it to the sampler by calling,

sampler.monitor("some.user.chosen.variable.name", m, monitor_eval);

When the function monitor_eval is called, iter is the iteration
number of the sampler, m is the same value as the one the user passed
to Sampler<T>::monitor; and thus one does not need global variable or
other similar techniques to access this value. The output pointer res
points to an $N \times m$ output array of row major order. That is, after the
calling of the function, res[i * dim + j] shall be $h_j(X_t^{(i)})$.

After each iteration of the sampler, the importance sampling estimate will be
computed automatically. See the reference manual for various ways to retrieve
the results. Usually it is sufficient to output the sampler by

std::ofstream output("file.name");
output << sampler << std::endl;

where the output file will contain the importance sampling estimates among
other things. Alternatively, one can use the Monitor<T>::record member
function to access specific historical results. See the reference manual for
details of various overloaded version of this member function.

A reference to the value collection T object can be retrieved through
the Particle<T>::value member function. Particle<T> objects
manage the weights through a weight set object, which by default is of type
WeightSet. The Particle<T>::weight_set member function returns a
reference to this weigh set object. A user defined weight set class that
abstracts $\{W_t^{(i)}\}_{i=1}^N$ can also be used. The details involve some
more advanced template techniques and are documented in the reference
manual. One possible reason for replacing the default is to provide special
memory management of the weight set. For example, the module provides a
special weight set class that manages weights across multiple nodes and
perform proper normalization, computation of , and other tasks.

The default WeightSet object provides some ways to retrieve weights.
Here we document some of the most commonly used. See the reference manual for
details of others. The weights can be accessed one by one, for example,

Particle<T>    particle = sampler.particle();
double w_i     = particle.weight_set().weight(i);
double log_w_i = particle.weight_set().log_weight(i);

One can also read all weights into a container, for example,

std::vector<double> w(particle.size());
double *lw = new double[particle.size()];

Note that these member functions accept general output iterators.

§.§.§ Implementing initialization and updating

So far we have only discussed how to add initialization and updating objects
to a sampler. To actually implement them, one writes callable objects that
operate on the Particle<T> object. For example, a move can be
implemented through the following function as mentioned before,

std::size_t move_func (std::size_t iter, Particle<T>    particle);

Inside the body of this function, one can change the value by manipulating the
object through the reference returned by particle.value().

When using the default weight set class, the weights can be updated through a
set of member functions of WeightSet. For example,

std::vector<double> weight(particle.size());

The set_equal_weight member function sets all weights to be equal. The
set_weight and set_log_weight member functions set the values of
weights and logarithm weights, respectively. The mul_weight and
add_log_weight member functions multiply the weights or add to the
logarithm weights by the given values, respectively. All these member functions
accept general input iterators as their arguments.

One important thing to note is that, whenever one of these member functions is
called, both the weights and logarithm weights will be re-calculated,
normalized, and the will be updated. The reason for not allowing changing
a single particle's weight is that, in a multi-threading environment, it is
possible for one to change one weight in one thread, and obtain another in
another thread without proper normalizing. Conceptually, changing one weight
actually changes all weights.

§.§.§ Generating random numbers

The Particle<T> object has a sub-object, a collection of engines
that can be used with  <random> or distributions. For
each particle i, one can obtain an engine that produces an stream
independent of others by

To generate distribution random variates, one can use the
<random> library. For example,

std::normal_distribution<double> rnorm(mean, sd);
double r = rnorm(particle.rng(i));

or use the library,

boost::random::normal_distribution<double> rnorm(mean, sd);
double r = rnorm(particle.rng(i));

itself also makes use of  <random> or depending
on the value of the configuration macro VSMC_HAS_CXX11LIB_RANDOM.
Though the user is free to choose which one to use in their own code, there is
a convenient alternative. For each class defined in  <random>,
it is imported to the vsmc::cxx11 namespace. Therefore one can use

cxx11::normal_distribution<double> rnorm(mean, sd);

while the underlying implementation of normal_distribution can be
either standard library or . The benefit is that if one needs to
develop on multiple platforms, and only some of them support and some
of them have the library, then only the configure macro
VSMC_HAS_CXX11LIB_RANDOM needs to be changed. This can be configured
through and other build systems. Of course, one can also use an
entirely different system than those provided by .

§.§ SMP module

§.§.§ The value collection

Many typical problems' value collections can be viewed as a matrix of certain
type. For example, a simple particle filter whose state is a vector of length
Dim and type double can be viewed as an $N$ by Dim matrix
where $N$ is the number of particles. A trans-dimensional problem can use an
$N$ by $1$ matrix whose type is a user defined class, say StateType.
For this kind of problems, provide a class template

template <MatrixOrder Order, std::size_t Dim, typename StateType>
class StateMatrix;

which provides the constructor and the copy member function required by
the core module interface, as well as methods for accessing individual values.
The first template parameter (possible value RowMajor or
ColMajor) specifies how the values are ordered in memory. Usually one
shall choose RowMajor to optimize data access. The second template
parameter is the number of variables, an integer value no less than $1$ or the
special value Dynamic, in which case StateMatrix provides a
member function resize_dim such that the number of variables can be
changed at runtime. The third template parameter is the type of the state

Each particle's state is thus a vector of length Dim, indexed from
0 to Dim - 1. To obtain the value at position pos of the
vector of particle i, one can use one of the following member

StateBase<RowMajor, Dim, StateType> value(N);
StateType val = value.state(i, pos);
StateType val = value.state(i, Position<Pos>());
StateType val = value.state<Pos>(i);

where Pos is a compile time constant expression whose value is the same
as pos, assuming the position is known at compile time. One can also
read all values. To read the variable at position pos,

std::vector<StateType> vec(value.size());
value.read_state(pos, vec.begin());

Or one can read all values through an iterator,

std::vector<StateType> mat(Dim * value.size());
value.read_state_matrix(ColMajor, mat.first());

Alternatively, one can also read all values through an iterator which points
to iterators,

std::vector<std::vector<StateType> > mat(Dim);
for (std::size_t i = 0; i != Dim; ++i)
std::vector<std::vector<StateType>::iterator> iter(Dim);
for (std::size_t i = 0; i != Dim; ++i)
iter[i] = mat[i].begin();

If the compiler support  <tuple>, vSMC also provides a
StateTuple class template, which is similar to StateMatrix
except that the types of values do not have to be the same for each variable.
This is similar to 's data.frame. For example, suppose each
particle's state is formed by two double's, an int and a user
defined type StateType, then the following constructs a value
collection using StateTuple,

StateTuple<ColMajor, double, double, int, StateType> value(N);

And there are a few ways to access the state values, similar to

double x0 = value.state(i, Position<0>());
int x2 = value.state<2>(i);
std::vector<StateType> vec(value.size());
state.read_state(Position<3>(), vec.begin());

See the reference manual for details.

§.§.§ A single particle

For a Particle<T> object, one can construct a SingleParticle<T>
object that abstracts one of the particle from the collection. For example,

Particle<T> particle(N);
SingleParticle<T> sp(i,    particle);

create a SingleParticle<T> object corresponding to the particle at
position i. There are a few member functions of
SingleParticle<T> that makes access to individual particles easier than
through the interface of Particle<T>. Firt sp.id() returns the
value of the argument i in the above code that created this
SingleParticle<T> object. In addition, sp.rng() is equivalent to
particle.rng(i). Also sp.particle() returns a constant reference
to the Particle<T> object. And sp.particle_ptr() returns a
pointer to such a constant Particle<T> object. Note that, one cannot
get write access to a Particle<T> object through interface of
SingleParticle<T>. Instead, one can only get write access to a single
particle. For example, If T is a StateMatrix instantiation or
its derived class, then sp.state(pos) is equivalent to
particle.value().state(i, pos) and the reference it returns is mutable.
See the reference manual for more informations on the interface of the
SingleParticle class template. SingleParticle<T> objects are
usually not constructed by users, but rather by the libraries' other classes
in the module, and passed to user defined functions, as we will see very

§.§.§ Sequential and parallel implementations

Once we have the SingleParticle<T> concept, we are ready to introduce
how to write implementations of algorithms that manipulate a single
particle and can be applied to all particles in parallel or sequentially. For
sequential implementations, this can be done through five base classes,

template <typename BaseState> class StateSEQ;
template <typename T, typename D = Virtual> class InitializeSEQ;
template <typename T, typename D = Virtual> class MoveSEQ;
template <typename T, typename D = Virtual> class MonitorEvalSEQ;
template <typename T, typename D = Virtual> class PathEvalSEQ;

The template parameter BaseState needs to satisfy the general value
collection requirements in addition to a copy_particle member function,
for example, StateMatrix. Other base classes expect T to satisfy
general value collection requirements. The details of all these class
templates can be found in the reference manual. Here we use the
MoveSEQ<T> class as an example to illustrate their usage. Recall that
Sampler<T> expect a callable object which has the following signature
as a move,

std::size_t move_func (std::size_t iter, Particle<T>    particle);

For the purpose of illustration, the type T is defined as,

typedef StateMatrix<RowMajor, 1, double> T;

Here is a typical example of implementation of such a function,

std::size_t move_func (std::size_t iter, Particle<T>    particle)

std::vector<double> inc_weight(particle.size());
for (Particle<T>::size_type i = 0; i != particle.size(); ++i)
cxx11::normal_distribution<double> rnorm(0, 1);
particle.value().state(i, 0) = rnorm(particle.rng(i));
inc_weight[i] = cal_inc_weight(particle.value().state(i, 0);

where cal_inc_weight is some function that calculates the logarithm
incremental weights. As we see, there are three main parts of a typical move.
First, we allocate a vector inc_weight. Second, we generate normal
random variates for each particle and calculate the incremental weights. This
is done through a for loop. Third, we add the logarithm incremental
weights. The first and the third are global operations while the second
is local. The first and the third are often optional and absent. The
local operation is also usually the most computational intensive part of algorithms and can benefit the most from parallelizations.

With MoveSEQ<T>, one can derive from this class, which defines
the operator() required by the core module interface,

std::size_t operator() (std::size_t iter, Particle<T>    particle);

and customize what this operator does by defining one or more of the following
three member functions, corresponding to the three parts, respectively,

void pre_processor (std::size_t iter, Particle<T>    particle);
std::size_t move_state (std::size_t iter, SingleParticle<T> sp);
void post_processor (std::size_t iter, Particle<T>    particle);

For example,

#include <vsmc/smp/backend_seq.hpp>

class move : public MoveSEQ<T>

public :

void pre_processor (std::size_t iter, Particle<T>    particle)

std::size_t move_state (std::size_t iter, SingleParticle<T> sp)

cxx11::normal_distribution<double> rnorm(0, 1);
sp.state(0) = rnorm(sp.rng());
inc_weight_[sp.id()] = cal_inc_weight(sp.state(0));

void post_processor (std::size_t iter, Particle<T>    particle)

private :

std::vector<double> inc_weight_;

The operator() of MoveSEQ<T> is equivalent to the single
function implementation as shown before.

In the simplest case, MoveSEQ only takes away the loop around Part 2.
However if one implement the move in such a way, and then replace
MoveSEQ with MoveOMP, the changing of the base class name causes
to use to parallelize the loop. For example, one can declare the
class as

#include <vsmc/smp/backend_omp.hpp>
class move : public MoveOMP<T>;

and use exactly the same implementation as before. Now when
move::operator() is called, it will be parallelized by . Other
backends are available in case is not available. Among them there are
and . In addition to these well known parallelization programming
models, also has its own implementation using  <thread>.
There are other backends documented in the reference manual. To use any of
these parallelization, all one need to do is to change a few base class names.
In practice, one can use conditional compilation, for example, to use a
sequential implementation or a parallelized one, we can write,

#ifdef USE_SEQ
#include <vsmc/smp/backend_seq.hpp>
#define BASE_MOVE MoveSEQ

#ifdef USE_OMP
#include <vsmc/smp/backend_omp.hpp>
#define BASE_MOVE MoveOMP

class move : public BASE_MOVE<T>;

And we can compile the same source into different samplers with
Makefile rules such as the following,

prog-seq : prog.cpp
$(CXX) $(CXXFLAGS) -DUSE_SEQ -o prog-seq prog.cpp
prog-omp : prog.cpp
$(CXX) $(CXXFLAGS) -DUSE_omp -o prog-omp prog.cpp

Or one can configure the source file through a build system such as ,
which can also determine which programming model is available on the system.

§.§.§ Adapter

Sometimes, the cumbersome task of writing a class to implement a move and
other operations can out weight the power we gain through object oriented
programming. For example, a simple move_state-like function is all that
we need to get the work done. In this case, one can create a move through the
MoveAdapter. For example, say we have implemented

std::size_t move_state (std::size_t iter, SingleParticle<T> sp);

as a function. Then one can create a callable object acceptable to
Sampler<T>::move through

MoveAdapter<T, MoveSEQ>  move_obj(move_state);
MoveAdapter<T, MoveSTD>  move_obj(move_state);
MoveAdapter<T, MoveTBB>  move_obj(move_state);
MoveAdapter<T, MoveCILK> move_obj(move_state);
MoveAdapter<T, MoveOMP>  move_obj(move_state);
sampler.move(move_obj, false);

These are respectively, sequential,  <thread>, , ,
and implementations. The first template parameter is the type of value
collection and the second is the name of the base class template. Actually,
the MoveAdapter's constructor accepts two more optional arguments, the
pre_processor and the post_processor, corresponding to the other
two aforementioned member functions. Similar adapters for the other three base
classes also exist.

Another scenario where an adapter is desired is that which backend to use
needs be decided at runtime. The above simple adapters can already be used for
this purpose. In addition, another form of the adapter is as the following,

class move;
MoveAdapter<T, MoveTBB, move> move_obj;
sampler.move(move_obj, false);

where the class move has the same definition as before but it no longer
derives from any base class. The class move is required to have a
default constructor, a copy constructor and an assignment operator.

§.§ Thread-safety and scalability considerations

When implementing parallelized algorithms, issues such as thread-safety
cannot be avoided even though the library hides most parallel constructs
from the user.

Classes in the library usually guarantee that their member functions are
thread-safe in the sense that calling the same member function on different
objects at the same time from different threads is safe. However, calling
mutable member functions on the same object from different threads is usually
not safe. Calling immutable member functions is generally safe. There are a
few exceptions,

* The constructors of Particle and Sampler are not
thread-safe. Therefore if one need to construct multiple Sampler
from different threads, a mutex protection is needed. However, subsequent
member function calls on the constructed objects are thread-safe according
to the above rules.

* Member functions that concern a single particle are generally
thread-safe in the sense that one can call them on the same object from
different threads as long as they are called for different particles. For
example Particle::rng and StateMatrix::state are

In general, one can safely manipulate different individual particles from
different threads, which is the minimal requirement for scalable
parallelization. But one cannot manipulate the whole particle collection from
different threads, for example WeightSet::set_log_weight.

User defined callbacks shall generally follow the same rules. For example, for
a MoveOMP subclass, pre_processor and post_processor does
not have be thread-safe, but move_state needs to be. In general, avoid
write access to memory locations shared by all particles from
move_state and other similar member functions. One needs to take some
extra care when using third-party libraries. For example, in our example
implementation of the move class, the rnorm object, which is
used to generate Normal random variates, is defined within move_state
instead of being a class member data even though it is created with the same
parameters for each particle. This is because the call rnorm(sp.rng())
is not thread-safe in some implementations, for example, when using the library.

For scalable performance, certain practices should be avoided when
implementing member functions such as move_state. For example, dynamic
memory allocation is usually lock-based and thus should be avoided.
Alternatively one can use a scalable memory allocator such as the one provided
by . In general, in functions such as move_state, one should avoid
using locks to guarantee thread-safety, which can be a bottleneck to parallel

§ EXAMPLE

§.§ A simple particle filter

§.§.§ The model and algorithm

This is an example used in . Through this example, we will show how to
re-implement a simple particle filter in . It shall walk one through the
basic features of the library. is the first general framework for
implementing algorithms in and serves as one of the most important
inspirations for the library. It is widely used by many researchers. One
of the goals of the current work is that users familiar with shall find
little difficulty in using the new library.

The state space model, known as the almost constant velocity model in the
tracking literature, provides a simple scenario. The state vector $X_t$
contains the position and velocity of an object moving in a plane. That is,
$X_t = (\xpos^t, \ypos^t, \xvel^t, \yvel^t)^T$. Imperfect observations $Y_t =
(\xobs^t, \yobs^t)^T$ of the positions are possible at each time instance. The
state and observation equations are linear with additive noises,

\begin{align*}
X_t &= AX_{t-1} + V_t \\
Y_t &= BX_t + \alpha W_t
\end{align*}

\begin{equation*}
A = \begin{pmatrix}
1 & \Delta & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{pmatrix} \qquad
B = \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
\end{pmatrix} \qquad
\alpha = 0.1
\end{equation*}

and we assume that the elements of the noise vector $V_t$ are independent
Normal with variance $0.02$ and $0.001$ for position and velocity,
respectively. The observation noise, $W_t$ comprises independent, identically
distributed $t$-distributed random variables with degree of freedom $\nu =
10$. The prior at time $0$ corresponds to an axis-aligned Gaussian with
variance $4$ for the position coordinates and $1$ for the velocity
coordinates. The particle filter algorithm is shown in Algorithm <ref>.

Set $t\leftarrow0$.
$\xpos^{(0,i)},\ypos^{(0,i)}\sim\rnorm(0,4)$ and
Weight $W_0^{(i)} \propto L(X_0^{(i)}|Y_0)$ where $L$ is
the likelihood function.

Set $t\leftarrow t + 1$.

\begin{align*}
\xpos^{(t,i)}&\sim\rnorm(\xpos^{(t-1,i)} + \Delta\xvel^{(t-1,i)}, 0.02) &
\xvel^{(t,i)}&\sim\rnorm(\xvel^{(t-1,i)}, 0.001) \\
\ypos^{(t,i)}&\sim\rnorm(\ypos^{(t-1,i)} + \Delta\yvel^{(t-1,i)}, 0.02) &
\yvel^{(t,i)}&\sim\rnorm(\yvel^{(t-1,i)}, 0.001)
\end{align*}

Weight $W_t^{(i)} \propto W_{t-1}^{(i)}L(X_t^{(i)}|Y_t)$.

Repeat the Iteration step until all data are processed.
Particle filter algorithm for the almost constant velocity model.

§.§.§ Implementations

We first introduce the body of the main function, showing the typical
work flow of a program.

#include <vsmc/core/sampler.hpp>
#include <vsmc/smp/adapter.hpp>
#include <vsmc/smp/state_matrix.hpp>
#include <vsmc/smp/backend_seq.hpp>

const std::size_t DataNum = 100;
const std::size_t ParticleNum = 1000;
const std::size_t Dim = 4;

int main ()

Sampler<cv> sampler(ParticleNum);
sampler.move(cv_move(), false);
MonitorEvalAdapter<cv, MonitorEvalSEQ> cv_est(cv_monitor_state);
sampler.monitor("pos", 2, cv_est);

sampler.initialize((void *)"pf.data");
sampler.iterate(DataNum - 1);

std::ofstream est("pf.est");
est << sampler << std::endl;

return 0;

In the main function, we constructed a sampler with ParticleNum
particles. And then we added the initialization function cv_init and
the move object of type cv_move. And then we added a monitor that will
record the importance sampling estimates of the two position parameters. Next,
we initialized the sampler with data file pf.data and iterate the
sampler until all data are processed. The last step is that we output the
results into a file called pf.est.

The class cv will be our value collection which is a subclass of
StateMatrix<RowMajor, Dim, double>. To illustrate both the core module
and the module, the initialization cv_init will be implemented
as a standalone function. The move cv_move will be implemented as a
derived class of MoveSEQ<cv>. To monitor the importance sampling
estimates of the two position parameters, we will implement a simple function
cv_monitor_state and use the adapter MonitorEvalAdapter<cv,

The value collection is an $N$ by Dim (in this case $\text{\code{Dim}}
= 4$) matrix of type double. We can simply use
StateMatrix<RowMajor, Dim, double> as our value collection. However, we
would like to enhance its functionality through inheritance. First, since the
data is shared by all particles, it is natural to bind it with the value
collection. Second, both the initialization and move will need to calculate
the log-likelihood. We can implement it as a standalone function, but since
the log-likelihood function will need to access the data, it is convenient to
implement it as a member function of the value collection cv. Here is
the full implementation of this simple value collection class

class cv : public StateMatrix<RowMajor, Dim, double>

public :

cv (size_type N) :
StateMatrix<RowMajor, Dim, double>(N),
x_obs_(DataNum), y_obs_(DataNum)

double log_likelihood (std::size_t iter, size_type id) const

const double scale = 10;
const double nu = 10;
double llh_x = scale * (state(id, 0) - x_obs_[iter]);
double llh_y = scale * (state(id, 1) - y_obs_[iter]);
llh_x = std::log(1 + llh_x * llh_x / nu);
llh_y = std::log(1 + llh_y * llh_y / nu);

return -0.5 * (nu + 1) * (llh_x + llh_y);

void read_data (const char *filename)

std::ifstream data(filename);
for (std::size_t i = 0; i != DataNum; ++i)
data >> x_obs_[i] >> y_obs_[i];

private :

std::vector<double> x_obs_;
std::vector<double> y_obs_;

The log_likelihood member function accepts the iteration number and the
particle's id as arguments. It returns the log-likelihood of the
id'th particle at iteration iter. The read_data member
function simply read the data from a file.

The initialization is implemented through the cv_init function,

std::size_t cv_init (Particle<cv>    particle, void *filename);

such that, it first checks if filename is NULL. If it is not,
then we use it to read the data. So the first initialization may look like

sampler.initialize((void *)filename);

And after that, if we want to re-initialize the sampler, we can simply call,

This will reset the sampler and initialize it again but without reading the
data. If the data set is large, repeated can be very expensive. After
reading the data, we will initialize each particle's value by Normal random
variates, and calculate its log-likelihood. The last step is to set the
logarithm weights of the particle collection. Since this is not a
accept-reject type algorithm, the returned acceptance count bares no meaning.
Here is the complete implementation,

std::size_t cv_init (Particle<cv>    particle, void *filename)

if (filename)
particle.value().read_data(static_cast<const char *>(filename));

const double sd_pos0 = 2;
const double sd_vel0 = 1;
cxx11::normal_distribution<double> norm_pos(0, sd_pos0);
cxx11::normal_distribution<double> norm_vel(0, sd_vel0);
std::vector<double> log_weight(particle.size());

for (Particle<cv>::size_type i = 0; i != particle.size(); ++i)
particle.value().state(i, 0) = norm_pos(particle.rng(i));
particle.value().state(i, 1) = norm_pos(particle.rng(i));
particle.value().state(i, 2) = norm_vel(particle.rng(i));
particle.value().state(i, 3) = norm_vel(particle.rng(i));
log_weight[i] = particle.value().log_likelihood(0,i);

return 0;

In this example, we read all data from a single file for simplicity. In a
realistic application, the data is often processed online – the filter is
applied when new data becomes available. In this case, the user can use the
optional argument of Sampler<T>::initialize to pass necessary
information to open a data connection instead of a file name.

In the above implementation we iterated over all particles. There are other
ways to iterate over a particle collection. First we can use
SingleParticle<cv> objects,

for (Particle<cv>::size_type i = 0; i != particle.size(); ++i)
SingleParticle<cv> sp(i,    particle);
sp.state(0) = norm_pos(sp.rng());

Second, provides style iterators, Particle<T>::iterator,
whose value_type is SingleParticle<T>. Therefore one can write
the following loop

for (Particle<cv>::iterator iter = particle.begin();
iter != particle.end(); ++iter)
iter->state(0) = norm_pos(iter->rng());

Third, if one has a compiler supporting  auto and range-based
for, the following is also supported,

for (auto    sp : particle) sp.state(0) = norm_pos(sp.rng());

There are little or no performance difference among all these forms of loops.
However, one can choose an appropriate form to work with the interfaces of
other libraries. For example, the iterator support allows to be used
with  <algorithm> and other similar libraries.

The updating method cv_move is similar to cv_init. It will
update the particle's value by adding Normal random variates. However, as we
see above, each call to cv_init causes a log_weight vector being
allocated. Its size does not change between iterations. So it can be viewed as
some resource of cv_move and it is natural to use a class object to
manage it. Here is the implementation of cv_move,

class cv_move : public MoveSEQ<cv>

public :

void pre_processor (std::size_t iter, Particle<cv>    particle)

std::size_t move_state (std::size_t iter, SingleParticle<cv> sp)

const double sd_pos = std::sqrt(0.02);
const double sd_vel = std::sqrt(0.001);
const double delta = 0.1;
cxx11::normal_distribution<double> norm_pos(0, sd_pos);
cxx11::normal_distribution<double> norm_vel(0, sd_vel);

sp.state(0) += norm_pos(sp.rng()) + delta * sp.state(2);
sp.state(1) += norm_pos(sp.rng()) + delta * sp.state(3);
sp.state(2) += norm_vel(sp.rng());
sp.state(3) += norm_vel(sp.rng());

incw_[sp.id()] = sp.particle().value().log_likelihood(iter, sp.id());

return 0;

void post_processor (std::size_t iter, Particle<cv>    particle)
particle.weight_set().add_log_weight(   incw_[0]);

private :

std::vector<double> incw_;

First before calling any move_state, the pre_processor will be
called, as described in Section <ref>. At this step, we will
resize the vector used for storing incremental weights. After the first
resizing, subsequent calls to resize will only cause reallocation if
the size changed. In our example, the size of the particle system is fixed, so
we don't need to worry about excessive dynamic memory allocations. The
move_state member function moves each particle according to our model.
And after move_state is called for each particle, the
post_processor will be called and we simply add the logarithm
incremental weights.

For each particle, we want to monitor the $\xpos$ and $\ypos$ parameters and
get the importance sampling estimates. To extract the two values from a
particle, we can implement the following function

void cv_monitor_state (std::size_t iter, std::size_t dim,
ConstSingleParticle<cv> csp, double *res)

assert(dim <= Dim);
for (std::size_t d = 0; d != dim; ++d)
res[d] = csp.state(d);

and in the main function we construct a monitor by

MonitorEvalAdapter<cv, MonitorEvalSEQ> cv_est(cv_monitor_state);

and add it to the sampler through

sampler.monitor("pos", 2, cv_est);

If later we decided to monitor all states, we only need to change the 2
in the above line to Dim.

After we implemented all the above, compiled and ran the program, a file
called pf.est was written by the following statement in the main

est << sampler << std::endl;

The output file contains the , resampling, importance sampling estimates
and other informations in a table. This file can be read by most statistical
softwares for further analysis. For instance, we can process this file with
, and get the plot of the estimates of positions against the
observations, as shown in Figure <ref>.

Observations and estimates of a simple particle filter.

§.§ Bayesian Modeling of Gaussian mixture model

§.§.§ The model and algorithm

Since Richardson and Green, 1997, the Gaussian mixture model () has
provided a canonical example of a model-order-determination problem. We use
the same model as in Del Moral et al., 2006 to illustrate the implementation
of this classical example in Monte Carlo literature. This model is also used
in Zhou et al., 2013 for demonstration of the use of in the context of
Bayesian model comparison, which provides more details of the following
setting. The model is as follows; data $\data = (y_1,\dots,y_n)$ are
independently and identically distributed as

\begin{equation*}
y_i|\theta_r \sim \sum_{j=1}^r \omega_j\rnorm(\mu_j,\lambda_j^{-1})
\end{equation*}

where $\rnorm(\mu_j,\lambda_j^{-1})$ denotes the Normal distribution with mean
$\mu_j$ and precision $\lambda_j$; $\theta_r =
(\mu_{1:r},\lambda_{1:r},\omega_{1:r})$ and $r$ is the number of components in
each model. The parameter space is thus $\Real^r\times\Real^{+r}\times S_r$
where $S_r = \{\omega_{1:r}:0\le\omega_j\le1; \sum_{j=1}^r\omega_j=1\}$ is the
standard $(r-1)$-simplex. The priors which are the same for each component are
taken to be $\mu_j\sim\rnorm(\xi,\kappa^{-1})$,
$\lambda_j\sim\rgamma(\nu,\chi)$ and $\omega_{1:r}\sim\rdir(\rho)$ where
$\rdir(\rho)$ is the symmetric Dirichlet distribution with parameter $\rho$
and $\rgamma(\nu,\chi)$ is the Gamma distribution with shape $\nu$ and scale
$\chi$. The prior parameters are set in the same manner as in
Richardson and Green, 1997; also see Zhou et al., 2013 for details. The data
is simulated from a four components model with $\mu_{1:4} = (-3, 0,3, 6)$, and
$\lambda_j =2$, $\omega_j = 0.25$, $j = 1,\dots,4$. Our interest is to
simulate the posterior distribution of models with $r$ components, denoted by
$M_r$ and obtaining the normalizing constant for the purpose of Bayesian model
comparison <cit.>.

Numerous strategies are possible to construct a sequence of distributions for
the purpose of sampling. One option is to use for each model $M_r$,
$r\in\{1,2,\dots\}$, the sequence $\{\pi_t\}_{t=0}^{T_r}$, defined by

\begin{equation}
\pi_t(\theta_r^t) \propto
\pi(\theta_r^t|M_r)p(\data|\theta_r^t,M_r)^{\alpha(t/T_r)}.
\label{eq:geometry}
\end{equation}

where the number of distribution, $T_r$, and the annealing schedule,
$\alpha:[0,1]\to[0,1]$ may be different for each model. This leads to
Algorithm <ref>.

For each model $M_r\in\calM$ perform the following algorithm.

Set $t\leftarrow0$.
Sample $\theta_r^{(i,t)}\sim\pi(\theta_r^{(i,t)}|M_r)$.
Weight $W_0^{(i)} \propto 1$.

Set $t\leftarrow t + 1$.
Weight $W_t^{(i)} \propto W_{t-1}^{(i)}
p(\data|\theta_r^{(i,t-1)},M_r)^{\alpha(t/T_r) - \alpha([t-1]/T_r)}$.
Apply resampling if necessary.
Sample $\theta_r^{(i,t)} \sim
K_t(\cdot|\theta_r^{(i,t-1)})$, a $\pi_t$-invariant kernel.

Repeat the Iteration step up to $t = T_r$.
algorithm for Bayesian modeling of Gaussian mixture

The kernel $K_t$ in Algorithm <ref> is constructed as a
three-blocks Metropolis random walk,

* Update $\mu_{1:r}$ through a Normal random walk.

* Update $\lambda_{1:r}$ through a Normal random walk on logarithm
scale, that is, on $\log\lambda_{j}$, $j = 1, \dots, r$.

* Update $\omega_{1:r}$ through a Normal random walk on logit scale,
that is, on $\omega_{j}/\omega_r$, $j = 1,\dots,r-1$.

The standard direct estimate of the normalizing constants
[Del Moral et al., 2006] can be obtained from the output of this algorithm,

\begin{equation}
\hat\lambda_{\text{\textsc{ds}}}^{T_r,N} =
\sum_{i=1}^N \frac{\pi(\theta_r^{(i,0)}|M_r)}{\nu(\theta_0^{(i,0)})} \times
\prod_{t=1}^{T_r} \sum_{i=1}^N W_{t-1}^{(i)}
p(\data|\theta_r^{(i,t)}M_r)^{\alpha(t/T_r) - \alpha([t-1]/T_r)}
\label{eq:smc-ds}
\end{equation}

where $W_{t-1}^{(i)}$ is the importance weight of sample $\theta_{t-1}^{(i)}$.

§.§.§ Path sampling for estimation of normalizing constants

As shown in Zhou et al., 2013 the estimation of the normalizing constant
associated with our sequence of distributions can also be achieved by a Monte
Carlo approximation to the path sampling formulation given by
Gelman and Meng, 1998, also known as thermodynamic integration or Ogata's
method. Given a parameter $\alpha$ which defines a family of distributions,
$\{p_{\alpha} = q_{\alpha} / Z_\alpha\}_{\alpha \in [0,1]}$ that move smoothly
from $p_0 = q_0 / Z_0$ to $p_1 = q_1 / Z_1$ as $\alpha$ increases from zero to
one, one can estimate the logarithm of the ratio of their normalizing
constants via a simple integral relationship,

\begin{equation}
\log\biggl( \frac{Z_1}{Z_0} \biggr) =
\int_{0}^{1} \Exp_\alpha \biggl[ \rnd{\log q_{\alpha}(\cdot)}{\alpha}
\biggr] \intd\alpha, \label{eq:path_identity}
\end{equation}

where $\Exp_\alpha$ denotes expectation under $p_\alpha$. The sequence of
distributions in the algorithm for this example can be interpreted as
belonging to such a family of distributions, with $\alpha = \alpha(t/T_r)$.

The sampler provides us with a set of weighted samples obtained from a
sequence of distributions suitable for approximating this integral. At each
time $t$ we can obtain an estimate of the expectation within the integral via
the usual importance sampling estimator, and this integral can then be
approximated via a trapezoidal integration. In summary, the path sampling
estimator of the ratio of normalizing constants $\lambda^{T_r} = \log(Z_1/Z_0)$
can be approximated by

\begin{equation}
\hat\lambda_{\text{\textsc{ps}}}^{T_r,N} = \sum_{t=1}^{T_r}
\frac{1}{2}(\alpha_t - \alpha_{t - 1})(U_t^N + U_{t-1}^N)
\label{eq:path_est}
\end{equation}

\begin{equation}
U_t^N = \sum_{i=1}^N
W_t^{(i)} \rnd{\log q_{\alpha}(X_t^{(i)})}{\alpha}\Bigm|_{\alpha = \alpha_t}
\label{eq:path_import}
\end{equation}

§.§.§ Implementations

In this example we will implement the following classes.

* gmm_param is a class that abstracts the parameters of the
model, $\theta_r = (\mu_{1:r},\lambda_{1:r},\omega_{1:r})$.

* gmm_state is the value collection class.

* gmm_init is a class that implements operations used to
initialize the sampler.

* gmm_move_smc is a class that implements operations used to
update the weights as well as selecting the random walk proposal scales
and distribution parameter $\alpha(t/T_r)$.

* gmm_move_mu, gmm_move_lambda and gmm_move_weight
are classes that implement the random walks, each for one of the three

* gmm_path is a class that implements monitors for the path
sampling estimator. This class is similar to the importance sampling
monitor introduced before. It is to be used with
Sampler<gmm_state>::path_sampling. Its interface requirement will
be documented later.

* gmm_alpha_linear and gmm_alpha_prior are classes that
implement two of the many possible annealing schemes, $\alpha(t/T_r) =
t/T_r$ (linear) and $\alpha(t/T_r) = (t/T_r)^p$, $p > 1$ (prior).

* And last, the main function, which configure, initialize and
iterate the samplers.

This example is considerably more complicated than the last one. Instead of
documenting all the implementation details, for many classes we will only
show the interfaces. In most cases, the implementations are straightforward as
they are either data member accessors or straight translation of mathematical
formulations. For member functions with more complex structures, detailed
explanation will be given. Interested readers can see the source code for more

Later we will build both sequential and parallelized samplers. A few
configuration macros will be defined at compile time. For example, the
sequential sampler is compiled with the following header and macros,

#include <vsmc/smp/backend_seq.hpp>
#define BASE_SATE StateSEQ
#define BASE_INIT InitializeSEQ
#define BASE_MOVE MoveSEQ
#define BASE_PATH PathEvalSEQ

The definitions of these macros will be changed at compile time to build
parallelized samplers. For example, when using parallelization, the
header backend_omp.hpp will be used instead of backend_seq.hpp;
and StateSEQ will be changed to StateOMP along with similar
changes to the other macros. In the distributed source code, this is
configured by the build system.

Again, we first introduce the main function. The required headers are
the same as the last particle filter example in addition to the backend
headers as described above. The following variables used in the main
function will be set by the user input.

int ParticleNum;
int AnnealingScheme;
int PriorPower;
int CompNum;
std::string DataFile;

In the main function, we will create objects that set the distribution
parameter $\alpha(t/T_r)$ at each iteration according to the user input of
AnnealingScheme. Below is the main function. Note that some
source code of operations which set the parameters above are omitted.

int main ()

gmm_move_smc::alpha_setter_type alpha_setter;
if (AnnealingScheme == 1)
alpha_setter = gmm_alpha_linear(IterNum);
if (AnnealingScheme == 2)
alpha_setter = gmm_alpha_prior(IterNum, PriorPower);

Sampler<gmm_state> sampler(ParticleNum);
.move(gmm_move_smc(alpha_setter), false);
.mcmc(gmm_move_mu(), false);
.mcmc(gmm_move_lambda(), true);
.mcmc(gmm_move_weight(), true);
.initialize((void *) DataFile.c_str());

double ds = sampler.particle().value().nc().log_zconst();
double ps = sampler.path().log_zconst();
std::cout << "Standard estimate :      " << ds << std::endl;
std::cout << "Path sampling estimate : " << ps << std::endl;

return 0;

The sampler first sets the number of components and allocate memory through
member function comp_num of gmm_state. Then it sets the
initialization and updating methods. Before possible resampling, a
gmm_move_smc object is added. After that, three Metropolis random walks
are appended. In addition, we add a gmm_path object to calculate the
path sampling integration. Then we initialize and iterate the sampler and get
the normalizing constant estimates.

It is obvious that the parameter class gmm_param need to store the
parameters $(\mu_{1:r},\lambda_{1:r},\omega_{1:r})$. We also associate with
each particle its log-likelihood and log-prior. Here is the definition of the
gmm_param class. We omitted definitions of some data access member

class gmm_param

public :

void comp_num (std::size_t num);
void save_old ();

double log_prior () const return log_prior_;
double    log_prior () return log_prior_;

double log_likelihood () const return log_likelihood_;
double    log_likelihood () return log_likelihood_;

int mh_reject_mu (double p, double u);
int mh_reject_lambda (double p, double u);
int mh_reject_weight (double p, double u);
int mh_reject_common (double p, double u);

double log_lambda_diff () const;
double logit_weight_diff () const;

void update_log_lambda ();

private :

std::size_t comp_num_;
double log_prior_, log_prior_old_, log_likelihood_, log_likelihood_old_;
std::vector<double> mu_, mu_old_;
std::vector<double> lambda_, lambda_old_;
std::vector<double> weight_, weight_old_;
std::vector<double> log_lambda_;

The comp_num member function allocate the memory for a given number of
components. The save_old member function save the current particle
states. It is used before the states are updated with the random walk
proposals, as we will see later when we implement the gmm_move_mu. The
mh_reject_mu member function accept the Metropolis acceptance
probability $p$ and a uniform $(0,1]$ random variate, say $u$; it rejects the
proposed change if $p < u$, and restore the particle state of the parameters
$\mu_{1:r}$ by those values saved by save_old. The member functions
mh_reject_lambda and mh_reject_weight do the same for the other
two set of parameters. All these three also call the mh_reject_common
which restore the stored log-likelihood and log-prior values. The use of these
member functions will be seen in the implementation of gmm_move_mu, in
the context of which their own implementation become obvious. Other member
functions provide some useful computations such as the logarithm of the
$\lambda_{1:r}$. They are used when compute the log-likelihood.

The class gmm_state contains some properties common to all particles,
such as the data and the distribution parameter $\alpha(t/T_r)$. Also, we will
have it to record the logarithm of the ratio of normalizing constants, using
the NormalizingConstant class. We will see how to update this variable
at each iteration in the implementation of gmm_move_smc. The prior
parameters are also stored in the value collection. Here is the definition of
this value collection class. Again, we omitted some data access member

class gmm_state : public BASE_STATE<StateMatrix<RowMajor, 1, gmm_param> >

public :

NormalizingConstant    nc () return nc_;
const NormalizingConstant    nc () const return nc_;

void alpha (double a)

a = a < 1 ? a : 1;
a = a > 0 ? a : 0;
if (a == 0)
alpha_inc_ = 0;
alpha_ = 0;
alpha_inc_ = a - alpha_;
alpha_ = a;

void comp_num (std::size_t num)

comp_num_ = num;
for (size_type i = 0; i != this->size(); ++i)
this->state(i, 0).comp_num(num);

double update_log_prior (gmm_param    param) const;

double update_log_likelihood (gmm_param    param) const

static const double log2pi = 1.8378770664093455; // log(2pi)
double ll = -0.5 * obs_.size() * log2pi;
for (std::size_t k = 0; k != obs_.size(); ++k)
double lli = 0;
for (std::size_t i = 0; i != param.comp_num(); ++i)
double resid = obs_[k] - param.mu(i);
lli += param.weight(i) * std::exp(
0.5 * param.log_lambda(i) -
0.5 * param.lambda(i) * resid * resid);

ll += std::log(lli);

return param.log_likelihood() = ll;

void read_data (const char *filename);

private :

NormalizingConstant nc_;
std::size_t comp_num_;
double alpha_, alpha_inc_;
double mu0_, sd0_, shape0_, scale0_;
double mu_sd_, lambda_sd_, weight_sd_;
std::vector<double> obs_;

The variable alpha_inc_ is $\Delta\alpha(t/T_r) = \alpha(t/T_r) -
\alpha((t-1)/T_r)$, which will be used when we update the weights. The
variable nc_ of type NormalizingConstant will be updated when
the weights are changed by gmm_move_smc and it will compute the
standard normalizing constant estimate
$\hat\lambda_{\text{\textsc{ds}}}^{T_r,N}$. The variables mu0_ and
sd0_ are the prior parameters of the means $\mu_{1:r}$. The variables
shape0_ and scale0_ are the prior parameters of the precisions
$\lambda_{1:r}$. The variables mu_sd_, lambda_sd_, and
weight_sd_ are the proposal scales of the three random walks,
respectively. The data access member functions of these variables are omitted
in the above source code snippet.

In the update_log_likelihood member function, the calculation is a
straightforward translation of the mathematical formulation. The
gmm_param::update_log_lambda member function is used before the loop,
which simply calculates $\log\lambda_j$ for $j = 1,\dots,r$, and stores their
values. The purpose is to avoid repeated computation of these quantities
inside the loop. When the function returns, it uses the mutable version of the
gmm_param::log_likelihood member function to update the log-likelihood
stored in the param object. This is the reason that the function is
named with a update prefix. As we will see later, whenever the
parameter values are updated, it will be followed by a call to
update_log_likelihood and update_log_prior, which is implemented
in a similar fashion. Therefore the value we get by calling
gmm_param::log_likelihood will always be “up-to-date” while no
repeated computation is involved. Surely there are other and possibly better
design choices. However, for this simple example, this design serves our
purpose well.

The initialization is implemented using the gmm_init class,

class gmm_init : public BASE_INIT<gmm_state, gmm_init>

public :

std::size_t initialize_state (vsmc::SingleParticle<gmm_state> sp)

const gmm_state    state = sp.particle().value();
gmm_param    param = sp.state(0);

vsmc::cxx11::normal_distribution<> rmu(
state.mu0(), state.sd0());
vsmc::cxx11::gamma_distribution<> rlambda(
state.shape0(), state.scale0());
vsmc::cxx11::gamma_distribution<> rweight(1, 1);

double sum = 0;
for (std::size_t i = 0; i != param.comp_num(); ++i)
param.mu(i)     = rmu(sp.rng());
param.lambda(i) = rlambda(sp.rng());
param.weight(i) = rweight(sp.rng());
sum += param.weight(i);

for (std::size_t i = 0; i != param.comp_num(); ++i)
param.weight(i) /= sum;

return 1;

void initialize_param (vsmc::Particle<gmm_state>    particle,
void *filename)

if (filename)
particle.value().read_data(static_cast<const char *>(filename));

The initialize_param member function is called before the
pre_processor, which is absent in this case and have the default
implementation which does nothing. And it processes the optional parameter of
Sampler::initialize, the file name of the data. The
initialize_state member function initialize the state values according
to the prior and update the log-prior and log-likelihood.

After initialization, at each iteration, gmm_move_smc class will
implement the updating of weights as well as selecting of the proposal scales
and the distribution parameter. For example, when using the linear annealing
scheme, we can implement a gmm_alpha_linear class as the following,

class gmm_alpha_linear

public :

gmm_alpha_linear (const std::size_t iter_num) : iter_num_(iter_num)

void operator() (std::size_t iter, Particle<gmm_state>    particle)
particle.value().alpha(static_cast<double>(iter) / iter_num_);

private :

std::size_t iter_num_;

It accepts the total number of iterations $T_r$ as an argument to its
constructor. And it implements an operator() that update the
distribution parameter $\alpha(t/T_r)$. The prior annealing scheme can be
implemented similarly. For simplicity and demonstration purpose, we only allow
gmm_move_smc to be configured with different annealing schemes, and
hard code the proposal scales. An industry strength design may make this class
a template with annealing scheme and proposal scales as policy template

class gmm_move_smc

public :

typedef cxx11::function<void (std::size_t, Particle<gmm_state>    )>

gmm_move_smc (const alpha_setter_type    alpha_setter) :

std::size_t operator() (std::size_t iter, Particle<gmm_state>    particle)

alpha_setter_(iter, particle);

double alpha = particle.value().alpha();
alpha = alpha < 0.02 ? 0.02 : alpha;
particle.value().mu_sd(0.15 / alpha);
particle.value().lambda_sd((1 + std::sqrt(1 / alpha)) * 0.15);
particle.value().weight_sd((1 + std::sqrt(1 / alpha)) * 0.2);

double coeff = particle.value().alpha_inc();
for (vsmc::Particle<gmm_state>::size_type i = 0;
i != particle.size(); ++i)
incw_[i] =
coeff * particle.value().state(i, 0).log_likelihood();

particle.value().nc().add_log_weight(   incw_[0],
particle.weight_set().add_log_weight(   incw_[0]);

return 0;

private :

alpha_setter_type alpha_setter_;
std::vector<double> incw_;
std::vector<double> weight_;

Note that, cxx11::function is an alias to either std::function
or boost::function, depending on the value of the configuration macro
VSMC_HAS_CXX11LIB_FUNCTIONAL. Objects of this class type can be added
to a sampler as a move. The operator() satisfies the interface
requirement of the core module. First it uses alpha_setter_ to set the
distribution parameter $\alpha(t/T_r)$. Second, it sets the proposal scales
for the three Metropolis random walks according to the current value of
$\alpha$. Then it computes the unnormalized incremental weights. The
NormalizingConstnat class has member function add_log_weight,
which is not unlike the one with the same name in WeightSet. It accepts
the logarithm of the incremental weights and a WeightSet object. The
standard normalizing constant estimates will be computed using these values.
The last, we also modify the WeightSet type object itself by adding the
logarithm of the incremental weights.

At each iteration, ramdom walks are also perfoemd. The implementations of the
random walks are straightforward. Below is the implementation of the random
walk on the mean parameters. The random walks on the other parameters are

class gmm_move_mu : public BASE_MOVE<gmm_state, gmm_move_mu>

public :

std::size_t move_state (std::size_t iter, SingleParticle<gmm_state> sp)

const gmm_state    state = sp.particle().value();
gmm_param    param = sp.state(0);

cxx11::normal_distribution<> rmu(0, state.mu_sd());
cxx11::uniform_real_distribution<> runif(0, 1);

double p =
param.log_prior() + state.alpha() * param.log_likelihood();
for (std::size_t i = 0; i != param.comp_num(); ++i)
param.mu(i) += rmu(sp.rng());
p = state.update_log_prior(param) +
state.alpha() * state.update_log_likelihood(param) - p;
double u = std::log(runif(sp.rng()));

return param.mh_reject_mu(p, u);

First we save the logarithm of the value of target density computed using the
old values in p, the acceptance probability. And then we call
gmm_param::save_old to save the old values. Next we update each
parameter with a proposed Normal random variates and compute the new log-prior
and the log-likelihood as well as the new value of the target density. Then we
reject it according to the Metropolis algorithm, as implemented in
gmm_param, which manages both the current states as well as the backup.

Last we need to monitor certain quantities for interference purpose. Recall
that, in the main function we used
sampler.path_sampling(gmm_path()) to set the monitoring of path
sampling integrands. The path_sampling member function requires a
callable objects with the following signature,

double path_eval (std::size_t iter, const Particle<T>    , double *res);

The input parameter iter is the iteration number, the value of $t$ in
Equation <ref>. The return value shall be the value of $\alpha_t$.
The output parameter res shall store the array of values of

\begin{equation*}
\rnd{\log q_{\alpha}(X_t^{(i)})}{\alpha}\Bigm|_{\alpha = \alpha_t}.
\end{equation*}

Our implementation of gmm_path is a sub-class of an module base
class, which provides an operator() that satisfies the above interface
requirement. Its usage is similar to the MoveSEQ template introduced in
Section <ref>.

The path sampling integrands under this geometry annealing scheme are simply
the log-likelihood. Therefore the implementation of gmm_path class is
rather simple,

class gmm_path : public BASE_PATH<gmm_state, gmm_path>

public :

double path_state (std::size_t, ConstSingleParticle<gmm_state> sp)
return sp.state(0).log_likelihood();

double path_grid (std::size_t, const Particle<gmm_state>    particle)
return particle.value().alpha();

§.§.§ Results

After compiling and running the algorithm, the results were consistent with
those reported in Del Moral et al., 2006. For a more in depth analyze of the
methodologies, extensions and the results see Zhou et al., 2013.

§.§.§ Extending the
implementation using

's module assumes that identical samplers are constructed on each
node, with possible different number of particles to accommodate the
difference in capacities among nodes. To extend the above implementation
for use with , first at the beginning of the main function, we add
the following,

MPIEnvironment env(argc, argv);

to initialize the environment. When the object env is destroyed at
the exit of the main function, the environment is finalized.
Second, we need to replace base value collection class template with
StateMPI. So now gmm_state is declared as the following,

class gmm_state :
public StateMPI<BASE_STATE<StateMatrix<RowMajor, 1, gmm_param> > >;

The implementation is exactly the same as before. Third, the gmm_param
class now needs to be transferable using . Unlike the situations, a
simple copy constructor is not enough. uses  library, and
thus one only needs to write a serialize member function for
gmm_param such that the data can be serialized into bytes. See
documents of  and serialization libraries for details. In
summary, the following member function accept an Archive object as
input, and it can perform a store or a load operation based on the
Archive type. In a load operation, the Archive object is like an
input stream and in a store operation, it is like an output stream.

template <typename Archive>
void serialize (Archive    ar, const unsigned)

int num = comp_num_;
ar     num;

ar     log_prior_;
ar     log_likelihood_;
for (std::size_t i = 0; i != comp_num_; ++i)
ar     mu_[i];
ar     lambda_[i];
ar     weight_[i];
ar     log_lambda_[i];

Fourth, after user input of the sampler parameters, we need to sync them with
all nodes. For example, for the ParticleNum parameter,

boost::mpi::communicator World;
boost::mpi::broadcast(World, ParticleNum, 0);

Last, any importance sampling estimates that are computed on each node, need
to be combined into final results. For example, the path sampling results are
now obtained through adding the results from each node together,

double ps_sum = 0;
boost::mpi::reduce(World, ps, ps_sum, std::plus<double>(), 0);
ps = ps_sum;

For the standard normalizing constant ratio estimator, we will replace
NormalizingConstant with NormalizingConstantMPI, which will
perform such and other tasks.

After these few lines of change, the sampler is now parallelized using and can be deployed to clusters and other distributed memory architecture. On
each node, the selected parallelization is used to perform
multi-threading parallelization locally. 's module will take care of
normalizing weights and other tasks.

§.§.§ Parallelization performance

One of the main motivation behind the creation of is to ease the
parallelization with different programming models. The same implementation can
be used to built different samplers based on what kind of parallel programming
model is supported on the users' platforms. In this section we compare the
performance of various parallel programming models and parallelization.

We consider five different implementations supported by 2013:
sequential, , , and  <thread>. The samplers
are compiled with

CXX=icpc -std=c++11 -gcc-name=gcc-4.7 -gxx-name=g++-4.7
CXXFLAGS=-O3 -xHost -fp-model precise           -DVSMC_HAS_CXX11LIB_FUNCTIONAL=1           -DVSMC_HAS_CXX11LIB_RANDOM=1

on a Ubuntu 12.10 workstation with an Xeon W3550 (3.06GHz, 4 cores, 8 hardware
threads through hyper-threading) . A four components model and $100$
iterations with a prior annealing scheme is used for all implementations. A
range of numbers of particles are tested, from $2^3$ to $2^{17}$.

For different number of particles, the wall clock time and speedup are shown
in Figure <ref>. For $10^4$ or more particles, the
differences are minimal among all the programming models. They all have
roughly 550% speedup. With smaller number of particles, 's parallelization is less efficient than other industry strength programming
models. However, with $1000$ or more particles, which is less than typical
applications, the difference is not very significant.

Performance of implementations of Bayesian modeling for
Gaussian mixture model (Linux; Xeon W3550, 3.06GHz, 4 cores, 8 threads).

implementations are also compared on the same workstation, which also
has an NVIDIA Quadro 2000 graphic card. programs can be compiled to
run on both s and s. For implementation, there are and
platforms. We use the implementation as a baseline for comparison.
The same implementation are used for all the and runtimes.
Therefore they are not particularly optimized for any of them. For the implementation, in addition to double precision, we also tested a single
precision configuration.  Unlike modern s, which have the same
performance for double and single precision floating point operations (unless
instructions are used, which can have at most a speedup by a factor of
2), s penalize double precision performance heavily.

For different number of particles, the wall clock time and speed up are
plotted in Figure <ref>. With smaller number of particles,
the implementations have a high overhead when compared to the implementation. With a large number of particles, has a similar
performance as the implementation. is about 40% faster than the
implementation. This is due to more efficient vectorization and compiler
optimizations. The double precision performance of the NVIDIA has a 220%
speedup and the single precision performance has near 1600% speedup. As a
rough reference for the expected performance gain, the has a theoretical
peak performance of 24.48 GFLOPS. The has a theoretical peak performance
of 60 GFLOPS in double precision and 480 GFLOPS in single precision. This
represents 245% and 1960% speedup compared to the , respectively.

Performance of implementations of Bayesian modeling for
Gaussian mixture model (Linux; Xeon W3550 , 3.06GHz, 4 cores, 8
threads; NVIDIA Quadro 2000).

It is widely believed that programming is tedious and hard. However,
provides facilities to manage platforms and devices as well as
common operations. Limited by the scope of this paper, the implementation (distributed with the source) is not documented in this
paper. Overall the implementation has about 800 lines including both
host and device code. It is not an enormous increase in effort when compared
to the 500 lines implementation. Less than doubling the code base but
gaining more than 15 times performance speedup, we consider the programming
effort is relatively small.

§ DISCUSSION

This paper introduced a template library intended for implementing
general algorithms and constructing parallel samplers with different
programming models. While it is possible to implement many realistic
application with the presented framework, some technical proficiency is still
required to implement some problem specific part of the algorithms. Some basic
knowledge of in general and how to use a template library are also

It is shown that with the presented framework it is possible to implement
parallelized, scalable samplers in an efficient and reusable way. The
performance of some common parallel programming models are compared using an

Some future work may worth the effort to ease the implementation of algorithms further. However, there is a balance between performance,
flexibility and the ease of use. aims to be developer-friendly and to
provide users as much control as possible for all performance related aspects.
For a -like interface, users may be interested in other software such as
. In addition provides a user friendly and high performance
alternative with a focus on state-space models. Compared with these recent
developments, is less accessible to those with little or no knowledge of
. However, for researchers with expertise in and template generic
programming in particular, provides a framework within which
potential superior performance can be obtained and greater flexibility and
extensibility are possible.

[Advanced Micro Devices, Inc., 2012]
Advanced Micro Devices, Inc (2012).
AMD Accelerated Parallel Processing OpenCL
Programming Guide.
Version 2.8,

[Cappé et al., 2007]
Cappé O, Godsill SJ, Moulines E (2007).
An Overview of Existing Methods and Recent Advances in
Sequential Monte Carlo.
Proceedings of the IEEE, 95(5), 899–924.

[Caron et al., 2012]
Caron F, Todeschini A, Legrand P, Moral PD (2012).
BiiPS: Bayesian Inference with Interacting Particle
Version 0.7.2,

[Dawes et al., 2013]
Dawes B, et al. (2013).
Boost – C++ Libraries.
Version 1.53, <http://www.boost.org/>.

[Del Moral et al., 2006]
Del Moral P, Doucet A, Jasra A (2006a).
Sequential Monte Carlo Methods for Bayesian Computation.
In Bayesian Statistics 8. Oxford University Press.

[Del Moral et al., 2006]
Del Moral P, Doucet A, Jasra A (2006b).
Sequential Monte Carlo Samplers.
Journal of Royal Statistical Society B, 68(3),

[Douc et al., 2005]
Douc R, Cappé O, Moulines E (2005).
Comparison of Resampling Schemes for Particle Filtering.
In Proceedings of the 4th International Symposium on Imange and
Signal Processing and Analysis, pp. 1–6.

[Doucet and Johansen, 2011]
Doucet A, Johansen AM (2011).
A Tutorial on Particle Filtering and Smoothing: Fifteen
Years Later.
In The Oxford Handbook of Non-linear Filtering. Oxford
University Press.

[Gelman and Meng, 1998]
Gelman A, Meng XL (1998).
Simulating Normalizing Constants: From Importance Sampling
to Bridge Sampling to Path Sampling.
Statistical Science, 13(2), 163–185.

[GNU Project, 2013]
GNU Project (2013).
GCC – GNU Compiler Collection.
Version 4.8.1, <http://gcc.gnu.org>.

[Intel Cooperation, 2011]
Intel Cooperation (2011).
Intel Cilk Plus Language Specification.
Version 1.1, <http://cilkplus.org/>.

[Intel Cooperation, 2013]
Intel Cooperation (2013a).
Intel C++ Composer XE.
Version 13.1,

[Intel Cooperation, 2013]
Intel Cooperation (2013b).
Intel SDK for OpenCL Applications 2013.

[Intel Cooperation, 2013]
Intel Cooperation (2013c).
Intel Threading Building Blocks.
Version 4.1, <http://threadingbuildingblocks.org/>.

[Johansen, 2009]
Johansen AM (2009).
SMCTC: Sequential Monte Carlo in C++.
Journal of Statistical Software, 30(6), 1–41.

[Kronos OpenCL Working Group, 2012]
Kronos OpenCL Working Group (2012).
The OpenCL Specification.
Version 1.2, <http://www.khronos.org/opencl/>.

[Lee et al., 2010]
Lee A, Yau C, Giles MB, Doucet A, Holmes CC (2010).
On the Utility of Graphics Cards to Perform Massively
Parallel Simulation of Advanced Monte Carlo Methods.
Journal of Computational and Graphical Statistics,
19(4), 769–789.

[Liu and Chen, 1998]
Liu JS, Chen R (1998).
Sequential Monte Carlo Methods for Dynamic Systems.
Journal of the American Statistical Association,
93(443), 1032–1044.

[Martin and Hoffman, 2010]
Martin K, Hoffman B (2010).
Mastering CMake.
Version 2.8, <http://www.cmake.org>.

[Message Passing Interface Forum, 2012]
Message Passing Interface Forum (2012).
MPI: A Message-Passing Interface Starndard.
Version 3.0, <http://www.mpi-forum.org/>.

[Microsoft Cooperation, 2012]
Microsoft Cooperation (2012).
Microsoft Visual C++ 2012.

[Murray, 2013]
Murray LM (2013).
Bayesain State-Space Modelling on High-Performance Hardware
Using LibBi.
arXiv, pp. 1–28.

[Neal, 2001]
Neal RM (2001).
Annealed Importance Sampling.
Statistics and Computing, 11(2), 125–139.

[OpenMP Architecture Review Board, 2011]
OpenMP Architecture Review Board (2011).
OpenMP Application Program Interface.
Version 3.1, <http://www.openmp.org>.

[Peters, 2005]
Peters GW (2005).
Topics in Sequential Monte Carlo Samplers.
Master's thesis, Department of Engineering, University of Cambridge.

[R Core Team, 2013]
R Core Team (2013).
R: A Language and Environment for Statistical
R Foundation for Statistical Computing, Vienna, Austria.
Version 2.15.2, <http://www.r-project.org>.

[Richardson and Green, 1997]
Richardson S, Green PJ (1997).
On Bayesian Analysis of Mixtures with an Unknown Number of
Journal of Royal Statistical Society B, 59(4),

[Robert, 2007]
Robert CP (2007).
The Bayesian Choice: From Decision-theoretic Foundations to
Computational Implementation.
2nd edition. Springer-Verlag, New York.

[Salmon et al., 2011]
Salmon JK, Moraes MA, Dror RO, Shaw DE (2011).
Parallel Random Numbers: As Easy as 1, 2, 3.

[Spiegelhalter et al., 1996]
Spiegelhalter D, Thomas A, Best N (1996).
BUGS – Bayesian Inference Using Gibbs Sampling.
Version 0.5, <http://www.mrc-bsu.cam.ac.uk/bugs/>.

[The LLVM Developer Group, 2013]
The LLVM Developer Group (2013a).
clang: A C Language Family Frontend for LLVM.
Version 3.2, <http://clang.llvm.org/>.

[The LLVM Developer Group, 2013]
The LLVM Developer Group (2013b).
“libc++” C++ Standard Library.
Version 1.1, <http://libcxx.llvm.org>.

[Torvalds et al., 2013]
Torvalds L, et al. (2013).
Git – Distributed Source Code Management.
Version 1.8.3.1, <http://git-scm>.

[van Heesch, 2013]
van Heesch D (2013).
Doxygen – Generating Documentation from Source Code.
Version 1.8.4,
<http://www.stack.nl/ dimitri/doxygen/index.html>.

[Zhou et al., 2013]
Zhou Y, Johansen AM, Aston JAD (2013).
Towards Automatic Model Comparison: Adaptive Sequential
Monte Carlo Approach.
arXiv, pp. 1–33.