Faster PET reconstruction with non-smooth priors by randomization and preconditioning

1.1.1. Subset acceleration with randomization

1.1.2. Preconditioning

1.1.3. Non-uniform sampling

1.1.4. Uncompressed data

Given the measured data vector $b \in {{\mathbb N}}^M$ and the projection model $\newcommand{\mP}{\mathbf P} \mP$, the PET reconstruction problem can be formulated as the solution to the optimization problem

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mP}{\mathbf P} \displaystyle \min_{u \geqslant 0} \Bigl\{D(\mP u) + \alpha R(u) \Bigr\} \label{EQU:PROB:OPT} \nonumber \end{align} \tag{ 1 }$$

where the data fidelity $\newcommand{\mP}{\mathbf P} D(\mP u)$ measures the match of the estimated image u with the data and the prior $\alpha R(u)$ penalizes features that are not desirable in the solution. In other words the prior can be used to avoid solutions which would fit the noisy data too closely. The data fidelity D is (up to constants independent of u) the negative log-likelihood of the multi-variate Poisson distribution

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle D(y) = \sum_{i=1}^M y_i + r_i - b_i + b_i \log\left(\frac{b_i}{y_i + r_i}\right) \, , \nonumber \end{align*}$$

with expected value being the sum of the projected image y  and the estimated background activity r. The latter is needed in order to model non-linear effects such as scatter and randoms. The data fidelity D measures the distance of the estimated data $\newcommand{\mP}{\mathbf P} \mP u + r$ to the measured data b in the sense that $\newcommand{\mP}{\mathbf P} D(\mP u) \geqslant 0$ and $\newcommand{\mP}{\mathbf P} D(\mP u) = 0$ if and only if $\newcommand{\mP}{\mathbf P} \mP u + r = b$. The operator $\newcommand{\mP}{\mathbf P} \mP$ performs the projection and includes geometric factors, attenuation and normalization.

While the main motivation is the efficient solution of non-smooth optimization problems, we first compare the method to ordered subsets expectation maximization (OSEM) (Hudson and Larkin 1994) for unregularized reconstruction. The 'ordered subsets' idea has subsequently been used for many algorithms related to non-smooth optimisation, see e.g. McGaffin and Fessler (2015). We would like to show in the next example (1) that the 'ordered subsets' idea is generally non-convergent and thus may be unstable and (2) that the proposed algorithm is as fast as OSEM—despite its proven convergence.

If there is no prior, i.e. $\alpha R = 0$, the most common algorithm to solve the optimization problem (1) is the maximum likelihood expectation maximization algorithm (MLEM) (Shepp and Vardi 1982) defined by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mP}{\mathbf P} \displaystyle u^{k+1} = \frac{u^k}{\mP^T 1} \mP^T\left(\frac{b}{\mP u^k + r} \right) \label{EQU:MLEM} \, , \nonumber \end{align} \tag{ 2 }$$

where all operations have to be understood element-wise. The computational bottleneck in the MLEM algorithm is the evaluation of the operator $\newcommand{\mP}{\mathbf P} \mP$ and its transpose $\newcommand{\mP}{\mathbf P} \mP^T$ in each iteration.

To overcome this hurdle, it has been proposed to change the update and evaluate the operator and its adjoint only on one out of m subsets of the data in each iteration. At every iteration k we choose $i = {\rm mod}(k, m)$ and change update formula (2) to

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mP}{\mathbf P} \displaystyle u^{k+1} = \frac{u^k}{\mP^T_i 1} \mP^T_i\left(\frac{b_i}{\mP_i u^k + r_i} \right) \label{EQU:OSEM} \, . \nonumber \end{align} \tag{ 3 }$$

This algorithm became known as OSEM. Here $\newcommand{\mP}{\mathbf P} \mP_i$ is the restriction of $\newcommand{\mP}{\mathbf P} \mP$ onto the ith subset, i.e. $\newcommand{\mP}{\mathbf P} \mP = (\mP_1^T, \ldots, \mP_m^T){}^T$. While this change of the update equation reduces the computational burden by 1/m, it is in general not guaranteed to converge to a solution of (1), illustrated in figures 1(a) and (b). A convergent version of OSEM, called complete-data OSEM (COSEM), has been developed (Hsiao et al 2002). While it comes with mathematical convergence guarantees, it is much slower than OSEM (see figure 1(c)) and therefore never became popular for the reconstruction of clinical PET data.

