Medical Image Retrieval Based On the Parallelization of the Cluster Sampling Algorithm

TL;DR

This paper introduces parallel cluster sampling algorithms, including a multi-chain version, which are highly parallelizable and effectively applied to medical image retrieval tasks.

Contribution

The paper presents novel parallel cluster sampling algorithms and demonstrates their efficiency in medical image retrieval applications.

Findings

Multi-chain version is embarrassingly parallel

Efficient for medical image retrieval

Applicable to other domains

Abstract

In this paper we develop parallel cluster sampling algorithms and show that a multi-chain version is embarrassingly parallel and can be used efficiently for medical image retrieval among other applications.

Equations

$${\mathbf{x}}^{\text{a}}=\underset{\mathbf{x}}{\min }T(\mathbf{x})={\Vert \mathcal{H}(\mathbf{x})-\mathbf{y}\Vert }_{R^{-1}}^2+\alpha \Vert \mathbf{x}{\Vert }_{\mathbf{C}}^2,$$

$${\Vert \mathbf{c}-\mathbf{d}\Vert }_M^2=(\mathbf{c}-\mathbf{d}{)}^{\mathbf{T}}\mathbf{M}(\mathbf{c}-\mathbf{d})$$

$${\nabla }_{\mathbf{x}}T(\mathcal{H}\mathbf{x})={\left[\partial \mathcal{H}(\mathbf{x})\right]}^{\ast }{\mathbf{R}}^{-1}(\mathcal{H}(\mathbf{x})-\mathbf{y})+\alpha \mathbf{Cx},$$

$$\mathcal{P}(\mathbf{x}|\mathbf{y})\propto \mathcal{P}(\mathbf{y}|\mathbf{x})\mathcal{P}(\mathbf{x}).$$

$$\mathcal{P}(\mathbf{y}|\mathbf{x})\propto \exp \left(-\frac{1}{2}\Vert \mathcal{H}(\mathbf{x})-\mathbf{y}{\Vert }_{{\mathbf{R}}^{-1}}^2\right).$$

$$\mathcal{P}(\mathbf{x})\propto \exp \left(-\frac{1}{2}\Vert \mathbf{x}{\Vert }_{\mathbf{B}}^2\right),$$

$$-\log \mathcal{P}(\mathbf{x}|\mathbf{y})\propto \frac{1}{2}\Vert \mathcal{H}(\mathbf{x})-\mathbf{y}{\Vert }_{{\mathbf{R}}^{-1}}^2+\frac{1}{2}\Vert \mathbf{x}{\Vert }_{\mathbf{B}}^2.$$

$${\mathcal{P}}^{\mathrm{a}}(\mathbf{x})=\mathcal{P}(\mathbf{x}|\mathbf{y})=\frac{\mathcal{P}(\mathbf{y}|\mathbf{x}){\mathcal{P}}^{\mathrm{b}}(\mathbf{x})}{\mathcal{P}(\mathbf{y})}\propto \mathcal{P}(\mathbf{y}|\mathbf{x}){\mathcal{P}}^{\mathrm{b}}(\mathbf{x}),$$

$${\mathcal{P}}^{\mathrm{b}}(\mathbf{x})=\sum _{i=1}^{\text{textsc}{n}_{\text{c}}}{\tau }_i\mathcal{N}({\mu }_i,{\text{Upsigma}}_i)=\sum _{i=1}^{\text{textsc}{n}_{\text{c}}}{\tau }_i\frac{(2\pi {)}^{-\frac{\text{textsc}{n}_{\text{var}}}{2}}}{\sqrt{|{\text{Upsigma}}_i|}}\exp \left(-\frac{1}{2}\Vert \mathbf{x}-{\mu }_i{\Vert }_{{{\text{Upsigma}}_{\mathbf{i}}}^{-1}}^2\right),$$

$$\mathcal{P}(\mathbf{y}|\mathbf{x})=\frac{(2\pi {)}^{-\frac{\text{m}}{2}}}{\sqrt{|\mathbf{R}|}}\exp \left(-\frac{1}{2}\Vert \mathcal{H}(\mathbf{x})-\mathbf{y}{\Vert }_{{\mathbf{R}}^{-1}}^2\right),$$

$${\mathcal{P}}^{\mathrm{a}}(\mathbf{x})$$

$$\sum _{i=1}^{\text{textsc}{n}_{\text{c}}}{\tau }_i\frac{(2\pi {)}^{-\frac{\text{textsc}{n}_{\text{var}}}{2}}}{\sqrt{|{\text{Upsigma}}_i|}}\exp \left(-\frac{1}{2}\Vert \mathbf{x}-{\mu }_i{\Vert }_{{{\text{Upsigma}}_{\mathbf{i}}}^{-1}}^2\right)$$

$$\propto \exp \left(-\frac{1}{2}\Vert \mathcal{H}(\mathbf{x})-\mathbf{y}{\Vert }_{{\mathbf{R}}^{-1}}^2\right)\sum _{i=1}^{\text{textsc}{n}_{\text{c}}}\frac{{\tau }_i}{\sqrt{|{\text{Upsigma}}_i|}}\exp \left(-\frac{1}{2}\Vert \mathbf{x}-{\mu }_i{\Vert }_{{{\text{Upsigma}}_{\mathbf{i}}}^{-1}}^2\right)$$

$$\mathcal{J}(\mathbf{x})=\frac{1}{2}\Vert \mathcal{H}(\mathbf{x})-\mathbf{y}{\Vert }_{{\mathbf{R}}^{-1}}^2-\log \left(\sum _{i=1}^{\text{textsc}{n}_{\text{c}}}\frac{{\tau }_i}{\sqrt{|{\text{Upsigma}}_i|}}\exp \left(-\frac{1}{2}\Vert \mathbf{x}-{\mu }_i{\Vert }_{{{\text{Upsigma}}_{\mathbf{i}}}^{-1}}^2\right)\right)$$

$${\nabla }_{\mathbf{x}}\mathcal{J}(\mathbf{x})$$

