Vector operations for accelerating expensive Bayesian computations -- a tutorial guide

TL;DR

This paper demonstrates how SIMD operations can significantly accelerate Bayesian computations on modern CPUs, providing practical techniques for leveraging parallel processing to improve performance.

Contribution

It introduces the use of SIMD for Bayesian applications, showing practical implementation and performance gains using standard programming libraries.

Findings

Up to 6x performance improvement with SIMD

SIMD benefits are multiplicative with multi-core processing

Provides practical techniques for SIMD exploitation in Bayesian computing

Abstract

Many applications in Bayesian statistics are extremely computationally intensive. However, they are often inherently parallel, making them prime targets for modern massively parallel processors. Multi-core and distributed computing is widely applied in the Bayesian community, however, very little attention has been given to fine-grain parallelisation using single instruction multiple data (SIMD) operations that are available on most modern commodity CPUs and is the basis of GPGPU computing. In this work, we practically demonstrate, using standard programming libraries, the utility of the SIMD approach for several topical Bayesian applications. We show that SIMD can improve the floating point arithmetic performance resulting in up to $6\times$ improvement in serial algorithm performance. Importantly, these improvements are multiplicative to any gains achieved through multi-core processing. We illustrate the potential of SIMD for accelerating Bayesian computations and provide the reader with techniques for exploiting modern massively parallel processing environments using standard tools.

Tables

Table 1. Table A.1 : Computational performance R and C implementations of prior predictive sampling for the toggle switch model.

Number of threads	E5-2680v3 Runtimes (Seconds)				6140 Runtimes (Seconds)
$P$	R Naïve	R Optimise	C Scalar	C AVX2	C Scalar	C AVX512
1	298,041	10,502	2,632	577	2,395	380
2	148,677	5,270	1,343	315	1,223	213
4	76,520	2,663	760	179	631	128
8	42,387	1,349	405	111	336	86
16	16,961	681	223	78	190	69

Table 2. Table A.2 : Computational performance for the weak informativity test.

Number of threads	E5-2680v3 Runtimes (Seconds)		6140 Runtimes (Seconds)
$P$	C Scalar	C AVX2	C Scalar	C AVX512
1	2,871	1,389	2,952	1,325
2	1,429	701	1,463	665
4	714	350	730	332
8	359	178	367	167
16	183	91	187	90

Table 3. Table A.3 : Computational performance for BEGE parameter inference.

Number of threads	E5-2680v3 Runtimes (Seconds)		6140 Runtimes (Seconds)
$P$	C Scalar	C AVX2	C Scalar	C AVX512
1	25,618	13,327	16,780	6,793
2	13,469	8,230	11,895	8,334
4	9,358	7,097	8,446	5,526
8	7,180	4,272	5,284	3,678
16	5,131	3,115	4,215	3,156

Table 4. Table B.1 : Computational performance of computing prior predictive samples under the tuberculosis model using an Intel Xeon E5-2680v3 (Haswell) processor. Parallel times are based on 4 4 4 -way multithreading and SIMD times using 256 256 256 bit vector operations.

Number of samples	Runtime (Seconds) [Speedup factor]
$N$	sequential	sequential+SIMD	parallel	parallel+SIMD
$500$	$46.7$ [$1.0\times$]	$25.5$ [$1.8\times$]	$64.5$ [$0.7\times$]	$36.0$ [$1.3\times$]
$1,000$	$230.2$ [$1.0\times$]	$127.7$ [$1.8\times$]	$72.1$ [$3.2\times$]	$40.1$ [$5.8\times$]
$2,000$	$385.0$ [$1.0\times$]	$212.7$ [$1.8\times$]	$91.8$ [$4.2\times$]	$50.8$ [$7.6\times$]
$4,000$	$657.4$ [$1.0\times$]	$363.1$ [$1.8\times$]	$229.5$ [$2.9\times$]	$127.7$ [$5.2\times$]
$8,000$	$1,175.6$ [$1.0\times$]	$649.2$ [$1.8\times$]	$386.7$ [$3.0\times$]	$212.5$ [$5.5\times$]
$16,000$	$2,601.7$ [$1.0\times$]	$1,449.4$ [$1.8\times$]	$697.3$ [$3.7\times$]	$389.6$ [$6.7\times$]

Table 5. Table B.2 : Computational performance of computing prior predictive samples under the tuberculosis model using an Intel Xeon Gold 6140 (Skylake) processor. Parallel times are based on 4 4 4 -way multithreading and SIMD times using 512 512 512 bit vector operations.

Number of samples	Runtime (Seconds) [Speedup factor]
$N$	sequential	sequential+SIMD	parallel	parallel+SIMD
$500$	$15.9$ [$1.0\times$]	$11.5$ [$1.4\times$]	$23.7$ [$0.7\times$]	$17.9$ [$0.9\times$]
$1,000$	$83.9$ [$1.0\times$]	$61.5$ [$1.4\times$]	$26.3$ [$3.2\times$]	$19.3$ [$4.4\times$]
$2,000$	$138.2$ [$1.0\times$]	$101.2$ [$1.4\times$]	$32.8$ [$4.2\times$]	$24.0$ [$5.8\times$]
$4,000$	$233.5$ [$1.0\times$]	$170.7$ [$1.4\times$]	$83.8$ [$2.8\times$]	$61.5$ [$3.8\times$]
$8,000$	$412.9$ [$1.0\times$]	$301.9$ [$1.4\times$]	$138.1$ [$3.0\times$]	$101.1$ [$4.1\times$]
$16,000$	$922.5$ [$1.0\times$]	$675.2$ [$1.4\times$]	$249.7$ [$3.7\times$]	$183.0$ [$5.0\times$]

Equations

$$s\le \frac{C_S+C_P}{C_S+C_P/P},$$

$$\pi (\mathcal{D})={\int }_{\mathbf{\Theta }}s(\mathcal{D};\mathbf{\theta })\pi (\mathbf{\theta })\text{d}\mathbf{\theta },$$

$$\pi (\mathbf{\theta }\mid {\mathcal{D}}_{\text{obs}})\approx \pi (\mathbf{\theta }\mid \rho (S({\mathcal{D}}_{\text{obs}}),S(\mathcal{D}))\le ϵ)=\frac{\mathbb{P}\left(\rho (S({\mathcal{D}}_{\text{obs}}),S(\mathcal{D}))\le ϵ|\mathbf{\theta }\right)\pi (\mathbf{\theta })}{\pi ({\mathcal{D}}_{\text{obs}})}$$

$$\begin{array}{rl}\text{d}{u}_i(t)& =\left(\frac{{\alpha }_u}{1+v_i(t{)}^{{\beta }_u}}-(1+0.03u_i(t))\right)\text{d}t+\frac{1}{2}\text{d}{W}_{i,u}(t),\\ \text{d}{v}_i(t)& =\left(\frac{{\alpha }_v}{1+u_i(t{)}^{{\beta }_v}}-(1+0.03v_i(t))\right)\text{d}t+\frac{1}{2}\text{d}{W}_{i,v}(t),\end{array}$$

$$y_i=u_i(T)+\mu +\mu \sigma \frac{{\eta }_i}{u_i(T{)}^{\gamma }},i=1,2,\dots ,C,$$

$$\begin{array}{rl}{u}_i(t+h)& =u_i(t)+\Delta t\frac{{\alpha }_u}{1+v_i(t{)}^{{\beta }_u}}-\Delta t(1+0.03u_i(t))+0.5\sqrt{\Delta t}{\xi }_{i,u}(t),\\ v_i(t+h)& =v_i(t)+\Delta t\frac{{\alpha }_v}{1+u_i(t{)}^{{\beta }_v}}-\Delta t(1+0.03v_i(t))+0.5\sqrt{\Delta t}{\xi }_{i,v}(t),\end{array}$$

$$p({\mathcal{D}}_{\text{obs}})=\mathbb{P}\left(\frac{1}{\pi (\mathcal{D})}\ge \frac{1}{\pi ({\mathcal{D}}_{\text{obs}})}\right)=\mathbb{P}\left(\pi (\mathcal{D})\le \pi ({\mathcal{D}}_{\text{obs}})\right),$$

$$\pi (\mathcal{D}\mid \mathbf{\lambda })={\int }_{\mathbf{\Theta }}\mathcal{L}(\mathbf{\theta };\mathcal{D})\pi (\mathbf{\theta }\mid \mathbf{\lambda })\text{d}\mathbf{\theta }$$

$$p_{\mathbf{\lambda }}({\mathcal{D}}_{\text{obs}})=\mathbb{P}\left(\pi (\mathcal{D}\mid \mathbf{\lambda })\le \pi ({\mathcal{D}}_{\text{obs}}\mid \mathbf{\lambda })\right).$$

$$\mathbb{P}\left(p_{\mathbf{\lambda }}({\mathcal{D}}_0)\le \gamma \right)<\mathbb{P}\left(p_{{\mathbf{\lambda }}_0}({\mathcal{D}}_0)\le \gamma \right).$$

