# Kaczmarz Iterative Projection and Nonuniform Sampling with Complexity Estimates

Authors: Tim Wallace; Ali Sekmen
Journal: Journal of Medical Engineering
 (2014)

---

## Abstract

Kaczmarz’s alternating projection method has been widely used for solving mostly over-determined linear system of equationsA
x
=
b 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 A 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.


---

## Body

## 1. Introduction

Kaczmarz (in [1]) introduced an iterative algorithm for solving a consistent linear system of equations A
x
=
b with A
∈
R
M
×
N. This method projects the estimate x
j onto a subspace normal to the row a
i at step j
+
1 cyclically with i
=
j
(
mod
⁡
M
)
+
1. The block Kaczmarz algorithm first groups the rows into matrices A
1
,
A
2
,
…
,
A
k and then it projects the estimate x
j onto the subspace normal to the subspace spanned by the rows of A
i at step j
+
1 cyclically with i
=
j
(
mod
⁡
k
)
+
1. Obviously, the block Kaczmarz is equivalent to the simple Kaczmarz for k
=
M. 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 l
2 norm of that row, then the algorithm converges in expectation exponentially with a rate that depends on a scaled condition number of A (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 A
x
=
b. 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 A 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 theO
(
N
2
) 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

Letx
* be the solution of consistent A
x
=
b
, where A
∈
R
M
×
M is full column rank. Let A be row-partitioned as {
A
1
,
…
,
A
k
} where A
i
∈
R
M
i
×
M. Then, the simple block Kaczmarz update is as follows:
(1)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
b
i
-
A
i
x
j
i
=
j
mod
⁡
k
+
1
,
where b
i is the section of b that corresponds to the rows of A
i. Note that since A
i is full row rank, A
i
T
(
A
i
A
i
T
)
-
1 is the right pseudo-inverse of A
i. This is equivalent to
(2)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
A
i
x
*
-
A
i
x
j
,
x
j
+
1
-
x
*
=
x
j
-
x
*
-
A
i
T
A
i
A
i
T
-
1
A
i
(
x
j
-
x
*
)
.
Note that A
i
T
(
A
i
A
i
T
)
-
1
A
i is the projection matrix for projection of the range of A
i
T:
(3)
x
j
+
1
-
x
*
=
x
j
-
x
*
-
P
S
p
A
i
T
x
j
-
x
*
,
(4)
x
j
+
1
-
x
*
=
I
-
P
S
p
A
i
T
x
j
-
x
*
,
x
j
+
1
-
x
*
=
P
S
p
⊥
(
A
i
T
)
(
x
j
-
x
*
)
.
For one cycle of the blocks,
(5)
x
k
-
x
*
=
P
S
p
⊥
A
k
T
P
S
p
⊥
A
k
-
1
T
⋯
P
S
p
⊥
(
A
1
T
)
(
x
0
-
x
*
)
.
Note that if A
∈
R
M
×
N is a full column rank with M
<
N, then the simple block Kaczmarz update is as follows:
(6)
x
j
+
1
=
x
j
+
A
i
†
b
i
-
A
i
x
j
=
x
j
+
A
i
†
A
i
x
*
-
x
j
i
=
j
mod
⁡
k
+
1
,
where A
i
† is the pseudo-inverse of A
i and A
i
†
A
i is the orthogonal projection onto S
p
(
A
i
T
). Then, we get the same equation as (3), and subsequently we get (5),
(7)
x
j
+
1
-
x
*
=
x
j
-
x
*
-
P
S
p
A
i
T
x
j
-
x
*
.


### 2.2. Exponential Convergence

Theorem 1.
Letx
* be the solution of consistent A
x
=
b
, where A
∈
R
M
×
M is full rank. Let A be row-partitioned as {
A
1
,
…
,
A
k
} where A
i
∈
R
M
i
×
M and M
1
+
M
2
,
…
,
+
M
k
=
M. Then, the simple block Kaczmarz converges exponentially and the convergence rate depends of the number of blocks.Proof.
By (3) and orthogonal projection,
(8)
x
j
+
1
-
x
*
2
2
=
x
j
-
x
*
2
2
-
P
S
p
(
A
i
T
)
(
x
j
-
x
*
)
2
2
.
So,
(9)
x
j
+
1
-
x
*
2
2
≤
x
j
-
x
*
2
2
;
x
j
-
x
* depends on the initial condition x
~
0
=
x
0
-
x
*, and this dependence is scale-invariant. To see this, let e
j
=
x
j
-
x
* and consider c
x
~
0 where c
∈
R. By (4),
(10)
e
j
+
1
c
x
~
0
=
P
S
p
⊥
A
j
+
1
T
e
j
c
x
~
0
=
P
S
p
⊥
A
j
+
1
T
P
S
p
⊥
A
j
T
⋯
P
S
p
⊥
A
1
T
e
0
c
x
~
0
=
P
S
p
⊥
A
j
+
1
T
P
S
p
⊥
A
j
T
⋯
P
S
p
⊥
A
1
T
c
x
~
0
=
c
P
S
p
⊥
A
j
+
1
T
P
S
p
⊥
A
j
T
⋯
P
S
p
⊥
A
1
T
e
0
x
~
0
=
c
e
j
+
1
(
x
~
0
)
.
We will first show that if x
0
≠
x
*, then x
k
-
x
*
2
<
x
0
-
x
*
2. By the way of contradiction, assume that x
0
≠
x
* and x
k
-
x
*
2
=
x
0
-
x
*
2. By (9),
(11)
x
k
-
x
*
2
≤
x
k
-
1
-
x
*
2
⋯
<
x
0
-
x
*
2
and therefore x
l
-
x
*
2
=
x
0
-
x
*
2 for all 1
≤
l
≤
k. By (3), P
S
p
(
A
l
T
)
(
x
l
-
1
-
x
*
)
=
0 for all 1
≤
l
≤
k. By (8), we get x
l
=
x
0 for all 1
≤
l
≤
k. This implies that P
S
p
(
A
l
T
)
(
x
0
-
x
*
)
=
0 for all 1
≤
l
≤
k. So,
(12)
P
S
p
⊥
A
k
T
∩
S
p
⊥
A
k
T
⋯
∩
S
p
⊥
A
1
T
x
0
-
x
*
=
0
,
P
S
p
⊥
(
A
T
)
(
x
0
-
x
*
)
=
0
.
Since A is full column rank, we get x
0
=
x
*, which is a contradiction. So we know that x
k
-
x
*
2
<
x
0
-
x
*
2 (for one full cycle of k-iterations).
By compactness, there exists anϵ
∈
(
0,1
) such that, for all x
~
0
=
x
0
-
x
*
∈
S
N
-
1,
(13)
x
k
-
x
*
2
≤
1
-
ϵ
.
By (10) and (13),
(14)
x
k
-
x
*
2
=
x
~
0
2
e
k
x
~
0
x
~
0
2
≤
1
-
ϵ
x
~
0
2
,
x
k
-
x
*
2
≤
(
1
-
ϵ
)
x
0
-
x
*
2
.
Now consider iteration for q cycles:
(15)
x
q
k
-
x
*
2
≤
1
-
ϵ
q
x
0
-
x
*
2
,
x
q
k
-
x
*
2
≤
1
-
ϵ
1
/
k
q
k
x
0
-
x
*
2
.
Therefore, we conclude that the exponential decay depends on the number of blocks k. Note that k
=
M 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 E
x
p
-
x
*
2
2
≤
(
1
-
κ
(
A
)
-
2
)
p
x
0
-
x
*
2
2, where κ
(
A
)
=
A
F
A
†
2 is the scaled condition number of matrix A with A
† is the pseudoinverse of A.


### 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.
LetM
1
,
M
2
,
…
,
M
k be closed subspaces of the real Hilbert space H. Let M
=
∩
i
=
1
k
M
i and P
M
i
(
i
=
1
,
…
,
k
) be orthogonal projection on M
i. Then, for each x
∈
H,
(16)
lim
⁡
q
→
∞
P
M
k
P
M
k
-
1
⋯
P
M
1
q
x
=
P
M
x
,
where P
M is the orthogonal intersection projection.The block Kaczmarz is an alternating projection method withM
1
=
S
p
⊥
(
A
1
T
)
,
…
,
M
k
=
S
p
⊥
(
A
k
T
). Also,
(17)
P
M
1
=
P
S
p
⊥
(
A
1
T
)
,
…
,
P
M
k
=
S
p
⊥
A
k
T
,
M
=
S
p
⊥
A
1
T
∩
⋯
∩
S
p
⊥
(
A
k
T
)
=
S
p
⊥
(
A
T
)
.
Since A is full column rank, S
p
⊥
(
A
T
)
=
{
0
} and P
M
=
{
0
}. After q cycles,
(18)
x
q
k
-
x
*
=
P
M
k
P
M
k
-
1
⋯
P
M
1
q
x
0
-
x
*
.
By Theorem 2, lim
⁡
q
→
∞
x
q
k
-
x
*
=
0 and lim
⁡
q
→
∞
x
q
k
=
x
*. Galántai in [9] gives a bound for x
q
k
-
x
*
2 in terms of principle angles between M
i’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.
LetM
1
,
M
2
,
…
,
M
k be closed subspaces of the real Hilbert space H. Let M
=
∩
i
=
1
k
M
i and P
M
i
(
i
=
1
,
…
,
k
) be orthogonal projection on M
i (P
M is the orthogonal intersection projection). Let θ
j
=
α
(
M
j
,
∩
i
=
j
+
1
k
M
i
); then, for each x
∈
H and integer q
≥
1,
(19)
P
M
k
P
M
k
-
1
⋯
P
M
1
q
x
-
P
M
x
2
2
≤
1
-
Π
j
=
1
k
-
1
sin
2
θ
j
q
x
-
P
M
x
2
2
,
where P
M is the orthogonal intersection projection.In the special case of the block Kaczmarz, we haveH
=
R
N, M
1
=
S
p
⊥
(
A
1
T
)
,
…
,
M
k
=
S
p
⊥
(
A
k
T
). Also, P
M
1
=
P
S
p
⊥
(
A
1
T
)
,
…
,
P
M
k
=
P
S
p
⊥
(
A
k
T
) and M
=
S
p
⊥
(
A
1
T
)
∩
⋯
∩
S
p
⊥
(
A
k
T
)
=
S
p
⊥
(
A
T
). Since A is full column rank, S
p
⊥
(
A
T
)
=
{
0
} and P
M
=
{
0
}. Therefore, after q cycles,
(20)
x
q
k
-
x
*
2
2
=
P
M
k
P
M
k
-
1
⋯
P
M
1
q
(
x
0
-
x
*
)
2
2
≤
1
-
Π
j
=
1
k
-
1
sin
2
θ
j
q
x
o
-
x
*
2
2
,
where θ
j is as defined in Theorem 3. Note that the exponential decay rate depends on the number of blocks k as shown below:
(21)
x
q
k
-
x
*
2
2
≤
1
-
Π
j
=
1
k
-
1
sin
2
θ
j
1
/
k
q
k
x
o
-
x
*
2
2
.
Galántai in [9] developed another bound (for A
∈
R
M
×
M) by defining a new matrix X
i for each block A
i as follows.Theorem 4.
Letx
* be the solution of A
x
=
b for a consistent linear system with A
∈
R
M
×
M. Let A be row-partitioned as {
A
1
,
…
,
A
k
} where A
i
∈
R
M
i
×
N. Let M
1
=
S
p
⊥
(
A
1
T
)
,
…
,
M
k
=
S
p
⊥
(
A
k
T
) and A
i
A
i
T
=
L
L
T be the Cholesky decomposition of A
i
A
i
T. Define X
i
=
A
i
T
L
-
T and X
=
[
X
1
,
…
,
X
k
]. Then, for each x
∈
R
N and integer q
≥
1,
(22)
x
q
k
-
x
*
2
2
≤
1
-
det
⁡
X
T
X
q
x
o
-
x
*
2
2
=
1
-
det
⁡
X
T
X
1
/
k
q
k
x
o
-
x
*
2
2
.


### 2.5. Special Case: Simple Kaczmarz forA
∈
R
M
×
M

Note that this section assumes thatA
∈
R
M
×
M. The block Kaczmarz algorithm is equivalent to the simple Kaczmarz algorithm if the number of blocks k is equal to the number of rows M. In this case, A
i
A
i
T
=
a
i
2
2
=
L
L
T. Therefore, L
=
a
i
2 and L
-
T
=
1
/
a
i
2. This implies that X
i
=
a
i
/
a
i
2. Then, X
∈
R
M
×
M is defined as
(23)
X
=
a
1
a
1
2
,
…
,
a
M
a
M
2
.
Assume the matrix A has normalized rows and we pick a row at each iteration uniformly randomly. Note that this assumption is feasible as scaling a row of A and the corresponding measurement in b does not change the solution x.X is the Gram matrix with 0
≤
det
⁡
(
X
T
X
)
≤
x
1
2
2
x
2
2
2
⋯
x
M
2
2. Since x
i
2
=
1 and X is full rank, we have 0
<
det
⁡
(
X
T
X
)
≤
1. Using Theorem 4, we get the following deterministic bound:
(24)
x
q
M
-
x
*
2
2
≤
1
-
det
⁡
X
T
X
1
/
M
q
M
x
0
-
x
*
2
2
.
Since A is normalized, we get X
=
A
T and therefore
(25)
x
q
M
-
x
*
2
2
≤
1
-
det
⁡
A
A
T
1
/
M
q
M
x
0
-
x
*
2
2
.
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 Rowl
2 Norms Method

Strohmer and Vershynin (in [4]) developed a randomized Kaczmarz algorithm that picks a row of A in a random fashion with probability proportional with l
2 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 l
2 norms of the rows of A, 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 l
2 norm does not change the geometry of the problem. Theorem 5 is from [4].Theorem 5.
Letx
* be the solution of A
x
=
b Then, Algorithm 1 converges to x
* in expectation, with the average error
(26)
E
x
p
-
x
*
2
2
≤
1
-
κ
A
-
2
p
x
0
-
x
*
2
2
,
where κ
(
A
)
=
A
F
A
†
2 is the scaled condition number of matrix A with A
† is the left pseudoinverse of A.Algorithm 1:Randomized Kaczmarz (of [4]).
Require: An over-determined linear set of consistent equations A
x
=
b, where A is M
×
N matrix and b
∈
R
M.
Leta
1
,
…
,
a
M be the rows of A and b
j be the jth element of b.
(1) Pick an arbitrary initial approximation x
0.
(2) Set p
=
0.
(3) while not converged do
(4)
Randomly choose r
(
i
) from {
1
,
…
,
M
} with probability proportional to ‖
a
r
(
i
)
‖
2
2.
(5)
x
p
+
1
=
x
p
+
b
r
(
i
)
-
〈
a
r
(
i
)
,
x
p
〉
‖
a
r
(
i
)
‖
2
2
a
r
(
i
)
(6)
Set p
=
p
+
1
(7) end whileNote thatA is a full column matrix (A
∈
R
M
×
N with rank
⁡
(
A
)
=
N) and therefore we define A
† as left pseudoinverse of A. 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 a
i and a
i
+
1 is small). Picking rows randomly (not necessarily based on the l
2 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 A 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.
LetA
x
=
b be a consistent linear system of equations (A
∈
R
M
×
N) and let x
0 be an arbitrary initial approximation to the solution of A
x
=
b. For k
=
1,2
,
…, compute
(27)
x
p
+
1
=
x
p
+
b
r
(
i
)
-
a
r
i
,
x
p
a
r
(
i
)
2
2
a
r
(
i
)
,
where r
(
i
) is chosen from the set {
1,2
,
…
,
M
} at random, with “any probability distribution.” Let x
* be the solution of A
x
=
b. Then,
(28)
E
x
p
-
x
*
2
2
≤
1
-
κ
B
-
2
p
x
0
-
x
*
2
2
,
where κ
(
B
)
=
B
F
B
†
2 is the scaled condition number of a matrix B that is obtained by some row-scaling of A.Proof.
This is due to the fact that row-scaling ofA (with scaling of the corresponding b) 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 B from A by scaling its rows in such a way that picking the rows of B based on the l
2 norms of the rows will be equivalent to picking the rows of A 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 B as follows and then pick the rows at random with probability proportional to the norms of the rows:
(29)
B
=
a
1
a
1
2
,
…
,
a
M
a
M
2
T
.


### 2.8. Randomization Based on Subspace Angles Method

Our approach iterates through the rows ofA 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) θ
i
,
j between the hyperplanes (represented by the rows of A
x
=
b). Then, it randomly picks two hyperplanes using this probability distribution. This is followed by a two-step projection on these hyperplanes (see Algorithm 2).Algorithm 2:Randomized Kaczmarz Hyperplane Angles.
Require: An over-determined linear set of consistent equations A
x
=
b, where A is M
×
N matrix and b
∈
R
M.
Leta
1
,
…
,
a
M be the rows of A and b
j be the jth element of b.
(1)  Pick an arbitrary initial approximationx
0.
(2)  Setk
=
0.
(3)  Randomly choosef
(
i
) from {
1,2
,
…
,
M
} with a uniform distribution.
(4)while not converged do
(5)Randomly choose g
(
i
) from {
1
,
…
,
M
} with probability proportional to 1
-
(
〈
a
f
(
i
)
,
a
g
(
i
)
〉
2
/
‖
a
f
(
i
)
‖
2
2
‖
a
g
(
i
)
‖
2
2
)
=
1
-
cos
2
⁡
θ
f
(
i
)
,
g
(
i
)
(6)Compute x
k
+
1
=
x
k
+
(
(
b
f
(
i
)
-
〈
a
f
(
i
)
,
x
k
〉
)
/
‖
a
f
(
i
)
‖
2
2
)
a
f
(
i
)
(7)Compute x
k
+
2
=
x
k
+
1
+
(
(
b
g
(
i
)
-
〈
a
g
(
i
)
,
x
k
〉
)
/
‖
a
g
(
i
)
‖
2
2
)
a
g
(
i
)
(8)Set f
(
i
)
=
g
(
i
)
(9)Set k
=
k
+
2
(10)end while


