# Pitman-Yor mixture based reconstruction for emission tomography

Eric Barat, Thomas Dautremer, Thierry Montagu

► To cite this version:

Eric Barat, Thomas Dautremer, Thierry Montagu. Pitman-Yor mixture based reconstruction for emission tomography. 25th International Conference on Machine Learning, Jul 2008, Helsinki, Finland. cea-05461122

HAL Id: cea-05461122

https://cea.hal.science/cea-05461122v1

Submitted on 17 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Éric Barat

Thomas Dautremer

Thierry Montagu

ERIC.BARAT@CEA.FR

THOMAS.DAUTREMER@CEA.FR

THIERRY.MONTAGU@CEA.FR

CEA, LIST, Stochastic Processes and Spectra Laboratory, F-91191 Gif-sur-Yvette, France.

Abstract

We introduce an emission tomography reconstruction algorithm following a nonparametric Bayesian approach. In contrast with the well known and used Expectation Maximization (EM), the proposed technique does not rely on any bulky (and rather *ad hoc*) space discretization. Namely, we cast the spatial emission density reconstruction problem in the more general Bayesian nonparametric point inverse problem framework:

$$F(\cdot) = \int_{\mathcal{X}} \mathcal{P}(\cdot|\mathbf{x}) G(dx) \\ Y_i|F \stackrel{\text{i.i.d}}{\sim} F, \text{ for } i=1, \dots, n \tag{1}$$

where $\mathcal{X}$ represents the object space, $\mathbf{Y} = \{Y_1, \dots, Y_n\}$ is the observed $F$-distributed dataset, $\mathcal{P}(\cdot)$ is a known probability distribution and $G(\cdot)$ is a random probability measure whose posterior distribution has to be estimated from observation $\mathbf{Y}$. In the context of emission tomography, $Y_i$ corresponds to the index of the line of response (LOR) of the $i^{th}$ observed event. For $l = 1, \dots, L$, $\mathcal{P}(l|\mathbf{x})$ is the projection distribution i.e. the probability of recording an event in the $l^{th}$ detection unit given an emission located in $\mathbf{x}$. Finally, $G(\cdot)$ is the spatial emission distribution (the so-called activity concentration).

1. Methodology

The classical MAP-EM approach for tomographic reconstruction requires in its basis the discretization of the object space. The framework of Bayesian nonparametrics allows to get free of any *ad hoc* space

discretization while taking into account the Poisson nature of observations in the inverse problem and provides estimates of the reconstruction error distribution.

1.1. Tomographic reconstruction problem

A characteristic of tomographic reconstruction problems is that only projections $\mathbf{Y}$ are observed, then the emission locations $\mathbf{X}$ are introduced as *latent variables*. We propose thus the following hierarchical model.

1.2. Generative hierarchical model for emission tomographic data.

$$\begin{align} (Y_i|X_i) &\stackrel{\text{ind}}{\sim} \mathcal{P}(Y_i|X_i) \nonumber \\ (X_i|Z_i) &\stackrel{\text{ind}}{\sim} \mathcal{N}(X_i|Z_i) \tag{2} \\ (Z_i|H) &\stackrel{\text{i.i.d.}}{\sim} H \nonumber \\ H &\sim \mathcal{PY}(\alpha, \theta, G_0) \nonumber \end{align}$$

where $H$ is a Pitman-Yor process (or two-parameter Poisson-Dirichlet process) and $H$ has the same distribution as $\tilde{H}$, see (Pitman, 1996) Definition 19,

$$\tilde{H}(\cdot) = \sum_{k=1}^{\infty} w_k \delta_{\tilde{Z}_k}(\cdot)$$

with

* $\tilde{Z} = \tilde{Z}_1, \tilde{Z}_2, \dots$ are i.i.d. samples generated from $G_0$ which is taken as a Normal-Inverse Wishart model ($\mathcal{NIW}_{\rho, n_0, \mu_0, \Sigma_0}$) s.t. $\tilde{Z}_k = (\mathbf{m}_k, \mathbf{\Sigma}_k)$ with