$$y_i\sim \text{Bin}\left(n_i,\frac{1}{1+e^{{\beta }_0+{\beta }_1{x}_i}}\right),i=1,2,\dots ,M,$$

$$\mathcal{L}(\mathbf{\theta };{\mathcal{D}}_{\text{obs}})=\prod _{i=1}^M\left(\genfrac{}{}{0ex}{}{n_i}{y_{\text{obs},i}}\right){\left(\frac{1}{1+e^{{\beta }_0+{\beta }_1{x}_i}}\right)}^{y_{\text{obs},i}}{\left(1-\frac{1}{1+e^{{\beta }_0+{\beta }_1{x}_i}}\right)}^{n_i-y_{\text{obs},i}}.$$

$$\begin{array}{rl}{r}_{t+1}=u_{t+1}+\mu ,& u_{t+1}={\sigma }_p{\omega }_{p,t+1}-{\sigma }_n{\omega }_{n,t+1},\\ {\omega }_{p,t+1}\sim \tilde{\Gamma }(p_t,1),& {\omega }_{n,t+1}\sim \tilde{\Gamma }(n_t,1),\end{array}$$

$$\begin{array}{rl}{p}_t& =p_0+{\rho }_p{p}_{t-1}+\frac{{\varphi }_p^{+}}{2{\sigma }_p^2}{u}_t^2\text{mathds}{1}_{[0,\infty )}\left(u_t\right)+\frac{{\varphi }_p^{-}}{2{\sigma }_p^2}{u}_t^2\text{mathds}{1}_{(-\infty ,0)}\left(u_t\right),\\ n_t& =n_0+{\rho }_n{n}_{t-1}+\frac{{\varphi }_n^{+}}{2{\sigma }_n^2}{u}_t^2\text{mathds}{1}_{[0,\infty )}\left(u_t\right)+\frac{{\varphi }_n^{-}}{2{\sigma }_n^2}{u}_t^2\text{mathds}{1}_{(-\infty ,0)}\left(u_t\right),\end{array}$$

$$\log \mathcal{L}(\mathbf{\theta };\mathcal{D})=\log \mathcal{L}(\mathbf{\theta };R_{\text{obs},T})=\sum _{t=1}^T\log \pi (r_{\text{obs},t}\mid r_{\text{obs},t-1},\mathbf{\theta }).$$

$$F_{\text{BEGE}}(u_t\mid {\sigma }_p,{\sigma }_n,p_{t-1},n_{t-1})={\int }_{-\infty }^{\infty }{G}_{\tilde{\Gamma }}({\omega }_{p,t}-u_t\mid n_{t-1},{\sigma }_n){\pi }_{\tilde{\Gamma }}({\omega }_{p,t}\mid p_{t-1},{\sigma }_p)\text{d}{\omega }_{p,t},$$

$$P(a,x)=\frac{1}{\Gamma (a)}{\int }_0^x{e}^{-t}{t}^{a-1}\text{d}t,$$

$$\begin{array}{r}{F}_{\text{BEGE}}(u_t\mid {\sigma }_p,{\sigma }_n,p_{t-1},n_{t-1})\approx \sum _{j=1}^{N_{\omega }}\{\left[1-P\left(n_{t-1},\frac{{\omega }_{p,t}^{j-1}-u_t+n_{t-1}{\sigma }_n}{{\sigma }_n}\right)\right]\\ \times \left[P\left(p_{t-1},\frac{{\omega }_{p,t}^j+p_{t-1}{\sigma }_p}{{\sigma }_p}\right)-P\left(p_{t-1},\frac{{\omega }_{p,t}^{j-1}+p_{t-1}{\sigma }_p}{{\sigma }_p}\right)\right]\},\end{array}$$

$$\begin{array}{rl}{P}^V(a,x_{1:V})=& \left[\left({\exp }^V\left(-x_{1:V}\right)\circ (x_{1:V}{)}^{\circ a}\right)\oslash \left(\Gamma (a){\mathbf{e}}_{1:V}\right)\right]\\ & \circ \left[\sum _{j=0}^{\infty }\left((x_{1:V}{)}^{\circ j}\right)\oslash \left((j+1){\mathbf{e}}_{1:V}\right)\right],\end{array}$$

$$N(t)=\sum _{i=1}^{G(t)}{X}_i(t),$$

$$a_b(X_i(t))=\alpha X_i(t),a_d(X_i(t))=\delta X_i(t),a_m(X_i(t))=\tau X_i(t),$$

Full text

1 Introduction

Practical applications in Bayesian statistics are computationally challenging since the posterior density is only known up to a normalising constant. Therefore, advanced Monte Carlo schemes are often required. Intractable likelihoods further compound these computational burdens [48]. In addition to posterior sampling, there are other computational challenges in Bayesian statistics. For example, it may be necessary to ensure priors are only weakly informative, and selected from a family of priors with a hyperparameter [16, 19].

The performance of computational methods and computer hardware continues to improve [26]. Most recent hardware improvements are due to increased parallel processing capacity. Implementation details for Monte Carlo algorithms are essential to fully exploit modern computational resources. Code optimisation for hardware acceleration is standard in the high performance computing (HPC) discipline, and Bayesian practitioners can benefit from these techniques [21].

Many Monte Carlo schemes are inherently parallel. For example, likelihood-free methods, such as approximate Bayesian computation (ABC; Sisson et al. 2018) require a large number of independent prior predictive samples that can be executed in parallel. However, more sophisticated samplers, such as Markov chain Monte Carlo (MCMC) [26, 38], sequential Monte Carlo (SMC) [14, 49] and multilevel Monte Carlo (MLMC) [30, 57] require more effort to efficiently parallelise [41] due to synchronisation and communication requirements. These overheads can inhibit scalability for multithreading on central processing units (CPUs).

General purpose graphics processing units (GPGPUs) are highly effective at accelerating advanced Monte Carlo schemes [31, 33]. GPGPUs make heavy use of single instruction multiple data (SIMD) [55], that is, instructions that operate on vectors element-wise. SIMD is also widely available in modern CPUs containing vector processing units (VPUs). By exploiting VPUs, CPU implementations can achieve performance boosts comparable with GPGPU implementations [34, 40]. Practical understanding of SIMD for modern CPUs is relevant to applied Bayesian analysis. We focus on code structures and algorithmic techniques for practitioners to harness VPUs available in most commodity CPUs, and provide example codes using both R and C222Example code is available from GitHub: https://github.com/davidwarne/Bayesian_SIMD_examples.

The paradigms of parallelism are distributed computing, multithreading, vectorisation, and pipelining [54, 55]. Each paradigm has a granularity that refers to the ratio of communication to computation for a parallel workload. We refer to distributed computing and multithreading as coarse-grained, whereas vectorisation and pipelining are fine-grained. Most exemplars of fine-grain parallelism in statistics are based on accelerators, such as, GPGPUs [33, 52], Intel Xeon Phis [29, 37], and custom co-processors [59]. However, accelerators are not practical for many practitioners, so we do not discuss them here. Instead we highlight parallelisation strategies for optimal utilisation of CPUs. In particular, we focus on the effective utilisation of vectorisation using SIMD.

In this paper, we demonstrate the utility of CPU-based SIMD for accelerating Bayesian inference. In Section 2 we introduce the principles of code optimisation and hardware acceleration using SIMD, then practically demonstrate these principles, in Section 3, through a tutorial that compares optimised R and C implementations of prior predictive sampling for ABC [48]. Finally, several topical case studies are presented in Section 4: the computation of Bayesian $p$-values for prior weak informativity tests [16], and parameter inference in econometrics [4]. We highlight algorithmic features that are suited for parallelism and demonstrate modifications for practical guidance. The clear demonstration of performance benefits using SIMD for CPUs is our main contribution. Optimised implementations are provided as supplementary material using the C language, OpenMP standard (version $\geq 4.5$), the Intel Math Kernel Library333https://software.intel.com/content/www/us/en/develop/tools/performance-libraries.html (MKL) (version $\geq 2018$), and the Intel C compiler 444https://software.intel.com/content/www/us/en/develop/tools/compilers/c-compilers.html icc (version $\geq 17.0.1$). Our guidelines are also directly applicable to the Julia language [6] and to external interfaces to pre-compiled C code, such as Matlab C-MEX and Rcpp combined with RcppXsimd (see Section 5).

2 An introduction to code optimisation

Here, we introduce CPU computing concepts essential for optimisation. We avoid technological details in favour of simple tools and rules-of-thumb for practitioners.

2.1 Vectorisation with SIMD

CPU cores are sequential555Although, they do not necessarily execute operations in the order defined by the users code. CPUs re-order operations as needed to optimise pipelining. This is called out-of-order execution and is a feature that is not readily available on accelerator architectures., performing a single operation, such as arithmetic or reading/writing data, per clock cycle. Without optimisations, the loop in Figure 1(a) requires four floating point addition operations. Modern CPU cores can perform scalar or vector arithmetic, thus Figure 1(a) takes only a single vector addition ($4\times$ speedup).