Zoom In Zoom Out Reset image size

MLEM has been extended to include smooth (Green 1990) and certain non-smooth (Sawatzky et al 2013) prior information, however, conceptually both algorithms intrinsically struggle with the ordered subset acceleration. Also other algorithms have been 'accelerated' based on the ordered subset idea, e.g. McGaffin and Fessler (2015) and Ross Schmidtlein et al (2017), but are similarly intrinsically unstable due to their non-convergence. See Cheng et al (2013) for a numerical comparison and Ahn et al (2015) and Teoh et al (2015) for a validation on clinical PET data. For differentiable priors, a surrogate based technique allows for stable subset acceleration (De Pierro and Yamagishi 2001, Ahn and Fessler 2003, Ahn et al 2006). In this work we propose the subset-accelerated algorithm SPDHG that is provably convergent and thus stable and robust, see figures 1(a) and (b). SPDHG is flexible enough to be applicable to a large variety of convex and non-smooth priors and is as efficient as OSEM if no explicit prior is being used, see figures 1(c) and (d).

As outlined above, PET reconstruction can be formulated in terms of the optimization problem (1). Computationally, it is convenient to rewrite (and solve) the optimization problem (1) in terms of subsets. We denote by M the number of projection bins. Let {S_i} be a partition of [M], in the sense that $\cup_{i=1}^m S_i = [M]$, where we used the notation $[M] := \{1, \ldots, M\}$. It is not necessary to assume that $\newcommand{\e}{{\rm e}} S_i \cap S_j = \emptyset$ for $i \neq j$. For notational simplicity we will restrict ourselves to the this case. We define

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle D_i(y) := \sum_{j \in S_i} \varphi_j(y_j) \label{EQU:DATATERM} \nonumber \end{align} \tag{ 4 }$$

with the distance function for every data point given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varphi(y) := \left\{\begin{array}{@{}ll@{}}{y + r - b \log(y + r) - b + b \log b }&{\rm if}~ y + r \geqslant 0 \nonumber \\ {\infty}&{{\rm else}}\end{array}\right. \, ,\label{EQU:LOCALDISTANCE} \nonumber \end{align} \tag{ 5 }$$

where we omitted the index j  at $\varphi, y, r$ and b for readability. Algorithms from convex optimization require the problem to be defined over an entire vector space which we satisfy by extending $\varphi$ to $\infty$ for non-positive estimated data y   +  r. The data and the background are photon counts and therefore have a natural non-negativity constraint. To allow for the concise notation in (5), we define $0 \log 0 := 0$ and $- \log 0 := \infty$.

We model the non-negativity constraint for the image u with the indicator function $\newcommand{\im}{{\rm Im}}\imath_+$, which is defined as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\im}{{\rm Im}} \displaystyle \imath_+(u) = \left\{\begin{array}{@{}ll@{}}{0}&{{\rm if}~ u \geqslant 0} \nonumber \\{\infty}&{{\rm else}}\end{array}\right. \, . \label{EQU:NONNEGATIVITY} \nonumber \end{align} \tag{ 6 }$$

Thus, this results in the unconstrained optimization problem

Problem 1 (PET reconstruction with subsets). 

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\im}{{\rm Im}} \newcommand{\mP}{\mathbf P} \newcommand{\bR}{{\mathbb R}} \displaystyle u^\sharp \in \arg \min_{u \in \bR^N} \left\{\sum_{i=1}^m D_i(\mP_i u) + \alpha R(u) + \imath_+(u) \right\} \label{EQU:PROB:OPT:SUBSETS} \, . \nonumber \end{align} \tag{ 7 }$$

We would like to stress that solving problem (7) is equivalent to solving the original problem (1) for any choice of subsets. In fact, the subset selection becomes a reconstruction parameter that may be varied to speed up the reconstruction procedure.

Often, our prior assumptions involve linear operators, too. One of the most prominent examples of this is the total variation (Rudin et al 1992)

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\TV}{{\rm TV}} \displaystyle R(u) = \TV(u) = \|\nabla u\|_{2,1} = \sum_i \|\nabla u_i\|_2 = \sum_i \left(\sum_{j=1}^3 (\partial_j u_i)^2\right)^{1/2} \, , \nonumber \end{align*}$$