$$\begin{align*} \mathbf{m}_k | \mathbf{\Sigma}_k &\stackrel{\text{ind.}}{\sim} \mathcal{N}(\mu_0, \mathbf{\Sigma}_k / \rho) \\ \mathbf{\Sigma}_k^{-1} &\stackrel{\text{i.i.d.}}{\sim} W(n_0, (n_0 \Sigma_0)^{-1}). \end{align*}$$

* $\mathbf{w} = w_1, w_2, \dots$ follow a Residual Allocation Model (RAM) – a.k.a stick-breaking – construction:

<footnote>Preliminary work. Under review by the International Conference on Machine Learning (ICML). Do not distribute.</footnote>
$V_1, V_2, \dots$ are i.i.d., s.t. $V_k = \text{Beta}(1 - \alpha, \theta + k\alpha)$ and $\mathbf{w} = w_1, w_2, \dots$ s.t. $w_1 = V_1$ and for all $k \ge 2$, $w_k = V_k \cdot \prod_{i=1}^{k-1} (1 - V_i)$.

Emission locations $\mathbf{X} = \{X_1, X_2, \dots, X_n\}$ are thus distributed following a Pitman-Yor mixture (PYM) (Pitman & Yor, 1997), (Ishwaran & James, 2001), with density

$$f_G(\cdot) = \sum_{k=1}^{\infty} w_k \mathcal{N}(\cdot | \tilde{Z}_k = (\mathbf{m}_k, \mathbf{\Sigma}_k)) \quad (3)$$

In (Pitman, 1996) Definition 19, the random distribution $H$ appears as the law governing an exchangeable sequence $\mathbf{Z}$ in a generalized Blackwell-MacQueen urn scheme while $\tilde{H}$ corresponds to a size-biased permutation of $H$. Both representation may be used equally for prior specification of the model. But for posterior inference, the use of stick-breaking representation implies working with an explicit labelling of clusters while atoms in the Pólya urn representation correspond to the equivalence classes of the random equivalence relation: $i \sim j$ iff $Z_i = Z_j$ which makes irrelevant any explicit cluster labelling. In the next section, inference will be based upon the urn characterization even if we make use of the stick-breaking construction for the data independent part of the posterior and follow a conditional approach for the sampler.

## 2. MCMC algorithm

Based on the hierarchical model (2), a Gibbs sampler is constructed in order to generate draws $f_G(\cdot) | \mathbf{Y}$.

Let first introduce $\mathbf{C} = C_1, C_2, \dots, C_n$ the latent classification variables of the emission locations to Pitman-Yor components, we have $Z_i = \tilde{Z}_{C_i}$ for all $i < n$. Note that at each iteration non empty clusters are relabelled according to their order of appearance in the sampling.

### 2.1. Sampling of Pitman-Yor mixtures using auxiliary variables

We propose a slice sampler inspired from (Walker, 2007) and (Kalli et al., 2008) slightly modified by the introduction of a threshold on the weights defining the domain of auxiliary variables. We summarize the main steps.

We introduce uniform auxiliary variables $\mathbf{u} = u_1, u_2, \dots, u_n$. The joint density for any $(X_i, u_i)$ is

defined as, for some $\zeta \in ]0, 1]$

$$f_{\zeta}(X_i, u_i | \mathbf{w}, \tilde{\mathbf{Z}}) = \sum_{k=1}^{\infty} w_k \mathcal{N}(X_i | \tilde{Z}_k) \mathcal{U}(u_i | 0, \min(\zeta, w_k)) \quad (4)$$

Obviously, $\int f_{\zeta}(X_i, u_i | \mathbf{w}, \tilde{\mathbf{Z}}) du_i = f_G(X_i)$ and

$$f_{\zeta}(X_i, u_i | \mathbf{w}, \tilde{\mathbf{Z}}) = \mathbf{1}(\zeta > u_i) \zeta^{-1} \sum_{w_k > \zeta} w_k \mathcal{N}(X_i | \tilde{Z}_k) \\ + \sum_{w_k \le \zeta} \mathbf{1}(w_k > u_i) \mathcal{N}(X_i | \tilde{Z}_k)$$

where both sums are finite since $\#\{j : w_j > \varepsilon\} < \infty$ for any $\varepsilon > 0$.