Vector arithmetic is implemented using VPUs that execute in SIMD. The VPU inputs are wider than the scalar arithmetic units. A VPU that accepts $256$ bit inputs can store eight single precision ($32$ bit) or four double precision ($64$ bit) floating point numbers. VPUs perform arithmetic element-wise on input vectors as shown in Figure 1(b). Acceleration through SIMD is multiplicative with gains achieved through multithreading. Using 18 threads and $512$ bit vectors of an Intel Xeon Gold 6140 CPU 666https://ark.intel.com/products/120485/Intel-Xeon-Gold-6140-Processor-24-75M-Cache-2-30-GHz-, could theoretically perform up to $18\times(512/64)=144$ times as many double precision operations per second than sequential processing with scalar operations [53].

2.2 Vectorisation and multithreading with OpenMP

VPUs may be accessed in a variety of ways. A simple approach is OpenMP, an open standard for automated parallel computing. OpenMP supports parallel computing through directives within standard C/C++ [55]. The directive #pragma omp simd (Figure 2(a)) guides the compiler to utilise VPUs.

The directive #pragma omp parallel tells the compiler to create a worker pool of threads within the next code block. Within a parallel block, the directive #pragma omp for before a loop causes iterations to be distributed across these threads. Figure 2(b) and (c) show the same loop as in Figure 2(a), but parallelised across cores. Determining when to use simd (Figure 2(a)), parallel and for (Figure 2(b)), and parallel for (Figure 2(c)) is key to performance. The following guidelines are useful:

Use simd (Figure 2(a)) when loop iterations are independent without conditionals. The loop body should only use arithmetic and standard mathematical functions (e.g. sin, cos, exp or pow). 2. 2.

Use parallel for (Figure 2(c)) when loop iterations are independent but conditionals cause iterations to perform different calculations. The loop body can be arbitrarily complex with calls to user defined functions. 3. 3.

A simd (Figure 2(a)) loop can be included within the body of a parallel for (Figure 2(c)) loop, but the converse is not true. 4. 4.

The full parallel and for (Figure 2(b)) construct can be useful to more explicitly control the behaviour of individual threads through the use of OpenMP functions.

The above guidelines are not all strict rules, for example, user defined functions may be called within a simd loop provided the function has been appropriately constructed and the declare simd pre-processor directive is used. Many other features are also available within OpenMP, and we refer the reader to the specifications777See https://www.openmp.org/specifications/ for full OpenMP specifications. for details. [54] and [55] also provide practical information on parallel computing for HPC. Alternatives to OpenMP are discussed in Section 5.

2.3 Memory access and alignment

To fully exploit SIMD, it is crucial to manage memory access. Commodity random access memory (RAM) bandwidth is around $20$ GB/sec, whereas the VPU can process floating point data at a rate of $600$ GB/sec. That means only $3\%$ VPU utilisation is possible for data retrieved from RAM. CPUs avoid this bottleneck with a hierarchy of memory caches, typically with three levels: L1, L2, and L3. Lower cache numbers correspond to higher bandwidth, but smaller capacity. L1 cache can keep the VPU close to $100\%$ utilised, but is typically less than $30$ KB. When data are requested for computation, the caches are tested in order. If the data are in L1 cache, then there is no memory access penalty. If not in L1 cache, then L2 is checked and so on. RAM is accessed only when the data are not in any cache, which is called a cache miss. For efficiency, data arrays should be accessed in patterns that minimise cache misses:

• •

Aim to reuse data soon after it was last accessed. If a partial calculation is pushed out of cache, then it can be faster to re-compute than access the result from RAM.

• •

Avoid random memory access patterns in favour of regular access a[i], then a[i+n], then a[i+2*n]888Access matrices according to the matrix storage format. For example, the row-major format is used in C/C++ (access row-by-row) and the column-major format is used in R (access column-by-column).. Ideally, n = 1 or is small to avoid cache churn.

Memory alignment is also important for SIMD. Memory is partitioned into blocks called cache lines. For many recent Intel CPUs, these are $64$ bytes in length999The CPU-Z utility (https://www.cpuid.com) is a useful tool to identify cache line sizes.. When the CPU loads data at an address, the entire cache line is loaded. When doing SIMD, the first byte of a data array should be aligned with the $64$byte cache line boundary, otherwise some vectors will cross cache lines. Performing SIMD on unaligned data incurs a performance penalty, or may fail completely for some older standards. In practice this means:

• •

Use either _mm_malloc (Intel) or memalign (GNU), in place of default malloc101010On most systems, malloc defaults is $8$byte or $16$byte..

• •

For static arrays, append __declspec(align(n)) to the array declaration111111For example, __declspec(align(64)) double A[100] is a $64$byte aligned array..

Efficient use of cache and memory access patterns is a large topic. See further discussion in [12] and [20].

2.4 Performance analysis

Code optimisation effort should be informed by performance analysis to ensure reductions in computation time are worth the increase in development time. It is common for this to be overlooked in practice [29, 34]. Some guidelines to assist are:

Use profiling and static code analysis software to obtain performance statistics on timing and memory bottlenecks, such as Intel VTune Amplifier121212https://software.intel.com/content/www/us/en/develop/tools/vtune-profiler.html and Intel Advisor131313https://software.intel.com/content/www/us/en/develop/tools/advisor.html. These tools are also essential to investigate and control compiler output. 2. 2.

Estimate the theoretical peak performance for your hardware. Compare with an estimate of the number of floating point operations your algorithm will perform, on average. Use this to estimate the minimum possible runtime [29]. 3. 3.

Most algorithms require some sequential operations, and not all calculations can be efficiently mapped to SIMD. Amdahl’s law [1] states the speedup factor, $s$, is bounded by

[TABLE]

where $C_{P}$ is the sequential runtime of code that can be parallelised, $C_{S}$ is the runtime of code that must remain sequential, and $P$ is the core count or vector length. Thus, $s=(C_{S}+C_{P})/C_{S}$ is the maximum speedup as $P\to\infty$.

2.5 A note on random number generation

When combining statistics with parallel computing it is important to consider initialisation and usage of random number generators (RNGs). Dealing with these aspects poorly can lead to invalid results. In particular, the use of different seed values (one per thread) for the same sequential RNG should be approached with caution. For example, sequences of Linear Congruential Generators can become correlated if a linear sequence of seeds is used [13].

Broadly there are three approaches to deal with this: (i) generate all required randomness serially and distribute among parallel processes; (ii) use a random number generator that can be split into independent sub-streams, for example a generator that supports Skip Ahead or Leapfrog methods; or (iii) use a sequence of random number generators that are statistically independent for the same seed value. We utilise the last option through the VSL_BRNG_MT2203 generator family (available within Intel MKL) that provides statistically independent Mersenne Twister generators. See [10] and [33] for details.

2.6 Summary

This introduction to vectorisation, multithreading, memory usage, code analysis, and RNGs is not specific to Bayesian statistics. However, the ideas are important to efficiently implement steps within the increasingly computationally expensive algorithms that form a core part of modern practical Bayesian techniques. In Sections 3 and 4 we demonstrate the application of these guidelines through a detailed tutorial and number of case studies in settings of direct interest to Bayesian practitioners.

3 A practical tutorial demonstration for R users

In this section, we provide a practical demonstration of the computational benefits of directly accessing CPU SIMD operations for ABC-based inference. We begin with R implementations, step through relevant optimisation within R, then demonstrate the computational advantages of direct SIMD access using C and OpenMP.

3.1 Prior predictive sampling for approximate Bayesian computation

ABC techniques are powerful for inference with intractable likelihood functions [48] and are routinely used to study complex stochastic models [46, 51]. The most accessible algorithm is ABC rejection sampling. Given observed data, $\mathcal{D}_{\text{obs}}$, the parameter prior distribution, $\pi(\boldsymbol{\theta})$, a discrepancy metric, $\rho$, and a vector of sufficient (or informative) summary statistics, $S(\mathcal{D}_{\text{obs}})$, then ABC rejection sampling generates approximate posterior samples by first generating artificial data from the prior predictive distribution,

[TABLE]

where $s(\mathcal{D};\boldsymbol{\theta})$ is the probability density of the data generation process for a fixed parameter vector $\boldsymbol{\theta}$ in parameter space $\boldsymbol{\Theta}$. Then a small proportion of samples are accepted to form an approximation to the posterior,

[TABLE]

with a sufficently small acceptance threshold $\epsilon$. This procedure is presented in algorithmic form in Algorithm 1.

In practice, ABC rejection sampling is rarely implemented in this direct manner (Algorithm 1 is serial, and produces a random number of candidate samples). Rather, it is common to generate a fixed number, $N$, of prior predictive joint samples, $(\boldsymbol{\theta}^{*},{\mathcal{D}}^{*})$. The acceptance threshold, $\epsilon$, is then selected a posteriori based on the empirical distribution of the discrepancy metric, $\rho$. See e.g. [17] and [58] for detailed reviews of Monte Carlo algorithms for ABC. ABC-based Monte Carlo estimators converge slowly in mean-square  [3], and therefore require substantial numbers of stochastic model simulations for reliable inference.

3.2 Example model: a genetic toggle switch

We consider a genetic toggle switch model [9]. Let $u_{i}(t)\geq 1$ and $v_{i}(t)\geq 1$ represent the expression levels of genes $u$ and $v$ at time $t$ for cells $i=1,2\ldots,C$. Gene expression evolves according to two coupled stochastic differential equations (SDEs)

[TABLE]

where parameters $\alpha_{u}$, $\beta_{u}$, $\alpha_{v}$ and $\beta_{v}$ define gene inhibition and $W_{i,u}(t)$ and $W_{i,v}(t)$ are independent Wiener processes. The observation process is

[TABLE]

where $T$ is the observation time, $\mu$, $\sigma$ and $\gamma$ control the error rate and $\eta_{i}$ are standard normal random variables. The data are the measurements for all cells, $\left\{y_{i}\right\}_{i=1}^{C}$. Sampling the prior predictive distribution is needed to infer the parameters $\boldsymbol{\theta}=(\mu,\sigma,\gamma,\alpha_{u},\alpha_{v},\beta_{u},\beta_{v})^{\prime}$ using ABC. We adopt independent uniform priors as chosen by Bonassi et al. [9], $\mu\sim\mathcal{U}(250,400)$, $\sigma\sim\mathcal{U}(0.05,0.5)$, $\gamma\sim\mathcal{U}(0.05,0.35)$, $\alpha_{u},\alpha_{v}\sim\mathcal{U}(0,50)$, and $\beta_{u},\beta_{v}\sim\mathcal{U}(0,7)$. We follow Vo et al. [56] in adopting 19 equally spaced quantiles of the empirical distribution of $\left\{y_{i}\right\}_{i=1}^{C}$ as the vector of summary statistics.

3.3 Implementation and optimisation using R

To simulate the toggle switch system (Equation (3)) and observation process (Equation (4)) we can use the Euler-Maruyama scheme [39],

[TABLE]

where $\Delta t$ is the discretisation step and $\xi_{i,u}(t)$ and $\xi_{i,v}(t)$ are independent normal random variables. For simplicity, we will assume $\Delta t=1$.

This can be implemented directly using R as shown in Listing 1. Here, the pair of SDEs for $u_{i}(t)$ and $v_{i}(t)$ are evolved one cell at a time, for $i=1,2,\ldots,C$. Such a naïve implementation is typical for initial prototypes, but it is profoundly inefficient with a single simulation with $T=600$ and $C=8000$ taking $35$ seconds. This is not practical for ABC sampling, even with individual simulations distributed across a HPC cluster.

One can achieve substantial computational improvements by re-writing the R code in terms of vector and matrix mathematics. This exploits optimised linear algebra libraries such as BLAS and LAPACK that often use SIMD. Listing 2 is an example of vectorised R code for the toggle switch model.

The key changes are: 1) genes of all cells are stored in vectors of length $C$; 2) all Gaussian random variates required for the entire simulation are generated in a single call of rnorm(); 3) all genes are simulated together using R’s vector maths functions; and 4) logical indexing is used instead of ifelse() to implement the boundary condition. The improvements are substantial, taking $1.3$ seconds to perform a single simulation with $T=600$ and $C=8000$. This is more than $25\times$ improvement, without any multithreading, explicit SIMD operations, or Rcpp code. See the excellent text by Gillespie and Lovelace [21] for details.

