Randomized Block Kaczmarz Method with Projection for Solving Least Squares

Introduction

The Kaczmarz method is a popular iterative solver of overdetermined systems of linear equations $\bm{A}\bm{x}=\bm{b}$. Because of its simplicity and performance, the method and its derivatives are used in a range of applications from image reconstruction to digital signal processing . The method performs a series of orthogonal projections and iteratively converges to the solution of the system of equations. It is therefore computationally feasible even for very large and overdetermined systems.

Given a vector $\bm{b}$ and an $n\times d$ full rank (real or complex) matrix $\bm{A}$ with rows $\bm{a}_{1},\bm{a}_{2},\ldots\bm{a}_{n}$, the algorithm begins with an initial estimate $\bm{x}_{0}$ and cyclically projects the estimation onto each of the solution spaces. This process can be described as follows:

When the system is perturbed by noise or no longer consistent, $\bm{A}\bm{x}_{\star}+\bm{e}=\bm{b}$, the randomized Kaczmarz method still provides expected linear convergence down to an error threshold ,

where $|e[i]|$ denotes the $i$th entry of $\bm{e}$. This result is sharp, and shows that the randomized Kaczmarz method converges with a radius proportional to the magnitude of the largest entry of the noise in the system. Since the iterates of the Kaczmarz method always lie in a single solution space, the method clearly will not converge to the least squares solution of an inconsistent system.

The bound (2) demonstrates that the randomized Kaczmarz method performs well when the noise in inconsistent systems is small. Zouzias and Freris introduced a variant of the method which utilizes a random projection to iteratively reduce the norm of the error . They show that the estimate of this Randomized Extended Kaczmarz (REK) method converges linearly in expectation to the least squares solution of the system, breaking the radius barrier of the standard method. The algorithm maintains not only an estimate $\bm{x}_{j}$ to the solution but also an approximation $\bm{z}_{j}$ to the projection of $\bm{b}$ onto the range of $\bm{A}$:

where in iteration $j$, $\bm{a}_{i}$ and $\bm{a}^{(k)}$ is the row and column of $\bm{A}$, respectively, each chosen randomly with probability proportional to their Euclidean norms. In this setting, we no longer require that the matrix $\bm{A}$ be full rank, and ask for the least squares solution,

where $\bm{A}^{\dagger}$ denotes the Moore-Penrose pseudoinverse of $\bm{A}$. Zouzias and Freris showed that the REK method converges linearly in expectation to the least squares solution ,

where $\sigma_{\min}(\bm{A})$ is the smallest non-zero singular value of $\bm{A}$ and $K(\bm{A})=\frac{\|\bm{A}\|_{F}}{\sigma_{\min}(\bm{A})}$ denotes its scaled condition number.

2 The block Kaczmarz method

Recently, Needell and Tropp analyzed a block version of the simple randomized Kaczmarz method . Like the traditional method, this version iteratively projects the current estimation onto the solution spaces. However, rather than using the solution space of a single equation, the block method projects onto the solution space of many equations simultaneously by selecting a block of rows rather than a single row. For a subset $\tau\subset\{1,2,\ldots,n\}$, denote by $\bm{A}_{\tau}$ the submatrix of $\bm{A}$ whose rows are indexed by $\tau$. We again begin with an arbitrary guess $\bm{x}_{0}$ for the solution of the system. Then for each iteration $j\geq 1$, select a block $\tau=\tau_{j}$ of rows. To obtain the next iterate, we project the current estimation onto the solution space of the equations listed in $\tau$ :

Here, the conditioning of the blocks $\bm{A_{\tau}}$ plays a crucial row in the behavior of the method. Indeed, if each block is well-conditioned, its pseudoinverse can be applied efficiently using an iterative method such as CGLS . To guarantee such properties, Needell and Tropp utilize a paving of the matrix $\bm{A}$.

A $(p,\alpha,\beta)$ row paving of a $n\times d$ matrix $\bm{A}$ is a partition $\mathcal{T}=\{\tau_{1},\dots,\tau_{p}\}$ of the rows such that

where again we denote by ${\bm{A}}_{\tau}$ the $|\tau|\times d$ submatrix of ${\bm{A}}$. We refer to the number $p$ as the size of the paving, and the numbers $\alpha$ and $\beta$ are called the lower and upper paving bounds.

