Journal of Medical Engineering

Volume 2014, Article ID 908984, 15 pages

http://dx.doi.org/10.1155/2014/908984

## Kaczmarz Iterative Projection and Nonuniform Sampling with Complexity Estimates

Department of Computer Science, Tennessee State University, 3500 John A. Merritt Boulevard, Nashville, TN 37209-1500, USA

Received 21 August 2014; Revised 26 October 2014; Accepted 27 October 2014; Published 15 December 2014

Academic Editor: Hengyong Yu

Copyright © 2014 Tim Wallace and Ali Sekmen. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Kaczmarz’s alternating projection method has been widely used for solving mostly over-determined linear system of equations in various fields of engineering, medical imaging, and computational science. Because of its simple iterative nature with light computation, this method was successfully applied in computerized tomography. Since tomography generates a matrix with highly coherent rows, randomized Kaczmarz algorithm is expected to provide faster convergence as it picks a row for each iteration at random, based on a certain probability distribution. Since Kaczmarz’s method is a subspace projection method, the convergence rate for simple Kaczmarz algorithm was developed in terms of subspace angles. This paper provides analyses of simple and randomized Kaczmarz algorithms and explains the link between them. New versions of randomization are proposed that may speed up convergence in the presence of nonuniform sampling, which is common in tomography applications. It is anticipated that proper understanding of sampling and coherence with respect to convergence and noise can improve future systems to reduce the cumulative radiation exposures to the patient. Quantitative simulations of convergence rates and relative algorithm benchmarks have been produced to illustrate the effects of measurement coherency and algorithm performance, respectively, under various conditions in a real-time kernel.

#### 1. Introduction

Kaczmarz (in [1]) introduced an iterative algorithm for solving a consistent linear system of equations with . This method projects the estimate onto a subspace normal to the row at step cyclically with . The block Kaczmarz algorithm first groups the rows into matrices and then it projects the estimate onto the subspace normal to the subspace spanned by the rows of at step cyclically with . Obviously, the block Kaczmarz is equivalent to the simple Kaczmarz for . The Kaczmarz method is a method of alternating projection (MAP) and it has been widely used in medical imaging as an algebraic reconstruction technique (ART) [2, 3] due to its simplicity and light computation. Strohmer and Vershynin [4] proved that if a row for each iteration is picked in a random fashion with probability proportional with norm of that row, then the algorithm converges in expectation exponentially with a rate that depends on a scaled condition number of (not on the number of equations). Needell (in [5]) extended the work of [4] for noisy linear systems and developed a bound for convergence to the least square solution for . Needell also developed a randomized Kaczmarz method that addresses coherence effects [6], and she analyzed the convergence of randomized block Kaczmarz method [7]. Chen and Powell (in [8]) consider a random measurement matrix instead of random selection of measurements. Galántai (in [9, 10]) provides convergence analysis for block Kaczmarz method by expanding the convergence analysis (based on subspace angles) of Deutsch [11]. Brezinski and Redivo-Zaglia (in [12]) utilizes the work of Galantai for accelerating convergence of regular Kaczmarz method.

The work of this paper endeavors to make the following contributions.(i)Research on regular and randomized Kaczmarz methods appears disconnected in the literature. Even though convergence rates have been studied separately, the link between them has not been explored sufficiently.(ii)A new randomization technique based on subspace angles has been developed which indicates an advantage with coherent data measurements.(iii)A further method is introduced which orthogonalizes the subspace blocks in order to mitigate the coherency. Convergence is consistent with statistical expectations from theory and simulations.(iv)The effects of measurement coherence are observed in the literature and illustrated in our simulations with norm and angle based iteration randomization.(v)A broader review and mathematical analysis of common methods is presented from both statistical and deterministic perspectives.(vi)Numerical simulations are provided to illustrate the typical effects of nonuniform sampling coherency upon convergence of Kaczmarz methods.(vii)Kaczmarz inversions versus matrix size were performed to allow comparison of the relative convergence rates of various well-known methods using typical hardware and software. The results show relative computational complexities of common methods for simple and randomized Kaczmarz, including the randomized Kaczmarz orthogonal block.

#### 2. Methods and Materials

Data inversion and reconstruction in computed tomography is most often based upon the iterative Kaczmarz algorithm due to the performance. First, in this section, given the number of methods currently in the literature, a broad but extensive overview of the mathematical theory for the more common methods is provided, such as simple, block, and randomized Kaczmarz. Where possible, the convergence results are compared from both random and deterministic perspectives to demonstrate similar results and convergence analysis methods. The concept of subspace projections is reviewed and the connection to iteration is noted.

Next, two new methods are proposed and analyzed in the context of coherent data measurements. These methods allow the algorithm to adapt to the changing environment of the sampling measurement system, in order to mitigate coherency.

Simulated methods for data acquisition under uniform and nonuniform X-ray beam measurements are included, and convergence results are computed comparing simple and random row selection methods.