The purpose of comparing Listings 1 and 2 is to highlight the importance of optimising serial code. Now that the R implementation is efficient, we can distribute the process across multiple cores using the doParallel package as shown in Listing 3. Distributing $N=8064$ draws from the prior predictive across $16$ cores (Intel Xeon E5-2680v3 processor141414https://ark.intel.com/products/81908/Intel-Xeon-Processor-E5-2680-v3-30M-Cache-2-50-GHz-) takes just over $11$ minutes!

For some applications, optimised R can be fast enough and we do not advocate the use of further optimisation for all Bayesian applications. However, the trade-off between development time and runtime is very application specific. For this example, even using 16 cores, the optimised R code will take around 24 hours to draw $N=1,000,000$ prior predictive samples, which may be insufficient for ABC inference with small $\epsilon$. To motivate progressing beyond R, note that the indirect access of SIMD through optimised BLAS and LAPACK libraries has limitations in terms of cache utilisation for this application. This is highlighted in the next section.

3.4 Optimisation using C and SIMD operations

Direct access to SIMD through C with OpenMP enables superior memory access patterns that ensure substantially higher utilisation of the VPUs (see Section 5 for alternatives). While the optimised C implementations presented here are more complex than the optimised R implementation (Listing 2), they provide a practical example that will help reduce the learning curve for practitioners.

Listing 4 provides the main structure of the C program. Note the inclusion of the header files for the MKL libraries and the OpenMP library; we will point out calls to these libraries as necessary. Definitions of the constants VECL and ALIGN refer to the width of the SIMD vectors and the memory alignment to ensure vectors never cross cache lines (see Section 2). The functions simulate_tsw_SDE and main represent the C equivalents of the R code Listings 2 and 3. C implementations are given in Listings 5 and 6, the details of which are discussed below.

The optimised C implementation of the toggle switch model is given in Listing 5. From an algorithmic perspective this code can be considered a hybrid between the two R implementations (Listings 1 and 2). That is, cells are evolved in small blocks of length VECL, as opposed to individually (Listing 1) or all together (Listing 2). The outer loop iterates over the cell index c in strides of VECL (the E5-2680v3 CPU supports 265 bit vectors, giving VECL = 4 doubles). As in the optimised R code (Listing 2), all the random variates required within each block are pre-computed using the vdRngGaussian function from MKL. In each block there are VECL noisy observations, y[c],y[c+1],$\ldots$,y[c+VECL-1], to be computed (one per cell). Using SIMD we can evolve the SDE pair associated with the y[c] observation and obtain the others in the same block at the same time. This is done with a second loop over the index c2 that represents the position in the SIMD vector.

Note the use of the aligned statement within the simd directive. This allows the compiler to assume the arrays zeta (random variates) and y (simulated data) are aligned to the cache boundary and enables more efficient machine code to be generated. The SIMD loop is not really iterating over values of c2, but rather synchronously executing each statement within the loop for all values of c2 concurrently. This leads to exceptional reuse of L1 cache: since the SDE state variable u_t and v_t are updated in-place, all computation is performed using the fast CPU memory.

It is possible to do even better, since computing all Gaussian random variates for each block represents a trade-off. To understand this, note that advanced interfaces, such as compiler intrinsic functions (see discussion in Section 5), can enable efficient RNG sampling within the rest of the Euler-Maruyama update. However, this cannot be done with OpenMP as the function vdRngGaussian has overheads that dominate for very small sample sizes. Conversely, generating all the Gaussian random variates at the start with a single vdRngGaussian call forces the CPU to access L2 cache after a few blocks have been processed, resulting in loss of performance. This trade-off to memory access presented in Listing 5 cannot be efficiently replicated in R. Using BLAS and LAPACK functions for small vectors results in similar performance issues as for the vdRngGaussian function. Therefore, the best that can be done with R is to use long vectors to evolve all cells together.

Armed with a highly efficient SIMD implementation of the model, we can generate prior predictive samples across multiple CPU cores. This uses the multithreading pre-processor directive, #pragma omp parallel, as shown in Listing 6. Each thread is provided with its own random number stream using the Intel MKL/VSL function vslNewStream with independence guarenteed through the VSL_BRNG_MT2203 generator family (see Section 2). The memory allocated to store prior samples, theta, and simulation outputs, obs_vals, are shared across threads, using the shared clause, and thread specific offsets are computed, using the OpenMP functions omp_get_thread_num() and omp_get_num_threads(). The loop over index k is performing the prior predictive sampling tasks assigned to each thread. Unlike the R implementation (Listing 3), we explicitly partition the work amongst the available threads to allow reuse of the memory allocated for the Gaussian random variates, zeta. OpenMP does support various forms of automated scheduling of thread tasks, but for this example explicit partitioning is more straight forward.

3.5 Benchmarks

We compare the R implementations (Listings 1,  2 and 3) with the C implementation (Listings 4, 5 and 6) in terms of the time taken to generate a fixed number, $N$, of prior predictive samples given $P$ cores. The fastest sampler would, in turn, produce a smaller Monte Carlo error for a fixed computational budget. For example, the standard central limit theorem suggests a $4\times$ speedup will result in roughly $1/2$ the Monte Carlo error for the same computational budget. In an ABC setting there is also bias to consider, so reducing the model simulation time is absolutely crucial since the acceptance threshold $\epsilon$ must be sufficiently small.

The codes are benchmarked using two Intel Xeon E5-2680v3 (Haswell) processors or a single Intel Xeon Gold 6140 (Skylake) processor. The E5-2680v3 supports the AVX2 instruction set (256 bit vector) and the 6140 supports the AVX512 instruction set (512 bit vectors). Both processor architectures have similar serial performance with cores clocked at around 2.5 GHz. We benchmark with different numbers of cores, $P=1,2,4,8,16$ for each implementation. For each value of $P$ the benchmark is repeated four times, to obtain mean computation times $\hat{C}_{T}$, with each replicate generating $N=8064$ prior predictive samples. For fairness, we compiled R (version 3.3.1) from source using the Intel compiler suite and linked against MKL for optimised BLAS and LAPACK routines. The Intel C compiler version is 17.0.1 (compatible with the GNU C complier version 6.3.0). To demonstrate the difference between speed-up obtained from the improved memory utilisation and the usage of SIMD, we also compile a version of the C implementation with SIMD disabled using the -no-vec compiler option.

Figure 3 summarises the results in terms of runtimes $\hat{C}_{T}$ (Figure 3(a),(c),(e)) and speed-up factor $\hat{C}_{T}/\hat{C}_{B}$ for various serial baselines $\hat{C}_{B}$ (Figure 3(b),(d),(f)). See the supplementary material for precise data tables of runtimes. Using the naïve R code as the baseline (Figure 3(a)–(b)), a speedup of $16\times$ is possible with $P=16$ as one would expect. However, this pales by comparison with the more than $430\times$ speedup using the optimised R code and $P=16$, demonstrating the value in optimising serial code before implementing parallelism. The improvement is even greater for the C code with speedups of more than $1000\times$ (no vectors, $P=16$), $3800\times$ (AVX2 vectors, $P=16$), and $4300\times$ (AVX512, $P=16$). These impressive numbers are against a naïve baseline and more sensible numbers are obtained using the optimised R code as a baseline (Figure 3(c)–(d)). Here optimised C code can achieve $130\times$ speedup (AVX2, $P=16$) and $150\times$ speedup (AVX512, $P=16$) compared with only $16\times$ speedup with optimised R code and $P=16$. Finally, even with the C code as a baseline (Figure 3(e)–(f)), SIMD enables $30\times$ speedup (AVX2, $P=16$) and almost $40\times$ speedup (AVX512, $P=16$). In this case, the speedup from multithreading alone is only $13\times$ for $P=16$. There are a few reasons that could cause this, but it is likely that each thread does not perform enough work to mask threading overheads. $P\leq 4$, seems to provide a large enough workload per thread leading to $15\times$ speedup (AVX2, $P=4$) and $20\times$ speedup (AVX512, $P=4$).

3.6 Summary

This tutorial has demonstrated the computational benefits of direct access to SIMD operations in comparison to indirect access through pre-compiled BLAS and LAPACK libraries available in languages such as R and Matlab. Our C implementation (Listings 5 and 6) is more that $20\times$ faster than our optimised R implementation (Listing 2 and 3) for the same CPU hardware and number of available cores, $P$. Restricting comparison to the C implementation, SIMD operations can improve the performance of prior predictive sampling for ABC by a factor of more than $6\times$ compared with standard scalar implementations. Given the enormous number of prior predictive samples required to obtain meaningful parameter estimates with ABC [3], we submit that these techniques are highly relevant for practitioners dealing with such applications. All codes presented here are available as supplementary material.

4 Case studies

In this section, we explore the benefit of SIMD operations in combination with mutlithreading using two case studies relevant to Bayesian practitioners.

4.1 Case study 1: weakly informative priors

Here, we consider the selection of weakly informative priors [19] within the framework proposed by Evans and Jang [16] using the implementation by Nott et al. [42]. Following Nott et al. [42], the prior predictive $p$-value,

[TABLE]

is a measure of Bayesian model criticism, where $\mathcal{D}_{\text{obs}}$ are the observational data, and $\pi(\mathcal{D})$ is the prior predictive distribution (Equation (2)). Equation (6) provides a $p$-value for prior-data conflict [16] associated with low prior density assigned to parameters with good model fit. If $\gamma$ is a pre-defined small cut-off, then $p(\mathcal{D}_{\text{obs}})\leq\gamma$ signifies prior-data conflict.

Consider a family of priors, $\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}})$ with the hyperparameter ${\boldsymbol{\lambda}}\in{\boldsymbol{\Lambda}}$. A base prior, $\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}}_{0})$ for ${\boldsymbol{\lambda}}_{0}\in{\boldsymbol{\Lambda}}$, represents the current best knowledge of the parameter $\boldsymbol{\theta}$. Let