Let denote by $K_n$ the (latent) number of distinct clusters after $n$ observations and $n_k = \#\{i : C_i = k\}$ for $k \le K_n$. As in (Kalli et al., 2008), we jointly sample $\mathbf{u}, \mathbf{w} | \mathbf{C}$ by first sampling $w_1, w_2, \dots, w_{K_n} | \mathbf{C}$, then sampling $\mathbf{u} | w_1, w_2, \dots, w_{K_n}, \mathbf{C}$, and finally sampling $w_{K_n+1}, w_{K_n+2}, \dots | \mathbf{u}$.

In this context, following (Pitman, 1996), Definition 19 and Corollary 20, we can express the posterior of the random Pitman-Yor distribution as the sum of a finite sum of atoms located at the distinct values whose weights are given by a Dirichlet distribution, and an unconditional Pitman-Yor process. The use of $\mathbf{u}$ allows to sample a finite number $k^*$ of weights and locations for the independent Pitman-Yor process.

* Sample $w_k$ for $k \le K_n$.

$$(w_1, w_2, \dots, w_{K_n}, r_{K_n} | \mathbf{C}) \sim \text{Dirichlet}(n_1 - \alpha, n_2 - \alpha, \dots, n_{K_n} - \alpha, \theta + K_n \alpha)$$

* Sample $u_i | w_1, w_2, \dots, w_{K_n}, \mathbf{C}$.

$$u_i | \mathbf{C}^{\text{ind.}} \sim \mathcal{U}(u_i | 0, \min(w_{C_i}, \zeta)) \\ u^* = \min(\mathbf{u})$$

* Sample $w_k$ for $k > K_n$. While $r_{k-1} > u^*$,

$$V_k \sim \text{Beta}(1 - \alpha, \theta + k \alpha) \\ w_k = V_k r_{k-1} \\ r_k = r_{k-1} (1 - V_k)$$

* Set $k^* = \min(\{k : r_k < u^*\})$.

Clearly, $w_k < u^*$ for all $k > k^*$, this is why we only have to sample a finite set of $w_k$. Note that we also drop in computations all $w_k$ s.t. $w_k < u^*$ for $k \le K_n$.
## 2.2. Sampling of classification variables

As underlined in (Walker, 2007), sampling of classification variables $C|\tilde{Z}, w, X, u$ requires the computation of a normalizing constant which becomes feasible using auxiliary variables since the choice of $C_i$ is from a finite set.

$$ (C_i|\tilde{Z}, w, X, u) \stackrel{\text{ind}}{\sim} \sum_{j=1}^{k^*} w_{j,i} \delta_j(\cdot) \quad (5) $$

with

$$ w_{j,i} \propto 1(w_j > u_i) \max(w_j, \zeta) f_N(X_i|\tilde{Z}_j) $$

and $\sum_{j=1}^{k^*} w_{j,i} = 1$

## 2.3. Sampling mixture parameters

Sampling from the conditional distribution of $(\tilde{Z}|C, X, u)$ is straightforward.

* for all $j \le K_n$, generate $(\tilde{Z}_j|C, X, u)$ from the Normal-Inverse Wishart model with updated parameters (conjugacy property)

$$ \mathbf{m}_j|\Sigma_j, C, X \stackrel{\text{ind}}{\sim} \mathcal{N}\left(\mu_j, \frac{\Sigma_j}{\rho_j}\right) \quad (6) $$

with

$$ \rho_j = \rho + n_j $$

$$ \mu_j = \frac{\rho \mu_0 + n_j \hat{\mu}_j}{\rho + n_j} $$

$$ \hat{\mu}_j = \frac{1}{n_j} \sum_{i:C_i=j} X_i $$

and

$$ \Sigma_j^{-1}|C, X \stackrel{\text{ind}}{\sim} W(\eta_j, (\eta_j \mathbf{S}_j)^{-1}) \quad (7) $$

with

$$ \eta_j = n_0 + n_j $$

$$ \mathbf{S}_j = \frac{1}{\eta_j} \left[ n_0 \mathbf{\Sigma}_0 + n_j \hat{\mathbf{\Sigma}}_j + \frac{(\mu_0 - \hat{\mu}_j)(\mu_0 - \hat{\mu}_j)^\top}{\frac{1}{\rho} + \frac{1}{n_j}} \right] $$