where we take the 2-norm locally, i.e. at every voxel i we take the 2-norm of the spatial gradient, and the 1-norm globally, i.e. we sum over all voxels. Forward difference discretization of the gradient operator $\nabla$ is used as in Chambolle and Pock (2011). Similarly, we use the directional total variation $\newcommand{\dTV}{{\rm dTV}} \newcommand{\mD}{\mathbf D} R(u) = \dTV(u) = \|\mD \nabla u\|_{2,1}$ to incorporate a-priori knowledge about the solution given by an anatomical prior image, see Ehrhardt and Betcke (2016), Ehrhardt et al (2016), Schramm et al (2017) and Bungert et al (2018) for details.

Solving problem (7) is challenging, even when the involved variables are small and matrix-vector products are easy to compute. The difficulty stems from its non-smoothness. The data term D_i is not finite everywhere and while it is differentiable on its effective domain $\newcommand{\dom}{{\rm dom}} \dom(D_i) := \{y \mid D_i(y) < \infty\}$, the gradient is not globally Lipschitz continuous. In addition, further non-smoothness comes from the constraint $\newcommand{\im}{{\rm Im}}\imath_+$ and the prior R may be non-smooth as well. All of this being said, in PET reconstruction, the variable sizes are actually very large and matrix-vector products expensive to compute.

To apply optimization algorithms to solve (7), we reformulate it as a generic optimization problem of the form

Problem 2 (Generic optimization problem). 

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mA}{\mathbf A} \displaystyle x^\sharp \in \arg \min_{x \in X} \left\{\Psi(x) := \sum_{i=1}^n f_i(\mA_i x) + g(x) \right\} \, . \label{EQU:PROB:OPT:GENERIC} \nonumber \end{align} \tag{ 8 }$$

For instance, for unregularized reconstructions, i.e. $\alpha R = 0$, we may make the association

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\im}{{\rm Im}} \newcommand{\mP}{\mathbf P} \newcommand{\mA}{\mathbf A} \displaystyle n = m, \quad g = \imath_+, \quad f_i = D_i, \quad \mA_i = \mP_i \nonumber \end{align*}$$

and reconstructions regularized by the total variation, i.e. $R(u) = \|\nabla u\|_{2,1}$, can be achieved by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\im}{{\rm Im}} \newcommand{\mP}{\mathbf P} \newcommand{\mA}{\mathbf A} \displaystyle n &= m + 1, \quad f_i = D_i, i \in [m], \,\,\,\quad f_n = \alpha \|\cdot\|_{2,1}\nonumber \\ g &= \imath_+, \quad \quad \,\,\mA_i = \mP_i, i \in [m], \quad\mA_n = \nabla \, . \label{EQU:SETTING:TV} \nonumber \end{align} \tag{ 9 }$$

Instead of solving problem (8) directly, it is more efficient to reformulate the minimization problem as a saddle point problem making use of the convex conjugate of a functional, see e.g. Bauschke and Combettes (2011).

Definition 1 (Convex conjugate). Let $\newcommand{\bR}{{\mathbb R}} \newcommand{\bRI}{\bR_\infty} f : Y \to \bRI := \bR \cup \{\infty\}$ be a functional with extended real values. Then we define the convex conjugate of f  as $\newcommand{\bR}{{\mathbb R}} \newcommand{\bRI}{\bR_\infty} f^\ast : Y \to \bRI$ with

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle f^\ast(y) = \sup_x \left\{\langle y, x\rangle - f(x) \right\} \, . \nonumber \end{align*}$$

For convex, proper and lower semi-continuous (lsc) functionals f  we have that $f^{\ast\ast} = f$, see e.g. Bauschke and Combettes (2011), and thus $f(x) = \sup_y \left\{\langle x, y\rangle - f^\ast(y) \right\}$. Then, with $Y = \prod_{i=1}^n Y_i$, problem (8) is equivalent to