[TABLE]

be the prior predictive distribution for the prior $\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}})$ and likelihood $\mathcal{L}(\boldsymbol{\theta};\mathcal{D})$, and let $p_{\boldsymbol{\lambda}}(\mathcal{D}_{\text{obs}})$ be the prior predictive $p$-value under $\pi(\mathcal{D}\mid{\boldsymbol{\lambda}})$ (equivalent to Equation (6)),

[TABLE]

Assume data are generated under the base prior predictive distribution, that is, $\mathcal{D}_{0}~{}\sim~{}\pi(\mathcal{D}\mid{\boldsymbol{\lambda}}_{0})$. The task is to find values of ${\boldsymbol{\lambda}}$ such that,

[TABLE]

Priors satisfying Equation (9) are weakly informative with respect to the base prior. Weak informativity indicates less prior-data conflict than the base. Equivalently, weak informativity at level $\alpha$ indicates $\gamma_{{\boldsymbol{\lambda}},\alpha}>\gamma_{{\boldsymbol{\lambda}}_{0},\alpha}$ where $\gamma_{{\boldsymbol{\lambda}},\alpha}$ and $\gamma_{{\boldsymbol{\lambda}}_{0},\alpha}$ are the prior-data conflict cut-offs such that $\mathbb{P}\left(p_{\boldsymbol{\lambda}}(\mathcal{D}_{0})\leq\gamma_{{\boldsymbol{\lambda}},\alpha}\right)=\mathbb{P}\left(p_{{\boldsymbol{\lambda}}_{0}}(\mathcal{D}_{0})\leq\gamma_{{\boldsymbol{\lambda}}_{0},\alpha}\right)=\alpha$.

Given a set of $K$ hyperparameter values, ${\boldsymbol{\lambda}}_{1},\ldots,{\boldsymbol{\lambda}}_{K}\in{\boldsymbol{\Lambda}}$, with ${\boldsymbol{\lambda}}_{0}$ corresponding to the base prior, the task is to compute $\gamma_{{\boldsymbol{\lambda}}_{k},\alpha}$ such that $\mathbb{P}\left(p_{{\boldsymbol{\lambda}}_{k}}(\mathcal{D}_{0})\leq\gamma_{{\boldsymbol{\lambda}}_{k},\alpha}\right)=\alpha$ for $k=0,1,\ldots,K$ and $\mathcal{D}_{0}\sim\pi(\mathcal{D}\mid{\boldsymbol{\lambda}}_{0})$. The process proceeds as in Algorithm 2. The output of Algorithm 2 is a set of $\alpha$ level cutoffs $\{\gamma_{{\boldsymbol{\lambda}}_{0},\alpha},\gamma_{{\boldsymbol{\lambda}}_{1},\alpha},\ldots,\gamma_{{\boldsymbol{\lambda}}_{K},\alpha}\}$, to derive a set of weakly informative prior distributions with respect to the base prior $\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}}_{0})$, that is, $\{\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}}_{k}):\gamma_{{\boldsymbol{\lambda}}_{k},\alpha}>\gamma_{{\boldsymbol{\lambda}}_{0},\alpha},k=1,2,\ldots,K\}$.

This process is extremely computationally intensive since the evidence term, $\pi(\mathcal{D}\mid{\boldsymbol{\lambda}})$, must be evaluated for many different $\mathcal{D}$. Adaptive sequential Monte Carlo (SMC) sampling (Algorithm 3) using likelihood annealing and Markov chain Monte Carlo (MCMC) proposals [5, 11] is used to estimate $\pi(\mathcal{D}\mid{\boldsymbol{\lambda}})$. Algorithm 2 requires $2KN$ executions of the adaptive SMC (Algorithm 3) with $N_{p}$ particles with tuning parameters $c$, related to the number of MCMC iterations, and $h$, related to the scaling of random walk MCMC proposals.

We consider weakly informative priors for the analysis of an acute toxicity test in which $M$ groups of animals are given different dosages of some toxin, and the number of deaths in each group are recorded [16, 42]. A logistic regression model is applied for the number of deaths $y_{i}$ in group $i\in[1,2,\ldots M]$,

[TABLE]

