# Continuous Analog of Accelerated OS-EM Algorithm for Computed Tomography

Authors: Kiyoko Tateishi; Yusaku Yamaguchi; Omar M. Abou Al-Ola; Tetsuya Yoshinaga
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

The maximum-likelihood expectation-maximization (ML-EM) algorithm is used for an iterative image reconstruction (IIR) method and performs well with respect to the inverse problem as cross-entropy minimization in computed tomography. For accelerating the convergence rate of the ML-EM, the ordered-subsets expectation-maximization (OS-EM) with a power factor is effective. In this paper, we propose a continuous analog to the power-based accelerated OS-EM algorithm. The continuous-time image reconstruction (CIR) system is described by nonlinear differential equations with piecewise smooth vector fields by a cyclic switching process. A numerical discretization of the differential equation by using the geometric multiplicative first-order expansion of the nonlinear vector field leads to an exact equivalent iterative formula of the power-based OS-EM. The convergence of nonnegatively constrained solutions to a globally stable equilibrium is guaranteed by the Lyapunov theorem for consistent inverse problems. We illustrate through numerical experiments that the convergence characteristics of the continuous system have the highest quality compared with that of discretization methods. We clarify how important the discretization method approximates the solution of the CIR to design a better IIR method.


---

## Body

## 1. Introduction

In computed tomography (CT), iterative image reconstruction (IIR) [1–4] gives better quality images compared with filtered backprojection, which is a transform reconstruction method based directly on the Radon inversion formula. Among IIR methods, the maximum-likelihood expectation-maximization (ML-EM) [4] algorithm shows a better performance with respect to cross-entropy minimization or likelihood maximization [5, 6] and is most suitable for the reconstruction of CT since its derivation from underlying the Poisson distribution. To accelerate the ML-EM algorithm, which has a drawback of slow convergence, the ordered-subset variation of the ML-EM, known as the ordered-subsets expectation-maximization (OS-EM) [7, 8], is an effective and popular method. In addition to the OS-EM method, a larger power factor that does not cause divergence in the iterative process was introduced [9–13] to further accelerate the convergence rate. It was asserted [12] that a power-based ML-EM algorithm with increased power (step size) resulted in not only accelerating the convergence rate but also maximizing the likelihood values. However, there is currently no general formula for optimizing the convergence rate by choosing appropriate values of the subsets number and the step size, which depend on the tomography system condition.In this paper, we propose a continuous analog to the power-based accelerated OS-EM, which is based on the approach of continuous-time dynamical optimization [14–19]. The system is described by a switched nonlinear differential equation with piecewise smooth vector fields. Reconstructed image pixels were obtained by using the initial value problem of differential equations describing the dynamical system, for example, a continuous-time image reconstruction (CIR) system. The proposed hybrid dynamical system has a different vector field from that of our previously presented continuous system [20–22]. We indicated that discretizing the differential equations using the geometric multiplicative first-order expansion of the nonlinear vector field leads to the exact same iterative formula of the power-based OS-EM with the scaling parameter as a step size of discretization. Iterative algorithms with the Runge-Kutta (RK), as well as Euler methods that are applied to multiplicatively and additively discretize the CIR system, can be used for image reconstruction. This paper shows that the discretization of a differential equation system can construct an IIR algorithm for solving bound-constrained inverse problems in CT.The proposed CIR system has an advantage in that the convergence of nonnegatively constrained solutions to a stable equilibrium is globally guaranteed by the Lyapunov stability theorem [23] for consistent inverse problems. Specifically, we prove theoretically that a weighted Kullback-Leibler (KL) divergence [24] between a solution and the equilibrium can become a common Lyapunov function for a continuous OS-EM system with an arbitrary number of subsets and under arbitrary switching signals. A cross-entropy function of the KL-divergence between projection and backprojection monotonically decreases along the solution to the continuous ML-EM system. This means the likelihood function monotonically increases in time.To study the cross-entropy convergence property in the case of inconsistent inverse problems, we applied numerical simulations using a dataset acquired from a clinical SPECT scanner for the CIR system and three kinds of discretization methods: the Euler discretization or equivalently the power-based OS-EM method, its rescaled method, and the third-order RK discretization method. We see that the convergence of solutions to the CIR system has the highest quality compared with that of other iterative methods. We clarify how important the discretization method approximates the solution of the CIR in designing a better IIR method.