We refer to a row paving of ${\bm{A}}^{*}$ as a column paving of ${\bm{A}}$. We thus seek pavings of ${\bm{A}}$ with small number of blocks $p$ and upper paving constant $\beta$. In the following, we will assume one has access to such a paving, and discuss in Section 3 how to construct such pavings for various types of matrices. When the matrix has unit-norm rows, equipped with such a row paving of $\bm{A}$, the main result of shows that the randomized block Kaczmarz algorithm (5) exhibits linear convergence in expectation:

where $C^{\prime}$ is an absolute constant and $\kappa(\bm{A})=\frac{\sigma_{\max}(\bm{A})}{\sigma_{\min}(\bm{A})}$ is the condition number of $\bm{A}$.

Since each iteration of the block method utilizes multiple rows, one can compare the rate of (6) and (1) by considering convergence per epoch (one cycle through the rows of $\bm{A}$). From this analysis, one finds the bounds to be comparable. However, the block method can utilize fast matrix multiplies and efficient implementation, yielding dramatic improvements in computational time. See for details and empirical results.

3 Contribution

The REK method breaks the so-called convergence horizon of standard Kaczmarz method, allowing convergence to the least squares solution of inconsistent systems. The block Kaczmarz method on the other hand, allows for significant computational speedup and accelerated convergence to within a fixed radius of the least squares solution. The main contribution of this paper analyzes a randomized block Kaczmarz method which also incorporates a blocked projection step, which provides accelerated convergence to the least squares solution. In this case we need a column partition for the projection step and a row partition for the Kaczmarz step. Our results show that this method offers both linear convergence to the least squares solution $\mathbf{x}_{\text{\tiny LS}}$ when a row paving can be obtained, and improved convergence speed due to the blocking of both the rows and columns. In addition, we present a block coordinate descent variant which utilizes only a column paving and also yields linear convergence. We will see that the desired column paving may be obtained easily for arbitrary matrices, so this variant is especially useful when a good row paving cannot be obtained.

4 Organization

In Section 2 we begin by presenting the double block randomized Kaczmarz method which utilizes both a row and column paving, and provides the strongest theoretical results overall. We discuss methods for obtaining the desired matrix pavings in Section 3. Section 4 introduces the variant of the block and extended Kaczmarz methods that requires only a column paving of the matrix, readily accessible for arbitrary matrices. In Section 6 we present some experimental results for the various algorithms. We conclude with a discussion of related work and open directions in Section 7. The appendix includes proofs of intermediate results used along the way.

The Randomized Double Block Kaczmarz Method

It is natural to ask whether one can consider both a row partition and a column partition in the Kaczmarz method, blocking both in the Kaczmarz update step and the projection step. Indeed, utilizing blocking in both steps yields Algorithm 1 below. We thus propose the following randomized block extended Kaczmarz method using double partitioning. We will see later that since this method requires a row paving, it works very well when each of the row norms of the matrix are relatively similar (see Sections 3 and 5).

Combining the theoretical approaches in we will prove the following result about the convergence of Algorithm 1. This result utilizes both a column paving and row paving.

where $\gamma=1-\frac{\sigma_{\min}^{2}({{\bm{A}}})}{{p}{\beta}}$ and $\overline{\gamma}=1-\frac{\sigma_{\min}^{2}({{\bm{A}}})}{\overline{p}\overline{\beta}}$.

Proof. We begin with the following lemma which is motivated by Lemma $2.2$ of , and shows that the iterates $\mathbf{z}_{k}$ converge linearly to the projection of ${\bm{b}}$ onto the kernel of ${\bm{A}}^{*}$.

Let ${\bm{b}}$ be a fixed vector and $\mathcal{T}$ be a $(\overline{p},\overline{\alpha},\overline{\beta})$ column paving of ${\bm{A}}$. Assuming the notation of Theorem 1, for every $k>0$ it holds that

where ${{\bm{b}}_{\mathcal{R}({\bm{A}})^{\bot}}}:=(\mathbf{I}-{\bm{A}}{{\bm{A}}}^{\dagger}){\bm{b}}$.

Proof. Let ${\bm{P}}_{\tau_{k}}={\bm{A}}_{\tau_{k}}{({\bm{A}}_{\tau_{k}})}^{\dagger}$ and notice $\mathbf{z}_{k}=({\bm{I}}-{\bm{P}}_{\tau_{k}})\mathbf{z}_{k-1}$. Define ${\mathbf{e}}_{k}=\mathbf{z}_{k}-{{\bm{b}}_{\mathcal{R}({\bm{A}})^{\bot}}}$ for $k\geq 0$. Then,