where $n_{i}$ and $x_{i}$ are respectively the number of animals and the toxin dose level in the $i$th group of the experiment, and $\boldsymbol{\theta}=(\beta_{0},\beta_{1})^{\prime}$ are the regression parameters. Given data $\mathcal{D}_{\text{obs}}=[y_{\text{obs},1},y_{\text{obs},2},\ldots,y_{\text{obs},M}]$, the likelihood function is

[TABLE]

From the data in [16] we have $M=4$, $n_{1}=n_{2}=n_{3}=n_{4}=5$, $x_{1}=-0.86$, $x_{2}=-0.3$, $x_{3}=-0.05$ and $x_{4}=-0.75$. We adopt bivariate Gaussian priors from [42], $\boldsymbol{\theta}~{}\sim~{}\mathcal{N}(\mathbf{0},\text{diag}({\boldsymbol{\lambda}})^{2})$ and ${\boldsymbol{\lambda}}_{0}=(10,2.5)^{\prime}$.

4.1.1 Parallelisation and SIMD opportunities

Parallel implementations of SMC samplers using multithreading require thread synchronisation for both the resampling and annealing steps [29, 33, 41]. Thus, it is more beneficial to distribute the $K$ hyperparameter values in Algorithm 2 across $P$ cores, with each thread computing the $\alpha$ quantile for $K/P$ hyperparameters, since these may be performed independently. For every hyperparameter, we sequentially process the $N$ $p$-value computations. The SMC sampler (Algorithm 3) can be accelerated through SIMD. SMC samplers are well suited to SIMD since sychronisation is maintained automatically. While there are many aspects of Algorithm 3 that use SIMD in the code example, we focus here on SIMD for the MCMC proposal kernel for diversification of particles (Steps 11–17 of Algorithm 3) as it has the greatest effect on performance.

The strategy is similar to the C implementation of the toggle switch model (Section 3). At SMC step $n$, each particle is perturbed via $R_{n}$ MCMC steps. The same operations are performed at each MCMC iteration, with the exception of a single branch operation that arises from the accept/reject step. We process the particle updates in blocks of length $V$ and evolve the $R_{n}$ MCMC steps for this block together using SIMD. All of the Gaussian and uniform random variates required for the block are generated together before the block is processed. The SIMD MCMC proposals are performed as in Algorithm 4.

Hadamard notation indicates element-wise division, $x\oslash y$, multiplication, $x\circ y$ and exponentiation, $x^{\circ a}$, for scalar $a$. Element-wise application of a function, $f$, over length $V$ vectors is denoted by $f^{V}(x_{i:i+V})=[f(x_{i}),f(x_{i+1}),\ldots,f(x_{i+V})]$. In practice, Algorithm 4 works with the likelihood and prior on the log scale to avoid numerical underflow. For optimal performance, SIMD forms for the likelihood $\mathcal{L}^{V}$, prior density $\pi^{V}$ and Gaussian proposal density $\phi^{V}$ are required. This is available for the logistic regression model, Gaussian priors and proposals used here. Other aspects of Algorithm 3 that can utilise SIMD include likelihood evaluations and the computation of ESS and weight updates.

The scalar bottleneck is the multinomial resampling step that utilises the look-up method. Parallel approximations for resampling have been proposed [41] and demonstrated to be very effective for large scale SMC samplers. Since we focus here on a SIMD implementation of SMC, we do not implement this, however, extending [41] to SIMD is an important piece of future work.

4.1.2 Performance

We test the performance improvement obtained through vectorisation and multithreading using the Xeon E5-2680v3 and Xeon Gold 6140 processors. The evaluation of the weak informativity test in Algorithm 2 is performed for $K=400$ hyperparameters and $N=400$ datasets. In all simulations, the $K$ hyperparameter values are generated using a bivariate uniform distribution ${\boldsymbol{\lambda}}_{k}\sim\mathcal{U}([0.1,10]\times[0.1,20])$ for $k=1,2,\ldots,K$. The SMC sampling is performed with $N_{p}=500$ particles, with the lower particle count enabling computation to remain largely within L1 and L2 cache. Tuning parameter values are specified as $c=0.01$ and $h=2.38/\sqrt{2}$. Results are provided in Figure 4 (see supplementary material for data tables).

Figure 4(b) shows almost perfect speedup from multithreading and consistent improvement due to SIMD when comparing scalar and SIMD performance for the same CPU and core count ($2.1\times$ for AVX2 and $2.3\times$ for AVX512). Diminishing returns are observed when stepping from $256$ bit vector operations to $512$ bit SIMD. In our bioassay example, the MCMC proposal kernel performs a small number of steps (usually no more than $R_{n}=30$) and the number of particles in the SMC sampler is $N_{p}=500$. The performance will improve for SIMD in cases where longer MCMC runs are required, since the MCMC step will dominate SMC iterations and more reuse of L1 cache will occur.

4.1.3 Summary

We have presented the more challenging problem of weak informativity tests over a family of priors $\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}})$ with respect to a base prior $\pi(\boldsymbol{\theta}\mid{\boldsymbol{\lambda}}_{0})$. This requires a large number of approximations of the posterior normalising constant for different datasets from the prior predictive distribution. The traditional approach of parallelisation of each SMC iteration would reduce the level of parallelism available across hyperparameters. Using the SIMD implementation of SMC we double our computational performance and can reserve multithreading across hyperparameters. A combination of parallelism and SIMD could be implemented to improve the performance of a single SMC step. However, this offers little benefit here since utilising threads for SMC forces hyperparameters to be processed serially. If HPC resources are available then hyperparameters could be distributed across servers: leaving both threading and SIMD for SMC steps.

4.2 Case study 2: Parameter inference for a non-Gaussian asymmetric volatility model

We now consider an adaptive SMC sampler to perform parameter inference in eleven dimensions for the “bad evironment – good environment” (BEGE) model of innovations on stock market returns [4, 50].

4.2.1 The BEGE model

The BEGE model [4] is a non-Gaussian generalisation of the Glosten-Jagannathan-Runkle (GJR) asymmetric volatility model [25] and is a generalised autoregressive conditional heteroskedasticity (GARCH) model [8]. The BEGE model describes the time-series of stock market returns, $\{r_{t}\}_{t\geq 0}$ using a model on innovation on returns, $\{u_{t}\}_{t\geq 0}$, that consists of a linear combination of positive “good environment” shocks, $\{\omega_{p,t}\}_{t\geq 0}$, and negative “bad environment” shocks, $\{\omega_{n,t}\}_{t\geq 0}$. The BEGE time-series evolves according to

[TABLE]

where $\mu$ is the conditional mean of returns, $\tilde{\Gamma}(k,1)$ is the centered (de-meaned) gamma distribution with shape $k$ and unit scale, $\{p_{t}\}_{t\geq 0}$ and $\{n_{t}\}_{t\geq 0}$ are respectively the shapes of the positive and negative shocks, and $\sigma_{p}$ and $\sigma_{n}$ are their constant scales. The shape parameters evolve according to

[TABLE]

where $p_{0}$, $n_{0}$ are initial conditions, and $\rho_{p}$, $\rho_{n}$, $\phi^{+}_{p}$, $\phi^{+}_{n}$, $\phi^{-}_{p}$, and $\phi^{-}_{n}$ are autoregression parameters.

We use S&P Composite Index returns over the period July 1926 to January 2018 (obtained from the Center of Research in Security Prices) consisting of $T=1099$ months of logged monthly divided-adjusted returns, $\mathcal{D}=R_{\text{obs},T}=\left\{r_{\text{obs},t}\right\}_{0\leq t\leq T}$. Using these data, Bayesian inference on the unknown parameters $\boldsymbol{\theta}=\left(p_{0},\sigma_{p},\rho_{p},\phi_{p}^{+},\phi_{p}^{-},n_{0},\sigma_{n},\rho_{n},\phi_{n}^{+},\phi_{n}^{-},\mu\right)$ is performed. Priors are adopted from [50] and are given by $p_{0}\sim\mathcal{U}(10^{-4},0.5)$, $\sigma_{p}\sim\mathcal{U}(10^{-4},0.3)$, $\rho_{p}\sim\mathcal{U}(10^{-4},0.99)$, $\phi_{p}^{+}\sim\mathcal{U}(10^{-4},0.5)$, $\phi_{p}^{-}\sim\mathcal{U}(10^{-4},0.5)$, $n_{0}\sim\mathcal{U}(10^{-4},1)$, $\sigma_{n}\sim\mathcal{U}(10^{-4},0.3)$, $\rho_{n}\sim\mathcal{U}(10^{-4},0.99)$, $\phi_{n}^{+}\sim\mathcal{U}(-0.2,0.1)$, $\phi_{n}^{-}\sim\mathcal{U}(10^{-4},0.75)$, and $\mu_{t}\sim\mathcal{U}(-0.9,0.9)$.

The major challenge in this inference problem is the computational cost associated with the evaluation of the log-likelihood function,

[TABLE]