$$\{\begin{array}{ll}O(\text{textsc}{n}_{\text{var}}^2)+O(\text{textsc}{n}_{\text{var}}^2)+O(1)=O(\text{textsc}{n}_{\text{var}}^2)& \text{full GMM covariances, and full Gaussian kernel}\\ O(\text{textsc}{n}_{\text{var}}^2)+O(\text{textsc}{n}_{\text{var}})+O(1)=O(\text{textsc}{n}_{\text{var}}^2)& \text{full GMM covariances, and diagonal Gaussian kernel}\\ O(\text{textsc}{n}_{\text{var}})+O(\text{textsc}{n}_{\text{var}}^2)+O(1)=O(\text{textsc}{n}_{\text{var}}^2)& \text{diagonal GMM covariances, and full Gaussian kernel}\\ O(\text{textsc}{n}_{\text{var}})+O(\text{textsc}{n}_{\text{var}})+O(1)=O(\text{textsc}{n}_{\text{var}})& \text{diagonal GMM covariances, and diagonal Gaussian kernel}\end{array}$$

$$T_{s}=\begin{cases}O\Bigl{(}\left(b_{s}+m_{s}\,\textsc{n}_{\textnormal{ens}}\right)\textsc{n}_{\textnormal{var}}\Bigr{)}&\text{\small diagonal prior covariances, and diagonal Gaussian kernel}\\ O\Bigl{(}\left(b_{s}+m_{s}\,\textsc{n}_{\textnormal{ens}}\right)\textsc{n}_{\textnormal{var}}^{2}\Bigr{)}&\text{\small otherwise}\end{cases}$$

$$\{\begin{array}{ll}O(\text{textsc}{n}_{\text{var}})+O(m\text{textsc}{n}_{\text{var}}^2)+O(\text{textsc}{n}_{\text{var}})+O(1)=O(m\text{textsc}{n}_{\text{var}}^2)& \text{diagonal prior covariances}\\ O(\text{textsc}{n}_{\text{var}})+O(m\text{textsc}{n}_{\text{var}}^2)+O(\text{textsc}{n}_{\text{var}}^2)+O(1)=O(m\text{textsc}{n}_{\text{var}}^2)& \text{non-diagonal prior covariances}\end{array}$$

$$T_{s}=\begin{cases}\left.\begin{aligned} O\left(\left(b_{s}+m_{s}\,\textsc{n}_{\textnormal{ens}}\right)\textsc{n}_{\textnormal{var}}\right)\,;\quad&\text{\begin{minipage}{150.00023pt} \small diagonal or spherical covariances of GMM, and proposal density \end{minipage} }\\ O\left(\left(b_{s}+m_{s}\,\textsc{n}_{\textnormal{ens}}\right)\textsc{n}_{\textnormal{var}}^{2}\right)\,;\quad&\text{ \small otherwise}\\ \end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small Gaussian\\ proposal\end{minipage}}&\\ \hskip 28.45274pt\\ \left.\begin{aligned} O\left(\left(b_{s}+m_{s}\,\textsc{n}_{\textnormal{ens}}\right)\,m\,\textsc{n}_{\textnormal{var}}^{2}\right)\,\quad\quad\end{aligned}\right\}\text{\begin{minipage}{100.00015pt}\small Hybrid Monte Carlo\end{minipage}}&\end{cases}$$

$$T_{p}=\begin{cases}\left.\begin{aligned} O\left(\frac{\textsc{n}_{\textnormal{c}}}{p}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\textsc{n}_{\textnormal{var}}\right)\,;\quad&\text{\begin{minipage}{150.00023pt} \small diagonal, or spherical covariances of GMM, and proposal density \end{minipage} }\\ O\left(\frac{\textsc{n}_{\textnormal{c}}}{p}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\textsc{n}_{\textnormal{var}}^{2}\right)\,;\quad&\text{ \small otherwise}\\ \end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small Gaussian\\ proposal\end{minipage}}&\\ \hskip 28.45274pt\\ \left.\begin{aligned} O\left(\frac{\textsc{n}_{\textnormal{c}}}{p}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\,m\,\textsc{n}_{\textnormal{var}}^{2}\right)\,\quad\quad\end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small Hybrid\\ Monte Carlo\end{minipage}}&\end{cases}$$

$$S=\frac{T_s}{T_p}=\frac{\left(b_s+m_s\text{textsc}{n}_{\text{ens}}\right)}{\frac{\text{textsc}{n}_{\text{c}}}{p}\left(b_s+m_s\frac{\text{textsc}{n}_{\text{ens}}}{\text{textsc}{n}_{\text{c}}}\right)}=\frac{p\left(b_s+m_s\text{textsc}{n}_{\text{ens}}\right)}{\text{textsc}{n}_{\text{c}}\left(b_s+m_s\frac{\text{textsc}{n}_{\text{ens}}}{\text{textsc}{n}_{\text{c}}}\right)}$$

$$E=\frac{S}{p}=\frac{\left(b_s+m_s\text{textsc}{n}_{\text{ens}}\right)}{\text{textsc}{n}_{\text{c}}\left(b_s+m_s\frac{\text{textsc}{n}_{\text{ens}}}{\text{textsc}{n}_{\text{c}}}\right)}$$