Problem 3 (Generic saddle point problem). 

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mA}{\mathbf A} \displaystyle \min_{x \in X} \sup_{y \in Y} \left\{\sum_{i=1}^n \langle \mA_i x, y_i \rangle - f_i^\ast (y_i) + g(x) \right\} \, . \label{EQU:SADDLE} \nonumber \end{align} \tag{ 10 }$$

We will refer to the variable x as the primal variable and to y  as the dual variable.

Example 1. The convex conjugate of the PET distance function (4) is given by $D_i^\ast(y) = \sum_{j \in S_i} \varphi^\ast_j(y_j)$ with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \varphi^\ast(y) = \left\{\begin{array}{@{}ll@{}}{-y r - b \log(1 - y)}&{{{\rm if}~ y \leqslant 1}} \nonumber \\{\infty}&{{\rm else}}\end{array}\right. \label{EQU:PET:CVX} \nonumber \end{align} \tag{ 11 }$$

where we omitted the index j  at $\varphi, y, b$ and r for readability.

The derivation of the formulas in this and the following example are omitted for brevity.

As some (or all) of the f _i and g in (8) are non-smooth, we make use of the proximal operator of these. Our definition varies slightly from the usual definition as we allow the step size parameter to be matrix-valued. For a symmetric and positive definite matrix $\newcommand{\mS}{\mathbf S} \mS$, we define the weighted norm $\newcommand{\mS}{\mathbf S} \|x\|_\mS$ as $\newcommand{\mS}{\mathbf S} \|x\|_\mS^2 := \|\mS^{-1/2} x\|^2 = \langle \mS^{-1} x, x \rangle$.

Definition 2 (Proximal operator). Let $\newcommand{\mS}{\mathbf S} \mS$ be a symmetric and positive definite matrix. Then we define the proximal operator of f  with metric (or step size) $\newcommand{\mS}{\mathbf S} \mS$ as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\mS}{\mathbf S} \newcommand{\Prox}{{\rm prox}} \displaystyle \Prox^\mS_f(x) := \arg \min_z \left\{\|z - x\|^2_\mS + f(z) \right\} \, . \nonumber \end{align*}$$

From here on, $\newcommand{\mS}{\mathbf S} \mS$ and $\newcommand{\mT}{\mathbf T} \mT$ will always be diagonal (and thus symmetric) and positive definite matrices.

Example 2. The proximity operator of the non-negativity constraint (6) is given element-wise by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\im}{{\rm Im}} \newcommand{\mT}{\mathbf T} \newcommand{\Prox}{{\rm prox}} \displaystyle \Prox^\mT_{\imath_+}(x) = \max(x, 0) \, . \nonumber \end{align*}$$

Example 3. Let $\newcommand{\Diag}{{\rm diag}} \newcommand{\mS}{\mathbf S} \mS_i = \Diag((\sigma_j)_{j \in S_i})$. The proximal operator of the convex conjugate of the PET distance (11) can be computed element-wise as $\newcommand{\Prox}{{\rm prox}} \newcommand{\mS}{\mathbf S} [\Prox^{\mS_i}_{D_i^\ast}(y)]_j = \Prox^{\sigma_j}_{\varphi_j^\ast}(y_j)$. For each element, the proximal operator is given by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Prox}{{\rm prox}} \displaystyle \Prox^\sigma_{\varphi^\ast}(y) = \frac12 \left[w + 1 - \Bigl((w - 1)^2 + 4\sigma b\Bigr)^{1/2} \right] \, , \nonumber \end{align*}$$

where we again omitted the indices j  for readability and denoted $w = y + \sigma r$.

SPDHG is guaranteed to converge for any f _i and g which are convex, proper and lsc. We now state a very general convergence result which can be derived from (Chambolle et al 2018, theorem 4.3). The actual proof is omitted here for brevity. For more details on convergence and convergence rates we refer the reader to Chambolle et al (2018).

Theorem 1 (Convergence). Assume that the sampling is proper, i.e. the probability p _i for an index $i \in [n]$ to be sampled is positive. Let the step length parameters $\newcommand{\mS}{\mathbf S} \newcommand{\mT}{\mathbf T} \mT = \min_{i \in [n]} \mT_i, \mS_i$ be chosen such that for all $i \in [n]$ the following bound on the operator norm

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mT}{\mathbf T} \newcommand{\mS}{\mathbf S} \newcommand{\mA}{\mathbf A} \displaystyle \left\|\mS^{1/2}_i \mA_i \mT^{1/2}_i\right\|^2 < p_i \label{EQU:STEPSIZE} \nonumber \end{align} \tag{ 12 }$$