Bekaert et al. [4] show that evaluation of the transitional densities, $\pi(r_{\text{obs},t}\mid r_{\text{obs},t-1},\boldsymbol{\theta})$, requires the so called BEGE density $\pi_{\text{BEGE}}(u_{t}\mid\sigma_{p},\sigma_{n},p_{t-1},n_{t-1})$. This can be approximated by computing the BEGE distribution function and taking finite differences [4]. The BEGE distribution function is given by

[TABLE]

with $G_{\tilde{\Gamma}}(\omega_{p,t}-u_{t}\mid n_{t-1},\sigma_{n})=1-F_{\tilde{\Gamma}}(\omega_{p,t}-u_{t}\mid n_{t-1},\sigma_{n})$, where $\pi_{\tilde{\Gamma}}(\cdot\mid k,s)$ and $F_{\tilde{\Gamma}}(\cdot\mid k,s)$ are, respectively, the probability density and distribution functions of a centered gamma distribution with shape $k$ and scale $s$. The integral in Equation (15) can be approximated numerically using quadrature. For each point in the discretisation of the $\omega_{p,t}$ parameter space, we require two evaluations of the incomplete gamma function,

[TABLE]

where $\Gamma(a)$ is the gamma function. Equation (16) can be computed using a power series [44].

We apply an adaptive SMC sampler to move $N_{p}$ particles from the prior $\pi(\boldsymbol{\theta})$ to the posterior $\pi(\boldsymbol{\theta}\mid R_{\text{obs},T})$ under the BEGE model. This is a very similar SMC sampler to that applied in Section 4.1 (Algorithm 3). The main difference is the scaling rule for the proposal kernel within the MCMC step, as the posterior is highly non-Gaussian. We apply the method of  [47] to evaluate a set of MCMC trials each with a random scale factor, $h\in[0.1,0.2,\ldots,1.0]$, at each SMC iteration. We then choose the scale factor, $h_{\text{opt}}$, which maximises the median expected squared jump distance across all particles. This continues until at least half of the particles have moved further than the median [47].

4.2.2 Parallelisation and vectorisation opportunities

In the BEGE model inference problem, we utilise both SIMD and multithreading to accelerate a single SMC sampler. Parallel implementations of SMC samplers and particle filters have been well studied in the literature [29, 33, 41]. Within a single SMC iteration, particles are completely independent of each other. However, as noted in Section 4.1, automatic synchronisation can occur within the MCMC proposal mechanism. Rather than distribute $N_{p}$ particles across $P$ cores, we ensure that the distribute occurs in contiguous blocks of length $V$. Each core will process $N_{p}/(PV)$ blocks of particles. All data associated with particles is processed in contiguous blocks of length $V$, allowing each thread to independently exploit SIMD within the MCMC proposal mechanism as per Algorithm 4.

Another way to exploit SIMD is in the evaluation of the BEGE log-likelihood. We extend the approximation of the integral in Equation (15). Consider a discretisation of $\omega_{p,t}$ of $N_{\omega}+1$ nodes with spacing $\Delta\omega$. Then Equation (15) can be approximated by

[TABLE]

where $\omega_{p,t}^{j}=\omega_{p,t}^{0}+j\Delta\omega$, for $j=0,1,\ldots,N_{\omega}$ and $\omega_{p,t}^{0}$ is the lower bound of the discretisation. For every point $\{\omega_{p,t}^{j}\}_{0\leq j\leq N_{\omega}}$, the incomplete gamma function is evaluated twice; once with $p_{t-1}$ and once with $n_{t-1}$. Therefore, we can consider SIMD for the incomplete gamma function for blocks of $\omega_{p,t}$ points of length $V$, that is, $P^{V}(a,x_{1:V})=\left[P(a,x_{1}),P(a,x_{2}),\ldots,P(a,x_{V})\right]$. We achieve this by extending the method of [44], resulting in the element-wise vector series expression

[TABLE]

where $\exp^{V}\left(-x_{1:V}\right)=\left[e^{-x_{1}},e^{-x_{2}},\ldots,e^{-x_{V}}\right]\in\mathbb{R}^{1\times V}$, and $\mathbf{e}_{1:V}~{}=~{}[1,1,1,\ldots,1]~{}\in~{}\mathbb{R}^{1\times V}$. Equation 18 allows efficient iteration that reuses previous steps, thus enabling good use of L1 and L2 cache. The series is truncated once it has converged under the $\infty$-norm.

We proceed to approximate $F_{\text{BEGE}}(u_{t}\mid\sigma_{p},\sigma_{n},p_{t-1},n_{t-1})$ by: 1) applying Equation (18) across the discretisation, $\{\omega_{p,t}^{j}\}_{0\leq j\leq N_{\omega}}$, in blocks of length $V$ with shape $p_{t-1}$; 2) applying Equation (18) across the discretisation, $\{\omega_{p,t}^{j}\}_{0\leq j\leq N_{\omega}}$, in blocks of length $V$ with shape $n_{t-1}$; and 3) accumulating the sum of products in Equation 17.

4.2.3 Performance

We implement the BEGE inference problem using adaptive SMC with likelihood annealing with $N_{p}=1024$ particles. We approximate $F_{\text{BEGE}}(u_{t}\mid\sigma_{p},\sigma_{n},p_{t-1},n_{t-1})$ using Equations (17) and (18) with the discretisation $\{\omega_{p,t}^{j}\}_{0\leq j\leq N_{\omega}}$, $N_{\omega}=100$, $\omega_{p,t}^{0}=10^{-4}-p_{t-1}\sigma_{p}$, and $\Delta\omega=(10\sigma_{p}\sqrt{p_{t-1}}-\omega_{p,t}^{0})/(N_{\omega}-1)$, as is performed by [4] and [50]. Results are provided in Figure 5 .

We observe an improvement of up to $2\times$ for AVX2 and almost $4\times$ for AVX512 regardless of the number of cores. This is improved VPU utilisation for AVX512 compared with the weak informativity test (Section 4.1). That is, the BEGE distribution function approximation is a sufficiently large proportion of the total computation cost, that efficient calculation can exploit the longer 512 bit vectors. However, the overall speedup factor is lower from multithreading and increases slowly with $P$. The SMC sampler requires threads to synchronise at the end of each iteration. This synchronisation is in the resampling and annealing steps, and in estimating the optimal MCMC proposal scaling could be causing bottlenecks. This can be seen in the diminishing returns on the parallel speed-up as the number cores increases.

4.2.4 Summary

We have demonstrated how SIMD can be used to further accelerate a parallel SMC sampler for a challenging inference problem from econometrics. Note that the log-likelihood approximation we apply, based on the work of Bekaert et al. [4] is biased. Recently, Li et al. [35] proposed an unbiased likelihood estimator, for which there are SIMD opportunities also. For the purposes of this manuscript, we find the biased approximation of [4] lends itself more direct discourse.

5 Discussion

Across each of our example applications there have been some common features, which allow several conclusions to be drawn regarding task suitability for vectorisation. Firstly, it is clear that the form of stochastic model under study can have a dramatic effect on the potential performance boost due to vectorisation. Simulation schemes such as Euler-Maruyama are ideal candidates for vectorisation, however, exact stochastic simulation methods like the Gillespie direct method are more challenging. We show this effect in the supplementary material using the application of ABC to models of Tuberculosis transmission [51] and observe only moderate improvements from SIMD. This highlights a difference between continuous time and discrete time Markov processes, and may motivate the practitioner to consider alternative algorithms that are more suited to SIMD, for example, the $\tau$-leap approximations to the Gillespie algorithm. Similarily, the success of the SIMD version of the MCMC proposal kernel in Section 4.1 relied on a SIMD implementation of the likelihood function.

Secondly, each application involved nested parallelisation. The utility of SIMD here is that within each parallel thread, the computational tasks may be further sub-divided through use of VPUs. Efficiency gains due to SIMD are multiplicative to those arising from the use of thread parallelisation. Finally, each application involved a mixture of tasks that can be performed completely independently and tasks that require synchronisation or communication. This is important, since independent parallel computation is ideal for multithreading and fine grain parallel operations with frequent synchronisation are well suited for SIMD. Our demonstrations are widely applicable, since these three features are common to many Bayesian applications.

Although we have focused on the acceleration of Monte Carlo methods for Bayesian applications, likelihood-based and frequentist applications can also benefit. This is especially true for complex optimisations for maximum likelihood parameter estimates [29]. Furthermore, other Monte Carlo schemes can benefit from our approach. Any application that has been demonstrated using GPGPUs could exploit SIMD on the CPU [33, 28]. This is rarely taken into consideration when comparison between GPGPUs and CPUs are made [28, 29, 34].

We have provided example C programs, using OpenMP and MKL libraries. We appreciate that, for very good reasons, higher level languages such as Matlab and R are often the preferred environments for many practitioners. The implementation techniques we present can be exploited through Matlab C-MEX or Rcpp interfaces. However, Matlab and R must be configured correctly to use the required compiler options. Of course, matrix operations performed within high level languages, such as Matlab or R, are likely already utilising SIMD via high-performance BLAS and LAPACK libraries.