where the first equality follows by the definition of $\mathbf{z}_{k}$, the second by orthogonality between the range of ${\bm{P}}_{\tau_{k}}$ and ${{\bm{b}}_{\mathcal{R}({\bm{A}})^{\bot}}}$, and the final equality by definition of ${\mathbf{e}}_{k-1}$. Next, we prove that

where the first equality follows since ${\bm{P}}_{\tau_{k}}={\bm{U}}_{\tau_{k}}{\bm{U}}_{\tau_{k}}^{*}$, the second equality from the unitary invariance property of the Euclidean norm, the third equality by replacing ${\bm{U}}_{\tau_{k}}^{*}$ with ${\bm{\Sigma}}_{\tau_{k}}^{-1}{\bm{V}}_{\tau_{k}}^{*}{\bm{A}}_{\tau_{k}}^{*}$, the next three lines follow since $\sigma_{\min}^{2}({\bm{\Sigma}}_{\tau_{k}}^{-1}{\bm{V}}_{\tau_{k}}^{*})=1/\sigma_{\max}^{2}({\bm{\Sigma}}_{\tau_{k}})=1/\sigma^{2}_{\max}({\bm{A}}_{\tau_{k}})$ and by the paving assumption, and the final inequality follows since ${\mathbf{e}}_{k}\in\mathcal{R}({\bm{A}})$ for all $k\geq 0$ (indeed, ${\mathbf{e}}_{0}={{\bm{b}}_{\mathcal{R}({\bm{A}})}}\in\mathcal{R}({\bm{A}})$ and it follows that ${\mathbf{e}}_{k}\in\mathcal{R}({\bm{A}})$ for every $k\geq 0$ by the recursive definition of ${\mathbf{e}}_{k}$). It follows that

Repeat the above inequality $k$ times and notice that ${\mathbf{e}}_{0}=\bm{b}-{{\bm{b}}_{\mathcal{R}({\bm{A}})^{\bot}}}={{\bm{b}}_{\mathcal{R}({\bm{A}})}}$ to conclude. $\hfill\Box$

Next, Lemma 2.2 of shows that for any vector $\bm{u}$,

Since the range of $\mathbf{I}-(\bm{A}_{\tau})^{\dagger}\bm{A}_{\tau}$ and $(\bm{A}_{\tau})^{\dagger}$ are orthogonal, we have

where for shorthand we will write $\bm{b}^{\perp}_{\tau}$ to mean $({{\bm{b}}_{\mathcal{R}({\bm{A}})^{\bot}}})_{\tau}$. Combining (9) with $\bm{u}=\bm{x}_{k-1}-\mathbf{x}_{\text{\tiny LS}}$ along with (10), we have

To apply this bound recursively, we will utilize an elementary lemma. It is essentially proved in [53, Theorem 8] but for completeness we recall its proof in the appendix.

Suppose that for some $\gamma,\overline{\gamma}<1$, the following bounds hold for all $k^{*}\geq 0$:

In Section 3 we will see that matrix row-pavings for standardized matrices, those whose rows have unit norm, can be obtained readily. One can thus use the column-normalized version of the matrix ${\bm{A}}$ in line 5 of Algorithm 1 and the corresponding column paving. Note that one need not have access to the complete column-standardized matrix, instead the columns of the submatrix in line 5 could be normalized on the fly. See Section 3 for more on obtaining row pavings and further details.

2 Comparison of Convergence Rates

The bound of Theorem 1 improves upon that of the randomized block Kaczmarz method because it demonstrates linear convergence to the least squares solution $\mathbf{x}_{\text{\tiny LS}}$, whereas (6) shows convergence only within a radius proportional to $\|{\mathbf{e}}\|_{2}^{2}$, which we call the convergence horizon. Algorithm 1 is able to break this barrier because it iteratively removes the component of ${\bm{b}}$ which is orthogonal to the range of ${\bm{A}}$. This of course is also true of the randomized Extended Kaczmarz method (3) as it also breaks this horizon barrier. To compare the rate of (4) to that of Theorem 5, we consider two important scenarios.

