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.
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 variable for pixel values satisfyingwhere 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 , with indicating the total number of divisions. We also definewhere is the element of and . To simplify the description, we denote vector-valued functions of each element in vector and of , respectively.
The objective function considered in this paper is expressed in the generalized KL-divergencefor the given nonnegative vectors and with and denoting the th elements of and , respectively, and indicating an all-ones vector. The divergence , known as Csiszár’s -divergence measure [25], for the vectors and in is nonnegative with 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 subsystemsfor , , and a series of times , where and . 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 th subsystem at every , for , and setting to be 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 value of the switched dynamical system in (4) is chosen in , its solution behaves in for all .
Proof. The vector field of the th element of the th subsystem can be rewritten as with the function . Therefore, holds for any on the subspace satisfying , 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 , 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 thatwhere . 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 is chosen in . It can be written asUsing the concavity of the log function and Jensen’s inequality, we then calculate its derivative with respect to the dynamical system in (4) aswhere , 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 .
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 and :
Proof. We show that the functionwhich is positive definite for , can be a Lyapunov function. Because for , we can obtain the derivative with respect to the system in (4) with as follows:where and indicate all-ones vectors. By puttingfor simplicity, (11) can be written asFor the first and third terms of the right-hand side in (13), one gets and obtains for the second and fourth terms Therefore,and the Lyapunov function 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 the th equation for the th subsystem in (4) as at for , and nonnegative integer . When applying the multiplicative Euler method to the th subsystem, we obtain a single step at a time:where denotes a step size depending on . For the block continuous-time hybrid system, by choosing and by connecting solutions to the subsystems at the discrete time , where andfor , with being the floor of , we have the block-iterative formThrough this discretization procedure, the discrete time is identical to the switching time with 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-EMwhere in (20) is substituted with . Here, the collection is a partition or a subset of the index set ; the definition is the same as the derivation of and for . 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.
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 7,500 (125 bins and 60 directions in 180 degrees) and ( pixels), respectively. The projection data from 60 directions were divided into nonoverlapping subsets as uniformly as possible with being 2, 4, and 8. When and 4, (5) holds with for ; 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 of IIR was set to 2.2 for any . 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.
(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 in (21) by the rescaling coefficientat every th iteration step with ; namely, the mappingwas applied as a postprocess to the OS-EM procedure. The multiplicative RK method with third-order is defined as follows.where with and , for and . 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).
Figure 2 plots the objective function in (9), which is a cross-entropy function and a Lyapunov function for the CIR system with , obtained by using CIR, the OS-EM, its rescaling, and the RK at time and the th discrete points for , 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 , the OS-EM algorithm with and the rescaled OS-EM with failed to reconstruct images corresponding to the divergence of the objective function . However, a further experiment indicated that the iterative sequence generated by the third-order RK algorithm with did not diverge nor oscillate for . 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.
(a)
(b)
(c)
(d)
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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was partially supported by JSPS KAKENHI Grant no. 15K06110.