On global randomized block Kaczmarz method for image reconstruction

1.   Introduction

Image reconstruction plays a vital role in magnetic resonance imaging ^[1], X-ray computed tomography ^[2] and radio astronomy imaging applications ^[3]. Reconstruction methods ^[4,5,6,7] are mainly divided into analytical reconstruction and iterative reconstruction, which are based on the Radon transform ^[8], the inverse Radon transform, and the projection slice theorem serving as mathematical foundations. The most commonly used analytical reconstruction methods are all in the form of the filtered back projection (FBP) method ^[9]. Iterative reconstruction methods formulate the final result as a solution of linear systems ^[10].

For example, image reconstruction, as a mathematical process, generates tomographic images from X-ray ^[11] projection data acquired at many different angles in computed tomography (CT). If there has a test image with $n \times n$ pixels, $p$ is the number of X-ray beams, and $t$ is the scanning angle, which ranges from ${0^ \circ}$ to ${179^ \circ }$, with intervals of ${1^ \circ }$. The X-ray beam $X_s$ passing through the test image is shown in the following figure.

In Figure 1, ${w_{X_s^{\left({t} \right)}}^{_{\left({i, j} \right)}}}$ represents the length of the X-ray beam $X_s (s = 1, 2, \cdots, p)$ when it passes through the pixel $\left\{ {i, j} \right\}$ along the line $l_{X_s}^{(t)} (t = 0, 2, \cdots, 179)$, and $\mu ^{{\left({i, j} \right)}}$ is the average attenuation coefficient for each pixel $\left\{ {i, j} \right\}$. The data ${g_{X_s^{\left({t} \right)}}^{_{\left({i, j} \right)}}}$ is measured when an X-ray beam $X_s^{(t)}$ crosses the pixel $\left\{ {i, j} \right\}$ and intensity is absorbed. Thus, the image reconstruction problem amounts to solving the linear systems with multiple right-hand sides of the form

where

The linear systems with multiple right-hand sides (1.1) can be written as

where $l = 1, 2, \cdots, n$. The numerical methods ^{[12,13,14,15]} for solving linear systems (1.2) have seen a significant maturation. Among them the Kaczmarz method ^[16] is a classic and effective row-action method. Out of its simplicity, the Kaczmarz method has been a preferable solution tool in many fields, such as image reconstruction ^[15,17], computerized tomography ^[9,18] and signal processing ^[19,20]. Based on the Kaczmarz method, a series of row iteration methods are proposed ^[21,22,23].

In order to reduce the computational costs, Censor proposed the block Kaczmarz method ^[18], Niu ^[24] presented the greedy block Kaczmarz method, Chen ^[25] declared the randomized double block Kaczmarz method and raised the upper bound of the error estimate in expectation. The approximate solution $x_k^{\left(l \right)} (l = 1, \cdots, n)$ can be obtained by solving the linear systems (1.2) with one of the above-mentioned solving methods, thereby getting the approximate solution $X_k = \left[ x_k^{\left(1 \right)}, \cdots, x_k^{\left(n \right)} \right]$ of the linear systems with multiple right-hand sides. In the above solving process, the calculation and storage requirements increase with the iteration.

Instead of solving each of the $n$ linear systems independently by using some iterative methods, it is more efficient to solve the linear systems with multiple right-hand sides globally. A great deal of research has been finished, including the block conjugate gradient method ^[26], the block generalized minimal residual method ^[27] and the GMRES seed projection method ^[28], see also ^{[29,30,31,32,33,34]} and the references therein.

In regard to the problem of image reconstruction, we transform it into solving a linear systems with multiple right-hand sides and propose the idea of global iteration. Based on the idea of K-means clustering, we put forward the global randomized block Kaczmarz method, and use this method to solve the linear systems with multiple right-hand sides effectively, so as to image reconstruction. Theoretical analysis proves the convergence of the global randomized block Kaczmarz method. Simulation results show that the reconstructed image quality of this method is superior to that of the filtered back projection method.