First, consider the case when ${\bm{A}}$ is nearly square, and each submatrix can be applied efficiently via a fast multiply. In this case, each iteration of Algorithm 1 incurs approximately the same computational cost as an iteration of the REK method. Thus, we may directly compare the convergence rates of Theorem 5 and (4) to find that Algorithm 1 is about $n/(p\beta)$ times faster than REK in this setting. Thus when $n$ is much larger than $p\beta$, this can result in a significant speedup.

Alternatively, if the matrix ${\bm{A}}$ does not admit a fast multiply, it is fair to only compare the convergence rate per epoch, since each iteration of Algorithm 1 may require more computational cost than those of REK. Since an epoch of Algorithm 1 and REK consist of $p$ and $n$ iterations, respectively, we see that the rate of the former is proportional to $\sigma^{2}_{\min}({\bm{A}})/\beta$ whereas that of REK is proportional to $\sigma^{2}_{\min}({\bm{A}})$. We see in this case that these bounds suggest REK exhibits faster convergence (assuming $\beta>1$). However, as observed in the randomized Block Kaczmarz method, the block methods still display faster convergence than their single counterparts because of implicit computational issues in the linear algebraic subroutines. See the discussion in and the experimental results below for further details.

Obtaining matrix pavings

We devote this section to a brief discussion about matrix pavings and how they may be obtained to utilize the results of Theorem 1. The results discussed here on matrix pavings stem from subset selection, the problem of selecting a large submatrix with desired geoemtric properties, whose origins come from the well-known Restricted Invertibility Principle of Bourgain Tzafriri . The literature now contains several results on subset selection, see e.g. and for recent advancements. We summarize here a few results useful for our purposes.

The simplest type of matrix to first consider is one which has unit-norm rows, which we call row-standardized (and a matrix with unit-norm columns is column-standardized). A surprising result shows that every row-standardized matrix admits a row paving with well-controlled paving parameters. Tropp proves the following result in [45, Thm. 1.2], whose origins are due to Bourgain and Tzafriri and Vershynin .

Fix a number $\delta\in(0,1)$ and row-standardized matrix ${\bm{A}}$ with $n$ rows. Then ${\bm{A}}$ admits a $(p,\alpha,\beta)$ row paving with

Proposition 4 shows the existence of such a paving, but the literature provides various efficient mechanisms for the construction of good pavings as well. In many cases one constructs such a paving simply by choosing a partition of an appropriate size at random, see for an efficient method to compute a paving satisfying (13). See also and the references therein for a thorough discussion of these types of results.

If the matrix ${\bm{A}}$ is row-standardized, one can thus construct a row paving satisfying (13). If the matrix also naturally admits a column paving (for example if it is symmetric or positive semi-definite), then both pavings will have such bounded parameters. If the latter does not hold, one can instead utilize the column-standardized version of ${\bm{A}}$, which we denote $\overline{{\bm{A}}}$, in line 5 of Algorithm 1. The standardization can either be done on the fly during the algorithm, or ahead of time. Either way, Proposition 4 can be combined with Theorem 1 to yeild the following corollary.

where $\gamma=1-\frac{C}{\kappa^{2}({\bm{A}})\log(1+n)}$, $\overline{\gamma}=1-\frac{C}{\kappa^{2}(\overline{{\bm{A}}})\log(1+d)}$and $\kappa(\bm{A})=\frac{\sigma_{\max}(\bm{A})}{\sigma_{\min}(\bm{A})}$ and $\kappa(\overline{{\bm{A}}})=\frac{\sigma_{\max}(\overline{{\bm{A}}})}{\sigma_{\min}(\overline{{\bm{A}}})}$ denote the condition numbers of $\bm{A}$ and $\overline{{\bm{A}}}$, respectively.

If the matrix ${\bm{A}}$ is not row-standardized, one can run Algorithm 1 on the standardized system, along with a paving guaranteed by Proposition 4. Clearly, the method will converge to the new least squares solution, which does not necessarily coincide with the original least squares solution in the inconsistent case.

If one uses this same paving for ${\bm{A}}$, one has (quite pessimistically) that the corresponding paving parameters $p$, $\alpha$, and $\beta$ for ${\bm{A}}$ satisfy

One can then directly apply Theorem 1 with this paving to obtain an analogous version of Corollary 5 for arbitrary matrices, which shows that