$$ \hat{\mathbf{\Sigma}}_j = \frac{1}{n_j} \sum_{i:C_i=j} (X_i - \hat{\mu}_j)(X_i - \hat{\mu}_j)^\top $$

* draw $\tilde{Z}_j$ i.i.d. $\mathcal{NIW}_{\rho,n_0,\mu_0,\Sigma_0}$, for $K_n < j \le k^*$.

## 2.4. MH sampling of emission locations

Sampling $(X|Y, w, \tilde{Z})$ requires a MH step (within Gibbs) and elicitation of a candidate distribution $Q$. We propose an independent MH scheme and to take for $Q$ a tube whose axis is the LOR of the observation. The proposal in the LOR direction is thus a re-weighted mixture of one-dimensional normal distributions. See (Barat et al., 2007) for details.

## 2.5. Remarks

The blocked Gibbs sampler structure allows easy implementation of the MCMC algorithm on a parallel computer since, at one iteration, the random distributions are retained for the whole dataset. In contrast with usual conditional blocked Gibbs sampler, we operate in the space of *equivalence classes over non empty cluster labels* which are thus *exchangeable* (Pitman, 1996), the stick-breaking prior is only used for independent clusters at a given iteration of the Gibbs sampler. As pointed out in (Porteous et al., 2006), this encourages good mixing over clusters in comparison to the usual stick-breaking sampling. In the same time, the sampler keeps the random distribution functions which is not accessible with collapsed Gibbs approach.

The introduction of threshold $\zeta$ in the slice sampler improves the mixing of the Markov chain whenever we choose $\zeta$ sufficiently small ($\zeta = 1$ in (Walker, 2007)). A typical value is the prior mean size of a cluster containing a single data allocated to ($\zeta \approx \frac{1-\alpha}{n+\theta}$). In this case, most weights of non empty clusters are used in (5) while just the kernel values take place in the sampling of classification variables in (Walker, 2007).

As an alternative to this sampler, one can resort to the approximation of taking $k^* = K_n + K_0$ with $K_0$ fixed. We have in this case a truncated Pitman-Yor process for the independent part of the posterior. This scheme is suggested in (Ishwaran & James, 2001), and it eliminates the need for auxiliary variables.

# 3. Application

## 3.1. Test of the method

A sensitive question is the ability of the nonparametric prior to capture the structure of spatial emission distribution. Since the number of observations and the data structure complexity can influence deeply the behavior of the nonparametric approach, a relevant test has to be performed on realistic data. We generated $3 \times 10^5$ decays in one slice of a voxelized 2D brain phantom (modified Zubal, see Fig. 1, image (a)). For sake of clarity, the projection distribution $P(\cdot|x)$ was assumed
purely geometrical, see (Barat et al., 2007). Attenuation and random coincidences were not simulated. We compare our results to MAP-EM algorithms, based on a Gibbs field with exponentially decreasing quadratic potentials ($\sigma = 1.7mm$), for a 128 $\times$ 128 image.

## 3.2. Algorithm parameterization

Parameters for the Pitman-Yor mixture were chosen as follows : $\theta = 30$, $\alpha = .1$, $\rho = 10^{-2}$, $\Sigma_0 = 3 \times I_2$, $n_0 = 3$, $\mu_0 = (0,0)^T$. The regularization of the ill-posed problem is entirely determined by this choice of parameters which are convenient to control. Actually, $\alpha$ and $\theta$ tune the number of distinct clusters, and the other parameters regulate the behavior of the *NTW* model. In particular, $\Sigma_0$ tunes the prior spread-out of clusters.

We computed 5000 burn-in iterations of the Gibbs sampler followed by 10000 iterations which are retained for Monte Carlo estimates.

## 3.3. Results

On Fig. 1, image (a) represents the 2D brain phantom, image (b) corresponds to the MAP-EM reconstruction after 500 iterations. The MCMC sampler allows straightforward estimation of any posterior functionals. We then computed and plotted on Fig. 1 the conditional mean *i.e.* the predictive density of activity concentration (image (c)) and the conditional standard deviations (image (d)).