In this paper, we adopt the following notations. Let $A^*$ be the conjugate transpose of matrix $A$, $A^\dagger$ the Moore-Penrose pseudoinverse, ${\lambda _{\min }}\left(A \right)$ the smallest nonzero eigenvalue, ${\left\| A \right\|_2}$ the spectral norm and ${\left\| A \right\|_F}$ the Frobenius norm respectively. $B^{(i)}$ represents the $i{\rm{th}}$ column of matrix $B$. Let $\left[ m \right]: = \left\{ {1, 2, \cdots, m} \right\}$, the index set $J = \left\{ {{J_1}, {J_2}, \cdots, {J_q}} \right\} \subseteq \left[ m \right]$ satisfies ${J_i} \cap {J_j} = \emptyset$ and $\cup _{i = 1}^q{J_i} = \left[ m \right]$. ${A_{{J_i}}}$ is the row submatrix indexed by ${J_i}$. ${X_ \star } = {A^ \dagger }B$ denotes the least Euclidean-norm solution of the linear systems with multiple right-hand sides (1.1). ${\mathbb{E}{_k}}$ indicates the expected value conditional on the first $k$ iterations, that is,

where ${J_l}\left({l = 0, 1, \cdots, k - 1} \right)$ is the ${J_l}{\rm{th}}$ submatrix chosen at the $l{\rm{th}}$ iteration, and then we can obtain that ${\mathbb{E}}\left[ {{{\mathbb{E}}_k}\left[ \cdot \right]} \right] = {\mathbb{E}}\left[ \cdot \right]$.

2.   Global randomized block Kaczmarz method

K-means clustering ^[35] is the simplest unsupervised learning method as well as the most widely used to solve clustering problems. Recently, several Kaczmarz-type methods that embed the clustering methods appeared, see e.g., ^[36,37], and the complete process of the K-means method is given in Method 2.1.

Using this method to divide the rows of coefficient matrix of linear systems with multiple right-hand sides (1.1), we can obtain

where ${A_{{J_i}}} \in {\mathbb{C}^{{m_{{J_i}}} \times n}}$, ${B_{{J_i}}} \in {\mathbb{C}^{{m_{{J_i}}}}}$ and $\sum\limits_{i = 1}^q {{m_{{J_i}}}} = m$. Starting from an initial guess ${X_0}$, the current iterate $X_k$ is orthogonally projected to the hyperplane ${A_{{J_{i_k}}}}{X_k} = {B_{{J_{i_k}}}}$ under the probability criterion

The global randomized block Kaczmarz (MGRBK) method for solving linear systems with multiple right-hand sides can be formulated as

where the index of working submatrix ${J_{i_k}}$ is chosen from the set $J = \left\{ {{J_1}, {J_2}, \cdots, {J_q}} \right\}$ at random, with probability proportional to $\left\| {{A_{{J_{i_k}}}}} \right\|_F^2$. The form of $A_{{J_{i_k}}}^\dagger$ in our simulation experiments can be seen in the reference ^[38]. Then the global randomized block Kaczmarz method can be described in the following.

Note that the index set $J$ is nonempty for all iteration index $k$. According to a large number of numerical experiments, we find that the number of blocks $q$ fluctuates around $\frac{m}{{2000}}$, where $m$ is the number of rows of coefficient matrix $A$. Hence, we have a rough estimate is that $q \approx \frac{m}{{2000}}$. For the convergence property of the global randomized block Kaczmarz method, we establish the following theorem.

Theorem 2.1. Let the linear systems with multiple right-hand sides (1.1), with the coefficient matrix $A \in {\mathbb{C}^{m \times n}}$ and the right-hand side $B \in \mathbb{C}^{m \times n}$, be consistent. The iteration sequence $\left\{ {{X_k}} \right\}_{k = 0}^\infty$ generated by the global randomized block Kaczmarz method starting from any initial approximation ${X_0} \in {\cal R}\left({{A^ * }} \right)$, converges to the unique least-norm solution ${X_ \star } = {A^\dagger }B$ in expectation. Moreover, the solution error in expectation for the iteration sequence $\left\{ {{X_k}} \right\}_{k = 0}^\infty$ obeys

Proof. From the iteration scheme of the global randomized block Kaczmarz method, we can straightforwardly obtain

Since ${B_{{J_{i_k}}}} = {A_{{J_{i_k}}}}{X_ \star }$, it holds that

From the properties of the Moore-Penrose pseudoinverse, we have

Therefore, we know that matrix ${X_k} - {X_{k + 1}}$ is orthogonal to ${X_{k + 1}} - {X_ \star }$ for any $k \ge 0$, in other words, ${\left({{X_k} - {X_{k + 1}}} \right)^*}\left({{X_{k + 1}} - {X_ \star }} \right) = 0$. It follows that

