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 .
| ||||||||||||||||||
Note that is a full column matrix ( with ) and therefore we define as left pseudoinverse of . We observe that the randomization should work better than the simple (cyclic) Kaczmarz algorithm for matrices with highly coherent rows (e.g., matrices generated by the computerized tomography). Since the Kaczmarz algorithm is based on projections, the convergence will be slow if the consecutive rows selected are highly coherent (i.e., the angle between and is small). Picking rows randomly (not necessarily based on the norms) makes picking more incoherent rows possible in each iteration. Therefore, the randomization may be useful for certain applications such as medical imaging. Note that matrix generated by computerized tomography has coherent and sparse rows due to physical nature of data collection. In fact, using Theorem 5, we can develop the following proposition.
Proposition 6. Let be a consistent linear system of equations () and let be an arbitrary initial approximation to the solution of . For , compute where is chosen from the set at random, with “any probability distribution.” Let be the solution of . Then, where is the scaled condition number of a matrix that is obtained by some row-scaling of .
Proof. This is due to the fact that row-scaling of (with scaling of the corresponding ) does not change the geometry of the problem and we can scale the rows to generate any probability distribution. In other words, we can obtain another matrix from by scaling its rows in such a way that picking the rows of based on the norms of the rows will be equivalent to picking the rows of based on the chosen probability distribution. Therefore, clearly, any randomization of the row selection will have exponential convergence; however, the rate will depend on the condition number of another matrix. For example, if we use uniform distribution, we can then normalize each row to have matrix as follows and then pick the rows at random with probability proportional to the norms of the rows:
2.8. Randomization Based on Subspace Angles Method
Our approach iterates through the rows of based on a probability distribution using the hyperplane (subspace) angles. Therefore, it is immune to scaling or normalization. This approach first generates a probability distribution based on the central angle(s) between the hyperplanes (represented by the rows of ). Then, it randomly picks two hyperplanes using this probability distribution. This is followed by a two-step projection on these hyperplanes (see Algorithm 2).
| ||||||||||||||||||||||||
2.9. -Subspaces Method
A new method has been developed which is intended to better accommodate the coherency of nonorthogonal data measurements. This next section makes contributions towards proving the statistical convergence of the randomized Kaczmarz orthogonal subspace (RKOS) algorithm. As described in [16], the RKOS initially uses -norm random hyperplane selection and subsequent projection into a constructed -dimensional orthogonal subspace comprised of an additional hyperplanes selected uniformly at random.
Algorithm 3 uses a recursive method to solve for the projections into the orthogonal subspace which is constructed using Gram-Schmidt (GS) procedure. However, a second approach demonstrates an alternate method of arriving at similar results, based upon an a closed form matrix for QR decomposition [17] of projection blocks.
| ||||||||||||||||||||||||||
In each of the above cases, vector operations inside the orthogonal subspace preserve the -norm and reduce errors that would normally be induced for coherent nonorthogonal projections which may be present in the simple Kaczmarz.
2.9.1. Orthogonal Subspaces
A statistical convergence analysis for Randomized Kaczmarz Orthogonal Subspace (RKOS) method is developed assuming identically and independently distributed (IID) random variables as vector components of each row of the measurement matrix .
(a) Orthogonal Construction. In many problems, and fast but optimal solutions are needed, often in noisy environments. In most cases, orthogonal data projection sampling is not feasible due to the constraints of the measurement system. The algorithm and procedure for the RKOS method are given in reference [16] and are intended to construct orthogonal measurements subspaces (see Algorithm 3).
The general technique is to solve using a constructed orthogonal basis from a full rank set of linearly independent measurements in for each subspace in Gram-Schmidt fashion [18, 19].
The subspace estimation may be computed as -dimensional subspace projection into the subspace orthonormal vector basis: where in subspace is the -dimensional solution approximation which becomes exact for for in the noiseless, self-consistent, case (the vector with the hat symbol indicates unit -norm).
(b) Modified Kaczmarz. The standard Kaczmarz equation is essentially iterative projections into a single subspace of dimension one; based upon the sampling hyperplanes, these projections are often oblique, especially in highly-coherent sampling.
The approach herein is motivated towards constructing an iterative algorithm based upon Kaczmarz which may be accelerated while controlling the potential projection errors and incurring reasonable computational penalty. The algorithm is simply to add subspaces of larger dimensions. Let It is convenient to make a substitution as follows: Using above substitution and orthonormal condition (It is worthwhile to note that in the problem setup, a fixed vector is projected into a randomized -dimensional subspace, where algebraic orthogonality was used to obtain (34). In the this statistical treatment of the same equation, the expectation of two random unit vectors vanishes for independent uncorrelated zero mean probability distribution functions, providing the statistical orthogonality on average satisfying (34).) , where the Kronecker find the -norm squared of : The ensemble average of the above equation (34) yields the convergence result, which is the main topic of this section.
2.9.2. Convergence for IID Measurement Matrix
Firstly, the expectation of a single random projection is computed. In the second step, the terms are summed for the -dimensional subspace. Experimental results are included in a latter section.
(a) Expectation of IID Projections. Consider the expectation of the -norm squared of the projection of fixed vector onto a random subspace basis of dimension : where the matrix basis is comprised of -columns of unit vectors in a constructed orthogonal basis for where the upper case components represent the th IID random variable component, and normalization constant is to be determined.
Further noting that complex conjugate reduces to transpose for real components, the -norm squared of the projection expands to In the next section, the goal is to find the expected value for outer product of the projection:
(b) Unit Vector. The deterministic identity for the magnitude of a unit vector is well known result for : The following statistical result must apply for the th column unit vector:
(c) Normalization of Random Unit Vector. Denote as the th random variable unit-norm vector associated with a set of column vectors comprising a random subspace matrix having IID random variable components . However, no additional assumptions on the distribution of the random variables are made at this time, other than IID. The expectation of both sides of (40) for random vector are found such that Solving above for each unit vector component in this treatment implies a random variable with zero mean and variance as follows: where is the associated IID probability distribution.
(d) P-Dimensional Random Projection. The next step is to compute the expectation of the magnitude of the projection of fixed vector onto random -dimensional orthonormal subspace projection term by term. Let be a column vector defined as and find the -norm squared: where Let upper case denote the th IID element random (this is not the same -variable as the Kaczmarz iteration variable) variable of the th column vector associated with column vector ; let vector denote a fixed point. Next, take the expectation of the term over the possible outcomes of random variables. Using the IID assumption, the expected value for a single projection component preserves terms squared as follows: It is now possible to determine the expectation for -terms of the projection as subject to IID constraint on where it is further noted that in (42).
(e) Error per Iteration. For a given th Kaczmarz iteration, the expectation of the projection of fixed vector onto the random -dimensional subspace is known from above. The total convergence expectation may then be computed, using a method similar to Strohmer’s, starting (recall that derivation of (47) requires orthogonality among the subspace basis vectors) with (47): We identify the term on the right as The results from the two equations ((49) and (48)) above may then be combined to obtain where the expectation on the right hand side includes accounting for the previous iteration. Next, apply induction to arrive at the expectation for the whole iterative sequence up to the th iteration given that :
(f) Asymptotic Convergence. The statistical ensemble average of the above equation (34) for the th iteration yields the convergence result given in (51). These results assume random variables identically and independently distributed but compare well to others in the literature, such as the convergence result in Strohmer and Vershynin [20].
The theoretical convergence iterative limit for uniform random IID sampling was compared to numerical simulations using random solution vector point on a unit sphere. Equation (52) has an asymptotic form: For comparison, recall the convergence for RK method of Strohmer for IID measurements with which is approximately Estimated noise bound convergence complexity to error is . Since the value of is given, the expectation is known to be the same.
(g) Theory and Simulation. Simulations in [16] compare theory to Gaussian IID with noise variance added to the measurements with magnitude (about five percent) and iteration termination at . In the first problem, the exact solution is chosen as a random point on the unit sphere—which is illustrated in Figure 1(a). In a second problem, a measurement of the standard phantom using parallel beam measurements is included, which contains coherent measurements.
(a) IID Gaussian Unit Vector Image
(b) CT Phantom Image
2.10. Regular Versus Randomized Kaczmarz Method
The randomized Kaczmarz’s algorithm developed by Strohmer and Vershynin in [4] has the following convergence in expectation: where is the scaled condition number of matrix with as the left pseudoinverse of . The bound for regular Kaczmarz is given in (25). Note that we assume . Now, we need to compare and to assess which bound is tighter. Let be ordered singular values of . Then, Also, note that where denotes the angles between the rows and of . Then, Note that therefore, Now, (54) and (25) become
2.11. Complexity Measurement Methods
In order to visualize and assess the relative performance of the well-known common methods of Kaczmarz, a simple routine was written in Octave [21] to record the total central processing unit (CPU) times for each method and variable matrix sizes of interest. The Linux (tm) kernel was modified to include the real-time patches for the x86/64 CPU architecture, and the process was run on a (see the cset in the cpuset package) shielded CPU set to single processor unit CPU7 to avoid process contention and interrupts. All system functions were locked to other CPU instances and not allowed to execute on CPU7.
Execution times were measured with rusage() calls before and after calls to each method in Octave. The elapsed time has microsecond resolution from the rusage() function. In addition, the interrupt balance kernel function was disabled. The status of the CPU cores and interrupts may be observed in /proc/interrupts. The CPU clock frequency was locked to a fixed value of 800 Mhz for the x86_64 system. Virtual memory was dropped prior to execution to ensure maximum physical memory availability.
In addition to the 64-bit platform, an embedded 32-bit ARM architecture was also configured for computational reference. Due to time constraints, no real-time kernel was implemented for the ARM. The CPU clock frequency was locked to a fixed value of 840 Mhz for the embedded ARM system.
In both of the above benchmark cases, all noncritical network and file system processes and daemons were stopped prior to code execution. In both cases, program execution was monitored and noted to reside in virtual (nonswaped) memory mapped to physical memory.
2.11.1. Methods
The kaczmarz(), randkaczmarz(), and cimmino() methods were called directly from the AIR toolkit [22]. The functions for the block method chenko() were implemented from Vasil’chenko and Svetlakov [23] of the form for .
The rkos() was implemented from [16] algorithm. The function for least squares ls() was computed from the equation The function for singular value decomposition svd() used the decomposition of the sampling matrix , and the solution was computed as follows: The resulting plots are intended to illustrate how the methods scale with increasing measurement matrix size and not directly used in absolute terms, since each method has various dependencies in hardware and software.
A simple baseline Kaczmarz code was written to ensure a baseline reference was available. The Octave source code for a typical method is included in the inset Listing 1.
| ||||||||||||||||||||||||||||||||||||||||||||||
2.11.2. Notes Regarding Simulations
The sizes of the matrix blocks equal matrix dimensions until constant block size of sixteen rows. The legend key is as follows:(1)least squares: LS(),(2)singular value decomposition: SVD(),(3)Kaczmarz(),(4)randomized Kaczmarz: RK(),(5)block method randomized Kaczmarz orthogonal subspaces: RKOS(),(6)Cimmino(),(7)block method Vasilchencko(),(8)reference P1: simple matrix multiply with known complexity,(9)reference P2: times P1 for known complexity,(10)reference Kaczmarz as shown in Listing 1.
3. Results and Discussion
Here, we compare our angle-based randomization with norm-based randomization of Strohmer and Vershynin [4] in the context of uniform and nonuniform measurement methods. In particular, Shepp-Logan and sinc2d() phantom images were used as the solutions in four different simulation experiments [3]. Figure 3(d) shows that our approach (angle-based randomization) provides a better convergence rate over the randomized Kaczmarz (norm-based randomization) in the case of fan-beam sampling of the sinc2d() phantom. However, our method is computationally more complex, and therefore we devised the -subspace algorithm (presented in the previous section).
In the experiments which follow, three factors of interest are the sampling angular distribution of the measurements (affects coherence), the algorithms’ method of iteration through the sampling hyperplanes, and the rate of convergence of the solution. The simulations were computed in parallel for each of the methods: Kaczmarz (K), randomized Kaczmarz hyperplane angles (RKHA), and randomized Kaczmarz (RK). Iteration is terminated at stopping point defined as the condition when at least one of the three methods attains a normalized error of 10% or less. Estimates for the signal to noise ratio are computed at same common stopping point.
3.1. Sampling Distributions
The following experiments compare Kaczmarz (K), randomized Kaczmarz (RK), and randomized Kaczmarz hyperplane angles (RKHA) via simulations. The objectives are to illustrate the effect of row randomization upon the convergence and observe the dependency upon the sampling methods.
3.1.1. Angular Distribution of Sampling Hyperplanes
A comparison of the distribution of hyperplane sampling angles in computed tomography (CT) was performed to investigate the convergence rate versus measurement strategy. Example results are presented for iterative convergence of methods K, RK, and RKHA under conditions of random, fan, and parallel beam sampling strategies using the Shepp-Logan phantom (see Figure 1(b)) [22], paralleltomo.m and fanbeamtomo.m from the AIRtools distribution [22] and randn() from the built-in function method [24]. The typical results for the angular distributions are provided in Figures 2(a) and 2(b) and discussed in more detail below.
(a) probability density distribution
(b) probability density distribution
(a) Random sampling of Shepp-Logan
(b) Random sampling of sinc2d()
(c) Fan sampling of Shepp-Logan
(d) Fan sampling of sinc2d()
3.1.2. Sampling Coherence
In linear algebra, the coherence or mutual coherence [25] of a row measurement matrix is defined as the maximum absolute value of the cross-correlations between the normalized rows of .
Formally, let be the set of row vectors of the matrix normalized such that where is the Hermitian conjugate and where . Let the mutual coherence of be defined as A lower bound was derived as in Welch [26].
In the distributions of Figures 2(a) and 2(b), it should be noted that the random sampling is concentrated near degrees probability but fan sampling is less concentrated across the interval degrees.
3.1.3. Observations on Effects of Sampling Distribution, Algorithm, and Convergence
Firstly, the convergence rates of K, RK, and RKHA are noted to be closely correlated for the case of random data sampling of the phantom in Figures 3(a) and 3(b). This is consistent with the mean values of coherence near zero for random sampling.
The cases for fan and parallel sampling have increasingly higher coherence and increasing numbers of iterations required to meet the 10% error stopping condition. Both test images for fan-beam sampling converge in example Figures 3(c) and 3(d) and show slight benefit from the methods which minimize the coherence, such as RK, RKHA, and RKOS. Example quantitative results for the two cases of fan sampling are shown in (a) Shepp-Logan Table 3 and Figure 4 and (b) sinc2d() Table 4 and Figure 5.
Comparison of convergence results to the estimated coherence for the three cases given in Table 5 suggest consistent interpretation. Comparison of percent error and SNR of Tables 1, 2, 3, and 4 indicate the randomization methods have slight advantage under coherent fan sampling, providing increased SNR and lower percent error.
3.1.4. Potential Applications
Since the iterative methods utilize projections, the angles between the optical lines-of-sight (LOS) forming the measurement hyperplanes are of considerable interest in terms of data acquisition and system design. Most methods of computed tomography do not reasonably allow random or orthogonal data sampling of the object of interest.
Therefore, these systems which acquire coherent data may benefit from use of the randomized methods RK, RKHA, or RKOS in data inversions. Typical examples of such systems may include computed tomography in medical and nonmedical X-ray, transmission ultrasound [27, 28], and resonance optical absorption and molecular florescence in-vivo imaging [29]. Each of the aforementioned systems are potentially feasible applications of RKHA or RKOS, since in each case, the measurements are path integrated along a mostly nonorthogonal set of electromagnetic LOS’s and generally require inversion to obtain the parameters of interest, such mole-fraction or species density.
3.2. Experimental Results for Convergence of K, RK, and RKHA
Iterative simulations were performed to estimate the relative convergence rates of methods K, RK, RKHA for the data examples above: random and fan beam sampling. Representative results are shown in Figures 3(a), 3(b), 3(c), and 3(d) for noiseless data measurement scenarios of the standard Shepp-Logan phantom and sinc2d() solutions under random and fan-beam sampling simulations. Parallel beam sampling and rect_sinc2d() solution vectors were also simulated. The results were consistent with observations reported herein but not included in this report.
3.2.1. Test Images for Algorithms
The Shepp Logan image was generated from the phantom code in reference [22], which is also the source for the fan-beam sampling algorithm. The sinc2d() function is defined as and was computed as the outer-product of two independent vectors constructed from the included code sinc.m in the signal processing package of Octave [21]. The rectangular rect_sinc2d.m function was used from source (https://engineering.purdue.edu/VISE/ee438/demos/2D_signals_systems/rect_sinc2d.m) which computes the two-dimensional sinc() function from the Fourier transform of the two-dimensional rectangle function.
3.2.2. Image Reconstruction Error
A comparison of the toy phantoms was performed based upon the estimated signal to noise ratio and total percent error versus iteration. Percent error and SNR were estimated for the Shepp-Logan phantom and the artifact based upon first method to obtain stopping error at percent normalized error.
The normalized error is defined as and estimate the signal to noise ratio (SNR) is estimated as with the vector set to the the zero vector.
3.3. Experimental Results for Complexity Estimates
Timings and plots of the methods were computed for a range of matrix sizes and plotted on log-log plots as shown in Figures 6 and 7, with measured run-times for 64-bit and 32-bit, respectively. For each method, the size of matrix was modified and average time to complete was measured. The data matrix was chosen to the be hadamard() function for ease of implementation and reference. The solution vector was chosen as a random point in the binomial distribution with using the binornd() function.
It is notable that the theoretical complexity of noise limited Kaczmarz is but the slope of the timings is closer to than . This is attributed to relatively small matrix sizes, and it is expected that the complexity asymptotically approaches ~2.0 for large .
The theoretical and numerical estimates for complexity are shown in Tables 6 and 7.
4. Conclusions
A new iterative selection rule based upon the relative central angle (RKHA) shows enhanced convergence in measurements which contain coherence. However, the method requires a computational penalty related to the dot-products of all to all rows, which may be overcome by a priori determination. A new block method using constructed orthogonal subspace projections provides enhanced tolerance to measurement coherence, but may be affected by noise at least as much as simple Kaczmarz. The exponential convergence is accelerated by the term and is computationally feasible for small relative to .
The theoretical convergence rates of above subspace methods were demonstrated using statistical IID assumptions or cyclical projections using the formalism of Galantai. Numerical results were presented from simulations of algorithm convergence under measurement distributions for fan and parallel X-ray beams, as well as random uniform sampling.
Relative performance benchmarks for complexity were obtained for typical hardware and software for various methods of well-known Kaczmarz algorithms. As expected, the Kaczmarz methods out-perform other methods, such as least-squares and singular value decomposition. The performance of embedded 32-bit ARM CPU architecture was sufficient to demonstrate functional capability over a range of low-power application environments such as mobile medicine platforms.
The new angle-based method (RKHA) and orthogonal block (RKOS) inversion method demonstrated herein showed quantitative convergence improvement consistent with increasing orthogonality and decreasing coherency of measurements. Future designs for tomography should consider optimization of angular sampling distributions in addition to other factors, such signal to noise ratio, as important system parameters, since these criteria ultimately affect the spatial-temporal resolution and uncertainty for a given number of samples per unit volume.
Conflict of Interests
The authors declare that there is no known conflict of interests regarding the publication of this paper.
Acknowledgments
The authors thank Akram Aldroubi and Alex Powell for their invaluable feedback. The research of Ali Sekmen is supported in part by NASA Grant NNX12AI14A.