where now $\gamma=1-\frac{Ca_{\text{min}}}{a_{\text{max}}\kappa^{2}({\bm{A}})\log(1+n)}=1-\frac{C}{a_{r}\kappa^{2}({\bm{A}})\log(1+n)}$ where we have written $a_{r}=a_{\text{max}}/a_{\text{min}}$ to denote the dynamic range of the row norms of ${\bm{A}}$ (and $\overline{\gamma}$ remains the same as in the corollary). This demonstrates that if the dynamic range is bounded, the convergence rate is the same as in the standardized case, up to constants. In either case, these results demonstrate linear convergence to the true least squares solution.

An alternative to these types of bounds can be obtained by simply only using a column paving for the matrix. As we elaborate below, the least squares solution can still be obtained from the column-normalized system. Since column-pavings can be easily attained in this case via Proposition 4, such a method offers a nice alternative in situations where the dynamic range of the row norms is unknown or unbounded. We propose such a method in the next section.

A Randomized Block Coordinate Descent Method

We next present a simple variant of the extended Kaczmarz method which utilizes only a column paving of the matrix. For that reason, one need not worry about whether the matrix ${\bm{A}}$ is row-standardized. Moreover, the column-standardized version can be used within the algorithm which guarantees bounded paving parameters, while still finding the true least squares solution.

Utilizing the benefits of both the block variant and the randomized extension, we propose the following randomized block coordinate descent method for the inconsistent case.

We may utilize some of the previous analysis to prove convergence of Algorithm 2. Observe that Step 6 of Algorithm 1 is identical to Steps 5 and 7 of Algorithm 2, therefore Lemma 2 implies that

To utilize this result, we aim to relate the iterates $\mathbf{z}_{k}$ to the estimation $\mathbf{x}_{k}$. The following claim quantifies precisely this relation.

For every $k\geq 0$, at the end of the $k$-th iteration, it holds that $\mathbf{z}_{k+1}={\bm{b}}-{\bm{A}}\mathbf{x}_{k+1}$.

Combining this lemma with (15) yields the following result which shows convergence of the estimation to the least squares solution under the map $\bm{A}$.

where the first equality follows by ${{\bm{b}}_{\mathcal{R}({\bm{A}})}}={\bm{A}}{{\bm{A}}}^{\dagger}{\bm{b}}={\bm{A}}\mathbf{x}_{\text{\tiny LS}}$, the second by orthogonality ${\bm{b}}={{\bm{b}}_{\mathcal{R}({\bm{A}})^{\bot}}}+{{\bm{b}}_{\mathcal{R}({\bm{A}})}}$ and the last equality from Lemma 6. Combined with inequality (7) this yields the desired result. $\hfill\Box$

When $\bm{A}$ has full column rank, we may bound the estimation error $\|\mathbf{x}_{\text{\tiny LS}}-\mathbf{x}_{T}\|_{2}$ by $\frac{1}{\sigma_{\min}(\bm{A})}\left\|{\bm{A}}(\mathbf{x}_{\text{\tiny LS}}-\mathbf{x}_{T})\right\|_{2}$ which combined with the fact that $\left\|{{\bm{b}}_{\mathcal{R}({\bm{A}})}}\right\|_{2}\leq\sigma_{\max}({\bm{A}})\left\|\mathbf{x}_{\text{\tiny LS}}\right\|_{2}$ implies the following corollary.

2 Implementation

The advantage to a single paving approach as in Algorithm 2 is that one can utilize the column-standardized version of ${\bm{A}}$ while maintaining the same convergence to the (scaled version of the) least squares solution. Utilizing Proposition 4, one is guaranteed a column-paving satisfying (13), so that paving parameters of Theorem 7 and Corollary 8 are bounded. Since re-normalizing the columns of the matrix only re-scales the entries of $\mathbf{x}_{\text{\tiny LS}}$, Theorem 7 and Corollary 8 grant one access to the original least squares solution $\mathbf{x}_{\text{\tiny LS}}$. To be precise, now let $\bm{D}$ be the diagonal matrix whose entries correspond to the reciprocals of the column norms of ${\bm{A}}$, so that $\overline{{\bm{A}}}={\bm{A}}\bm{D}$ has unit-norm columns. Let ${\overline{\mathbf{x}}_{\text{\tiny LS}}}=\overline{{\bm{A}}}^{\dagger}{\bm{b}}$ denote the least squares solution of the re-normalized system, so that one has ${\overline{\mathbf{x}}_{\text{\tiny LS}}}=\bm{D}^{-1}\mathbf{x}_{\text{\tiny LS}}$. Then