The number $K_n$ of PYM clusters affected with at least one data is of the order of 400 during iterations. The nonparametric behavior of the method is characterized by the fact that this number can increase with the amount of observations, and then automatically improves the activity distribution resolution. Note that the discretization of the spatial distribution of PYM steps in only for visualization purposes and has been chosen a *posteriori* to 256 $\times$ 256 corresponding to the phantom image size.

As a reconstruction quality indicator, we give in table 1 the Kullback-Leibler divergence (KL) and Hellinger distance (DH) w.r.t. the brain phantom for the PYM based reconstruction approach and MAP-EM. Datasets of various sizes have been considered. These results highlight the good performances of the Bayesian nonparametric reconstruction scheme in terms of both signal to noise ratio and spatial resolution. Indeed, the conditional expectation reconstruc- tion exhibits smooth regions (white matter) and sharp edges (grey matter, skull). This is due to the ability of the *NTW* model to estimate a covariance matrix for

<figure><img src="image_2.png" /><figcaption>Figure 1. (a) 2D brain phantom, (b) MAP EM (500 iter.), (c) PYM Conditional expectation (10000 draws), (d) PYM Conditional standard deviation (10000 draws), scale in (d) is a quarter of scale in (a), (b), (c).</figcaption></figure>

each component of the mixture. By this versatility, the Pitman-Yor mixture of normal components turns out to be a convenient regularization model for the brain spatial activity distribution.

Finally, on Figs. 3 and 4, we present results of 3D tomographic reconstructions based on $10^7$ events in coincidence from a 3D phantom (Fig.2).

## References

* Barat, E., Comtat, C., Dautremer, T., Montagu, T., & Trebossen, R. (2007). A Bayesian nonparametric approach for PET reconstruction. *IEEE Medical Imaging Conference Record, NSS/MIC '07*. (pp. 4155–4162).

* Ishwaran, H., & James, L. F. (2001). Gibbs sampling methods for stick-breaking priors. *J. Am. Stat. Assoc.*, 96, 161–173.

* Kalli, M., Griffin, J. E., & Walker, S. G. (2008). *Slice sampling the mixture models* (Technical Report UKC/IMS/08/024). IMSAS, University of Kent.

* Pitman, J. (1996). Some developments of the Blackwell-MacQueen urn scheme. In et T. F. al. (Ed.), *Statistics, probability and game theory; papers*
<table><caption>Table 1. Kullback-Leibler divergence (%) and Hellinger distance (%) for MAP-EM and PYM approaches.</caption><thead><tr><th></th><th>n</th><th>300k</th><th>100k</th><th>30k</th></tr></thead><tbody><tr><td rowspan="3">DH</td><td>PYM</td><td>26.23</td><td>28.94</td><td>32.4</td></tr><tr><td>MAP-EM, σ = 1.7</td><td>28.22</td><td>31.0</td><td>35.9</td></tr><tr><td>MAP-EM, σ = 2.5</td><td>30.69</td><td>32.01</td><td>34.34</td></tr><tr><td rowspan="3">KL</td><td>PYM</td><td>19.41</td><td>23.64</td><td>31.42</td></tr><tr><td>MAP-EM, σ = 1.7</td><td>20.99</td><td>27.45</td><td>42.96</td></tr><tr><td>MAP-EM, σ = 2.5</td><td>23.50</td><td>26.39</td><td>32.90</td></tr></tbody></table>

in honor of david blackwell, vol. 30 of Lecture Notes-Monograph Series, 245–267. Institute of Mathematical Statistics.

Pitman, J., & Yor, M. (1997). The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator. *Ann. Probab.*, 25, 855–900.

Porteous, I. R., Ihler, A., Smyth, P., & Welling, M. (2006). Gibbs sampling for (coupled) infinite mixture models in the stick breaking representation. *UAI*. AUAI Press.

Walker, S. G. (2007). Sampling the Dirichlet mixture model with slices. *Comm. Statist.*, 36, 45–54.

<figure><img src="image_3.png" /><figcaption>Figure 2. 3D brain phantom</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. PYM Conditional expectation (10000 draws)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. PYM Conditional standard dev. (10000 draws)</figcaption></figure>