## 2. System Description

Image reconstruction in CT is a problem to obtain unknown variablex∈R+J for pixel values satisfying(1)y=Ax,where y∈R++I and A∈R+I×J, respectively, denote the projection and a projection operator representing the discretization of the Radon transform (R+ and R++, resp., indicate the set of nonnegative and positive real numbers). If the system in (1) has a nonnegative solution, it is consistent; otherwise, it is inconsistent. The problem can be reduced to finding x using an optimization scheme that minimizes an objective function with respect to the system of (1).Before describing the proposed system, we will prepare a number of definitions and notations. Letym∈R+Im and Am∈R+Im×J be, respectively, a subvector consisting of Im partial elements of y and a submatrix of A with the same corresponding rows of ym for m=1,2,…,M, with M indicating the total number of divisions. We also define(2)λj≔∑i=1IAij-1,λjm≔∑i=1ImAijm-1,j=1,2,…,J,where Aij is the (i,j) element of A and λ≔(λ1,λ2,…,λJ)⊤. To simplify the description, we denote vector-valued functions Log⁡(α)≔(log⁡α1,log⁡α2,…,log⁡(αL))⊤ of each element in vector α=(α1,α2,…,αL)⊤ and Exp(β)≔(exp⁡β1,exp⁡β2,…,exp⁡(βL))⊤ of β=(β1,β2,…,βL)⊤, respectively.The objective function considered in this paper is expressed in the generalized KL-divergence(3)KLα,β=∑l=1Lαllog⁡αlβl+βl-αl=α⊤Log⁡α-Log⁡β+v⊤β-α,for the given nonnegative vectors α and β with αl and βl denoting the lth elements of α and β, respectively, and v indicating an all-ones vector. The divergence KL(α,β), known as Csiszár’s I-divergence measure [25], for the vectors α and β in R+L is nonnegative with KL(α,β)=0 if and only if α=β.The proposed method for solving the tomographic inverse problem is based on the use of an initial value problem of the switched nonlinear dynamical system [26, 27], which consists of a family of M subsystems(4)dxtdt=XtLog⁡ΛmAm⊤Exp⁡Log⁡ym-Log⁡Amxt,t-kτ∈tm-1,tm,t≥0,x0=x0∈R++Jfor m=1,2,…,M, k=0,1,…, and a series of times 0=t0<t1<t2<⋯<tM=τ, where X(t)≔diag(x(t)) and Λ≔diag(λ). Note that each subsystem is described by autonomous differential equations with a sufficiently smooth vector field. The solutions to the hybrid dynamical system [27] are constructed through the connection of the last state of the previous mth subsystem and the initial state of the next (m+1)th subsystem at every tm, for m=1,2,…,M, and setting M+1 to be 1 for cyclic switching process.


## 3. Theoretical Results