By Eq (2.4), Eq (2.6) can be rewritten as

With this, we obtain

Here the last inequality is achieved with the use of the estimate

where $C \in {\mathbb{C}^{n \times c}}$. Thus, we can get the estimate

Finally, by taking the full expectation on both sides of (2.8), we obtain

By induction on the iteration index $k$, we have

Consequently, the convergence property of the global randomized block Kaczmarz method is proved.

3.   Simulation results

In this section, we investigate the performance of the global randomized block Kaczmarz (MGRBK) method in image reconstruction through practical experiments. We also compare the image reconstruction results of the MGRBK method and filtered back projection (FBP) method, in which the FBP method is realized by using the MATLAB function $\bf{iradon}$. All computations are started from the initial matrix $X_0 = 0$, and the row vectors of coefficient matrix $A$ are divided into $q$ blocks.

Example 3.1. The MGRBK method is used to solve the linear systems with multiple right-hand sides to reconstruct the test image. The test image is obtained through the image processing toolbox function ${\bf{phantom}}$ in MATLAB. The Shepp-Logan phantom is shown in the following Figure.

For the Shepp-Logan phantom with $100 \times 100$ pixels, $p = 4$, the corresponding linear systems with multiple right-hand sides is $AX = B$, where the dimension of coefficient matrix $A$ is $72000 \times 100$. Then we solve this linear systems with multiple right-hand sides by using the MGRBK method, $q = 60$ and the number of iterations is 30. The computational process of image reconstruction is shown in the following figure.

In Figure 3, we can see that at the 18th iteration, the test image can be reconstructed with the global randomized block Kaczmarz method. Through observation, the reconstructed image can retain the details of the original image.

In Figure 4, the image reconstruction is realized by using the global randomized block Kaczmarz method and the filtered back-projection method respectively, and the reconstruction results are compared. From the visual inspection of the reconstructed image, the MGRBK method outperforms the FBP method.

Example 3.2. Select the test image from the photo gallery in MATLAB, and then solve the linear systems with multiple right-hand sides through the MGRBK method to reconstruct the test image. The test image is shown in the following figure.

For the cameraman image with $200 \times 200$ pixels, $p = 5$, the corresponding linear systems with multiple right-hand sides is $AX = B$, where the dimension of coefficient matrix $A$ is $180000 \times 200$. Then we solve this linear systems with multiple right-hand sides by using the MGRBK method. Let $q = 70$ and the number of iterations is 40. The computational process of image reconstruction is shown in the following figure.

In Figure 6, we can see that test image can be reconstructed with the global randomized block Kaczmarz method as reach the 38th iteration. Through observation, the reconstructed image can retain the details of the original image.

In Figure 7, the reconstruction results of the MGRBK method and the FBP method are compared. From the visual inspection of the reconstructed image, the former method is superior to the latter method.

Example 3.3. The test image is obtained through the image processing toolbox function ${\bf{phantom}}$ in MATLAB, which is the same as that in Example 3.1. We compare the reconstructed image of the global randomized block Kaczmarz (MGRBK) method with that of a typical iterative method, that is Kaczmarz method. The reconstructed images of these two iterative methods are shown in the following figure.

From the visual inspection, the reconstructed image obtained by using the MGRBK method has higher resolution. In the MGRBK method, the relative residual (Rres) is defined by

and Rres of the MGRBK method is 1.6049e-04. In ^[39], the relative residual of Kaczmarz method is 0.0015. This means that the approximate solution of MGRBK method is approach to exact solution. Hence, the reconstructed image quality of the MGRBK method is much better.

4.   Conclusions

Iterative reconstruction refers to iterative methods used to reconstruct 2D and 3D images in certain imaging techniques. A variety of reconstruction methods have been investigated to improve the quality of reconstructed images. In this paper, we transform the problem of image reconstruction into solving a linear systems with multiple right-hand sides. In order to solve the linear systems with multiple right-hand sides effectively, we propose the global randomized block Kaczmarz method based on the idea of K-means clustering, and use this method to image reconstruction. Theoretical analysis proves the convergence of this method. Simulation results show that the reconstructed image quality of this method is better than that of the FBP method.

Acknowledgments

This research was supported in part by NSFC under grant 11401305 and 11571171.

Conflict of interest

The authors declare there is no conflicts of interest.