$$S=\{\begin{array}{ll}\mathbf{p}& ;p\le \text{textsc}{n}_{\text{c}}\\ \text{textsc}{n}_{\text{c}}& ;p>\text{textsc}{n}_{\text{c}}\end{array},\text{and,}E=\frac{S}{p}=\{\begin{array}{ll}1& ;p\le \text{textsc}{n}_{\text{c}}\\ \frac{nc}{p}& ;p>\text{textsc}{n}_{\text{c}}\end{array}.$$

$$\begin{cases}\Bigl{(}t_{s}+t_{w}\left((2\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right)\Bigr{)}\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small diagonal, or spherical GMM covariances \end{minipage} }\\ \Bigl{(}t_{s}+t_{w}\left((\frac{\textsc{n}_{\textnormal{var}}^{2}}{\textsc{n}_{\textnormal{c}}}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right)\Bigr{)}\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small tied GMM covariances \end{minipage} }\\ \Bigl{(}t_{s}+t_{w}\left((\textsc{n}_{\textnormal{var}}^{2}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right)\Bigr{)}\log{(p)}\,;&\text{\begin{minipage}{100.00015pt} \small full GMM covariances \end{minipage} }\end{cases}$$

$$\begin{cases}\Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}+2\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small diagonal, or spherical GMM covariances \end{minipage} }\\ \Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}+\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+1\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small tied GMM covariances \end{minipage} }\\ \Bigl{(}2t_{s}+t_{w}\left[(\frac{\textsc{n}_{\textnormal{ens}}\times\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+\textsc{n}_{\textnormal{var}}^{2}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}\log{(p)}\,;&\text{\begin{minipage}{100.00015pt} \small full GMM covariances \end{minipage} }\end{cases}$$

$$pT_{p}=\begin{cases}&\left.\begin{aligned} \begin{cases}\Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}+2\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}+O\left(\textsc{n}_{\textnormal{c}}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\textsc{n}_{\textnormal{var}}\right)\,;\quad&\text{\begin{minipage}{150.00023pt} \small diagonal, or spherical covariances of GMM, and proposal density \end{minipage} }\\ \hskip 28.45274pt&\\ \Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}+\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+1\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}+O\left(\textsc{n}_{\textnormal{c}}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\textsc{n}_{\textnormal{var}}^{2}\right)\,;\quad&\text{\begin{minipage}{120.00018pt} \small tied covariances of GMM \end{minipage} }\\ \hskip 28.45274pt&\\ \Bigl{(}2t_{s}+t_{w}\left[(\frac{\textsc{n}_{\textnormal{ens}}\times\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+\textsc{n}_{\textnormal{var}}^{2}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}+O\left(\textsc{n}_{\textnormal{c}}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\textsc{n}_{\textnormal{var}}^{2}\right)\,;\quad&\text{ \small otherwise}\\ \end{cases}\end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small Gaussian\\ proposal\end{minipage}}\\ &\vspace{1cm}\textstyle\\ &\left.\begin{aligned} O\left(\textsc{n}_{\textnormal{c}}\left(b_{s}+m_{s}\,\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}\right)\,m\,\textsc{n}_{\textnormal{var}}^{2}\right)\,+\,\begin{cases}\Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}+2\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}\,;&\text{\begin{minipage}{130.0002pt} \small diagonal, or spherical GMM covariances \end{minipage} }\\ \hskip 28.45274pt\\ \Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}+\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+1\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small tied GMM covariances \end{minipage} }\\ \hskip 28.45274pt\\ \Bigl{(}2t_{s}+t_{w}\left[(\frac{\textsc{n}_{\textnormal{ens}}\times\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+\textsc{n}_{\textnormal{var}}^{2}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small full GMM covariances \end{minipage} }\end{cases}\quad\end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small HMC \end{minipage}}\end{cases}$$

$$T_{o}=\begin{cases}&\left.\begin{aligned} \begin{cases}\Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}+2\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}+O\left(\left(\textsc{n}_{\textnormal{c}}-1\right)b_{s}\,\textsc{n}_{\textnormal{var}}\right)\,;\quad&\text{\begin{minipage}{150.00023pt} \small diagonal, or spherical covariances of GMM, and proposal density \end{minipage} }\\ \hskip 28.45274pt\\ \Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}+\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+1\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}+O\left(\left(\textsc{n}_{\textnormal{c}}-1\right)b_{s}\,\textsc{n}_{\textnormal{var}}^{2}\right)\,;\quad&\text{\begin{minipage}{120.00018pt} \small tied covariances of GMM \end{minipage} }\\ \hskip 28.45274pt\\ \Bigl{(}2t_{s}+t_{w}\left[(\frac{\textsc{n}_{\textnormal{ens}}\times\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+\textsc{n}_{\textnormal{var}}^{2}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}+O\left(\left(\textsc{n}_{\textnormal{c}}-1\right)b_{s}\,\textsc{n}_{\textnormal{var}}^{2}\right)\,;\quad&\text{\small otherwise}\\ \end{cases}\end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small Gaussian\\ proposal\end{minipage}}\\ \vspace{1cm}\\ &\left.\begin{aligned} O\left(\left(\textsc{n}_{\textnormal{c}}-1\right)\,b_{s}\,m\,\textsc{n}_{\textnormal{var}}^{2}\right)\,+\,\begin{cases}\Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}+2\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}\,;&\text{\begin{minipage}{140.00021pt} \small diagonal, or spherical GMM covariances \end{minipage} }\\ \hskip 28.45274pt\\ \Bigl{(}2t_{s}+t_{w}\left[\left(\left(\frac{\textsc{n}_{\textnormal{ens}}+\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+1\right)\textsc{n}_{\textnormal{var}}\,+1\right)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}\,;&\text{\begin{minipage}{140.00021pt} \small tied GMM covariances \end{minipage} }\\ \hskip 28.45274pt\\ \Bigl{(}2t_{s}+t_{w}\left[(\frac{\textsc{n}_{\textnormal{ens}}\times\textsc{n}_{\textnormal{var}}}{\textsc{n}_{\textnormal{c}}}+\textsc{n}_{\textnormal{var}}^{2}+\textsc{n}_{\textnormal{var}}+1)\,\textsc{n}_{\textnormal{c}}\right]\Bigr{)}p\log{(p)}\,;&\text{\begin{minipage}{120.00018pt} \small full GMM covariances \end{minipage} }\end{cases}\quad\end{aligned}\right\}\text{ \begin{minipage}{150.00023pt}\small HMC \end{minipage}}\end{cases}$$

$$W(W,p)=\begin{cases}k\,t_{w}\,\textsc{n}_{\textnormal{ens}}\,\textsc{n}_{\textnormal{var}}\,p\,\log{(p)}\,;\quad&\text{\begin{minipage}{180.00027pt} \small diagonal, or spherical GMM covariances, and {HMC, or Gaussian proposal with diagonal covariance} \end{minipage} }\\ \hskip 28.45274pt&\\ k\,t_{w}\,\left(\textsc{n}_{\textnormal{ens}}+\textsc{n}_{\textnormal{var}}\right)\,\textsc{n}_{\textnormal{var}}\,p\,\log{(p)}\,;\quad&\text{\begin{minipage}{160.00024pt} \small tied covariances of GMM, and {HMC, or Gaussian proposal} \end{minipage} }\\ \hskip 28.45274pt&\\ k\,t_{w}\,\left(\textsc{n}_{\textnormal{ens}}\,+\textsc{n}_{\textnormal{var}}\,\textsc{n}_{\textnormal{c}}\right)\,\textsc{n}_{\textnormal{var}}\,p\,\log{(p)}\,;\quad&\text{\begin{minipage}{160.00024pt} \small full covariances of GMM, and {HMC, or Gaussian proposal} \end{minipage} }\\ \end{cases}$$

$$\begin{array}{r}\{({\tau }_i;{\mu }_i,{\sigma }_i^2){\}}_{i=1,\dots ,5}=\{(0.09;-6.0,0.20),(0.19;-2.5,0.28),(0.09;0.0,0.08),\\ (0.28;2,5,0.24),(0.15;6.0,0.28),(0.15;6.5,0.08),\\ (0.03;7.5,0.12),(0.02;8.0,0.04)\}.\end{array}$$

$$\centering\mathcal{P}(\mathbf{y}|\mathbf{x})=\frac{1}{\sqrt{2.2}\,\sqrt{2\pi}}\,\,\exp{\Bigl{(}-\frac{1}{2}\,\frac{\left(\mathbf{x}-\mathbf{y}\right)^{2}}{2.2}\Bigr{)}}\,.\@add@centering$$

$$\begin{array}{r}\{({\tau }_i;{\mu }_i,{\sigma }_i^2){\}}_{i=1,\dots ,7}=\{(0.111;-5.78,0.123),(0.177;-2.49,0.223),\\ (0.045;-1.49,0.001),(0.065;0.12,0.061),(0.146;2.05,0.032),\\ (0.225;2.78,0.148),(0.231;6.12,0.164)\}.\end{array}$$

Full text

The parallel nature of the prposed algorithm makes the it a strong candidate for practical and large scale applications.

**Keywords: ** Parallel programming, Medical image reconstruction, Inverse problems, Bayes’ theorem, Markov chain Monte-Carlo, Hamiltonian Monte-Carlo

1 Introduction

Signal retrieval from noisy measurements (observations) involves solving an inverse problem. Inverse problems are essential in many fields such as image reconstruction or retrieval, tomography, weather prediction, and other predictions based on space-time models. The solution of inverse problems in case of space-time models usually employs a data assimilation (DA)[6, 2, 4] methodology. DA refers to the process of fusing information about a physical system obtained from different sources in order to produces more accurate conclusions about the physical system of concern.

Two approaches are widely employed to solve an inverse problem. The first approach is a variational approach that involves solving an optimization problem with a regularized solution. The second approach is the statistical formulation of the DA problem which incorporates a prior distribution that encapsulates the knowledge about the system produced by the model, prior to the incorporation of any other source of information. Given the prior information, a likelihood function, the posterior is formulated as best estimate of the truth.

Markov Chain Monte-Carlo (MCMC) is one of the most powerful simulation techniques for sampling a high-dimensional probability distribution, given only it’s shape function, without the intrinsic need to the associated scaling factor.

HMC sampling filter [3] is an accelerated Markov chain Monte-Carlo (MCMC) algorithm for solving the non-Gaussian sequential DA “filtering problem”. This algorithm works by sampling the posterior distribution to produce description of the system state along with associated uncertainty. Specifically, these algorithms follow a Hamiltonian Monte-Carlo (HMC) approach to sample the posterior.

Cluster sampling filters ($\mathcal{C}\mathrm{\ell}$HMC, and MC-$\mathcal{C}\mathrm{\ell}$HMC) [1] are developed as extension of the Hamiltonian Monte-Carlo (HMC) sampling filter presented in [3] where the true (unknown) prior distribution is approximate using a Gaussian mixture model (GMM).

Given the current computational power, it is natural to try to run Monte-Carlo simulations in parallel. However, Markov chains in general have to satisfy the so-called the “Markovian” property which makes the chain generation an inherently sequential problem. This restriction is mainly posed by the transition density function used to generate a proposal state given the current state of the chain.

Two approaches have gained wide popularity to parallelize MCMC samplers. The parallel-chain approach proceeds by running several chains in parallel from different initial states. The main disadvantage of this approach is that the burn-in stage has to be carried out by all chains independently, which limits the efficiency gained by running the chains on different processors. Another difficulty with this approach is the aggregation of samples generated on different processors such that the combined ensemble correctly represents the mass of the target distribution. The second approach is to parallelize a single chain.

The parallel chain approach turns out to be surprisingly effective in practice. Moreover, if sufficient information about the geometry of the target distribution is available, we can guide the parallel chains to sample effectively from the target distribution.

The accuracy of $\mathcal{C}\mathrm{\ell}$HMC filters to handle nonlinearity in both model dynamics, and observational mapping operator, puts it on the right direction of applicability to practical problems. The cost of serial $\mathcal{C}\mathrm{\ell}$HMC is nearly similar to the cost of the original HMC sampling filter, however the MC-$\mathcal{C}\mathrm{\ell}$HMC algorithm is naturally parallelizable.

Following a Bayesian approach, $\mathcal{C}\mathrm{\ell}$HMC algorithm can be easily modified and applied for image retrieval given noisy image and a probabilistic representation of prior knowledge. This can be very useful in settings where several medical snapshots are collected for the same object, e.g. a tumor, of different resolution or uncertainty levels.

Mathematical regularization is amongst the most popular methods for image reconstruction from noisy sources [10]. Among the regularization methods, the Tikhonov scheme is most popular due to the Gaussianity assumption about data noise, and the easiness to incorporate prior information. Disbite simplicity, the perfomace of this approach is highly influenced by the choice of the regularization parameter.

Widely used methodologies for solving the Bayesian image retrieval problem include the algorithms discussed in [9, 5]. In [5], the authors investigate statistical image reconstruction (SIR) with regularization based on the Markov random field (MRF) model.

While, regularization approach is popular, it is sensitive to the choice of the regularization parameter.

Our main interest here is to develop highly accurate parallel Bayesian sampling algorithms that can be efficiently used for solving large-scale inverse problems, and show that they are suitable for a wide spectrum of applications including medical image retrieval.

In this work, we develop parallel cluster sampling algorithms, and show that a multi-chain version is embarrassingly parallel, and can be used efficiently for medical image retrieval amongst other applications. The approach discussed in this work does not require regularization, and is designed to work in both Gaussian, and non-Gaussian case, where the computationl expense is minimized via parallelization. Specifically, in this paper, we focus on describing the complexity analysis of a specific scenario where the MC-$\mathcal{C}\mathrm{\ell}$HMC is parallelized by running several chains in parallel to sample the posterior distribution. The algorithm proceeds by running several Markov chains in parallel such that the number of chains is specified by the the number of components in the mixture model. We will focus on the case where an ensemble of states, generated from an unknown prior distribution is available, and the likelihood function relating observations to target states is either a linear or a nonlinear map. The prior distribution is approximated using a Gaussian mixture distribution which parameters are approximated based on the given prior ensemble by running an expectation maximization (EM) algorithm.

In Section 2 we review the general iterative and Bayesian frameworks for inverse problems and image reconstruction. Section 3 formulates the problem, and reviews the $\mathcal{C}\mathrm{\ell}$HMC filter formulation. In Section 4 we discuss opportunities for parallelization of $\mathcal{C}\mathrm{\ell}$HMC. Section 5 presents a detailed complexity analysis of the proposed parallel version of $\mathcal{C}\mathrm{\ell}$HMC filter. Numerical results are presented in Section 6. Conclusions and future works are drawn in Section 7

2 Iterative and Bayesian Image Reconstruction

As mentioned in Section 1, one of the most popular iterative reconstruction algorithms is Regularization-based algorithms. For the sake of completeness, we review the Tikhonov regularization approach [10, 8] next, then we present the Bayesian formulation.

The Tikhonov regularization approach involves solving the following optimization problem:

[TABLE]

where $\alpha$ is the regularization parameter, and $\mathbf{C}$ is the regularization matrix and it can be chosen in many clever ways. Here $\mathcal{H}$ is an observation operator that maps the model space to the observation space. If the target state is directly observed then $\mathcal{H}=\mathbf{I}$, where $\mathbf{I}$ is the identity operator. The weighted norm in Equation (1) is described as follows:

[TABLE]

The traditional approach to regularization is the variational formulation in which equation (1) is minimized w.r.t $\mathbf{x}$. Usually, derivative-based iterative minimizaiton algorithms are employed to solve the problem described by (1). The derivative of the objective function $T(\mathcal{H}(\mathbf{x}))$ w.r.t the parameter $]x$ is given by:

[TABLE]

where ${\left[\partial\mathcal{H}\right]}^{*}$ is the adjoint of the derivative, e.g. the Jacobian-transpose, of the observation operator $\mathcal{H}$. In the case of a linear observaiton operator this is simply the transpose of the observation operator.

In the statistical approach we infer the underlying state $\mathbf{x}$ based on a formulation constructed using Bayesian theory, where the goal is to represent the state as a random variable which distribution is of interest. Assume $\mathcal{P}^{\rm b}(\mathbf{x})$ is a prior probability density function (PDF) representing prior knowladge about the state $\mathbf{x}$. Assume also that $\mathcal{P}(\mathbf{y}|\mathbf{x})$ is the data likelihood function that describes the observational error distribution. Using Bayes’ theorem, the probability of the state $\mathbf{x}$ given the collected meauserments is characterized by the posterior distribution:

[TABLE]

A common practice is to assume that the prior distribution is, generally speaking, a multivariate normal (Gaussian) distribution centered around a forecasted or a first-guess state $\mathbf{x}^{\textnormal{b}}$, and have a predefined covariance structure $\mathbf{C}$, i.e. $\mathbf{x}\sim\mathcal{N}(\mathbf{x}^{\textnormal{b}},C)$.

For Gaussian priors and consequently Gaussian posteriors, the variational approach corresponds to finding the maximum aposterior estimate (MAP) of the posterior PDF. The MAP is the maximizer of the posterior PDF, or equivalently, the minimizer of its negative logarithm $-\log{\left(\mathcal{P}(\mathbf{x}|\mathbf{y})\right)}$. Following Tikhonov regularization approach (1), and assuming Gaussian noise, the likelihood function reads:

[TABLE]

Without loss of generality, if we assume $\mathbf{x}\sim\mathbf{N}(0,\mathbf{C})$, the prior would be on the form

[TABLE]

where $\mathbf{B}$ is the precision matrix, i.e. the inverse of the covariance $\mathbf{B}=\mathbf{C}^{-1}$. The MAP estimator in this formulation is the minimizer of

[TABLE]

This shows the equivalence between Tikhonov regularization approach with the Bayesian formulation in the Gaussian linear settings.

In the Bayesian approach, once the posteiror is constructed, a sampling mechanism is usually employed to estimate all the desired statistics of the posterior PDF, such as the posteiror mean $\mathbb{E}(\mathbf{x}|\mathbf{y})$ that can be used as a reliable estimate of the state given the data. Moreover, the generated ensemble can be used to estimate the posteiror covariance that can be used as a proxy prior error covariance for future applications of the inverse problem. Sampling the posterior PDF is usually carried out following a Monte-Carlo approach. The most powerful Monte-Carlo sampling methodology is the general family Markov-Chain Monte Carlo (MCMC) samplers. Sampling high dimensional distribution however is a very expensive process, and requires parallel efficient implementation to be considered practical. As explained in the next Section, MCMC is not limited to Gaussian or linear settings, and can be very efficient if implemented in parallel.

3 Cluster Sampling Filter

Let $\mathbf{x}\in\mathbb{R}^{\textsc{n}_{\textnormal{var}}}$ is a discretized approximation of the true state of the model, for example the entensities of an image pixels.

The prior distribution $\mathcal{P}^{\rm b}(\mathbf{x})$ encapsulates the knowledge about the system state before additional information is incorporated. The likelihood function $\mathcal{P}(\mathbf{\mathbf{y}}|\mathbf{x})$ quantifies the deviation of the prediction of model observations from the collected measurements $\mathbf{y}\in\mathbb{R}^{\textsc{n}_{\textnormal{obs}}}$, where $\textsc{n}_{\textnormal{obs}}\leq\textsc{n}_{\textnormal{var}}$.

From Bayes’ theorem, the posterior distribution $\mathcal{P}^{\rm a}(\mathbf{x})$ reads:

[TABLE]

where $\mathcal{P}^{\rm b}(\mathbf{x})$ is the prior distribution, $\mathcal{P}(\mathbf{y}|\mathbf{x})$ is the likelihood function. $\mathcal{P}(\mathbf{y})$ acts as a scaling factor and is ignored in in the MCMC context.

Assuming the prior distribution is approximated by a Gaussian Mixture distribution, the prior takes the form:

[TABLE]

where the weight $\tau_{i}$ quantifies the probability that an ensemble member $\mathbf{x}[e]$ belongs to the $i^{th}$ component, and $(\mu_{i},\,\Upsigma_{i})$ are the mean and the covariance matrix associated with the $i^{th}$ component of the mixture model. Here $\mathbf{x}\in\mathtt{R}^{\textsc{n}_{\textnormal{var}}}$, where $\textsc{n}_{\textnormal{var}}$ the dimension of the target state space.

Assuming the observation errors are characterized by a Gaussian distribution $\mathcal{N}(0,\mathbf{R})$, the likelihood reads:

[TABLE]

where $\mathbf{R}\,\in\mathtt{R}^{\textnormal{m}\times\textnormal{m}}$ is the observation error covariance matrix, and $\mathcal{H}:\mathtt{R}^{\textsc{n}_{\textnormal{var}}}\rightarrow\mathtt{R}^{\textnormal{m}}$ is the observation operator that maps the state space to the observation space. Here $\mathbf{y}\in\mathtt{R}^{\textnormal{m}}$ is the observation vector.

From Equations (8), (9), (10), the posterior takes the form:

[TABLE]

This mixture distribution may not correspond to a Gaussian mixture in general if the observation operator is a nonlinear map.

The negative-logarithm (negative-log) of the posterior distribution kernel (11) is required by the HMC sampling algorithm. Specifically, the posterior negative-log is viewed as the potential energy function in the extended Hamiltonian phase space. The posterior negative-log is given by:

[TABLE]

The derivative of the posterior negative-log reads:

[TABLE]