Unfortunately, many of the SIMD and memory optimisations we presented cannot be directly exploited using a high level language alone. This is demonstrated practically in Section 3. The only reliable way to access SIMD level parallelism from within a high level language is to use built-in matrix or vector functions that are already optimised (for example, by compiling R to use MKL). Memory access is still a significant challenge here. Between successive high level functions calls, it is unlikely that the caches are preserved, and as a result, our SIMD versions of the Euler-Maruyama scheme and MCMC proposal cannot be directly replicated in Matlab or R alone. A notable exception to this rule is the Julia language [6] using the @simd macro to achieve high performance.

OpenMP is also not the only way to access SIMD within C/C++. For example, while OpenCL is primarily applied to GPCPUs and other accelerators, OpenCL kernels may be compiled for CPUs that support SIMD units [29, 36]. Instruction level intrinsic functions (also available in R via the RcppXsimd package) allow advanced features such as efficient random variates for small vectors, but this approach is very challenging and akin to machine code.

Many other Monte Carlo and Bayesian applications could benefit from these approaches. One possibility is large scale particle filters, perhaps operating within a pseudo-marginal scheme [2], which could be implemented as a hybrid algorithm in which proposals and weight updates are performed in SIMD blocks that are distributed across multiple threads. Such a scheme would enable prefix summations, that is, parallel computation of cumulative sums, to be performed [27] in the resampling step, which may improve the performance. Another possibility is within the class of ABC validation or post-processing procedures. For example, in the recalibration post-processing technique of [45], each of the $N$ samples $\boldsymbol{\theta}^{*}$ from the approximate posterior distribution are individually recalibrated by the construction of a further approximate (ABC) posterior for each $\boldsymbol{\theta}^{*}$ using all previously generated $(\boldsymbol{\theta}^{*},\mathcal{D}^{*})$ pairs, but based on observing the associated (simulated) dataset ${\mathcal{D}}^{*}$, and constructing all univariate marginal posterior distribution functions. This re-use of ABC algorithm and previously generated parameter values and datasets is very common in ABC [[, e.g.]]Blum2013, and makes them particularly suited for performance gains through parallelism and SIMD.

We have demonstrated that by following a few simple guidelines to maximise the utilisation of modern CPUs, advanced Monte Carlo methods may be relatively straightforwardly accelerated by a factor of up to $6\times$ in addition to further speedups obtained through multithreading. These techniques will only become more relevant in the future as CPUs architectures are released with wider VPUs and statisticians develop more complex and sophisticated inferential algorithms.

Acknowledgements

C.D. was supported by the Australian Research Council (ARC) under the Discovery Project scheme (DP200102101). S.A.S. was supported by the ARC under the Discovery Project scheme (DP160102544). C.D., S.A.S., and D.J.W. are supported by the Australian Centre of Excellence for Mathematical and Statistical Frontiers (ACEMS; CE140100049). C.D. and D.J.W. also acknowledge support frm the Centre for Data Science, Queensland University of Technology. Computational resources were provided by the eResearch Office, Queensland University of Technology.

Software availability

Source code for the tutorial and case studies are freely available on GitHub https://github.com/davidwarne/Bayesian_SIMD_examples

Appendix A Benchmark Tables

Data tables of mean runtimes used for Figures 3–5 in the main text are provided in Tables A.1,A.2 and A.3

Appendix B Tuberculosis transmission example

The second model describes tuberculosis transmission dynamics [51, 49, 18]. The evolution of the number of tuberculosis cases over time, $N(t)$, is given by

[TABLE]

where $G(t)$ is the number of unique genotypes of the tuberculosis bacterium at time $t$ and $X_{i}(t)$ is the number of cases caused by the $i$th genotype. The time-varying evolution of the $X_{i}(t)$’s is driven by a linear birth-death process with birth rate $\alpha$ and death rate $\delta$. Mutation events occur at rate $\tau$ resulting in new genotypes and increases $G(t)$.

Given values for $\alpha$, $\delta$ and $\tau$, exact realisations can be generated using the Gillespie Direct Method [22] with propensity functions,

[TABLE]

where $a_{b}$, $a_{d}$ and $a_{m}$ are the birth, death and mutation propensities respectively.

[51] demonstrate the application of ABC methods to infer the parameters $\theta=(\alpha,\delta,\tau)^{\prime}$ for a real dataset of tuberculosis cases in San Francisco during the early 1990s that consists of $326$ distinct genotypes. The number of distinct genotypes and the genetic diversity are used as summary statistics (see [51]).

B.1 Parallelisation and vectorisation opportunities

Given $P$ cores, the process of i.i.d. sampling $(\theta^{*},{\mathcal{D}}^{*})$ from the joint prior predictive can be trivially divided between cores. For the toggle switch model (Equation (3)) it is reasonable to divide the workload evenly since the number of floating point operations per realisation is constant. This is called a static schedule for OpenMP. Each core will generate $N_{P}=N/P$ prior predictive samples. We assume $N_{P}$ is an integer.

For the tuberculosis model, generation of the prior predictive samples involves a discrete-state, continuous-time Markov process that is implemented with the Gillespie algorithm. The simulation times can vary wildly and a static schedule can result in a load imbalances. Therefore, guided or dynamic scheduling could perform better. In practice, we find the overhead of load imbalance is relatively low for large $N$. Furthermore, an early termination rule, as used in the Lazy ABC approach of [43] would further reduce imbalances. We conclude that static distribution of prior predictive samples is appropriate as a general strategy for ABC applications. Reproducibility of RNG sequences is also possible with static scheduling.

No such general strategy exists for vectorisation. The simulation algorithms for the two models are very different. The Euler-Maruyama scheme is applied to toggle switch model and the Gillespie Direct method is applied to the tuberculosis model. The performance improvement opportunity naturally differs.

Since each Euler-Maruyama simulation performs identical operations for each cell $i=~{}1,2\ldots,C$, this is ideal for vectorisation. We evolve cells in blocks of length $V$ where $64V$ is the bit width of the VPU inputs. Aside from the generation of Gaussian random variates, every operation is replaced with SIMD operation. We generate $2n_{t}V$ Gaussian variates using the Intel MKL RNG ahead of processing each block. This is efficient when $n_{t}V$ is within L1 cache capacity, but large enough to take advantage of the Intel routines. Typically, the Intel routines require at least $1,000$ variates to be generated for peak performance.

Fewer opportunities exist for vectorisation of the Gillespie method used for the tuberculosis model. The only suitable vectorisable step for updating genotype weights for the lookup method to select the genotype to apply the next event. The discrepancy measure components, $g$ and $H$, can also be vectorised. The length of the state vector in the tuberculosis model is large enough for some performance boost from the vectorisation to be noticeable. The genetic diversity $H$ calculation is particularly efficient as the sum of squares operation can exploit fused-multiply-add (FMA) SIMD operations that performs $A\leftarrow A+(B\circ C)$.

B.2 Performance

We report performance improvements obtained through vectorisation and multithreading for two CPU architectures, the Xeon E5-2680v3 (Haswell)151515https://ark.intel.com/products/81908/Intel-Xeon-Processor-E5-2680-v3-30M-Cache-2-50-GHz- and Xeon Gold 6140 (Skylake)161616https://ark.intel.com/products/120485/Intel-Xeon-Gold-6140-Processor-24-75M-Cache-2-30-GHz-. The Xeon E5-2680v3 supports Intel’s AVX2 instruction set ($256$ bit vectors), and the Xeon Gold 6140 supports Intel’s AVX512 instruction set ($512$ bit vectors).

Tables B.1 and B.2 show the results for the tuberculosis model for a range of prior predictive sample counts, $N$. The speedup factor from vectorisation is unsurprisingly not as significant as for the toggle switch model. This demonstrates the limitations of SIMD operations (including GPGPU-based accelerators). There is still an improvement from vectorisation of between $1.4\times$ to $1.8\times$. A good result in the context of the HPC literature [15, 32, 40].

The speedup factor due to vectorisation is larger, at $1.8\times$, for the Xeon E5-2680v3 compared with $1.4\times$ for the Xeon Gold 6140. However, the Xeon Gold 6140 results demonstrate almost $3\times$ improvement in the scalar performance over the Xeon E5-2680v3. Therefore, we conclude that the tuberculosis model simulation is accessing L3 cache and main memory more often, and thus taking advant memory bandwidth and extra memory channels on the Xeon Gold 6140. If the tuberculosis model is not utilising L1 and L2 cache as intensively, then this explains the reduced speedup for wider vectors as the VPUs is not efficiently utilised.

Differences in performance of the toggle switch model and the tuberculosis model suggest the choice of stochastic simulation algorithm is crucial. Introducing approximations [24, 23] for the tuberculosis model could significantly improve performance.