Lastly, after the algorithm methods section, a brief overview of the methods used to obtain the complexity estimates is presented. Methods include using common software and hardware under dedicated kernel conditions. Simulations for Kaczmarz convergence and complexity were performed using Octave software.

##### 2.1. Convergence of Regular Block Kaczmarz Method

Let be the solution of consistent **,** where is full column rank. Let be row-partitioned as where . Then, the simple block Kaczmarz update is as follows:
where is the section of that corresponds to the rows of . Note that since is full row rank, is the right pseudo-inverse of . This is equivalent to
Note that is the projection matrix for projection of the range of :
For one cycle of the blocks,
Note that if is a full column rank with , then the simple block Kaczmarz update is as follows:
where is the pseudo-inverse of and is the orthogonal projection onto . Then, we get the same equation as (3), and subsequently we get (5),

##### 2.2. Exponential Convergence

Theorem 1. *Let be the solution of consistent , where is full rank. Let be row-partitioned as where and . Then, the simple block Kaczmarz converges exponentially and the convergence rate depends of the number of blocks.*

*Proof. *By (3) and orthogonal projection,
So,
depends on the initial condition , and this dependence is scale-invariant. To see this, let and consider where . By (4),
We will first show that if , then . By the way of contradiction, assume that and . By (9),
and therefore for all . By (3), for all . By (8), we get for all . This implies that for all . So,
Since is full column rank, we get , which is a contradiction. So we know that (for one full cycle of -iterations).

By compactness, there exists an such that, for all ,
By (10) and (13),
Now consider iteration for cycles:
Therefore, we conclude that the exponential decay depends on the number of blocks . Note that for regular simple Kaczmarz and the exponential decay depends on the number of rows in this case. The randomized Kaczmarz algorithm proposed by Strohmer and Vershynin [4] avoids this and it converges in expectation as , where is the scaled condition number of matrix with is the pseudoinverse of .

##### 2.3. Iterative Subspace Projection Approach

We can use the following theorem (in [9, 11]) to show the convergence of regular block Kaczmarz method.

Theorem 2. *Let be closed subspaces of the real Hilbert space . Let and be orthogonal projection on . Then, for each ,
**
where is the orthogonal intersection projection.*

The block Kaczmarz is an alternating projection method with . Also, Since is full column rank, and . After cycles, By Theorem 2, and . Galántai in [9] gives a bound for in terms of principle angles between ’s.

##### 2.4. Bound for Block Kaczmarz in Terms of Principle Angles

Smith et al. established the following convergence theorem for applying the alternating projection method in tomography [9, 13].

Theorem 3. *Let be closed subspaces of the real Hilbert space . Let and be orthogonal projection on ( is the orthogonal intersection projection). Let ; then, for each and integer ,
**
where is the orthogonal intersection projection.*

In the special case of the block Kaczmarz, we have , . Also, and . Since is full column rank, and . Therefore, after cycles, where is as defined in Theorem 3. Note that the exponential decay rate depends on the number of blocks as shown below: Galántai in [9] developed another bound (for ) by defining a new matrix for each block as follows.

Theorem 4. *Let be the solution of for a consistent linear system with . Let be row-partitioned as where . Let and be the Cholesky decomposition of . Define and . Then, for each and integer ,
*

*2.5. Special Case: Simple Kaczmarz for *

*Note that this section assumes that . The block Kaczmarz algorithm is equivalent to the simple Kaczmarz algorithm if the number of blocks is equal to the number of rows . In this case, . Therefore, and . This implies that . Then, is defined as
Assume the matrix has normalized rows and we pick a row at each iteration uniformly randomly. Note that this assumption is feasible as scaling a row of and the corresponding measurement in does not change the solution .*

* is the Gram matrix with . Since and is full rank, we have . Using Theorem 4, we get the following deterministic bound:
Since is normalized, we get and therefore
Bai and Liu (in [14]) uses the Meany Inequality to develop a general form of this inequality.*

*2.6. Randomized Kaczmarz Method*

*Several methods of randomized Kaczmarz are discussed in this section.*

*2.7. Randomization Based on Row Norms Method*

*
Strohmer and Vershynin (in [4]) developed a randomized Kaczmarz algorithm that picks a row of in a random fashion with probability proportional with norm of that row. They proved that this method has exponential expected convergence rate. Since the rows are picked based on a probability distribution generated by the norms of the rows of , it is clear that scaling some of the equations does not change the solution set. However, it may drastically change the order of the rows picked at each iteration. Censor et al. discuss (in [15]) that this should not be better than the simple Kaczmarz as picking a row based on its norm does not change the geometry of the problem. Theorem 5 is from [4].*

*Theorem 5. Let be the solution of Then, Algorithm 1 converges to in expectation, with the average error
where is the scaled condition number of matrix with is the left pseudoinverse of .*