In this work, we sample the posterior distribution (11) following a parallel chains approach given only ensemble of states generated from the prior distribution.

4 Parallelization of $\mathcal{C}\mathrm{\ell}$HMC Filter

In this section, we present a brief review of the MCMC sampling algorithm, and discuss opportunities of running $\mathcal{C}\mathrm{\ell}$HMC filter on parallel architecture. We discuss the parallelism of $\mathcal{C}\mathrm{\ell}$HMC filter even if the Hamiltonian system is replaced with a Gaussian proposal kernel build around the forecast.

4.1 Markov Chain Monte-Carlo (MCMC)

MCMC is a sampling scheme capable of producing ensembles from an arbitrary distribution given it’s shape function, without the need for the associated scaling factor.

The choice of the proposal kernel has the greatest influence on the performance of the sampler. Here, we chose two proposal kernels; a) a Gaussian density function centered around the current state of the chain, b) Hamiltonian Monte Carlo (HMC). The standard MCMC sampler is described in Algorithm 1.

The standard MCMC algorithm1 generally suffers from random walk behaviour, slow convergence to the target density, low acceptance rate, and slow space exploration. Moreover, the generated samples are highly correlated when the vanilla MCMC algorithm is used. Many of these problems can be addressed by using Hybrid Monte Carlo (HMC). HMC uses a Hamiltonian system which plays the role of the proposal density.

4.2 The multi-chain MCMC algorithm (MC-MCMC)

Although the traditional MCMC algorithm is inherently serial, there are several modifications that can be made to allow for parallelization. In our approach, instead of constructing a single long Markov chain to produces $\textsc{n}_{\textnormal{ens}}$ samples, we generate several shorter Markov chains and divide the ensemble size $\textsc{n}_{\textnormal{ens}}$ over these chains.

The constructed chains can run in parallel to sample different regions of the target distribution independently.

The parallel (MC-MCMC) sampler starts by running an Expectation Maximization step to build a Gaussian Mixture Model (GMM) approximation of the prior distribution.

Once the GMM is constructed on the root processor, the GMM information is broadcasted to all the working nodes. Of course, if the number of processors is exactly the same as the number of components, each node is assigned one chain. If the number of processors is less than the number of components/chains, we can assign several chains to each processor, e.g. based on the local ensemble sizes or simply in a round-robin fashion. Caution has to be exercised to maintain a balanced load. Once all chains have generated their assigned local samples, the ensembles are gathered to the root processor and returned as output. Of course parallel output can be considered to reduce the communication overhead. The steps of the proposed MC-MCMC scheme is detailed in Algorithm 2.

By running a Markov chain starting at each component of the mixture distribution, we ensure that the proposed algorithm navigates all modes of the posterior distribution, and covers all regions of high probability.

While parallelization of HMC itself can be considered for further increase in the performance, it has been avoided for clarity and to simplify the idea. We elected to use the Master-Slave parallel pattern, where a master core sends the information required for creating the individual chains.

5 Complexity Analysis

In this section we provide a detailed theoretical discussion of the computational cost of the proposed parallel algorithm. Since sampling from a Gaussian distribution is essential for the two flavors tested here, we start with discussion the computational cost of sampling from a Gaussian distribution

5.1 Cost of sampling a Gaussian $\mathcal{N}(\mu,\,\Sigma)$

A scalar normal distribution can be sampled using many accurate algorithms such as the Mersenne Twister, Box-Muller transform, Marsaglia polar method, and Ziggurat algorithm. The least expensive is Ziggurat algorithm, where a typical value produced only requires the generation of one random floating-point value and one random table index, followed by one table lookup, one multiply operation and one comparison. The cost of generating an $\textsc{n}_{\textnormal{var}}-$dimensional standard normal random vector $\mathbf{x}\in\mathbb{R}^{\textsc{n}_{\textnormal{var}}}$ is $O(\textsc{n}_{\textnormal{var}})$. To generate a multivariate normal (MVN) random vector Ziggurat algorithm from a general Gaussian distribution $\mathcal{N}(\mu,\,\Sigma)$ the following steps are required:

Factorization of the covariance matrix $\Sigma$, to generate $\Sigma^{\frac{1}{2}}$, e.g. using Cholesky decomposition 2. 2.

Draw a standard normal random vector $\mathbf{y}\in\mathcal{N}\mathbf{(0\,I})$, 3. 3.

Scaling: $\mathbf{x}=\Sigma^{\frac{1}{2}}\mathbf{y}+\mu$

If the covariance matrix $\Sigma$ is diagonal, the factorization costs $O(\textsc{n}_{\textnormal{var}})$, while the cost of Cholesky decomposition in general is $O(\textsc{n}_{\textnormal{var}}^{3})$. If the covariance matrix $\Sigma$ is diagonal, the cost of the scaling step is $O(\textsc{n}_{\textnormal{var}})$, otherwise the scaling cost is $O(\textsc{n}_{\textnormal{var}}^{2})$.

5.2 Cost of MCMC with a Gaussian proposal density

Assuming independent observations, evaluating the posterior PDF at a given state requires the evaluation of $\textsc{n}_{\textnormal{c}}+1$ matrix-vector products. The cost of evaluating the posterior PDF is $O\left((\textsc{n}_{\textnormal{c}}+1)\,\textsc{n}_{\textnormal{var}}\right)=O(\textsc{n}_{\textnormal{var}})$ if the covariance matrices of the components the GMM prior are diagonal, otherwise, the cost is $O(\textsc{n}_{\textnormal{var}}^{2})$.

Cost of one MCMC step:

In the presence of a Gaussian kernel, each step of the chain construction requires the following:

One MVN random vector drawn from the Gaussian kernel, 2. 2.

Two function (posterior PDF) evaluations, 3. 3.

One draw from a uniform random distribution $\mathcal{U}(0,1)$, 4. 4.

One scalar comparison in the Metropolis-Hastings step.

The cost of one MCMC step can be summarized as follows:

[TABLE]

To find the total cost of MCMC sampling algorithm, we need to evaluate the total number of steps in the Markov chain. If the ensemble size is $\textsc{n}_{\textnormal{ens}}$, we need to construct a chain of length $b_{s}+m_{s}\times\textsc{n}_{\textnormal{ens}}$, where $b_{s}$ are burn-in steps carried out to achieve convergence to the stationary distribution, and $m_{s}$ are mixing steps introduced to improve the sampler mixing and reduce the correlation between selected ensembles

Total cost of MCMC sampling:

The total cost of a serial MCMC with with a Gaussian proposal density reads:

[TABLE]

5.3 Cost of HMC