holds. Then for any initialization, the iterates $(x, y)$ of SPDHG (algorithm 2) converge to a saddle point of (10) almost surely in a Bregman distance.

Algorithm 1. Primal-dual hybrid gradient (PDHG) to solve (10). Default values given in brackets.

Input: iterates $x (= 0)$, $y (= 0)$, step parameters $\newcommand{\mS}{\mathbf S} \mS = \{\mS_i\}$, $\newcommand{\mT}{\mathbf T} \mT$

1: $\newcommand{\mA}{\mathbf A} \overline{z} = z = \mA^T y \; (= 0)$

2: for $k = 1, \ldots$ do

3:   $\newcommand{\Prox}{{\rm prox}} \newcommand{\mT}{\mathbf T} x = \Prox^\mT_g\left(x - \mT \overline{z}\right)$

4:   $\newcommand{\Prox}{{\rm prox}} \newcommand{\mA}{\mathbf A} \newcommand{\mS}{\mathbf S} y_i^+ = \Prox^{\mS_i}_{f^\ast_i}\left(y_i + \mS_i \mA_i x\right)$   for $i = 1, \ldots, n$

5:   $\newcommand{\mA}{\mathbf A} \Delta z = \sum_{i=1}^n \mA_i^T \left(y_i^+ - y_i\right)$

6:  $z = z + \Delta z, \quad y = y^+$

7:   $\overline{z} = z + \Delta z$

Algorithm 2. Stochastic primal-dual hybrid gradient (SPDHG) to solve (10). Default values given in brackets.

Input: iterates $x (= 0)$, $y (= 0)$, step parameters $\newcommand{\mS}{\mathbf S} \mS = \{\mS_i\}$, $\newcommand{\mT}{\mathbf T} \mT$

1: $\newcommand{\mA}{\mathbf A} \overline{z} = z = \mA^T y \; (= 0)$

2: for $k = 1, \ldots$ do

3:   $\newcommand{\Prox}{{\rm prox}} \newcommand{\mT}{\mathbf T} x = \Prox^\mT_g\left(x - \mT \overline{z}\right)$

4:   Select $i \in [n]$ at random with probability p i

5:   $\newcommand{\Prox}{{\rm prox}} \newcommand{\mA}{\mathbf A} \newcommand{\mS}{\mathbf S} y_i^+ = \Prox^{\mS_i}_{f^\ast_i}\left(y_i + \mS_i \mA_i x\right)$

6:   $\newcommand{\mA}{\mathbf A} \Delta z = \mA_i^T \left(y_i^+ - y_i\right)$

7:   $z = z + \Delta z, \quad y_i = y_i^+$

8:   $\overline{z} = z + \frac{1}{p_i} \Delta z$

Remark 1 (Computational efficiency). Each iteration of algorithm 2 is computationally efficient as only projections and backprojections corresponding to the randomly selected subset i of the data are required. However, the algorithm maintains the whole backprojected dual variable $\newcommand{\mP}{\mathbf P} z = \mP^T y = \sum_{i=1}^m \mP_i^T y_i$ and in each iteration updates the primal variable with it.

Remark 2 (Memory requirements). The memory requirement of algorithm 2 is higher compared to OSEM or gradient descent but still reasonably low. It requires memory equivalent to two images $(z, \overline{z})$ and up to twice the binned sinogram data (y ,y ⁺) in addition to the necessary memory consumption (output image, sinogram data, background and normalization).

Remark 3 (Sampling). SPDHG allows any kind of random selection as long as the draws are independent and the probability that block i is being selected with positive probability p _i  >  0. We will investigate two choices of sampling in the numerical section of this paper. A more thorough numerical and theoretical investigation will be subject of future work.

We will now discuss two different choices of step sizes under which SPDHG is guaranteed to converge. The proof of the following theorem uses arguments from Chambolle et al (2018) and Pock and Chambolle (2011) and is omitted here for brevity.