In this section, the theoretical results of the solutions to the hybrid dynamical system in (4) are given.We first show that a solution can be made to stay positive. This property is preserved whether or not the inverse problem is consistent and indicates that the CIR system does not produce images with unphysical negative pixel values.Proposition 1.
If the initial valuex0 of the switched dynamical system in (4) is chosen in R++J, its solution ϕ(t,x0) behaves in R++J for all t≥0.Proof.
The vector field of thejth element of the mth subsystem can be rewritten as dxj/dt=xjfjm(x) with the function fjm. Therefore, dϕj/dt≡0 holds for any j on the subspace satisfying xj=0, which means the subspace becomes invariant, and, on the basis of the uniqueness criteria of solutions to the initial value problem, any flow cannot pass through every invariant subspace. This leads to the proof.Assuming that the individual subsystems have the common equilibriume satisfying ym=Ame for any m=1,2,…,M, we can prove the existence of a Lyapunov function for all subsystems in (4), which guarantees that the corresponding switched system has a stable equilibrium e.The OS-EM [7, 8] was introduced to accelerate the computation time of ML-EM. However, even in the consistent case, the convergence proof in [7] requires the values of λjm to be independent of the subset m, which is referred to as the “subset balance.” This condition is restrictive in practical applications. Byrne [6] presented a more general sufficient condition called “a generalized subset balance,” meaning that λjm is separable; positive values δm exist such that(5)λjm=δmλj,j=1,2,…,J,where m=1,2,…,M. Our convergence proof on the stability of an equilibrium observed in the switched system in (4) needs the same condition as (5).One of the main results is a stability theorem for the continuous OS-EM system as follows.Theorem 2.
If there exists a unique solutione∈R++J to the system y=Ax with the matrix A satisfying the condition of (5), the equilibrium e of the dynamical system in (4) is asymptotically stable.Proof.
This follows from Lyapunov’s stability theorem. We define a possible candidate for a Lyapunov function as a weighted KL-divergence,(6)Vx=∑j=1Jλj-1KLej,xj=∑j=1Jλj-1ejlog⁡ejxj+xj-ej,which is positive definite and is well-defined via Proposition 1 when initial value x0 is chosen in R++J. It can be written as(7)Vx=∑j=1Jλj-1∫ejxjs-ejsds.Using the concavity of the log function and Jensen’s inequality, we then calculate its derivative with respect to the dynamical system in (4) as(8)dVdtx4=∑j=1Jxj-ejλjxjdxjdt=-∑j=1Jejλjlog⁡λjmAjm⊤Exp⁡Log⁡ym-Log⁡Amx+∑j=1Jxjλjlog⁡λjmAjm⊤Exp⁡Log⁡ym-Log⁡Amx≤-∑j=1JejλjλjmAjm⊤Log⁡ym-Log⁡Amx+∑j=1JxjλjλjmAjm⊤Exp⁡Log⁡ym-Log⁡Amx-1=-δm∑i=1Imyimlog⁡yim-log⁡Amxi+δm∑i=1ImAmxiexp⁡log⁡yim-log⁡Amxi-1≤-δm∑i=1Imyim1-exp⁡log⁡Amxi-log⁡yim+δm∑i=1ImAmxiexp⁡log⁡yim-log⁡Amxi-1=0,where m=1,2,…,M, and x in R++J. The derivative equals zero at x=e∈R++J. Consequently, the hybrid system consisting of the family of subsystems in (4) has a common Lyapunov function defined by (6), so the equilibrium e of the system is asymptotically stable under arbitrary switching signals.The proof of Theorem2 can be rephrased as follows. If the system y=Ax has a unique solution e∈R++J, the objective function V(x) in (6) decreases monotonically in time for the solution to the system in (4) with x0∈R++J.We use another Lyapunov function as an objective function that has to be minimized.Theorem 3.
If the systemy=Ax possesses a unique solution e∈R++J, the following objective function W(x) decreases monotonically in time for solutions of the system in (4) with M=1 and x0∈R++J:(9)Wx=KLy,Ax.Proof.
We show that the function(10)Wx=y⊤Log⁡y-Log⁡Ax+v⊤Ax-y,which is positive definite for x∈R++J, can be a Lyapunov function. Because (s-1)/s≤log⁡s≤s-1 for s>0, we can obtain the derivative with respect to the system in (4) with M=1 as follows:(11)dWdtx4=v-Exp⁡Log⁡y-Log⁡Ax⊤AXLog⁡ΛA⊤Exp⁡Log⁡y-Log⁡Ax≤-Exp⁡Log⁡y-Log⁡Ax⊤AXu-Exp⁡-Log⁡ΛA⊤Exp⁡Log⁡y-Log⁡Ax+v⊤AXΛA⊤Exp⁡Log⁡y-Log⁡Ax-u,where u∈RJ and v∈RI indicate all-ones vectors. By putting(12)p≔Exp⁡Log⁡y-Log⁡Axfor simplicity, (11) can be written as(13)dWdtx4≤-p⊤AXu+p⊤AXExp⁡-Log⁡ΛA⊤p+v⊤AXΛA⊤p-v⊤AXu.For the first and third terms of the right-hand side in (13), one gets (14)-p⊤AXu+v⊤AXΛA⊤p=-p⊤AXu+p⊤AΛXA⊤v=p⊤AX-u+ΛA⊤v=0and obtains for the second and fourth terms (15)p⊤AXExp⁡-Log⁡ΛA⊤p-v⊤AXu=u⊤XΛ-1u-v⊤AXu=u⊤XA⊤v-v⊤AXu=0.Therefore,(16)dWdtx4≤0and the Lyapunov function W(x) decreases along the flow. This concludes the proof.