The HMC sampler requires the evaluation of the gradient of the negative-log of the target distribution.

Cost of one HMC step:

In the case of HMC, each step of the chain construction requires the following:

One draw of a momentum $\mathbf{p}\sim\mathcal{N}(0,\mathbf{M})$ costs $O(\textsc{n}_{\textnormal{var}})$ with diagonal $M$. 2. 2.

Forward propagation of the pair $(\mathbf{p},\,\mathbf{x})$.

The cost of the second step is dominated by the cost of evaluating a Jacobian-vector product $O(\textsc{n}_{\textnormal{var}}^{2})$. With step parameters $T=m\times h$, the cost is $O(m\,\textsc{n}_{\textnormal{var}}^{2})$, where $m$ is the number of steps in the Hamiltonian trajectory. Evaluating loss of energy requires evaluating the negative-log of the posterior shape function. Again, the cost of evaluating the posterior PDF is $O\left((\textsc{n}_{\textnormal{c}}+1)\,\textsc{n}_{\textnormal{var}}\right)=O(\textsc{n}_{\textnormal{var}})$ if the covariance matrices of the components the GMM prior are diagonal, otherwise, the cost is $O(\textsc{n}_{\textnormal{var}}^{2})$.

For Hamiltonian Monte Carlo (HMC), Assuming a diagonal mass matrix $\mathbf{M}$ the cost of one step of the chain is

[TABLE]

Total cost of HMC sampling:

Following the discussion above, the cost of serial HMC sampler is T_{s}=O\Bigl{(}\left(b_{s}+m_{s}\,\textsc{n}_{\textnormal{ens}}\right)\,m\,\textsc{n}_{\textnormal{var}}^{2}\Bigr{)}.

5.4 Cost of MC-MCMC sampling

The serial complexity $T_{s}$ of MC-MCMC is summarized as follows:

[TABLE]

Parallel cost:

The parallel complexity can be studied clearly under a simplified (ideal) assumption, where each chain is to sample $\frac{\textsc{n}_{\textnormal{ens}}}{\textsc{n}_{\textnormal{c}}}$ ensemble points. In this case, since we have $\textsc{n}_{\textnormal{c}}$ chains the parallel cost (discarding the communication cost) of MC-MCMC is given by:

[TABLE]

Speedup:

The speedup of MC-MCMC is given by:

[TABLE]

Parallel efficiency:

The parallel efficiency of MC-MCMC is given by:

[TABLE]

If we discard the burn-in stage, i.e. set $b_{s}=0$, the speedup, and the parallel efficiency simplify to:

[TABLE]

It follows that both speedup, and parallel efficiency are independent from the state space dimension $\textsc{n}_{\textnormal{var}}$.

Communication overhead:

Assuming serial GMM run on the root node, the cost of broadcasting GMM information to all nodes is the cost of broadcasting, the means, the covariance and/or precision matrices, and the weights.

Assuming a linear communication model [7], and assuming that $t_{s}$, and $t_{w}$ are the startup, and the per-word transfer times respectively, the cost of broadcasting GMM information to $p$ nodes is given by:

[TABLE]

After sampling in parallel, the collected ensembles are gathered on the root node, at a cost $\left(t_{s}+t_{w}\left(\textsc{n}_{\textnormal{ens}}\times\textsc{n}_{\textnormal{var}}\right)\right)\log{(p)}\,.$ Consequently, the total communication cost reads:

[TABLE]

Total parallel cost:

It follows immediately from the discussion above, that the total parallel cost of MC-MCMC sampling algorithm simplifies to:

[TABLE]

Total overhead:

The total overhead function ($T_{o}=PT_{p}-T_{s}$) of MC-MCMC reads:

[TABLE]

Isoefficiency:

Assuming the burn-in stage is discarded, i.e. set $b_{s}=0$. With $\textsc{n}_{\textnormal{ens}}\geq 2\textsc{n}_{\textnormal{c}}$, the isoefficiency function $W(p)=\frac{E}{1-E}T_{o}=kT_{o}$ simplifies to (the dominant terms):

[TABLE]

6 Numerical Experiments and Performance

In this section we present numerical experiments to assess the complexity analysis provided in Section 5. As mentioned above, the speedup, and parallel efficiency are independent of problem dimensionality. We discuss the computational cost and the performance using one dimensional examples.

To execute the Markov chains in parallel, we used the MPI4Py package, which provides bindings of the Message Passing Interface (MPI) standard for the Python programming language, allowing any Pythonprogram to exploit multiple processors.

we present numerical experiments to assess the complexity analysis provided in Section 5. Specifically, we show numerical results for the parallel $\mathcal{C}\mathrm{\ell}$HMC flavors discussed in this paper using one dimensional synthetic example, and a two dimensional image retrieval experiment.

6.1 One-dimensional example:

Following the strategy described in [1], we start with a synthetic prior ensemble generated from a GMM with $\textsc{n}_{\textnormal{c}}=5$. A GMM approximation of the true prior probability distribution is constructed using the EM algorithm . The model selection criterion used here is AIC. The parameters of the true GMM prior are:

[TABLE]

Assume the observation errors follow Gaussian distribution with zero mean, and variance $2.2$. Assuming a synthetic observation $\mathbf{y}=-1.0$, the observation likelihood function given by:

[TABLE]

The generated GMM approximation of the prior has $\textsc{n}_{\textnormal{c}}=7$ and the following parameters:

[TABLE]

Figure 1 shows the results of sampling a one dimensional posterior (11) with seven components. Since we are mainly interested in asymptotic behaviour of the sampler, and to avoid the effect of sampling error, the ensemble size here is set to $1000$. Despite the good mass distribution resulting from the use of a Gaussian density with the serial MCMC sampler 1(a), the acceptance rate is $\approx 45\%$ resulting in a large amount of wasted calculations. On the other hand, the parallel MC-MCMC with Gaussian proposal kernel 1(b) improves the acceptance rate to $\approx 81\%$. The acceptance rate is increased due to the local adjustment of the sampler hyperparameters based on the local ensemble under the corresponding prior component in the mixture.

While the Gaussian-based MCMC sampler represents an acceptable mass distribution, it suffers from random walk behaviour leading to the demonstrated low acceptance rate. One the other hand HMC sampling results in general in high acceptance rate. Unfortunately, as explained by the results in Figure 1(c) is unable to sample all probability modes. The acceptance rate in the serial $\mathcal{C}\mathrm{\ell}$HMC sampler here is $96\%$. By running $\mathcal{C}\mathrm{\ell}$HMC sampling methodology in parallel, the acceptance rate drops to $94\%$ (which is still very high). However, the mass distribution is much better than the serial case. This is supported by results in Figure 1(d) compared to Figure 1(c).