### 2.9.P-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 l
2-norm random hyperplane selection and subsequent projection into a constructed P-dimensional orthogonal subspace S
P comprised of an additional P
-
1 hyperplanes selected uniformly at random.Algorithm3 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.Algorithm 3:
P-Subspace Kaczmarz Projections.
Require: Matrix A
∈
R
M
×
N full-rank consistent measurements subject to A
x
=
b, for b
∈
R
M.
(1) Set x
0 to initial approximation, i
=
1
(2) while not converged do
(3)   Select dim
⁡
(
S
P
)
=
P
<
N distinct linearly independent rows of A relative to random rule. Construct block
matrix A
i
∈
R
P
×
N comprised of rows {
a
i
,
1
,
…
,
a
i
,
P
}.
(4)   Perform Gram-Schmidt procedure on A
i to obtain the orthonormal set of columns {
u
i
,
1
,
…
,
u
i
,
P
}.
Let Q
i
=
{
u
i
,
1
,
…
,
u
i
,
P
}
∈
R
N
×
P
(5)   Update x
i as follows:
x
i
=
x
i
-
1
+
Proj
S
p
(
Q
i
)
(
x
i
-
1
-
x
),
x
i
=
x
i
-
1
-
Q
i
Q
i
T
(
x
-
x
i
-
1
),
(6)   Compute Q
i
T
x iteratively using {
a
i
,
1
,
…
,
a
i
,
P
}, {
b
i
,
1
,
…
,
b
i
,
P
}, {
u
i
,
1
,
…
,
u
i
,
P
}
(7)   Update i
=
i
+
1
(8) end whileIn each of the above cases, vector operations inside the orthogonal subspace preserve thel
2-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 matrixA.(a) Orthogonal Construction. In many problems, M
≫
N 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 asP-dimensional subspace projection into the subspace orthonormal vector basis:
(30)
x
S
P
=
∑
l
=
1
P
‍
u
^
l
,
x
u
^
l
,
where x
S
P in S
P
⊆
S
N subspace is the P-dimensional solution approximation which becomes exact for S
P
=
N for x
S
P
=
N
∈
R
N in the noiseless, self-consistent, case (the u vector with the hat symbol u
^ indicates unit l
2-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(31)
x
-
x
k
+
1
=
x
-
x
k
-
∑
l
=
1
P
‍
u
^
l
,
x
-
x
k
u
^
l
.
It is convenient to make a substitution as follows:
(32)
z
k
+
1
=
x
-
x
k
+
1
.
Using above substitution and orthonormal condition (It is worthwhile to note that in the problem setup, a fixed vector is projected into a randomized P-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).)  〈
u
^
j
,
u
^
k
〉
=
δ
j
,
k, where the Kronecker
(33)
δ
j
,
k
=
0
,
if
j
≠
k
,
1
,
if
j
=
k
,
find the l
2-norm squared of z
k
+
1:
(34)
z
k
+
1
2
2
=
z
k
2
2
-
∑
l
=
1
P
‍
u
^
l
,
z
k
2
.
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 theP-dimensional subspace. Experimental results are included in a latter section.(a) Expectation of IID Projections. Consider the expectation of the l
2-norm squared of the projection of fixed vector x
∈
R
N
×
1 onto a random subspace basis U
P
∈ of dimension P:
(35)
E
U
P
T
x
2
2
,
where the matrix basis U
P
∈
R
N
x
P is comprised of P-columns of unit vectors u
^
j
∈
R
N in a constructed orthogonal basis for
(36)
u
^
j
⟶
U
^
j
=
U
j
,
1
,
…
,
U
j
,
N
1
C
σ
=
U
j
U
j
2
2
∀
j
∈
1
,
…
,
P
,
where the upper case components U
j
,
i represent the (
j
,
i
)th IID random variable component, and normalization constant C
σ is to be determined.Further noting that complex conjugate(
·
)
* reduces to transpose (
·
)
T for real components, the l
2-norm squared of the projection expands to
(37)
U
P
T
x
2
2
=
x
T
U
P
U
P
T
x
.
In the next section, the goal is to find the expected value for outer product of the projection:
(38)
E
x
T
U
^
j
U
^
j
T
x
∀
j
∈
[
1
,
…
,
P
]
.(b) Unit Vector. The deterministic identity for the magnitude of a unit vector is well known result for u
^
∈
R
N:
(39)
u
^
2
2
=
∑
i
=
1
N
‍
u
i
2
u
2
2
=
1
.
The following statistical result must apply for the jth column unit vector:
(40)
E
U
^
j
2
2
=
E
U
^
j
T
U
^
j
=
1
=
E
U
j
,
1
2
+
⋯
+
U
j
,
N
2
1
C
σ
2
.(c) Normalization of Random Unit Vector. Denote U
^
j as the jth random variable unit-norm vector associated with a set of column vectors U
j
j
∈
1
,
…
,
P comprising a random subspace matrix U
N
×
P having IID random variable components U
j
,
i. 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 U
j are found such that
(41)
E
∑
i
=
1
N
‍
U
j
,
i
2
C
σ
2
=
∑
i
=
1
N
‍
E
U
j
,
i
2
C
σ
2
=
1
,
N
×
E
U
j
,
i
2
C
σ
2
=
1
.
Solving above for each unit vector component in this treatment implies a random variable U
j
,
i with zero mean and variance as follows:
(42)
E
U
j
,
i
2
=
σ
j
,
i
2
=
C
σ
2
N
∀
U
j
,
i
∈
1
,
…
,
N
∈
f
(
U
j
,
i
)
,
where f
(
U
i
,
j
) 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 x onto random P-dimensional orthonormal subspace U
P projection term by term. Let α
∈
R
P be a column vector defined as α
=
U
P
T
x and find the l
2-norm squared:
(43)
α
2
2
=
α
1
2
+
α
2
2
+
⋯
+
α
P
2
=
U
P
T
x
2
2
=
x
T
U
P
U
P
T
x
,
where
(44)
α
j
2
=
u
^
j
,
x
2
=
u
j
,
1
x
1
+
⋯
+
u
j
,
N
x
N
2
=
∑
i
,
k
N
,
N
‍
u
j
,
k
u
j
,
i
x
k
x
i
u
j
2
2
.
Let upper case U
j
,
k denote the kth IID element random (this is not the same k-variable as the Kaczmarz iteration variable) variable of the jth column vector U
j associated with column vector u
j; let x vector denote a fixed point. Next, take the expectation of the term over the possible outcomes of U
j
,
k random variables. Using the IID assumption, the expected value for a single projection component preserves terms squared as follows:
(45)
E
α
j
2
=
E
∑
i
,
k
N
,
N
‍
U
j
,
k
U
j
,
i
x
k
x
i
C
σ
2
=
∑
i
,
k
N
,
N
‍
E
U
j
,
k
U
j
,
i
x
k
x
i
C
σ
2
=
∑
k
=
1
N
‍
E
U
k
2
x
k
2
C
σ
2
=
∑
k
N
‍
E
U
j
,
k
2
C
σ
2
x
k
2
=
E
U
j
,
k
2
C
σ
2
∑
k
N
‍
x
k
2
=
E
U
j
,
k
2
C
σ
2
x
2
2
=
1
C
σ
2
C
σ
2
N
x
2
2
=
1
N
x
2
2
.
It is now possible to determine the expectation for P-terms of the projection as
(46)
E
α
2
2
=
E
∑
j
=
1
P
‍
α
j
2
=
P
N
x
2
2
subject to IID constraint on U
^
j where it is further noted that σ
2
N
=
C
σ
2 in (42).(e) Error per Iteration. For a given kth Kaczmarz iteration, the expectation of the projection of fixed vector x onto the random P-dimensional subspace U
P 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 u
^
l subspace basis vectors) with (47):
(47)
z
k
+
1
2
2
=
z
k
2
2
-
∑
l
=
1
P
‍
z
k
,
u
^
l
2
(48)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
+
1
2
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
-
∑
l
=
1
P
‍
z
k
,
u
^
l
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
-
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
∑
l
=
1
P
‍
z
k
,
u
^
l
2
.
We identify the term on the right as
(49)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
∑
l
=
1
P
‍
z
k
,
u
^
l
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
U
P
z
k
2
2
=
P
N
×
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
.
The results from the two equations ((49) and (48)) above may then be combined to obtain
(50)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
+
1
2
2
=
1
-
P
N
×
E
k
∣
z
0
,
z
1
,
…
,
z
k
-
1
z
k
2
2
,
where the expectation on the right hand side includes k
+
1
→
k 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 z
0
≡
x
-
x
0:
(51)
E
β
+
1
∣
z
0
z
β
+
1
2
2
=
1
-
P
N
β
z
0
2
2
∀
β
∈
1,2
,
3
,
…
.(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:
(52)
E
{
β
+
1
|
z
0
}
z
β
+
1
2
2
z
0
2
2
=
lim
⁡
β
→
∞
1
-
P
N
β
≃
e
-
β
P
N
,
P
=
dim
⁡
(
S
P
)
,
β
≫
1,2
,
3
,
…
⟶
k
∈
P
,
2
P
,
3
P
,
…
.
For comparison, recall the convergence for RK method of Strohmer for IID measurements with R
=
N which is approximately
(53)
E
{
k
+
1
∣
z
0
}
z
k
+
1
2
2
z
0
2
2
=
1
-
1
N
k
lim
⁡
k
→
∞
1
-
1
N
k
≃
e
-
k
/
N
∀
k
≫
1,2
,
3
,
…
.
Estimated noise bound convergence complexity to ϵ error is O
(
N
2
). Since the value of  z
0 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 β
=
0.05 (about five percent) and iteration termination at β
=
0.05
/
4
=
0.0125. In the first problem, the exact solution x 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.Representative test data.
(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:
(54)
E
x
q
M
-
x
*
2
2
≤
1
-
1
κ
A
2
q
M
x
0
-
x
*
2
2
,
where κ
(
A
)
=
A
F
A
†
2 is the scaled condition number of matrix A with A
† as the left pseudoinverse of A. The bound for regular Kaczmarz is given in (25). Note that we assume A
∈
R
M
×
M. Now, we need to compare (
1
-
(
1
/
A
F
2
A
†
2
2
)
) and (
1
-
det
⁡
(
A
A
T
)
)
1
/
M to assess which bound is tighter. Let σ
1
≥
σ
2
≥
⋯
≥
σ
M
>
0 be ordered singular values of A. Then,
(55)
A
†
2
2
=
1
σ
N
2
,
A
F
2
=
∑
i
=
1
M
‍
σ
i
2
.
Also, note that
(56)
A
A
T
=
1
cos
⁡
θ
12
⋯
cos
⁡
θ
1
M
cos
⁡
θ
21
1
⋯
cos
⁡
θ
2
M
⋮
⋮
⋯
⋮
cos
⁡
θ
M
1
cos
⁡
θ
M
2
⋯
1
,
where θ
i
j denotes the angles between the rows a
i and a
j of A. Then,
(57)
det
⁡
(
A
A
T
)
≤
∏
i
=
1
M
‍
∑
j
=
1
M
‍
cos
⁡
2
θ
i
j
.
Note that
(58)
∏
i
=
1
M
‍
σ
i
2
A
=
∏
i
=
1
M
‍
λ
i
A
T
A
=
det
⁡
A
T
A
=
det
⁡
A
A
T
;
therefore,
(59)
1
-
det
⁡
A
A
T
1
/
M
=
1
-
∏
i
=
1
M
‍
σ
i
2
1
/
M
.
Now, (54) and (25) become
(60)
E
x
q
M
-
x
*
2
2
≤
1
-
σ
M
2
∑
i
=
1
M
σ
i
2
q
M
x
0
-
x
*
2
2
,
x
q
M
-
x
*
2
2
≤
1
-
∏
i
=
1
M
‍
σ
i
2
1
/
M
q
M
x
0
-
x
*
2
2
.


### 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 thecset 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 withrusage() calls before and after calls to each method in Octave. The elapsed time has microsecond resolution from therusage() 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

Thekaczmarz(),randkaczmarz(), andcimmino() methods were called directly from the AIR toolkit [22]. The functions for the block methodchenko() were implemented from Vasil’chenko and Svetlakov [23] of the form
(61)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
b
i
-
A
i
x
j
i
=
j
(
mod
⁡
k
)
+
1
for A
i
=
E
i
T
A.Therkos() was implemented from [16] algorithm. The function for least squaresls() was computed from the equation
(62)
x
*
=
A
T
A
-
1
A
T
b
.
The function for singular value decompositionsvd() used the decomposition of the sampling matrix A
=
U
Σ
V
T, and the solution was computed as follows:
(63)
x
*
=
V
Σ
-
1
U
T
A
x
*
=
V
Σ
-
1
U
T
b
.
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 Listing1.Listing 1:Reference Kaczmarz Baseline.
(1)  function [X, info] = kazbase (A, b)
% note that matrix A is assumed to have normalized rows in below.
load  x_solution.data  x_sol
x_exact  =  x_sol
load  convergence limit.data
(6)
convergence_limit
myKmax  =  ceil(−log(convergence_limit))
ntimes  =  myKmax
[M,  N]  =  size  (A)
jj  =  0
(11)
xk0  =  zeros  (M,  1);
x0  =  xk0;
for  n  =  1:  ntimes
for  kn  =  1: M
jj  =  jj  +  1;
(16)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delta_xk_plus_10  =  (eye  (N)  −  A(kn,  :)′  ∗ A(kn,  :))∗(x_exact −  xk0);
xk0 = x_exact  − delta_xk_plus_10;
%error0  (jj)=norm(delta_xk_plus_10,  2)/norm(x_exact  −  x0,  2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(21)end
end
X  =  xk0


#### 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 methodVasilchencko(),
(8)
reference P1: simple matrix multiply with knownO
(
N
2
) complexity,
(9)
reference P2:N times P1 for known O
(
N
3
) complexity,
(10)
reference Kaczmarz as shown in Listing1.


## 2.1. Convergence of Regular Block Kaczmarz Method

Letx
* be the solution of consistent A
x
=
b
, where A
∈
R
M
×
M is full column rank. Let A be row-partitioned as {
A
1
,
…
,
A
k
} where A
i
∈
R
M
i
×
M. Then, the simple block Kaczmarz update is as follows:
(1)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
b
i
-
A
i
x
j
i
=
j
mod
⁡
k
+
1
,
where b
i is the section of b that corresponds to the rows of A
i. Note that since A
i is full row rank, A
i
T
(
A
i
A
i
T
)
-
1 is the right pseudo-inverse of A
i. This is equivalent to
(2)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
A
i
x
*
-
A
i
x
j
,
x
j
+
1
-
x
*
=
x
j
-
x
*
-
A
i
T
A
i
A
i
T
-
1
A
i
(
x
j
-
x
*
)
.
Note that A
i
T
(
A
i
A
i
T
)
-
1
A
i is the projection matrix for projection of the range of A
i
T:
(3)
x
j
+
1
-
x
*
=
x
j
-
x
*
-
P
S
p
A
i
T
x
j
-
x
*
,
(4)
x
j
+
1
-
x
*
=
I
-
P
S
p
A
i
T
x
j
-
x
*
,
x
j
+
1
-
x
*
=
P
S
p
⊥
(
A
i
T
)
(
x
j
-
x
*
)
.
For one cycle of the blocks,
(5)
x
k
-
x
*
=
P
S
p
⊥
A
k
T
P
S
p
⊥
A
k
-
1
T
⋯
P
S
p
⊥
(
A
1
T
)
(
x
0
-
x
*
)
.
Note that if A
∈
R
M
×
N is a full column rank with M
<
N, then the simple block Kaczmarz update is as follows:
(6)
x
j
+
1
=
x
j
+
A
i
†
b
i
-
A
i
x
j
=
x
j
+
A
i
†
A
i
x
*
-
x
j
i
=
j
mod
⁡
k
+
1
,
where A
i
† is the pseudo-inverse of A
i and A
i
†
A
i is the orthogonal projection onto S
p
(
A
i
T
). Then, we get the same equation as (3), and subsequently we get (5),
(7)
x
j
+
1
-
x
*
=
x
j
-
x
*
-
P
S
p
A
i
T
x
j
-
x
*
.


## 2.2. Exponential Convergence

Theorem 1.
Letx
* be the solution of consistent A
x
=
b
, where A
∈
R
M
×
M is full rank. Let A be row-partitioned as {
A
1
,
…
,
A
k
} where A
i
∈
R
M
i
×
M and M
1
+
M
2
,
…
,
+
M
k
=
M. Then, the simple block Kaczmarz converges exponentially and the convergence rate depends of the number of blocks.Proof.
By (3) and orthogonal projection,
(8)
x
j
+
1
-
x
*
2
2
=
x
j
-
x
*
2
2
-
P
S
p
(
A
i
T
)
(
x
j
-
x
*
)
2
2
.
So,
(9)
x
j
+
1
-
x
*
2
2
≤
x
j
-
x
*
2
2
;
x
j
-
x
* depends on the initial condition x
~
0
=
x
0
-
x
*, and this dependence is scale-invariant. To see this, let e
j
=
x
j
-
x
* and consider c
x
~
0 where c
∈
R. By (4),
(10)
e
j
+
1
c
x
~
0
=
P
S
p
⊥
A
j
+
1
T
e
j
c
x
~
0
=
P
S
p
⊥
A
j
+
1
T
P
S
p
⊥
A
j
T
⋯
P
S
p
⊥
A
1
T
e
0
c
x
~
0
=
P
S
p
⊥
A
j
+
1
T
P
S
p
⊥
A
j
T
⋯
P
S
p
⊥
A
1
T
c
x
~
0
=
c
P
S
p
⊥
A
j
+
1
T
P
S
p
⊥
A
j
T
⋯
P
S
p
⊥
A
1
T
e
0
x
~
0
=
c
e
j
+
1
(
x
~
0
)
.
We will first show that if x
0
≠
x
*, then x
k
-
x
*
2
<
x
0
-
x
*
2. By the way of contradiction, assume that x
0
≠
x
* and x
k
-
x
*
2
=
x
0
-
x
*
2. By (9),
(11)
x
k
-
x
*
2
≤
x
k
-
1
-
x
*
2
⋯
<
x
0
-
x
*
2
and therefore x
l
-
x
*
2
=
x
0
-
x
*
2 for all 1
≤
l
≤
k. By (3), P
S
p
(
A
l
T
)
(
x
l
-
1
-
x
*
)
=
0 for all 1
≤
l
≤
k. By (8), we get x
l
=
x
0 for all 1
≤
l
≤
k. This implies that P
S
p
(
A
l
T
)
(
x
0
-
x
*
)
=
0 for all 1
≤
l
≤
k. So,
(12)
P
S
p
⊥
A
k
T
∩
S
p
⊥
A
k
T
⋯
∩
S
p
⊥
A
1
T
x
0
-
x
*
=
0
,
P
S
p
⊥
(
A
T
)
(
x
0
-
x
*
)
=
0
.
Since A is full column rank, we get x
0
=
x
*, which is a contradiction. So we know that x
k
-
x
*
2
<
x
0
-
x
*
2 (for one full cycle of k-iterations).
By compactness, there exists anϵ
∈
(
0,1
) such that, for all x
~
0
=
x
0
-
x
*
∈
S
N
-
1,
(13)
x
k
-
x
*
2
≤
1
-
ϵ
.
By (10) and (13),
(14)
x
k
-
x
*
2
=
x
~
0
2
e
k
x
~
0
x
~
0
2
≤
1
-
ϵ
x
~
0
2
,
x
k
-
x
*
2
≤
(
1
-
ϵ
)
x
0
-
x
*
2
.
Now consider iteration for q cycles:
(15)
x
q
k
-
x
*
2
≤
1
-
ϵ
q
x
0
-
x
*
2
,
x
q
k
-
x
*
2
≤
1
-
ϵ
1
/
k
q
k
x
0
-
x
*
2
.
Therefore, we conclude that the exponential decay depends on the number of blocks k. Note that k
=
M 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 E
x
p
-
x
*
2
2
≤
(
1
-
κ
(
A
)
-
2
)
p
x
0
-
x
*
2
2, where κ
(
A
)
=
A
F
A
†
2 is the scaled condition number of matrix A with A
† is the pseudoinverse of A.


## 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.
LetM
1
,
M
2
,
…
,
M
k be closed subspaces of the real Hilbert space H. Let M
=
∩
i
=
1
k
M
i and P
M
i
(
i
=
1
,
…
,
k
) be orthogonal projection on M
i. Then, for each x
∈
H,
(16)
lim
⁡
q
→
∞
P
M
k
P
M
k
-
1
⋯
P
M
1
q
x
=
P
M
x
,
where P
M is the orthogonal intersection projection.The block Kaczmarz is an alternating projection method withM
1
=
S
p
⊥
(
A
1
T
)
,
…
,
M
k
=
S
p
⊥
(
A
k
T
). Also,
(17)
P
M
1
=
P
S
p
⊥
(
A
1
T
)
,
…
,
P
M
k
=
S
p
⊥
A
k
T
,
M
=
S
p
⊥
A
1
T
∩
⋯
∩
S
p
⊥
(
A
k
T
)
=
S
p
⊥
(
A
T
)
.
Since A is full column rank, S
p
⊥
(
A
T
)
=
{
0
} and P
M
=
{
0
}. After q cycles,
(18)
x
q
k
-
x
*
=
P
M
k
P
M
k
-
1
⋯
P
M
1
q
x
0
-
x
*
.
By Theorem 2, lim
⁡
q
→
∞
x
q
k
-
x
*
=
0 and lim
⁡
q
→
∞
x
q
k
=
x
*. Galántai in [9] gives a bound for x
q
k
-
x
*
2 in terms of principle angles between M
i’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.
LetM
1
,
M
2
,
…
,
M
k be closed subspaces of the real Hilbert space H. Let M
=
∩
i
=
1
k
M
i and P
M
i
(
i
=
1
,
…
,
k
) be orthogonal projection on M
i (P
M is the orthogonal intersection projection). Let θ
j
=
α
(
M
j
,
∩
i
=
j
+
1
k
M
i
); then, for each x
∈
H and integer q
≥
1,
(19)
P
M
k
P
M
k
-
1
⋯
P
M
1
q
x
-
P
M
x
2
2
≤
1
-
Π
j
=
1
k
-
1
sin
2
θ
j
q
x
-
P
M
x
2
2
,
where P
M is the orthogonal intersection projection.In the special case of the block Kaczmarz, we haveH
=
R
N, M
1
=
S
p
⊥
(
A
1
T
)
,
…
,
M
k
=
S
p
⊥
(
A
k
T
). Also, P
M
1
=
P
S
p
⊥
(
A
1
T
)
,
…
,
P
M
k
=
P
S
p
⊥
(
A
k
T
) and M
=
S
p
⊥
(
A
1
T
)
∩
⋯
∩
S
p
⊥
(
A
k
T
)
=
S
p
⊥
(
A
T
). Since A is full column rank, S
p
⊥
(
A
T
)
=
{
0
} and P
M
=
{
0
}. Therefore, after q cycles,
(20)
x
q
k
-
x
*
2
2
=
P
M
k
P
M
k
-
1
⋯
P
M
1
q
(
x
0
-
x
*
)
2
2
≤
1
-
Π
j
=
1
k
-
1
sin
2
θ
j
q
x
o
-
x
*
2
2
,
where θ
j is as defined in Theorem 3. Note that the exponential decay rate depends on the number of blocks k as shown below:
(21)
x
q
k
-
x
*
2
2
≤
1
-
Π
j
=
1
k
-
1
sin
2
θ
j
1
/
k
q
k
x
o
-
x
*
2
2
.
Galántai in [9] developed another bound (for A
∈
R
M
×
M) by defining a new matrix X
i for each block A
i as follows.Theorem 4.
Letx
* be the solution of A
x
=
b for a consistent linear system with A
∈
R
M
×
M. Let A be row-partitioned as {
A
1
,
…
,
A
k
} where A
i
∈
R
M
i
×
N. Let M
1
=
S
p
⊥
(
A
1
T
)
,
…
,
M
k
=
S
p
⊥
(
A
k
T
) and A
i
A
i
T
=
L
L
T be the Cholesky decomposition of A
i
A
i
T. Define X
i
=
A
i
T
L
-
T and X
=
[
X
1
,
…
,
X
k
]. Then, for each x
∈
R
N and integer q
≥
1,
(22)
x
q
k
-
x
*
2
2
≤
1
-
det
⁡
X
T
X
q
x
o
-
x
*
2
2
=
1
-
det
⁡
X
T
X
1
/
k
q
k
x
o
-
x
*
2
2
.


## 2.5. Special Case: Simple Kaczmarz forA
∈
R
M
×
M

Note that this section assumes thatA
∈
R
M
×
M. The block Kaczmarz algorithm is equivalent to the simple Kaczmarz algorithm if the number of blocks k is equal to the number of rows M. In this case, A
i
A
i
T
=
a
i
2
2
=
L
L
T. Therefore, L
=
a
i
2 and L
-
T
=
1
/
a
i
2. This implies that X
i
=
a
i
/
a
i
2. Then, X
∈
R
M
×
M is defined as
(23)
X
=
a
1
a
1
2
,
…
,
a
M
a
M
2
.
Assume the matrix A has normalized rows and we pick a row at each iteration uniformly randomly. Note that this assumption is feasible as scaling a row of A and the corresponding measurement in b does not change the solution x.X is the Gram matrix with 0
≤
det
⁡
(
X
T
X
)
≤
x
1
2
2
x
2
2
2
⋯
x
M
2
2. Since x
i
2
=
1 and X is full rank, we have 0
<
det
⁡
(
X
T
X
)
≤
1. Using Theorem 4, we get the following deterministic bound:
(24)
x
q
M
-
x
*
2
2
≤
1
-
det
⁡
X
T
X
1
/
M
q
M
x
0
-
x
*
2
2
.
Since A is normalized, we get X
=
A
T and therefore
(25)
x
q
M
-
x
*
2
2
≤
1
-
det
⁡
A
A
T
1
/
M
q
M
x
0
-
x
*
2
2
.
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 Rowl
2 Norms Method

Strohmer and Vershynin (in [4]) developed a randomized Kaczmarz algorithm that picks a row of A in a random fashion with probability proportional with l
2 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 l
2 norms of the rows of A, 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 l
2 norm does not change the geometry of the problem. Theorem 5 is from [4].Theorem 5.
Letx
* be the solution of A
x
=
b Then, Algorithm 1 converges to x
* in expectation, with the average error
(26)
E
x
p
-
x
*
2
2
≤
1
-
κ
A
-
2
p
x
0
-
x
*
2
2
,
where κ
(
A
)
=
A
F
A
†
2 is the scaled condition number of matrix A with A
† is the left pseudoinverse of A.Algorithm 1:Randomized Kaczmarz (of [4]).
Require: An over-determined linear set of consistent equations A
x
=
b, where A is M
×
N matrix and b
∈
R
M.
Leta
1
,
…
,
a
M be the rows of A and b
j be the jth element of b.
(1) Pick an arbitrary initial approximation x
0.
(2) Set p
=
0.
(3) while not converged do
(4)
Randomly choose r
(
i
) from {
1
,
…
,
M
} with probability proportional to ‖
a
r
(
i
)
‖
2
2.
(5)
x
p
+
1
=
x
p
+
b
r
(
i
)
-
〈
a
r
(
i
)
,
x
p
〉
‖
a
r
(
i
)
‖
2
2
a
r
(
i
)
(6)
Set p
=
p
+
1
(7) end whileNote thatA is a full column matrix (A
∈
R
M
×
N with rank
⁡
(
A
)
=
N) and therefore we define A
† as left pseudoinverse of A. 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 a
i and a
i
+
1 is small). Picking rows randomly (not necessarily based on the l
2 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 A 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.
LetA
x
=
b be a consistent linear system of equations (A
∈
R
M
×
N) and let x
0 be an arbitrary initial approximation to the solution of A
x
=
b. For k
=
1,2
,
…, compute
(27)
x
p
+
1
=
x
p
+
b
r
(
i
)
-
a
r
i
,
x
p
a
r
(
i
)
2
2
a
r
(
i
)
,
where r
(
i
) is chosen from the set {
1,2
,
…
,
M
} at random, with “any probability distribution.” Let x
* be the solution of A
x
=
b. Then,
(28)
E
x
p
-
x
*
2
2
≤
1
-
κ
B
-
2
p
x
0
-
x
*
2
2
,
where κ
(
B
)
=
B
F
B
†
2 is the scaled condition number of a matrix B that is obtained by some row-scaling of A.Proof.
This is due to the fact that row-scaling ofA (with scaling of the corresponding b) 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 B from A by scaling its rows in such a way that picking the rows of B based on the l
2 norms of the rows will be equivalent to picking the rows of A 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 B as follows and then pick the rows at random with probability proportional to the norms of the rows:
(29)
B
=
a
1
a
1
2
,
…
,
a
M
a
M
2
T
.


## 2.8. Randomization Based on Subspace Angles Method

Our approach iterates through the rows ofA 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) θ
i
,
j between the hyperplanes (represented by the rows of A
x
=
b). Then, it randomly picks two hyperplanes using this probability distribution. This is followed by a two-step projection on these hyperplanes (see Algorithm 2).Algorithm 2:Randomized Kaczmarz Hyperplane Angles.
Require: An over-determined linear set of consistent equations A
x
=
b, where A is M
×
N matrix and b
∈
R
M.
Leta
1
,
…
,
a
M be the rows of A and b
j be the jth element of b.
(1)  Pick an arbitrary initial approximationx
0.
(2)  Setk
=
0.
(3)  Randomly choosef
(
i
) from {
1,2
,
…
,
M
} with a uniform distribution.
(4)while not converged do
(5)Randomly choose g
(
i
) from {
1
,
…
,
M
} with probability proportional to 1
-
(
〈
a
f
(
i
)
,
a
g
(
i
)
〉
2
/
‖
a
f
(
i
)
‖
2
2
‖
a
g
(
i
)
‖
2
2
)
=
1
-
cos
2
⁡
θ
f
(
i
)
,
g
(
i
)
(6)Compute x
k
+
1
=
x
k
+
(
(
b
f
(
i
)
-
〈
a
f
(
i
)
,
x
k
〉
)
/
‖
a
f
(
i
)
‖
2
2
)
a
f
(
i
)
(7)Compute x
k
+
2
=
x
k
+
1
+
(
(
b
g
(
i
)
-
〈
a
g
(
i
)
,
x
k
〉
)
/
‖
a
g
(
i
)
‖
2
2
)
a
g
(
i
)
(8)Set f
(
i
)
=
g
(
i
)
(9)Set k
=
k
+
2
(10)end while


## 2.9.P-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 l
2-norm random hyperplane selection and subsequent projection into a constructed P-dimensional orthogonal subspace S
P comprised of an additional P
-
1 hyperplanes selected uniformly at random.Algorithm3 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.Algorithm 3:
P-Subspace Kaczmarz Projections.
Require: Matrix A
∈
R
M
×
N full-rank consistent measurements subject to A
x
=
b, for b
∈
R
M.
(1) Set x
0 to initial approximation, i
=
1
(2) while not converged do
(3)   Select dim
⁡
(
S
P
)
=
P
<
N distinct linearly independent rows of A relative to random rule. Construct block
matrix A
i
∈
R
P
×
N comprised of rows {
a
i
,
1
,
…
,
a
i
,
P
}.
(4)   Perform Gram-Schmidt procedure on A
i to obtain the orthonormal set of columns {
u
i
,
1
,
…
,
u
i
,
P
}.
Let Q
i
=
{
u
i
,
1
,
…
,
u
i
,
P
}
∈
R
N
×
P
(5)   Update x
i as follows:
x
i
=
x
i
-
1
+
Proj
S
p
(
Q
i
)
(
x
i
-
1
-
x
),
x
i
=
x
i
-
1
-
Q
i
Q
i
T
(
x
-
x
i
-
1
),
(6)   Compute Q
i
T
x iteratively using {
a
i
,
1
,
…
,
a
i
,
P
}, {
b
i
,
1
,
…
,
b
i
,
P
}, {
u
i
,
1
,
…
,
u
i
,
P
}
(7)   Update i
=
i
+
1
(8) end whileIn each of the above cases, vector operations inside the orthogonal subspace preserve thel
2-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 matrixA.(a) Orthogonal Construction. In many problems, M
≫
N 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 asP-dimensional subspace projection into the subspace orthonormal vector basis:
(30)
x
S
P
=
∑
l
=
1
P
‍
u
^
l
,
x
u
^
l
,
where x
S
P in S
P
⊆
S
N subspace is the P-dimensional solution approximation which becomes exact for S
P
=
N for x
S
P
=
N
∈
R
N in the noiseless, self-consistent, case (the u vector with the hat symbol u
^ indicates unit l
2-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(31)
x
-
x
k
+
1
=
x
-
x
k
-
∑
l
=
1
P
‍
u
^
l
,
x
-
x
k
u
^
l
.
It is convenient to make a substitution as follows:
(32)
z
k
+
1
=
x
-
x
k
+
1
.
Using above substitution and orthonormal condition (It is worthwhile to note that in the problem setup, a fixed vector is projected into a randomized P-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).)  〈
u
^
j
,
u
^
k
〉
=
δ
j
,
k, where the Kronecker
(33)
δ
j
,
k
=
0
,
if
j
≠
k
,
1
,
if
j
=
k
,
find the l
2-norm squared of z
k
+
1:
(34)
z
k
+
1
2
2
=
z
k
2
2
-
∑
l
=
1
P
‍
u
^
l
,
z
k
2
.
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 theP-dimensional subspace. Experimental results are included in a latter section.(a) Expectation of IID Projections. Consider the expectation of the l
2-norm squared of the projection of fixed vector x
∈
R
N
×
1 onto a random subspace basis U
P
∈ of dimension P:
(35)
E
U
P
T
x
2
2
,
where the matrix basis U
P
∈
R
N
x
P is comprised of P-columns of unit vectors u
^
j
∈
R
N in a constructed orthogonal basis for
(36)
u
^
j
⟶
U
^
j
=
U
j
,
1
,
…
,
U
j
,
N
1
C
σ
=
U
j
U
j
2
2
∀
j
∈
1
,
…
,
P
,
where the upper case components U
j
,
i represent the (
j
,
i
)th IID random variable component, and normalization constant C
σ is to be determined.Further noting that complex conjugate(
·
)
* reduces to transpose (
·
)
T for real components, the l
2-norm squared of the projection expands to
(37)
U
P
T
x
2
2
=
x
T
U
P
U
P
T
x
.
In the next section, the goal is to find the expected value for outer product of the projection:
(38)
E
x
T
U
^
j
U
^
j
T
x
∀
j
∈
[
1
,
…
,
P
]
.(b) Unit Vector. The deterministic identity for the magnitude of a unit vector is well known result for u
^
∈
R
N:
(39)
u
^
2
2
=
∑
i
=
1
N
‍
u
i
2
u
2
2
=
1
.
The following statistical result must apply for the jth column unit vector:
(40)
E
U
^
j
2
2
=
E
U
^
j
T
U
^
j
=
1
=
E
U
j
,
1
2
+
⋯
+
U
j
,
N
2
1
C
σ
2
.(c) Normalization of Random Unit Vector. Denote U
^
j as the jth random variable unit-norm vector associated with a set of column vectors U
j
j
∈
1
,
…
,
P comprising a random subspace matrix U
N
×
P having IID random variable components U
j
,
i. 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 U
j are found such that
(41)
E
∑
i
=
1
N
‍
U
j
,
i
2
C
σ
2
=
∑
i
=
1
N
‍
E
U
j
,
i
2
C
σ
2
=
1
,
N
×
E
U
j
,
i
2
C
σ
2
=
1
.
Solving above for each unit vector component in this treatment implies a random variable U
j
,
i with zero mean and variance as follows:
(42)
E
U
j
,
i
2
=
σ
j
,
i
2
=
C
σ
2
N
∀
U
j
,
i
∈
1
,
…
,
N
∈
f
(
U
j
,
i
)
,
where f
(
U
i
,
j
) 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 x onto random P-dimensional orthonormal subspace U
P projection term by term. Let α
∈
R
P be a column vector defined as α
=
U
P
T
x and find the l
2-norm squared:
(43)
α
2
2
=
α
1
2
+
α
2
2
+
⋯
+
α
P
2
=
U
P
T
x
2
2
=
x
T
U
P
U
P
T
x
,
where
(44)
α
j
2
=
u
^
j
,
x
2
=
u
j
,
1
x
1
+
⋯
+
u
j
,
N
x
N
2
=
∑
i
,
k
N
,
N
‍
u
j
,
k
u
j
,
i
x
k
x
i
u
j
2
2
.
Let upper case U
j
,
k denote the kth IID element random (this is not the same k-variable as the Kaczmarz iteration variable) variable of the jth column vector U
j associated with column vector u
j; let x vector denote a fixed point. Next, take the expectation of the term over the possible outcomes of U
j
,
k random variables. Using the IID assumption, the expected value for a single projection component preserves terms squared as follows:
(45)
E
α
j
2
=
E
∑
i
,
k
N
,
N
‍
U
j
,
k
U
j
,
i
x
k
x
i
C
σ
2
=
∑
i
,
k
N
,
N
‍
E
U
j
,
k
U
j
,
i
x
k
x
i
C
σ
2
=
∑
k
=
1
N
‍
E
U
k
2
x
k
2
C
σ
2
=
∑
k
N
‍
E
U
j
,
k
2
C
σ
2
x
k
2
=
E
U
j
,
k
2
C
σ
2
∑
k
N
‍
x
k
2
=
E
U
j
,
k
2
C
σ
2
x
2
2
=
1
C
σ
2
C
σ
2
N
x
2
2
=
1
N
x
2
2
.
It is now possible to determine the expectation for P-terms of the projection as
(46)
E
α
2
2
=
E
∑
j
=
1
P
‍
α
j
2
=
P
N
x
2
2
subject to IID constraint on U
^
j where it is further noted that σ
2
N
=
C
σ
2 in (42).(e) Error per Iteration. For a given kth Kaczmarz iteration, the expectation of the projection of fixed vector x onto the random P-dimensional subspace U
P 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 u
^
l subspace basis vectors) with (47):
(47)
z
k
+
1
2
2
=
z
k
2
2
-
∑
l
=
1
P
‍
z
k
,
u
^
l
2
(48)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
+
1
2
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
-
∑
l
=
1
P
‍
z
k
,
u
^
l
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
-
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
∑
l
=
1
P
‍
z
k
,
u
^
l
2
.
We identify the term on the right as
(49)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
∑
l
=
1
P
‍
z
k
,
u
^
l
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
U
P
z
k
2
2
=
P
N
×
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
.
The results from the two equations ((49) and (48)) above may then be combined to obtain
(50)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
+
1
2
2
=
1
-
P
N
×
E
k
∣
z
0
,
z
1
,
…
,
z
k
-
1
z
k
2
2
,
where the expectation on the right hand side includes k
+
1
→
k 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 z
0
≡
x
-
x
0:
(51)
E
β
+
1
∣
z
0
z
β
+
1
2
2
=
1
-
P
N
β
z
0
2
2
∀
β
∈
1,2
,
3
,
…
.(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:
(52)
E
{
β
+
1
|
z
0
}
z
β
+
1
2
2
z
0
2
2
=
lim
⁡
β
→
∞
1
-
P
N
β
≃
e
-
β
P
N
,
P
=
dim
⁡
(
S
P
)
,
β
≫
1,2
,
3
,
…
⟶
k
∈
P
,
2
P
,
3
P
,
…
.
For comparison, recall the convergence for RK method of Strohmer for IID measurements with R
=
N which is approximately
(53)
E
{
k
+
1
∣
z
0
}
z
k
+
1
2
2
z
0
2
2
=
1
-
1
N
k
lim
⁡
k
→
∞
1
-
1
N
k
≃
e
-
k
/
N
∀
k
≫
1,2
,
3
,
…
.
Estimated noise bound convergence complexity to ϵ error is O
(
N
2
). Since the value of  z
0 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 β
=
0.05 (about five percent) and iteration termination at β
=
0.05
/
4
=
0.0125. In the first problem, the exact solution x 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.Representative test data.
(a)
IID Gaussian Unit Vector Image
(b)
CT Phantom Image


## 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 matrixA.(a) Orthogonal Construction. In many problems, M
≫
N 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 asP-dimensional subspace projection into the subspace orthonormal vector basis:
(30)
x
S
P
=
∑
l
=
1
P
‍
u
^
l
,
x
u
^
l
,
where x
S
P in S
P
⊆
S
N subspace is the P-dimensional solution approximation which becomes exact for S
P
=
N for x
S
P
=
N
∈
R
N in the noiseless, self-consistent, case (the u vector with the hat symbol u
^ indicates unit l
2-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(31)
x
-
x
k
+
1
=
x
-
x
k
-
∑
l
=
1
P
‍
u
^
l
,
x
-
x
k
u
^
l
.
It is convenient to make a substitution as follows:
(32)
z
k
+
1
=
x
-
x
k
+
1
.
Using above substitution and orthonormal condition (It is worthwhile to note that in the problem setup, a fixed vector is projected into a randomized P-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).)  〈
u
^
j
,
u
^
k
〉
=
δ
j
,
k, where the Kronecker
(33)
δ
j
,
k
=
0
,
if
j
≠
k
,
1
,
if
j
=
k
,
find the l
2-norm squared of z
k
+
1:
(34)
z
k
+
1
2
2
=
z
k
2
2
-
∑
l
=
1
P
‍
u
^
l
,
z
k
2
.
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 theP-dimensional subspace. Experimental results are included in a latter section.(a) Expectation of IID Projections. Consider the expectation of the l
2-norm squared of the projection of fixed vector x
∈
R
N
×
1 onto a random subspace basis U
P
∈ of dimension P:
(35)
E
U
P
T
x
2
2
,
where the matrix basis U
P
∈
R
N
x
P is comprised of P-columns of unit vectors u
^
j
∈
R
N in a constructed orthogonal basis for
(36)
u
^
j
⟶
U
^
j
=
U
j
,
1
,
…
,
U
j
,
N
1
C
σ
=
U
j
U
j
2
2
∀
j
∈
1
,
…
,
P
,
where the upper case components U
j
,
i represent the (
j
,
i
)th IID random variable component, and normalization constant C
σ is to be determined.Further noting that complex conjugate(
·
)
* reduces to transpose (
·
)
T for real components, the l
2-norm squared of the projection expands to
(37)
U
P
T
x
2
2
=
x
T
U
P
U
P
T
x
.
In the next section, the goal is to find the expected value for outer product of the projection:
(38)
E
x
T
U
^
j
U
^
j
T
x
∀
j
∈
[
1
,
…
,
P
]
.(b) Unit Vector. The deterministic identity for the magnitude of a unit vector is well known result for u
^
∈
R
N:
(39)
u
^
2
2
=
∑
i
=
1
N
‍
u
i
2
u
2
2
=
1
.
The following statistical result must apply for the jth column unit vector:
(40)
E
U
^
j
2
2
=
E
U
^
j
T
U
^
j
=
1
=
E
U
j
,
1
2
+
⋯
+
U
j
,
N
2
1
C
σ
2
.(c) Normalization of Random Unit Vector. Denote U
^
j as the jth random variable unit-norm vector associated with a set of column vectors U
j
j
∈
1
,
…
,
P comprising a random subspace matrix U
N
×
P having IID random variable components U
j
,
i. 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 U
j are found such that
(41)
E
∑
i
=
1
N
‍
U
j
,
i
2
C
σ
2
=
∑
i
=
1
N
‍
E
U
j
,
i
2
C
σ
2
=
1
,
N
×
E
U
j
,
i
2
C
σ
2
=
1
.
Solving above for each unit vector component in this treatment implies a random variable U
j
,
i with zero mean and variance as follows:
(42)
E
U
j
,
i
2
=
σ
j
,
i
2
=
C
σ
2
N
∀
U
j
,
i
∈
1
,
…
,
N
∈
f
(
U
j
,
i
)
,
where f
(
U
i
,
j
) 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 x onto random P-dimensional orthonormal subspace U
P projection term by term. Let α
∈
R
P be a column vector defined as α
=
U
P
T
x and find the l
2-norm squared:
(43)
α
2
2
=
α
1
2
+
α
2
2
+
⋯
+
α
P
2
=
U
P
T
x
2
2
=
x
T
U
P
U
P
T
x
,
where
(44)
α
j
2
=
u
^
j
,
x
2
=
u
j
,
1
x
1
+
⋯
+
u
j
,
N
x
N
2
=
∑
i
,
k
N
,
N
‍
u
j
,
k
u
j
,
i
x
k
x
i
u
j
2
2
.
Let upper case U
j
,
k denote the kth IID element random (this is not the same k-variable as the Kaczmarz iteration variable) variable of the jth column vector U
j associated with column vector u
j; let x vector denote a fixed point. Next, take the expectation of the term over the possible outcomes of U
j
,
k random variables. Using the IID assumption, the expected value for a single projection component preserves terms squared as follows:
(45)
E
α
j
2
=
E
∑
i
,
k
N
,
N
‍
U
j
,
k
U
j
,
i
x
k
x
i
C
σ
2
=
∑
i
,
k
N
,
N
‍
E
U
j
,
k
U
j
,
i
x
k
x
i
C
σ
2
=
∑
k
=
1
N
‍
E
U
k
2
x
k
2
C
σ
2
=
∑
k
N
‍
E
U
j
,
k
2
C
σ
2
x
k
2
=
E
U
j
,
k
2
C
σ
2
∑
k
N
‍
x
k
2
=
E
U
j
,
k
2
C
σ
2
x
2
2
=
1
C
σ
2
C
σ
2
N
x
2
2
=
1
N
x
2
2
.
It is now possible to determine the expectation for P-terms of the projection as
(46)
E
α
2
2
=
E
∑
j
=
1
P
‍
α
j
2
=
P
N
x
2
2
subject to IID constraint on U
^
j where it is further noted that σ
2
N
=
C
σ
2 in (42).(e) Error per Iteration. For a given kth Kaczmarz iteration, the expectation of the projection of fixed vector x onto the random P-dimensional subspace U
P 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 u
^
l subspace basis vectors) with (47):
(47)
z
k
+
1
2
2
=
z
k
2
2
-
∑
l
=
1
P
‍
z
k
,
u
^
l
2
(48)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
+
1
2
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
-
∑
l
=
1
P
‍
z
k
,
u
^
l
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
-
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
∑
l
=
1
P
‍
z
k
,
u
^
l
2
.
We identify the term on the right as
(49)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
∑
l
=
1
P
‍
z
k
,
u
^
l
2
=
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
U
P
z
k
2
2
=
P
N
×
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
2
2
.
The results from the two equations ((49) and (48)) above may then be combined to obtain
(50)
E
k
+
1
∣
z
0
,
z
1
,
…
,
z
k
z
k
+
1
2
2
=
1
-
P
N
×
E
k
∣
z
0
,
z
1
,
…
,
z
k
-
1
z
k
2
2
,
where the expectation on the right hand side includes k
+
1
→
k 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 z
0
≡
x
-
x
0:
(51)
E
β
+
1
∣
z
0
z
β
+
1
2
2
=
1
-
P
N
β
z
0
2
2
∀
β
∈
1,2
,
3
,
…
.(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:
(52)
E
{
β
+
1
|
z
0
}
z
β
+
1
2
2
z
0
2
2
=
lim
⁡
β
→
∞
1
-
P
N
β
≃
e
-
β
P
N
,
P
=
dim
⁡
(
S
P
)
,
β
≫
1,2
,
3
,
…
⟶
k
∈
P
,
2
P
,
3
P
,
…
.
For comparison, recall the convergence for RK method of Strohmer for IID measurements with R
=
N which is approximately
(53)
E
{
k
+
1
∣
z
0
}
z
k
+
1
2
2
z
0
2
2
=
1
-
1
N
k
lim
⁡
k
→
∞
1
-
1
N
k
≃
e
-
k
/
N
∀
k
≫
1,2
,
3
,
…
.
Estimated noise bound convergence complexity to ϵ error is O
(
N
2
). Since the value of  z
0 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 β
=
0.05 (about five percent) and iteration termination at β
=
0.05
/
4
=
0.0125. In the first problem, the exact solution x 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.Representative test data.
(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:
(54)
E
x
q
M
-
x
*
2
2
≤
1
-
1
κ
A
2
q
M
x
0
-
x
*
2
2
,
where κ
(
A
)
=
A
F
A
†
2 is the scaled condition number of matrix A with A
† as the left pseudoinverse of A. The bound for regular Kaczmarz is given in (25). Note that we assume A
∈
R
M
×
M. Now, we need to compare (
1
-
(
1
/
A
F
2
A
†
2
2
)
) and (
1
-
det
⁡
(
A
A
T
)
)
1
/
M to assess which bound is tighter. Let σ
1
≥
σ
2
≥
⋯
≥
σ
M
>
0 be ordered singular values of A. Then,
(55)
A
†
2
2
=
1
σ
N
2
,
A
F
2
=
∑
i
=
1
M
‍
σ
i
2
.
Also, note that
(56)
A
A
T
=
1
cos
⁡
θ
12
⋯
cos
⁡
θ
1
M
cos
⁡
θ
21
1
⋯
cos
⁡
θ
2
M
⋮
⋮
⋯
⋮
cos
⁡
θ
M
1
cos
⁡
θ
M
2
⋯
1
,
where θ
i
j denotes the angles between the rows a
i and a
j of A. Then,
(57)
det
⁡
(
A
A
T
)
≤
∏
i
=
1
M
‍
∑
j
=
1
M
‍
cos
⁡
2
θ
i
j
.
Note that
(58)
∏
i
=
1
M
‍
σ
i
2
A
=
∏
i
=
1
M
‍
λ
i
A
T
A
=
det
⁡
A
T
A
=
det
⁡
A
A
T
;
therefore,
(59)
1
-
det
⁡
A
A
T
1
/
M
=
1
-
∏
i
=
1
M
‍
σ
i
2
1
/
M
.
Now, (54) and (25) become
(60)
E
x
q
M
-
x
*
2
2
≤
1
-
σ
M
2
∑
i
=
1
M
σ
i
2
q
M
x
0
-
x
*
2
2
,
x
q
M
-
x
*
2
2
≤
1
-
∏
i
=
1
M
‍
σ
i
2
1
/
M
q
M
x
0
-
x
*
2
2
.


## 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 thecset 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 withrusage() calls before and after calls to each method in Octave. The elapsed time has microsecond resolution from therusage() 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

Thekaczmarz(),randkaczmarz(), andcimmino() methods were called directly from the AIR toolkit [22]. The functions for the block methodchenko() were implemented from Vasil’chenko and Svetlakov [23] of the form
(61)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
b
i
-
A
i
x
j
i
=
j
(
mod
⁡
k
)
+
1
for A
i
=
E
i
T
A.Therkos() was implemented from [16] algorithm. The function for least squaresls() was computed from the equation
(62)
x
*
=
A
T
A
-
1
A
T
b
.
The function for singular value decompositionsvd() used the decomposition of the sampling matrix A
=
U
Σ
V
T, and the solution was computed as follows:
(63)
x
*
=
V
Σ
-
1
U
T
A
x
*
=
V
Σ
-
1
U
T
b
.
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 Listing1.Listing 1:Reference Kaczmarz Baseline.
(1)  function [X, info] = kazbase (A, b)
% note that matrix A is assumed to have normalized rows in below.
load  x_solution.data  x_sol
x_exact  =  x_sol
load  convergence limit.data
(6)
convergence_limit
myKmax  =  ceil(−log(convergence_limit))
ntimes  =  myKmax
[M,  N]  =  size  (A)
jj  =  0
(11)
xk0  =  zeros  (M,  1);
x0  =  xk0;
for  n  =  1:  ntimes
for  kn  =  1: M
jj  =  jj  +  1;
(16)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delta_xk_plus_10  =  (eye  (N)  −  A(kn,  :)′  ∗ A(kn,  :))∗(x_exact −  xk0);
xk0 = x_exact  − delta_xk_plus_10;
%error0  (jj)=norm(delta_xk_plus_10,  2)/norm(x_exact  −  x0,  2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(21)end
end
X  =  xk0


### 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 methodVasilchencko(),
(8)
reference P1: simple matrix multiply with knownO
(
N
2
) complexity,
(9)
reference P2:N times P1 for known O
(
N
3
) complexity,
(10)
reference Kaczmarz as shown in Listing1.


## 2.11.1. Methods

Thekaczmarz(),randkaczmarz(), andcimmino() methods were called directly from the AIR toolkit [22]. The functions for the block methodchenko() were implemented from Vasil’chenko and Svetlakov [23] of the form
(61)
x
j
+
1
=
x
j
+
A
i
T
A
i
A
i
T
-
1
b
i
-
A
i
x
j
i
=
j
(
mod
⁡
k
)
+
1
for A
i
=
E
i
T
A.Therkos() was implemented from [16] algorithm. The function for least squaresls() was computed from the equation
(62)
x
*
=
A
T
A
-
1
A
T
b
.
The function for singular value decompositionsvd() used the decomposition of the sampling matrix A
=
U
Σ
V
T, and the solution was computed as follows:
(63)
x
*
=
V
Σ
-
1
U
T
A
x
*
=
V
Σ
-
1
U
T
b
.
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 Listing1.Listing 1:Reference Kaczmarz Baseline.
(1)  function [X, info] = kazbase (A, b)
% note that matrix A is assumed to have normalized rows in below.
load  x_solution.data  x_sol
x_exact  =  x_sol
load  convergence limit.data
(6)
convergence_limit
myKmax  =  ceil(−log(convergence_limit))
ntimes  =  myKmax
[M,  N]  =  size  (A)
jj  =  0
(11)
xk0  =  zeros  (M,  1);
x0  =  xk0;
for  n  =  1:  ntimes
for  kn  =  1: M
jj  =  jj  +  1;
(16)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
delta_xk_plus_10  =  (eye  (N)  −  A(kn,  :)′  ∗ A(kn,  :))∗(x_exact −  xk0);
xk0 = x_exact  − delta_xk_plus_10;
%error0  (jj)=norm(delta_xk_plus_10,  2)/norm(x_exact  −  x0,  2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
(21)end
end
X  =  xk0


## 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 methodVasilchencko(),
(8)
reference P1: simple matrix multiply with knownO
(
N
2
) complexity,
(9)
reference P2:N times P1 for known O
(
N
3
) complexity,
(10)
reference Kaczmarz as shown in Listing1.


## 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 andsinc2d() 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 thesinc2d() phantom. However, our method is computationally more complex, and therefore we devised the P-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 Figure1(b)) [22],paralleltomo.m andfanbeamtomo.m from the AIRtools distribution [22] andrandn() 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.Example angles distribution (y-axis) from A
A
T where θ
i
,
j
=
cos
⁡
-
1
(
〈
a
^
i
,
a
^
j
〉
)  ∀
i
,
j
∈
{
1
,
…
,
M
} versus angles (x-axis) degrees using (a) random and (b) fan-beam data acquisition strategy where the unit norm vectors a
^
i, a
^
j are selected rows of matrix A.
(a)
θ
(
i
,
j
) probability density distribution
(b)
θ
(
i
,
j
) probability density distributionSemilog (y-axis) plot example of normalized convergence error results for K, RK, and RKHA on Shepp-Logan phantom using random and fan tomographic data acquisition sampling and stopping at 10 percent error.
(a)
Random sampling of Shepp-Logan
(b)
Random sampling ofsinc2d()
(c)
Fan sampling of Shepp-Logan
(d)
Fan sampling ofsinc2d()


#### 3.1.2. Sampling Coherence

In linear algebra, the coherence or mutual coherence [25] of a row measurement matrix A is defined as the maximum absolute value of the cross-correlations between the normalized rows of A.Formally, leta
1
,
…
,
a
M
∈
R
N be the set of row vectors of the matrix A
∈
R
M
×
N normalized such that 〈
a
i
,
a
i
〉
=
a
i
H
a
i
=
1 where (
·
)
H is the Hermitian conjugate and where M
>
N. Let the mutual coherence of A be defined as
(64)
ϕ
i
,
j
=
max
⁡
1
≤
i
≠
j
≤
M
a
i
H
a
j
.
A lower bound was derived as ϕ
≥
(
M
-
N
)
/
N
(
M
-
1
) in Welch [26].In the distributions of Figures2(a) and 2(b), it should be noted that the random sampling is concentrated near 90 degrees probability but fan sampling is less concentrated across the interval [
0,90
] 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 Figures3(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 Figures3(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.Figure 4
Reconstructed test images for Shepp-Logan with fan sampling with error from left to right: exact, K (9.9% error), RKHA, and RK Shepp Logan.Figure 5
Reconstructed test images forsinc2d() with fan sampling with error from left to right: exact, K, RKHA (9.9%), and RK.Comparison of convergence results to the estimated coherence for the three cases given in Table5 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.Table 1
Random sampling of Shepp-Logan SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.1
5.60
4.5
% error
9.9
19.9
22.1Table 2
Random sampling ofsinc2d() SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.0
5.6
5.2
% error
9.9
17.6
19.1Table 3
Fan beam sampling of Shepp Logan, SNR, and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.0
9.95
9.5
% error
9.9
10.0
10.5Table 4
Fan beam sampling ofsinc2d() SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
7.1
10.0
9.2
% error
14.0
9.9
10.7Table 5
Typical coherence estimates forN
=
100, M
=
200 for random randn() and N
=
100, M
=
222 for fan fanbeamtomo() and parallel paralleltomo().
Coherence versus measurement method
Random
Fan
Parallel
Coherence (64)
.4
1.0
1.0
Average value ofG
i
,
j
=
〈
a
^
i
,
a
^
j
〉
(
1
≤
i
≠
j
≤
M
)
−.0013
.06
.18
Median value ofG
i
,
j
=
〈
a
^
i
,
a
^
j
〉
(
1
≤
i
≠
j
≤
M
)
−0.0009
0
.12


#### 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 florescencein-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 Figures3(a), 3(b), 3(c), and 3(d) for noiseless data measurement scenarios of the standard Shepp-Logan phantom andsinc2d() solutions under random and fan-beam sampling simulations. Parallel beam sampling andrect_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. Thesinc2d() function is defined as
(65)
s
i
n
c
2
d
(
x
,
y
)
=
sin
(
π
x
)
π
x
sin
(
π
y
)
π
y
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-dimensionalsinc() 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 at10
% percent normalized error.The normalized error is defined as(66)
ϵ
k
=
x
-
x
k
2
x
-
x
0
2
,
and estimate the signal to noise ratio (SNR) is estimated as 1
/
ϵ
k with the x
0 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 Figures6 and 7, with measured run-times for 64-bit and 32-bit, respectively. For each method, the size of A matrix was modified and average time to complete was measured. The data matrix was chosen to the behadamard() function for ease of implementation and reference. The solution vector was chosen as a random point in the binomial distribution with p
=
0.5 using thebinornd() function.Figure 6
Relative computational time using dedicated custom real-time Linux kernel 3.12.5 patched versus square matrix dimensionsN
=
2
q  ∀
q
∈
[
4
,
…
,
10
]; CPU is shielded and frequency locked for single process and single core (one of eight) 800 MHz Intel 64-bit i7.Figure 7
Relative computational time for embedded microprocessor dedicated task custom Linux kernel 3.12.7 ARMv7 relative computational time versus matrix dimensionsN
=
2
q  ∀
q
∈
[
4
,
…
,
10
]; CPU ARM frequency locked kernel only process stack 840 MHz 32-bit.It is notable that the theoretical complexity of noise limited Kaczmarz isO
(
N
2
) but the slope of the timings is closer to 1.1 than 2.0. This is attributed to relatively small matrix sizes, and it is expected that the complexity asymptotically approaches ~2.0 for large N.The theoretical and numerical estimates for complexity are shown in Tables6 and 7.Table 6
Numerical estimates for computational complexity using naive codes.
Computational complexity
LS
SVD
Kaczmarz
RK
RKOS
Theoretical
O
(
N
3
)
O
(
N
3
)
O
(
N
2
)
O
(
N
2
)
O
(
N
2
)
Numerical estimates
2.7
3.1
1.1
1.1
1.1Table 7
Numerical estimates for computational complexity using naive codes.
Computational complexity
Cimmino
Vasilchenko
Ref.P1
Ref.P2
Ref. Kaczmarz
Theoretical
O
(
N
2
)
O
(
N
2
)
O
(
N
2
)
O
(
N
3
)
O
(
N
2
)
Numerical estimates
1.6
2.1
2
3.2
1.1


## 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 Figure1(b)) [22],paralleltomo.m andfanbeamtomo.m from the AIRtools distribution [22] andrandn() 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.Example angles distribution (y-axis) from A
A
T where θ
i
,
j
=
cos
⁡
-
1
(
〈
a
^
i
,
a
^
j
〉
)  ∀
i
,
j
∈
{
1
,
…
,
M
} versus angles (x-axis) degrees using (a) random and (b) fan-beam data acquisition strategy where the unit norm vectors a
^
i, a
^
j are selected rows of matrix A.
(a)
θ
(
i
,
j
) probability density distribution
(b)
θ
(
i
,
j
) probability density distributionSemilog (y-axis) plot example of normalized convergence error results for K, RK, and RKHA on Shepp-Logan phantom using random and fan tomographic data acquisition sampling and stopping at 10 percent error.
(a)
Random sampling of Shepp-Logan
(b)
Random sampling ofsinc2d()
(c)
Fan sampling of Shepp-Logan
(d)
Fan sampling ofsinc2d()


### 3.1.2. Sampling Coherence

In linear algebra, the coherence or mutual coherence [25] of a row measurement matrix A is defined as the maximum absolute value of the cross-correlations between the normalized rows of A.Formally, leta
1
,
…
,
a
M
∈
R
N be the set of row vectors of the matrix A
∈
R
M
×
N normalized such that 〈
a
i
,
a
i
〉
=
a
i
H
a
i
=
1 where (
·
)
H is the Hermitian conjugate and where M
>
N. Let the mutual coherence of A be defined as
(64)
ϕ
i
,
j
=
max
⁡
1
≤
i
≠
j
≤
M
a
i
H
a
j
.
A lower bound was derived as ϕ
≥
(
M
-
N
)
/
N
(
M
-
1
) in Welch [26].In the distributions of Figures2(a) and 2(b), it should be noted that the random sampling is concentrated near 90 degrees probability but fan sampling is less concentrated across the interval [
0,90
] 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 Figures3(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 Figures3(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.Figure 4
Reconstructed test images for Shepp-Logan with fan sampling with error from left to right: exact, K (9.9% error), RKHA, and RK Shepp Logan.Figure 5
Reconstructed test images forsinc2d() with fan sampling with error from left to right: exact, K, RKHA (9.9%), and RK.Comparison of convergence results to the estimated coherence for the three cases given in Table5 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.Table 1
Random sampling of Shepp-Logan SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.1
5.60
4.5
% error
9.9
19.9
22.1Table 2
Random sampling ofsinc2d() SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.0
5.6
5.2
% error
9.9
17.6
19.1Table 3
Fan beam sampling of Shepp Logan, SNR, and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.0
9.95
9.5
% error
9.9
10.0
10.5Table 4
Fan beam sampling ofsinc2d() SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
7.1
10.0
9.2
% error
14.0
9.9
10.7Table 5
Typical coherence estimates forN
=
100, M
=
200 for random randn() and N
=
100, M
=
222 for fan fanbeamtomo() and parallel paralleltomo().
Coherence versus measurement method
Random
Fan
Parallel
Coherence (64)
.4
1.0
1.0
Average value ofG
i
,
j
=
〈
a
^
i
,
a
^
j
〉
(
1
≤
i
≠
j
≤
M
)
−.0013
.06
.18
Median value ofG
i
,
j
=
〈
a
^
i
,
a
^
j
〉
(
1
≤
i
≠
j
≤
M
)
−0.0009
0
.12


### 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 florescencein-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.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 Figure1(b)) [22],paralleltomo.m andfanbeamtomo.m from the AIRtools distribution [22] andrandn() 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.Example angles distribution (y-axis) from A
A
T where θ
i
,
j
=
cos
⁡
-
1
(
〈
a
^
i
,
a
^
j
〉
)  ∀
i
,
j
∈
{
1
,
…
,
M
} versus angles (x-axis) degrees using (a) random and (b) fan-beam data acquisition strategy where the unit norm vectors a
^
i, a
^
j are selected rows of matrix A.
(a)
θ
(
i
,
j
) probability density distribution
(b)
θ
(
i
,
j
) probability density distributionSemilog (y-axis) plot example of normalized convergence error results for K, RK, and RKHA on Shepp-Logan phantom using random and fan tomographic data acquisition sampling and stopping at 10 percent error.
(a)
Random sampling of Shepp-Logan
(b)
Random sampling ofsinc2d()
(c)
Fan sampling of Shepp-Logan
(d)
Fan sampling ofsinc2d()


## 3.1.2. Sampling Coherence

In linear algebra, the coherence or mutual coherence [25] of a row measurement matrix A is defined as the maximum absolute value of the cross-correlations between the normalized rows of A.Formally, leta
1
,
…
,
a
M
∈
R
N be the set of row vectors of the matrix A
∈
R
M
×
N normalized such that 〈
a
i
,
a
i
〉
=
a
i
H
a
i
=
1 where (
·
)
H is the Hermitian conjugate and where M
>
N. Let the mutual coherence of A be defined as
(64)
ϕ
i
,
j
=
max
⁡
1
≤
i
≠
j
≤
M
a
i
H
a
j
.
A lower bound was derived as ϕ
≥
(
M
-
N
)
/
N
(
M
-
1
) in Welch [26].In the distributions of Figures2(a) and 2(b), it should be noted that the random sampling is concentrated near 90 degrees probability but fan sampling is less concentrated across the interval [
0,90
] 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 Figures3(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 Figures3(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.Figure 4
Reconstructed test images for Shepp-Logan with fan sampling with error from left to right: exact, K (9.9% error), RKHA, and RK Shepp Logan.Figure 5
Reconstructed test images forsinc2d() with fan sampling with error from left to right: exact, K, RKHA (9.9%), and RK.Comparison of convergence results to the estimated coherence for the three cases given in Table5 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.Table 1
Random sampling of Shepp-Logan SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.1
5.60
4.5
% error
9.9
19.9
22.1Table 2
Random sampling ofsinc2d() SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.0
5.6
5.2
% error
9.9
17.6
19.1Table 3
Fan beam sampling of Shepp Logan, SNR, and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
10.0
9.95
9.5
% error
9.9
10.0
10.5Table 4
Fan beam sampling ofsinc2d() SNR and percent error at stopping.
Performance estimates
Kaczmarz
RKHA
RK
SNR
7.1
10.0
9.2
% error
14.0
9.9
10.7Table 5
Typical coherence estimates forN
=
100, M
=
200 for random randn() and N
=
100, M
=
222 for fan fanbeamtomo() and parallel paralleltomo().
Coherence versus measurement method
Random
Fan
Parallel
Coherence (64)
.4
1.0
1.0
Average value ofG
i
,
j
=
〈
a
^
i
,
a
^
j
〉
(
1
≤
i
≠
j
≤
M
)
−.0013
.06
.18
Median value ofG
i
,
j
=
〈
a
^
i
,
a
^
j
〉
(
1
≤
i
≠
j
≤
M
)
−0.0009
0
.12


## 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 florescencein-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 Figures3(a), 3(b), 3(c), and 3(d) for noiseless data measurement scenarios of the standard Shepp-Logan phantom andsinc2d() solutions under random and fan-beam sampling simulations. Parallel beam sampling andrect_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. Thesinc2d() function is defined as
(65)
s
i
n
c
2
d
(
x
,
y
)
=
sin
(
π
x
)
π
x
sin
(
π
y
)
π
y
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-dimensionalsinc() 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 at10
% percent normalized error.The normalized error is defined as(66)
ϵ
k
=
x
-
x
k
2
x
-
x
0
2
,
and estimate the signal to noise ratio (SNR) is estimated as 1
/
ϵ
k with the x
0 vector set to the the zero vector.


## 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. Thesinc2d() function is defined as
(65)
s
i
n
c
2
d
(
x
,
y
)
=
sin
(
π
x
)
π
x
sin
(
π
y
)
π
y
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-dimensionalsinc() 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 at10
% percent normalized error.The normalized error is defined as(66)
ϵ
k
=
x
-
x
k
2
x
-
x
0
2
,
and estimate the signal to noise ratio (SNR) is estimated as 1
/
ϵ
k with the x
0 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 Figures6 and 7, with measured run-times for 64-bit and 32-bit, respectively. For each method, the size of A matrix was modified and average time to complete was measured. The data matrix was chosen to the behadamard() function for ease of implementation and reference. The solution vector was chosen as a random point in the binomial distribution with p
=
0.5 using thebinornd() function.Figure 6
Relative computational time using dedicated custom real-time Linux kernel 3.12.5 patched versus square matrix dimensionsN
=
2
q  ∀
q
∈
[
4
,
…
,
10
]; CPU is shielded and frequency locked for single process and single core (one of eight) 800 MHz Intel 64-bit i7.Figure 7
Relative computational time for embedded microprocessor dedicated task custom Linux kernel 3.12.7 ARMv7 relative computational time versus matrix dimensionsN
=
2
q  ∀
q
∈
[
4
,
…
,
10
]; CPU ARM frequency locked kernel only process stack 840 MHz 32-bit.It is notable that the theoretical complexity of noise limited Kaczmarz isO
(
N
2
) but the slope of the timings is closer to 1.1 than 2.0. This is attributed to relatively small matrix sizes, and it is expected that the complexity asymptotically approaches ~2.0 for large N.The theoretical and numerical estimates for complexity are shown in Tables6 and 7.Table 6
Numerical estimates for computational complexity using naive codes.
Computational complexity
LS
SVD
Kaczmarz
RK
RKOS
Theoretical
O
(
N
3
)
O
(
N
3
)
O
(
N
2
)
O
(
N
2
)
O
(
N
2
)
Numerical estimates
2.7
3.1
1.1
1.1
1.1Table 7
Numerical estimates for computational complexity using naive codes.
Computational complexity
Cimmino
Vasilchenko
Ref.P1
Ref.P2
Ref. Kaczmarz
Theoretical
O
(
N
2
)
O
(
N
2
)
O
(
N
2
)
O
(
N
3
)
O
(
N
2
)
Numerical estimates
1.6
2.1
2
3.2
1.1


## 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 bya 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 P
/
N term and is computationally feasible for small P relative to N.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.



---

*Source: 908984-2014-12-15.xml*