## 4. Discretization

Consider the numerical discretization of the piecewise differential equations describing the CIR hybrid system. The integration of differential equations by discretization via the multiplicative calculus is derived as a multiplicative counterpart of classical Euler and RK methods [28, 29]. The following shows the geometric multiplicative first-order expansion of a nonlinear function that appeared in the vector field. See [29] for details on constructing a multiplicative RK method.We describe thejth equation for the mth subsystem in (4) as (17)dxjdtt=xjtfjmxtat t∈[tm-1+kτ,tm+kτ) for m=1,2,…,M, and nonnegative integer k. When applying the multiplicative Euler method to the mth subsystem, we obtain a single step at a time:(18)xjt+hm=xjtexp⁡hmfjmxt,where hm>0 denotes a step size depending on m. For the block continuous-time hybrid system, by choosing tm-tm-1≕hm and by connecting solutions to the subsystems at the discrete time tn, where t0=0 and(19)tn=∑l=1n-kMhl+k∑l=1Mhlfor n=1,2,…, with k being the floor of n/M, we have the block-iterative form(20)xjtn+1=xjtnexp⁡hmfjmxtn,n=0,1,2,…,m=nmod⁡M+1,xt0=x0.Through this discretization procedure, the discrete time tn is identical to the switching time t=tm+kτ with τ=∑l=1Mhl in the hybrid dynamical system under the corresponding special switching signals.The iterative formula for the continuous-time system in (4) is equivalent to the power-based OS-EM(21)zjn+1=zjnλjmAjm⊤Exp⁡Log⁡ym-Log⁡Amznhm=zjnλjm∑i∈BmAijyiAznihm,n=0,1,2,…,m=nmod⁡M+1,z0=x0,where x(tn) in (20) is substituted with z(n). Here, the collection {B1,B2,…,BM} is a partition or a subset of the index set {i=1,2,…,I}; the definition is the same as the derivation of Am and ym for m=1,2,…,M. The step size hm in (20) and (21) corresponds to the scaling parameters in the OS-EM algorithm. Note that because the step size is derived in the discretization procedure of the hybrid dynamical system, its value does not affect the theoretical results of the solutions for the CIR system defined in Section 2.


## 5. Numerical Example