Figure 2 shows the CPU time and speedup results of the clustering sampling algorithms with Gaussian and HMC proposal mechanisms. As suggested by the analysis in Section 5, the CPU time becomes flat once the number of processors reaches $7$, the number of components in the mixture.

The parallel efficiency results are shown in Figure 3.

The numerical results shown here suggest that, running the clustering filters in parallel not only results in computational saving, but also can potentially increase the sampling accuracy. While the discussed parallelization of the sampler reduces the computational time, the numerical results suggest that more parallelization effort is needed in order to achieve higher efficiency. In the current settings, the chains are assigned to processes in a round-robin fashion. The performance of the sampler can be greatly enhanced if the a smart scheduler is used such that the parallel chains are assigned to processors based on the sample size per chain.

6.2 Two-dimensional example

Here we test the the paralle $\mathcal{C}\mathrm{\ell}$HMC sampling algorithm as a tool for statistical medical image retrieval. We employ a non-linear Gaussian convolution filter as a forward operator $\mathcal{H}$. The convolution filter is applied to a two-dimensional image resulting in a blurred image, then Gaussian noise is added to collect synthetic measured/observed image. Here, the vector $\mathbf{x}$ represents intensities of the image pixels arranged in a column vector. The observation nosie level is set to be $9\%$ of the average intensity of the original image. This formulation clearly results in a nonlinear inverse problem that can be challanging for traditional approaches such as Tikhonov regularization. For , the Jacobian of the convoluted image with respect to the intensities of the image pixels is found to be a Toeplitz matrix.

The goal here is to retrieve the original image given the noisy measurement, and a sample drawn from the probability distribution from which the blurred image is drawn. This is is a relevant problem description in many cases, where several low resolution or blurrd images are taken along with the collected measurement.

Figure 4 shows the original (true) image, the blurred image constructed by the convolution filter, and the noisy image, i.e. the blurred image with additive noise.

To create a synthetic non-Gaussian sample, we sampled $50$ images from a Gaussian distribution centered around the blurred image with variances equal to $8\%$ of the average entensity of the original image. To create a synthetic prior sample, we have selected $\textsc{n}_{\textnormal{ens}}=30$ uniformly random distributed images from the generated Gaussin sample. This procedure is guaranteed to result in a non-Gaussian prior, and is powerful to test the GMM prior assumption.

The $\mathcal{C}\mathrm{\ell}$HMC sampling algorithm with both Gaussian and Hamiltonian kernels performed similarily. The acceptance rates however were quite different. Specifically, the rejection rate in the case of Gaussian proposal was $56\%$ resulting in wasted computations, while HMC rejection rate was $7\%$.

The initial state of each chain is chosen as the mean of the corresponding component in the prior mixture, and the parameters of the Hamiltonian system are tuned empirically to give acceptable acceptance rate. $30$ sample members are collected from the posterior distribution.

The mean and median of posterior samples collectecd using the parallel $\mathcal{C}\mathrm{\ell}$HMC sampling algorithm are shown in Figure 5 (panels 5(b), and 5(c)). For the sake of comparison to one of the most popular and widely used approaches, we show the results obtained using Tikhonov regularization approach with regularization parameter optimally tuned following an L-curve approach. The Tikhonov-regularized solution is shown in panel 5(a) of Figure 5.

By comparing results in Figure 5, to the blurred (prior) image 4(b) and the noisy image 4(c), one can see that the posterior samples produce statistics those are closer to the original image. Moreover, the retrieved results using the non-Gaussian $\mathcal{C}\mathrm{\ell}$HMC (serial or parallel) algorithm is much better than the solution obtined by the traditional Tikhonov regularization approach. The relative error of the mean obtained using the parallel $\mathcal{C}\mathrm{\ell}$HMC sampling algorithm is $0.01663$, while the relative error of the regularized solution is $0.021447$. The relative error is defined as $\frac{\lVert\mathbf{x}-\mathbf{x}^{\rm true}\rVert}{\lVert\mathbf{x}^{\rm true}\rVert}\,,$ where $\mathbf{x}$ is the retrieved image, and $\mathbf{x}^{\rm true}$ is the true image. The proposed parallel algorithms achieve an improvement of $29\%$ over the optimally tuned Tikhonov-based solution. These results show the capability and accuracy of the $\mathcal{C}\mathrm{\ell}$HMC sampling algorithm in non-Gaussian settings.

While we have shown the blurred image in Figure 4, it is of utmost importance to highlight the fact that the formulation and the sampler is unaware of this blurred image. This is mainly due to the fact that we have uniformly sampled the Gaussian sample centered around this blurred image. To further enhance the retrieved image, one can follow the inverse problem solution by deblurring filter which is out of the scope of this work.

7 Conclusions and Future Work

In this work, we have developed parallel cluster sampling algorithms for solving Bayesian inverse problems. Specifically, we have proposed parallel sampling algorithms based on the cluster sampling filter ($\mathcal{C}\mathrm{\ell}$HMC) developed mainly for non-Gaussian data assimilation. The proposed algorithms can be efficiently used for solving various large-scale problems including medical image retrieval from noisy observations.

We have introduced a detailed complexity analysis of the proposed parallel clustering sampling ($\mathcal{C}\mathrm{\ell}$HMC) algorithms with mixture model represenation of the prior information. Generally speaking, aside from parallelization, the parallel versions of the algorithm result in higher acceptance rates. Specifically, the parallel $\mathcal{C}\mathrm{\ell}$HMC increases the acceptance rate of the sampler from $44\%$ to $93\%$ with Gaussian proposal kernel, leading to massive saving of computations. The proposed sampling algorithms achieve an improvement of $29\%$ over the optimally-tuned Tikhonov-based solution for image retrieval. The algorithm can run significantly faster than the serial sampler in ideal settings. However, the algorithm can be slower than the serial sampler if too many outliers exist where some chains are assigned much smaller ensemble size than the others. The $\mathcal{C}\mathrm{\ell}$HMC sampling algorithm, in addition to desired parallelization features, has proved powerful in the context of Bayesian image retrieval.

In future work, we will investigate the possibility of parallelizing other components of the sampling algorithm such as the likelihood function and the proposal mechanisms. For example a parallel version of EM can be considered to construct the GMM approximation to the prior distribution. In the case of HMC, the symplectic integrator can be parallelized. Also, matrix-vector products can be parallelized. Methods for parallelizing a single chain can be considered for a second level of parallelization. While our implementation is not a direct parallelization of the MCMC algorithm, it still provides acceptable sampling with better performance.