Theorem 2 (Step size parameters). Let $\rho < 1$ and $\gamma > 0$. Then, condition (12) of theorem 1 is satisfied by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mT}{\mathbf T} \newcommand{\mS}{\mathbf S} \newcommand{\mI}{\mathbf I} \newcommand{\mA}{\mathbf A} \displaystyle \mS_i = \gamma \frac{\rho}{\|\mA_i\|} \mI \; , \quad \mT_i = \gamma^{-1} \frac{\rho p_i}{\|\mA_i\|} \mI \, . \label{EQU:STEPSIZE:SCALAR} \nonumber \end{align} \tag{ 13 }$$

Moreover, if $\newcommand{\mA}{\mathbf A} \mA_i$ has only non-negative elements, then condition (12) is also satisfied by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\mT}{\mathbf T} \newcommand{\mS}{\mathbf S} \newcommand{\mA}{\mathbf A} \newcommand{\Diag}{{\rm diag}} \displaystyle \mS_i = \gamma \Diag\left(\frac{\rho}{\mA_i 1}\right) \; , \quad \mT_i = \gamma^{-1} \Diag\left(\frac{\rho p_i}{\mA_i^T 1}\right) \, . \label{EQU:STEPSIZE:PRECOND} \nonumber \end{align} \tag{ 14 }$$

An example of preconditioned step sizes (14) is shown in figure 2.

Zoom In Zoom Out Reset image size

Remark 4. If n  =  1 and p _i  =  1, then the step sizes (13) can be identified with the scalar step sizes $\newcommand{\mA}{\mathbf A} \sigma_i = \gamma \rho / \|\mA_i\|$ and $\newcommand{\mA}{\mathbf A} \tau = \gamma^{-1} \rho / \|\mA_i\|$ which are commonly chosen for PDHG.

Remark 5. Note that the non-negativity condition holds for the PET projection operator (and any other ray tracing based operator). Moreover, the step size $\newcommand{\mT}{\mathbf T} \mT$ in (14) resembles the sensitivities used in the update of MLEM (2) and OSEM (3). In addition, a similar preconditioning is performed for the dual variable in the data space.

We validate the numerical performance of the proposed algorithm on two clinical PET data sets which we refer to as $\newcommand{\dataFdg}{\texttt{FDG}} \dataFdg$ and $\newcommand{\dataAmy}{\texttt{florbetapir}} \dataAmy$. The two separate PET brain datasets each use a distinct radiotracer: [¹⁸F]FDG for epilepsy and [¹⁸F]florbetapir for the neuroscience sub-study Insight'46 of the Medical Research Council National Survey of Health and Development (Lane et al 2017). The epileptic patient was injected with 250 Mbq of FDG, one hour before the 15 min PET acquisition. The neuroscience volunteer was injected with 370 MBq of florbetapir and scanned dynamically for one hour, starting at the injection time. The last ten minutes were used as a measurement of amyloid deposition, which for the participant was negative.

In this section we analyze the impact of various choices within SPDHG on its performance, from randomness over sampling to preconditioning. The test case is total variation prior as defined in (9).

4.2.1. Randomness

Figure 3 shows the effect of randomness where we compare the deterministic PDHG to SPDHG with uniform sampling and scalar step sizes (13) for two different number of subsets. The horizontal axis reflects the number of projections in each algorithm, we call one full projection for the whole data one 'epoch'. Here and in the following dashed lines represent deterministic and solid lines randomized algorithms. We can easily see that both random variants are faster than then deterministic PDHG. Moreover, the randomized SPDHG becomes faster by choosing a larger number of subsets.

Zoom In Zoom Out Reset image size

4.2.2. Sampling

The effect of different choices of sampling is shown in figure 4. We compare two different samplings: uniform sampling and balanced sampling. The uniform sampling chooses all indices $i \in [n]$ with equal probability p _i  =  1/n. In contrast, for balanced sampling we choose with uniform probability either data or prior. If we choose data, then we select a subset again randomly with uniform probability. Thus, the probability for each subset of the data to be selected is $p_i = 1 / (2m), i \in [m]$ and for the prior to be selected p _n  =  1/2.

Zoom In Zoom Out Reset image size

We make two observations. First, balanced sampling is always faster than uniform sampling. This shows the importance of updating the dual variable associated to the prior. Second, for either sampling choosing a larger number of subsets again improves the performance.