Thus applying Theorem 7 for $\overline{{\bm{A}}}$, and utilizing the fact that the range of ${\bm{A}}$ is the same as that of $\overline{{\bm{A}}}$, one has

This then implies that when the matrix is full rank,

Since $\overline{p}$ and $\overline{\beta}$ are the paving parameters of $\overline{A}$, one has by utilizing Proposition 4 and substituting the bounds of (13) into (17) that

Although this bound does depend on the conditioning of both ${\bm{A}}$ and $\overline{{\bm{A}}}$ (note the rate itself only depends on the conditioning of the latter), it is the first to guarantee linear convergence to the true least squares solution utilizing a paving which is guaranteed by Proposition 4 while not placing any restrictions on the matrix ${\bm{A}}$ itself (such as standardization).

In considering the improvements offered by both the REK method and the block Kaczmarz method, one may ask whether it is advantageous to run a traditional REK projection step as in (3) along with a traditional block Kaczmarz update step as in (5). However, empirically we have observed that such a combination actually leads to a degradation in performance and requires far more epochs to converge than the algorithms discussed above. We conjecture that it is important to run both the projection update and the Kaczmarz update “at the same speed”; if the Kaczmarz update utilizes many rows at once, so should the projection update, and vice versa.

Summary of Approaches

Here we detail the various approaches proposed, and summarize the practical implementation in several frameworks. Since the block variants of the methods require the blocks be well conditioned, it is natural to rely on matrix pavings for the analysis. Unfortunately, such pavings are only readily available when the matrix itself is properly normalized. For that reason, we have proposed several practical alternatives for the important setting in which normalization is not feasible. Both of our proposed methods still offer computational advantages over the standard approaches (see the next section). We summarize these approaches here, which cover all possible settings.

When the system is consistent, one does not lose any convergence properties by re-scaling the system to be standardized, since the solution remains the same. In this simple setting, one can utilize the standardized system (either standardizing a priori or on the fly), and benefit from the convergence guaranteed by Corollary 5.

When the matrix is already standardized (as naturally occurs for example in Vandermonde matrices used in trigonometric approximation), as in the above case one can immediately utilize Corollary 5 to guarantee convergence to the least squares solution.

If the system is not standardized, but the ratio of the largest to smallest row norm is bounded (as is the case in random matrices for example, which have tightly concentrated row norms), one can still utilize this result. Indeed, by utilizing the same paving one would use if the matrix was actually standardized – but not actually standardizing the system, Corollary 5 can be used to obtain the convergence rate given in (14). One sees that if the dynamic range $a_{r}$ is bounded (by say, a constant or even $\log n$), the method convergences to the original least squares solution with approximately the same convergence rate.

If the matrix has a large variety of row norms, it may clearly be challenging to guarantee the desired row paving. For that reason, it will be advantageous to use a column paving. To that end, we propose Algorithm 2 which is designed for systems that cannot be separated into well-conditioned row blocks. As discussed in Section 4.2, Theorem 7 can be used via the column-standardized paving to obtain the convergence rate given in (18). This bound guarantees linear convergence to the true least squares solution, even for matrices which are far from standardized. The disadvantage of course is that the rate depends on the conditioning of the column-standardized version, which may be hard to explicitly bound (as is true for general large matrices anyway).

In any of these cases, our results may be used to guarantee linear convergence in expectation to the true least squares solution of the system. Moreover, because matrix blocks can be utilized, we often see a significant speedup in runtime due to practical considerations. See the next section for examples of such behavior.

Experimental Results

Here we present some experiments using simple examples to illustrate the benefits of block methods. We do not claim optimized implementations of the method, and only run on small problem sizes; our purpose is only to demonstrate that even in these simple examples, the block method offers advantages to the standard method. We refer the reader to for more empirical results for both REK and block methods.

