Continuous Analog of Accelerated OSEM Algorithm for Computed Tomography

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 continuoustime 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-basedOS-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.


Introduction
In computed tomography (CT), iterative image reconstruction (IIR) [1][2][3][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 expectationmaximization (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 orderedsubset variation of the ML-EM, known as the orderedsubsets 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][10][11][12][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][15][16][17][18][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][21][22].We indicated that discretizing the differential equations using the geometric multiplicative 2 Mathematical Problems in Engineering 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 powerbased OS-EM method, its rescaled method, and the thirdorder 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.

System Description
Image reconstruction in CT is a problem to obtain unknown variable  ∈   + for pixel values satisfying  = , where  ∈   ++ and  ∈  × + , respectively, denote the projection and a projection operator representing the discretization of the Radon transform ( + and  ++ , 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  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.Let   ∈    + and   ∈    × + be, respectively, a subvector consisting of   partial elements of  and a submatrix of  with the same corresponding rows of   for  = 1, 2, . . ., , with  indicating the total number of divisions.We also define where   is the (, ) element of  and  fl ( for the given nonnegative vectors  and  with  ℓ and  ℓ denoting the ℓth elements of  and , respectively, and V indicating an all-ones vector.The divergence KL(, ), known as Csiszár's -divergence measure [25], for the vectors  and  in   + 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  subsystems for  = 1, 2, . . ., ,  = 0, 1, . .., and a series of times 0 =  0 <  1 <  2 < ⋅ ⋅ ⋅ <   = , where () fl diag(()) and Λ fl 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 th subsystem and the initial state of the next ( + 1)th subsystem at every   , for  = 1, 2, . . ., , and setting  + 1 to be 1 for cyclic switching process.

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 value  0 of the switched dynamical system in ( 4) is chosen in   ++ , its solution (,  0 ) behaves in    ++ for all  ≥ 0.
Proof.The vector field of the th element of the th subsystem can be rewritten as   / =      () with the function    .Therefore,   / ≡ 0 holds for any  on the subspace satisfying   = 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 equilibrium  satisfying   =    for any  = 1, 2, . . ., , 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 .
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    to be independent of the subset , 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    is separable; positive values   exist such that where  = 1, 2, . . ., .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 solution  ∈   ++ to the system  =  with the matrix  satisfying the condition of ( 5), the equilibrium  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, which is positive definite and is well-defined via Proposition 1 when initial value  0 is chosen in   ++ .It can be written as 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 where  = 1, 2, . . ., , and  in   ++ .The derivative equals zero at  =  ∈   ++ .Consequently, the hybrid system consisting of the family of subsystems in (4) has a common Lyapunov function defined by (6), so the equilibrium  of the system is asymptotically stable under arbitrary switching signals.
The proof of Theorem 2 can be rephrased as follows.If the system  =  has a unique solution  ∈   ++ , the objective function () in (6) decreases monotonically in time for the solution to the system in (4) with  0 ∈   ++ .We use another Lyapunov function as an objective function that has to be minimized.

Theorem 3. If the system 𝑦 = 𝐴𝑥 possesses a unique solution 𝑒 ∈ 𝑅 𝐽
++ , the following objective function () decreases monotonically in time for solutions of the system in ( 4) with  = 1 and  0 ∈   ++ : Proof.We show that the function which is positive definite for  ∈   ++ , can be a Lyapunov function.Because ( − 1)/ ≤ log() ≤  − 1 for  > 0, we can obtain the derivative with respect to the system in (4) with  = 1 as follows: where  ∈   and V ∈   indicate all-ones vectors.By putting for simplicity, (11) can be written as For the first and third terms of the right-hand side in (13), one gets and obtains for the second and fourth terms Therefore,   ()        (4) ≤ 0 (16) and the Lyapunov function () decreases along the flow.This concludes the proof.

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 the th equation for the th subsystem in (4) as at  ∈ [ −1 +,   +) for  = 1, 2, . . ., , and nonnegative integer .When applying the multiplicative Euler method to the th subsystem, we obtain a single step at a time: where ℎ  > 0 denotes a step size depending on .For the block continuous-time hybrid system, by choosing   −  −1 š ℎ  and by connecting solutions to the subsystems at the discrete time   , where  0 = 0 and for  = 1, 2, . .., with  being the floor of /, we have the block-iterative form Through this discretization procedure, the discrete time   is identical to the switching time  =   + with  = ∑  ℓ=1 ℎ ℓ 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 where (  ) in ( 20) is substituted with ().Here, the collection { 1 ,  2 , . . .,   } is a partition or a subset of the index set { = 1, 2, . . ., }; the definition is the same as the derivation of   and   for  = 1, 2, . . ., .The step size ℎ  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.

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  = 7,500 (125 bins and 60 directions in 180 degrees) and  = 7,569 (87 × 87 pixels), respectively.The projection data from 60 directions were divided into  nonoverlapping subsets as uniformly as possible with  being 2, 4, and 8.When  = 2 and 4, (5) holds with   =  for  = 1, 2, . . ., ; 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 ℎ fl ℎ  of IIR was set to 2.2 for any  = 1, 2, . . ., .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 ℎ was chosen to characterize the difference between lower-and higher-order discretization methods (i.e., OS-EM and RK, resp.) in the convergence process.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 (+1) in ( 21) by the rescaling coefficient  fl ∑ ∈    ∑ ∈  ( ( + 1))  (22) at every th iteration step with  = 1, 2, . . ., ; namely, the mapping R :   ++ →   ++ ;  ( + 1)  →  ( + 1) was applied as a postprocess to the OS-EM procedure.The multiplicative RK method with third-order is defined as follows.
Figure 2 plots the objective function  in ( 9), which is a cross-entropy function and a Lyapunov function for the CIR system with  = 1, obtained by using CIR, the OS-EM, its rescaling, and the RK at time  ∈ [0,10] and the th discrete points for  = 0, 1, . . ., ⌊10/ℎ⌋, with the number of subsets  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 .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  = ⌊10/ℎ⌋ℎ, the OS-EM algorithm with  ≥ 4 and the rescaled OS-EM with  = 8 failed to reconstruct images corresponding to the divergence of the objective function .However, a further experiment indicated that the iterative sequence generated by the thirdorder RK algorithm with  ≤ 8 did not diverge nor oscillate for ℎ ≤ 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.

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.

Figure 1 :
Figure 1: (a) Sinogram of brain scan and (b) profile of counts along vertical line in middle of sinogram.