4.2.3. Preconditioning

As shown in theorem 2, the step size parameters $\newcommand{\mT}{\mathbf T} \mT$ and $\newcommand{\mS}{\mathbf S} \mS_i$ can be chosen either as scalars (13) or as vectors (14), the latter can be seen as a form of preconditioning. Results are shown in figure 5, where we see that preconditioning may accelerate the convergence of either the deterministic PDHG or the randomized SPDHG. Moreover, combining randomization and preconditioning yields an even faster algorithm.

Zoom In Zoom Out Reset image size

4.2.4. Performance of proposed algorithm

Based on the previous three examples, we propose to combine randomization, balanced sampling and preconditioning, which we refer to as SPDHG+. Figure 6 shows the visual performance of PDHG and SPDHG+. In contrast to the deterministic PDHG, the proposed SPDHG+  yields a good approximation of the optimal solution after only 10 epochs.

Zoom In Zoom Out Reset image size

4.3.1. Anisotropic total variation

Anisotropic total variation decouples the penalization of the derivatives. The mathematical model is similar to the isotropic TV model (9), the only difference being the norm how the total variation is measured: $f_n = \alpha \|\cdot\|_{1,1}$. It can be seen in figure 7 for $\newcommand{\dataAmy}{\texttt{florbetapir}} \dataAmy$ that with randomization and preconditioning only a few epochs are needed to obtain a good approximation of the optimal solution.

Zoom In Zoom Out Reset image size

4.3.2. Directional total variation

Anatomical information from a co-registered MRI is available on combined PET-MR scanners. The structural information of the anatomy can be utilized by the directional total variation prior, see Ehrhardt and Betcke (2016), Ehrhardt et al (2016), Schramm et al (2017) and Bungert et al (2018) for details. The mathematical model is similar to the total variation model (9), except for an additional matrix $\newcommand{\mD}{\mathbf D} \mD$. Thus, the only difference is $\newcommand{\mA}{\mathbf A} \newcommand{\mD}{\mathbf D} \mA_n = \mD \nabla$. A numerical example is shown in figure 8 for the data set $\newcommand{\dataAmy}{\texttt{florbetapir}} \dataAmy$.

Zoom In Zoom Out Reset image size

4.3.3. Total generalized variation

More sophisticated regularization can be achieved by the total generalized variation (TGV) (Bredies et al 2010, Bredies and Holler 2015)

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\TGV}{{\rm TGV}} \displaystyle \TGV_{\alpha_0, \alpha_1}(u) = \inf_w \left\{\alpha_0 \|\nabla u - w\|_{2,1} + \alpha_1 \|{\mathcal E} w\|_{2,1} \right\} \nonumber \end{align*}$$

which can balance first and second order regularization and achieves edge-preserved reconstruction while avoiding the stair-casing artifact. We can solve the TGV regularized PET reconstruction problem by solving problem (8) with the assignment $x = (u, w)$ and

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\im}{{\rm Im}} \newcommand{\mP}{\mathbf P} \newcommand{\mI}{\mathbf I} \newcommand{\mA}{\mathbf A} \displaystyle n = m + 2, & \quad\mA_i = (\mP_i, 0), i \in [m], \quad\mA_{n-1} = (\nabla, - \mI), \quad\mA_n = (0, {\mathcal E})\nonumber \\ g(x) = \imath_+(u), & \quad \,\,\,f_i = D_i, i \in [m], \quad \quad \quad f_{n-1} = \alpha_0 \|\cdot\|_{2,1}, \quad f_n = \alpha_1 \|\cdot\|_{2,1} \, , \nonumber \end{align*}$$

where ${\mathcal E}$ is a symmetrized gradient operator, see Bredies et al (2010) and Bredies and Holler (2015) for more details.

The numerical results shown in figure 9 are in line with the previous findings indicating that randomization and preconditioning can significantly speed up the reconstruction. However, we notice a significant increase in performance by increasing the number of subsets from 21 to 252.

Zoom In Zoom Out Reset image size

4.3.4. Comparison of mathematical models

We conclude this section by a comparison of various methods on both data sets in figures 10 and 11. While we leave the detailed visual comparisons to the reader, we would like to note that all these images use the same number of projections so have basically the same computational cost.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size