Lastly, we tested the methods on tomography problems, generated using the Matlab Regularization Toolbox by P.C. Hansen (http://www.imm.dtu.dk/~pcha/Regutools/) . In particular we present a 2D tomography problem $\bm{A}\bm{x}=\bm{b}$ for an $n\times d$ matrix with $n=fN^{2}$ and $d=N^{2}$. Here $\bm{A}$ corresponds to the absorption along a random line through an $N\times N$ grid. In our experiments we set $N=20$ and the oversampling factor $f=3$. This yielded a matrix $\bm{A}$ with condition number $\kappa(\bm{A})=2.08$. Since for this matrix it may be difficult to obtain a row paving, we instead use Algorithm 2 and obtain a random column paving from the standardized version and use that for the matrix $\bm{A}$. The results for various choices of paving size are displayed in Figure 5, which are in line with previous experiments.

Related Work and Discussion

The Kaczmarz method was first introduced in the 1937 work of Kaczmarz himself . Since then, the method has been revitalized by researchers in computer tomography, under the name Algebraic Reconstruction Technique (ART) . Deterministic convergence results for the method often depend on properties of the matrix that are difficult to compute or analyze . Moreover, it has been well observed that random choice of row selection often speeds up the convergence .

Recently, Strohmer and Vershynin derived the first provable convergence rate of the Kaczmarz method, showing that when each row is selected with probability proportional to its norm the method exhibits the expected linear convergence of (1). This work was extended to the inconsistent case in , which shows linear convergence to within some fixed radius of the least squares solution. The almost-sure guarantees were recently derived by Chen and Powell . To break the convergence barrier, relaxation parameters can be introduced, so that each iterate is over or under projected onto each solution space. Whitney and Meany prove that if the relaxation parameters tend to zero that the iterates converge to the least squares solution . Further results using relaxation have also been obtained, see for example . An alternative to relaxation parameters was recently proposed by Zouzias and Freris as the REK method described by (3). Rather than alter the projection step, motivated by ideas of Popa they introduce a secondary step which aims to reduce the residual.

The Kaczmarz method has been extended beyond linear systems as well. For example, Leventhal and Lewis analyze the method for systems with polyhedral constraints and inequalities, which was also extended to the block case , and Richtárik and Takávc build on these results for general optimization problems.

Another important aspect of research in this area focuses on accelerating the convergence of the methods. Geometric brute force methods can be used , additional row directions may be added , or instead one can select blocks of rows rather than a single row in each iteration. The block version of the Kaczmarz method is originally due to work of Elfving and Eggermont et al. . Its convergence rates were recently studied in and analyzed via pavings by Needell and Tropp . The block Kaczmarz method is of course a special instance in a broader class of block projection algorithms, see for example for a more general analysis and for a presentation of other block variants.

To use block methods effectively, one needs to obtain a suitable partition of the rows (and/or columns). Popa constructs such partitions by creating orthogonal blocks , whereas Needell and Tropp promote the use of row pavings to construct the partition .

Construction of pavings has been studied for quite some time now, and most early results rely on random selection. The guarantee of lower and upper paving bounds has been derived by Bourgain and Tzafriri and Kashin and Tzafriri , respectively. Simultaneous guarantees were later derived by Bourgain and Tzafriri with suboptimal dependence on the matrix norm. Recently, Spielman and Srivastava and Youssef provided simple proofs of the results from and , respectively. Vershynin and Srivastava extend the paving results to general matrices with arbitrary row norms; see also . Proposition 4 follows from the work of Vershynin and Tropp , and is attributed to the seminal work of Bourgain and Tzafriri . For particular classes of matrices, the paving can even be obtained from a random partition of the rows with high probability. This is proved by Tropp using ideas from , and is refined in .

Acknowledgments

D.N. is thankful to the Simons Foundation Collaboration grant and the Alfred P. Sloan Fellowship. A.Z. has received funding from the European Research Council under the European Union’s Seventh Framework Program (FP7/2007-2013) / ERC grant agreement $n^{o}$ 259569. We would also like to thank Anna Ma for thoughtful discussions, and the reviewers for useful comments which significantly improved the manuscript.

Appendix A Proof of intermediate results

Assume the bounds (12) hold. Applying the first bound in (12) recursively yields

where the second inequality holds by the assumption that $r_{k}\leq\gamma^{k}B\leq B$, and the last by the properties of the geometric summation. Similarly, observe that for any $k$ and $k^{*}$ we have

Now we choose $k$ and $k^{*}$ such that $T=k+k^{*}$ and $k=k^{*}$ if $T$ is even, or $k=k^{*}+1$ if $T$ is odd. Combining the two inequalities above, we have

References