This section numerically illustrates the convergence property of continuous solutions to the CIR system and iterative solutions to certain IIR procedures derived from its discretizations. We used a set of SPECT projection data available from a publicly accessible database of Monte Carlo simulated datasets for Emission Tomography (the MC-ET database) [30]. The sinogram of a brain scan is shown in Figure 1. The numbers of projections and pixels of the reconstructed images were I= 7,500 (125 bins and 60 directions in 180 degrees) and J=7,569 (87×87 pixels), respectively. The projection data from 60 directions were divided into M nonoverlapping subsets as uniformly as possible with M being 2, 4, and 8. When M=2 and 4, (5) holds with δm=M for m=1,2,…,M; namely, the subsets into which the angles are uniformly split satisfy the subset balance [7] condition. Whereas, nonuniform division into eight subsets gives an unbalanced dataset. The step size h≔hm of IIR was set to 2.2 for any m=1,2,…,M. Although bigger step size is effective for accelerating IIR method, the algorithm diverges or oscillates if it is too big [12, 13]. The value of h was chosen to characterize the difference between lower- and higher-order discretization methods (i.e., OS-EM and RK, resp.) in the convergence process.Figure 1
(a) Sinogram of brain scan and (b) profile of counts along vertical line in middle of sinogram.
(a)
(b)The IIR procedures considered here are the OS-EM (or equivalently the multiplicative Euler method) in (21), a rescaled OS-EM, and the multiplicative third-order RK method of the CIR. The rescaled OS-EM was achieved by multiplying the state z(n+1) in (21) by the rescaling coefficient(22)ζ≔∑i∈Bmyi∑i∈BmAzn+1iat every nth iteration step with m=1,2,…,M; namely, the mapping(23)R:R++J⟶R++J;zn+1⟼ζzn+1was applied as a postprocess to the OS-EM procedure. The multiplicative RK method with third-order is defined as follows.(24)zjn+1=zjnpjmznhm/6qjmzn2hm/3rjmznhm/6,where (25)pjmz≔λjmAjm⊤Exp⁡Log⁡ym-Log⁡Amz,qjmz≔λjmAjm⊤Exp⁡Log⁡ym-Log⁡Amw,rjmz≔λjmAjm⊤Exp⁡Log⁡ym-Log⁡Amvwith wj≔zjpjm(z)hm/2 and vj≔zjpjm(z)-hmqjm(z)2hm, for j=1,2,…,J and m=(nmod⁡M)+1. To solve the differential equation in (4) for CIR, we used the solver ode45, which is a variable-step RK method in MATLAB (MathWorks, Natick, USA).Figure2 plots the objective function W in (9), which is a cross-entropy function and a Lyapunov function for the CIR system with M=1, obtained by using CIR, the OS-EM, its rescaling, and the RK at time t∈[0,10M] and the nth discrete points for n=0,1,…,⌊10M/h⌋, with the number of subsets M being equal to 1, 2, 4, and 8. Since there is no theory available to guarantee the convergence of iterative solutions to IIR with a step size greater than one, the iterative points diverge or oscillate depending on the conditions set, including the selection of M. We see that the OS-EM, having the relatively large step size, accumulates an error per iteration because of the insufficient accuracy of the computation by using first-order approximation. As shown in Figure 3, which represents the reconstructed images at t=⌊10M/h⌋h, the OS-EM algorithm with M≥4 and the rescaled OS-EM with M=8 failed to reconstruct images corresponding to the divergence of the objective function W. However, a further experiment indicated that the iterative sequence generated by the third-order RK algorithm with M≤8 did not diverge nor oscillate for h≤2.4. One can confirm that the OS-EM and the other iterative procedures originate from the meaning of convergence from the CIR system and the third-order RK has the best robust performance among three IIR algorithms. It is important when designing a better IIR method to choose a discretization method that approximates the solution of the CIR.Figure 2
Values ofW for CIR, OS-EM, rescaled OS-EM, and third-order RK (abbreviated as CIR, OSEM, ROSEM, and MRK3, resp.) at time t, which is the independent variable of the continuous system (CIR), and at discretized time series defined by tn=nh in t, with h=2.2 and n=0,1,2,…,⌊10M/h⌋, for the discrete systems (OSEM, ROSEM, and MRK3).
(a)
M
=
1
(b)
M
=
2
(c)
M
=
4
(d)
M
=
8Figure 3
Images reconstructed by CIR, OS-EM, rescaled OS-EM, and third-order RK (abbreviated as CIR, OSEM, ROSEM, and MRK3, resp.) withM=1, 2, 4, and 8. Color scale is common to all four images for each M.


## 6. Concluding Remarks

We proposed a hybrid dynamical system (CIR system), a multiplicative Euler discretization that is exactly the same as a power-based accelerated OS-EM. We have theoretically shown the positiveness of solutions and the stability of a common equilibrium corresponding to the exact solution of the consistent CT inverse problem. The proof of the convergence to the stable equilibrium was made using the Lyapunov stability theorem. As a result, the common Lyapunov function or the KL-divergence measure monotonically decreases along the time course.We confirmed through numerical experiments that the convergence characteristics of the CIR system had the highest quality compared with that of any discretization method. This fact illustrates that the OS-EM method and its rescaled method are considered to be derived from discretization of the CIR system.



---

*Source: 1564123-2017-08-06.xml*