Two-subspace Projection Method for Coherent Overdetermined Systems (Technical Report)

Deanna Needell, Rachel Ward

Introduction

We consider a consistent system of linear equations of the form

Theoretical results about the rate of convergence of the Kaczmarz method have been difficult to obtain, and most are based on quantities which are themselves hard to compute . Even more importantly, the method as we have just described depends heavily on the ordering of the rows of $\bm{A}$ . A malicious or unlucky ordering may therefore lead to extremely slow convergence. To overcome this, one can select the rows of $\bm{A}$ in a random fashion rather than cyclically . Strohmer and Vershynin analyzed a randomized version of the Kaczmarz method that in each iteration selects a row of $\bm{A}$ with probability proportional to its Euclidean norm . Thus in the standardized case we consider here, a row of $\bm{A}$ is chosen uniformly at random. This randomized Kaczmarz method is described by the following pseudocode.

Note that this method as stated selects each row with replacement, see for a discussion on the differences in performance when selecting with and without replacement. Strohmer and Vershynin show that this method exhibits exponential convergence in expectation ,

Leventhal and Lewis show that for certain probability distributions, the expected rate of convergence can be bounded in terms of other natural linear-algebraic quantities. They propose generalizations to other convex systems . Recently, Chen and Powell proved that for certain classes of random matrices $\bm{A}$ , the randomized Kaczmarz method convergences exponentially to the solution not only in expectation but also almost surely .

In the presence of noise, one considers the possibly inconsistent system $\bm{A}\bm{x}+\bm{w}\approx\bm{b}$ for some error vector $\bm{w}$ . In this case the randomized Kaczmarz method converges exponentially fast to the solution within an error threshold ,

where $R$ the the scaled condition number as in (1.1) and $\|\cdot\|_{\infty}$ denotes the largest entry in magnitude of its argument. This error is sharp in general . Modified Kaczmarz algorithms can also be used to solve the least squares version of this problem, see for example and the references therein.

Although the convergence results for the randomized Kaczmarz method hold for any consistent system, the factor $\frac{1}{R}$ in the convergence rate may be quite small for matrices with many correlated rows. Consider for example the reconstruction of a bandlimited function from nonuniformly spaced samples, as often arises in geophysics as it can be physically challenging to take uniform samples. Expressed as a system of linear equations, the sampling points form the rows of a matrix $\bm{A}$ ; for points that are close together, the corresponding rows will be highly correlated.

To be precise, we examine the coherence of a standardized matrix $\bm{A}$ by defining the quantities

Note that because $\bm{A}$ is standardized, $0\leq\delta\leq\Delta\leq 1$ . It is clear that when $\bm{A}$ has high coherence parameters, $\|\bm{A}^{-1}\|$ is very small and thus the factor $R$ in (1.1) is also small, leading to a weak bound on the convergence. Indeed, when the matrix has highly correlated rows, the angles between successive orthogonal projections are small and convergence is stunted. We can explore a wider range of orthogonal directions by looking towards solution hyperplanes spanned by pairs of rows of $\bm{A}$ . We thus propose a modification to the randomized Kaczmarz method where each iteration performs an orthogonal projection onto a two-dimensional subspace spanned by a randomly-selected pair of rows. We point out that the idea of projecting in each iteration onto a subspace obtained from multiple rows rather than a single row has been previously investigated numerically, see e.g. .

With this as our goal, a single iteration of the modified algorithm will consist of the following steps. Let ${\bm{x_{k}}}$ denote the current estimation in the $k$ th iteration.

Select two distinct rows $\bm{a_{r}}$ and $\bm{a_{s}}$ of the matrix $\bm{A}$ at random

Compute the translation parameter $\varepsilon$

Perform an intermediate projection: $\bm{y}\leftarrow\bm{x_{k}}+\varepsilon(b_{r}-\langle\bm{x_{k}},\bm{a_{r}}\rangle)\bm{a_{r}}$

Perform the final projection to update the estimation: $\bm{x_{k+1}}\leftarrow\bm{y}+(b_{s}-\langle\bm{y},\bm{a_{s}}\rangle)\bm{a_{s}}$

In general, the optimal choice of $\varepsilon$ at each iteration of the two-step procedure corresponds to subtracting from $\bm{x_{k}}$ its orthogonal projection onto the solution space $\{\bm{x}:\langle\bm{a_{r}},\bm{x}\rangle=b_{r}\text{ and }\langle\bm{a_{s}},\bm{x}\rangle=b_{s}\}$ , which motivates the name two-subspace Kaczmarz method. By optimal choice of $\varepsilon$ , we mean the value $\varepsilon_{opt}$ minimizing the residual $\|\bm{x}-\bm{x_{k+1}}\|_{2}^{2}$ . Expanded, this reads

Using that the minimizer of $\|\gamma\bm{w}+\bm{z}\|_{2}^{2}$ is $\gamma=-\frac{\left\langle\bm{w},\bm{z}\right\rangle}{\|\bm{w}\|_{2}^{2}}$ , we see that

Note that the unknown vector $\bm{x}$ appears in this expression only through its observable inner products, and so $\varepsilon_{opt}$ is computable. After some algebra, one finds that the two-step procedure with this choice of $\varepsilon_{opt}$ can be re-written as follows .

Our main result shows that the two-subspace Kaczmarz algorithm provides the same exponential convergence rate as the standard method in general, and substantially improved convergence when the rows of $\bm{A}$ are coherent . Figure 1 plots two iterations of the one-subspace random Kaczmarz and compares this to a single iteration of the two-subspace Kaczmarz algorithm.

Let $\bm{A}$ be a full-rank standardized matrix with $n$ columns and $m>n$ rows and suppose $\bm{A}\bm{x}=\bm{b}$ . Let $\bm{x_{k}}$ denote the estimation to the solution $\bm{x}$ in the $k$ th iteration of the two-subspace Kaczmarz method. Then

where $D=\min\Big\{\frac{\delta^{2}(1-\delta)}{1+\delta},\frac{\Delta^{2}(1-\Delta)}{1+\Delta}\Big\}$ , $\Delta$ and $\delta$ are the coherence parameters (1.3), and $R=\|\bm{A}\|_{F}^{2}\|\bm{A}^{-1}\|^{2}$ denotes the scaled condition number.

Remarks. 1. When $\Delta=1$ or $\delta=0$ we recover the same convergence rate as provided for the standard Kaczmarz method (1.1) since the two-subspace method utilizes two projections per iteration.

2. The bound presented in Theorem 1.1 is a pessimistic bound. Even when $\Delta=1$ or $\delta=0$ , the two-subspace method improves on the standard method if any rows of $\bm{A}$ are highly correlated (but not equal). This is evident in the proof of Theorem 1.1 in Section 2 but we present this bound for simplicity. See also Section 4 for more details on improved convergence bounds.

Figure 2 shows the value of $D$ of Theorem 1.1 for various values of $\Delta$ and $\delta$ . This demonstrates that in the best case (when $\delta\approx\Delta\approx 0.62$ ), the convergence rate is improved by at least a factor of 0.1.

2. Organization

The remainder of the report is organized as follows. In Section 2 we state and prove the main lemmas which serve as the proof of Theorem 1.1. Section 3 discusses the two-subspace Kaczmarz method in the presence of noise and shows that in this case the method exhibits exponential convergence to an error threshold. Section 4 presents further modifications of the two-subspace Kaczmarz method which provide even more improvements on the provable convergence bounds. A discussion of these methods is provided in Section 6. We conclude with numerical experiments demonstrating the improvements from our method in Section 5.

Main Results

We now present the proof of Theorem 1.1. We first derive a bound for the expected progress made in a single iteration. Since the two row indices are chosen independently at each iteration, we will be able to apply the bound recursively to obtain the desired overall expected convergence rate.

Our first lemma shows that the expected estimation error in a single iteration of the two-subspace Kaczmarz method is decreased by a factor strictly less than that of the standard randomized method.

Let $\bm{x_{k}}$ denote the estimation to the solution of $\bm{A}\bm{x}=\bm{b}$ in the $k$ th iteration of the two-subspace Kaczmarz method. Denote the rows of $\bm{A}$ by $\bm{a}_{1},\bm{a}_{2},\ldots\bm{a}_{m}$ . Then we have the following bound,

where $C_{r,s}=\frac{|\mu_{r,s}|-\mu_{r,s}^{2}}{\sqrt{1-\mu_{r,s}^{2}}}$ , $\mu_{r,s}=\left\langle\bm{a_{r}},\bm{a_{s}}\right\rangle$ , and $R=\|\bm{A}^{-1}\|^{2}\|\bm{A}\|_{F}^{2}$ denotes the scaled condition number.

We fix an iteration $k$ and for convenience refer to $\bm{v}_{k}$ , $\mu_{k}$ , and $\bm{y}_{k}$ as $\bm{v}$ , $\mu$ , and $\bm{y}$ , respectively. We will also denote $\gamma=\langle\bm{a_{r}},\bm{v}\rangle$ .

First, observe that by the definitions of $\bm{v}$ and $\bm{x_{k}}$ we have

Since $\bm{a_{s}}$ and $\bm{v}$ are orthonormal, this gives the estimate

We wish to compare this error with the error from the standard randomized Kaczmarz method. Since we utilize two rows per iteration in the two-subspace Kaczmarz method, we compare its error with the error from two iterations of the standard method. Let $\bm{z}$ and $\bm{z^{\prime}}$ be two subsequent estimates in the standard method following the estimate $\bm{x_{k-1}}$ , and assume $\bm{z}\neq\bm{z^{\prime}}$ . That is,

Recalling the definitions of $\bm{v}$ , $\mu$ and $\gamma$ , we have

Now substituting this into (2.2) and taking the orthogonality of $\bm{a_{s}}$ and $\bm{v}$ into account,

For convenience, let $\bm{e_{k-1}}=\bm{x}-\bm{x_{k-1}}$ denote the error in the $(k-1)$ st iteration of two-subspace Kaczmarz. Then we have

The third equality follows from the orthonormality of $\bm{a_{s}}$ and $\bm{v}$ . We now expand the last term,

It thus remains to analyze the last term. Since we select the two rows $r$ and $s$ independently from the uniform distribution over pairs of distinct rows, the expected error is just the average of the error over all $m^{2}-m$ ordered choices ${r,s}$ . To this end we introduce the notation $\mu_{r,s}=\langle\bm{a_{r}},\bm{a_{s}}\rangle$ . Then by definitions of $\bm{v}$ , $\mu$ and $\gamma$ ,

We now recall that for any $\theta,\pi,u,$ and $v$ ,

Setting $\theta_{r,s}=\frac{\mu_{r,s}^{2}}{\sqrt{1-\mu_{r,s}^{2}}}$ and $\pi_{r,s}=\frac{\mu_{r,s}}{\sqrt{1-\mu_{r,s}^{2}}}$ , we have by rearranging terms in the symmetric sum,

Since selecting two rows without replacement (i.e. guaranteeing not to select the same row back to back) can only speed the convergence, we have from (1.1) that the error from the standard randomized Kaczmarz method satisfies

Combining this with (2.4) and (2) yields the desired result.

Although the result of Lemma 2.1 is tighter, using the coherence parameters $\delta$ and $\Delta$ of (1.3) allows us to present the following looser but simpler result.

Let $\bm{x_{k}}$ denote the estimation to $\bm{A}\bm{x}=\bm{b}$ in the $k$ th iteration of the two-subspace Kaczmarz method. Denote the rows of $\bm{A}$ by $\bm{a}_{1},\bm{a}_{2},\ldots\bm{a}_{m}$ . Then

where $D=\min\Big\{\frac{\delta^{2}(1-\delta)}{1+\delta},\frac{\Delta^{2}(1-\Delta)}{1+\Delta}\Big\}$ , $\delta$ and $\Delta$ are the coherence parameters as in (1.3), and $R=\|\bm{A}^{-1}\|^{2}\|\bm{A}\|_{F}^{2}$ denotes the scaled condition number.

By the assumption that $\delta\leq|\langle\bm{a_{r}},\bm{a_{s}}\rangle|\leq\Delta$ , we have

In the last inequality we have employed the fact that for any $\bm{z}$ ,

Combining (2) and (2.6) along with the definition of $R$ yields the claim.

Applying Lemma 2.2 recursively and using the fact that the selection of rows in each iteration is independent yields our main result Theorem 1.1.

Noisy Systems

Next we consider systems which have been perturbed by noise. The inconsistent system $\bm{b}=\bm{A}\bm{x}$ now becomes (the possibly inconsistent system) $\bm{b}=\bm{A}\bm{x}+\bm{w}$ for some error vector $\bm{w}$ . As evident from (1.2), the standard method with noise exhibits exponential convergence down to an error threshold, which is proportional to $\|\bm{w}\|_{\infty}$ . Our main result in the noisy case is that the two-subspace version again exhibits even faster exponential convergence, down to a threshold also proportional to $\|\bm{w}\|_{\infty}$ .

Let $\bm{A}$ be a full rank matrix with $m$ rows and suppose $\bm{b}=\bm{A}\bm{x}+\bm{w}$ is a noisy system of equations. Let $\bm{x_{k}}$ denote the estimation to the solution $\bm{x}$ in the $k$ th iteration of the two-subspace Kaczmarz method. Then

where $\eta=\left(1-\frac{1}{R}\right)^{2}-\frac{D}{R}$ , $D=\min\Big\{\frac{\delta^{2}(1-\delta)}{1+\delta},\frac{\Delta^{2}(1-\Delta)}{1+\Delta}\Big\}$ , $\Delta$ and $\delta$ are the coherence parameters (1.3), and $R=\|\bm{A}^{-1}\|^{2}\|\bm{A}\|_{F}^{2}$ denotes the scaled condition number.

As in the case of our main result Theorem 1.1, this bound is not tight. The same improvements mentioned in the remarks about Theorem 1.1 can also be applied here. In particular, the dependence on $\Delta$ seems to be only an artifact of the proof (see Section 5. Nonetheless, this result still shows that the two-subspace Kaczmarz method provides expected exponential convergence down to an error threshold which is analagous to that of the standard method. The convergence factors are again substantially better than the standard method for coherent systems.

Fix an iteration $k$ and denote by $\bm{y}$ , $\bm{v}$ , $\mu$ and $\beta$ the values of $\bm{y_{k}}$ , $\bm{v_{k}}$ , $\mu_{k}$ and $\beta_{k}$ for convenience. Let $\bm{y^{\prime}}$ and $\beta^{\prime}$ be the values of $\bm{y_{k}}$ , and $\beta_{k}$ as if there were noise (i.e. $\bm{w}=0$ ). In other words, we have

Then by the definition of $\bm{x_{k}}$ , we have

where $\bm{x^{*}_{k}}=\bm{y^{\prime}}+(\beta^{\prime}-\left\langle\bm{x_{k-1}},\bm{v}\right\rangle)\bm{v}$ denotes the next estimation from $\bm{x_{k-1}}$ if there were no noise. Therefore, we have that

Further improvements

Next we state and prove a lemma which demonstrates even more improvements on the convergence rate from the standard method in the case where the correlations between the rows are non-negative. If this is not the case, we may alter one step of the two-subspace method to generalize the result to matrices with arbitrary correlations. This modification will decrease the factor yet again in the exponential convergence rate of the two-subspace method. We consider the noiseless case here, although results analagous to those in Section 3 can easily be obtained using the same methods.

We first define an $m^{2}\times n$ matrix $\bm{\Omega}$ whose rows $\omega$ are differnces of the rows of $\bm{A}$ :

Let $\bm{x_{k}}$ denote the estimation to $\bm{A}\bm{x}=\bm{b}$ in the $k$ th iteration of the two-subspace Kaczmarz method. For indices $r$ and $s$ , set $\mu_{r,s}=\langle\bm{a_{r}},\bm{a_{s}}\rangle$ . We have the following bound,

where $C_{r,s}=\frac{\mu_{r,s}^{2}(1-\mu_{r,s})}{1+\mu_{r,s}}$ and $E_{i,j}=4\mu_{i,j}^{3}$ .

If the correlations $\mu_{r,s}=\langle\bm{a_{r}},\bm{a_{s}}\rangle$ between the rows of $\bm{A}$ are non-negative, then the constants $E_{i,j}$ are all non-negative and thus this lemma offers a strict improvement over Lemma 2.1. However, if this is not the case, this bound may actually be worse than that of Lemma 2.1. To overcome this, we may simply modify the two-subspace Kaczmarz algorithm so that in each iteration $\mu_{r,s}$ is non-negative (by possibly using $-\bm{a_{r}}$ instead of $\bm{a_{r}}$ when needed for example). This modification seems necessary only for the proof, and empirical results for the modified and unmodified methods remain the same. From this point on, we will assume this modification is in place.

We again fix an iteration $k$ and for convenience refer to $\bm{v}_{k}$ , $\mu_{k}$ , and $\bm{y}_{k}$ as $\bm{v}$ , $\mu$ , and $\bm{y}$ , respectively. We will also let $\bm{e_{k-1}}=\bm{x}-\bm{x_{k-1}}$ be the error in the $(k-1)$ st iteration.

We now analyze the last term carefully. To take expectation we must look over all combinations of choices ${r,s}$ , so to that end denote $\langle\bm{a_{r}},\bm{a_{s}}\rangle$ by $\mu_{r,s}$ . Since we select two rows uniformly at random (with replacement), using the definitions of $\bm{c}$ , $\mu$ and $\gamma$ , we have

Thus taking advantage of the symmetry in the sum we have,

We may now use the coherence parameters $\delta$ and $\Delta$ from (1.3) to obtain the following simplified result.

Let $\bm{x_{k}}$ denote the estimation to $\bm{A}\bm{x}=\bm{b}$ in the $k$ th iteration of the two-subspace Kaczmarz method. Then,

where $D=\min\Big\{\frac{\delta^{2}(1-\delta)}{1+\delta},\frac{\Delta^{2}(1-\Delta)}{1+\Delta}\Big\}$ , $E=4\delta^{3}$ and $R=\|\bm{A}^{-1}\|^{2}\|\bm{A}\|_{F}^{2}$ and $Q=\|\bm{\Omega}^{-1}\|^{2}\|\bm{\Omega}\|_{F}^{2}$ denote the scaled condition numbers of $\bm{A}$ and $\bm{\Omega}$ (from (4.1), respectively.

In light of Lemma 4.1 and the proof of Lemma 2.2, it suffices to show that

By the definition (1.3) of $\delta$ and $\|\Omega^{-1}\|$ , we have

The last equality follows since the rows of $\Omega$ are unit norm and equal to zero for $i=j$ .

Applying Lemma 4.2 recursively yields our main theorem.

Let $\bm{x_{k}}$ denote the estimation to $\bm{A}\bm{x}=\bm{b}$ in the $k$ th iteration of the two-subspace Kaczmarz method. Then,

Numerical Results

Next we perform several experiments to compare the convergence rate of the two-subspace randomized Kaczmarz with that of the standard randomized Kaczmarz method. As discussed, both methods exhibit exponential convergence in expectation, but in many regimes the constant factor in the exponential bound of the two-subspace method is much smaller, yielding much faster convergence.

To test these methods, we construct various types of $500\times 50$ matrices $\bm{A}$ . To get a range of $\delta$ and $\Delta$ , we set the entries of $\bm{A}$ to be independent indentically distributed uniform random variables on some interval $[c,1]$ . Changing the value of $c$ will appropriately change the values of $\delta$ and $\Delta$ . Note that there is nothing special about this interval, other intervals (both negative and positive or both) of varying widths yield the same results. For each matrix construction, both the randomized Kaczmarz and two-subspace randomized methods are run with the same initial (randomly selected) estimate. The estimation errors are computed at each iteration. Since each iteration of the two-subspace method utilizes two rows of the matrix $\bm{A}$ , we call a single iteration of the standard method two iterations in the Algorithm 1.1 for fair comparison.

Figure 3 demonstrates the regime where the two-subspace method offers the most improvement over the standard method. Here the matrix $A$ has highly coherent rows, with $\delta\approx\Delta\approx 1$ .

Our result Theorem 1.1 suggests that as $\delta$ becomes smaller the two-subspace method should offer less and less improvements over the standard method. When $\delta=0$ the convergence rate bound of Theorem 1.1 is precisely the same as that of the standard method (1.1). Indeed, we see this precise behavior as is depicted in Figure 4.

Next we performed experiments on noisy systems. We used the same dimensions and construction of the matrix $\bm{A}$ as well as the signal type. Then we added i.i.d. Gaussian noise with norm $0.1$ to the measurements $\bm{b}$ . Figure 5 demonstrates the exponential convergence of the methods in the presence of noise for various values of $\delta$ and $\Delta$ .

Discussion

As is evident from Theorems 1.1 and 4.3, the two-subspace Kaczmarz method provides exponential convergence in expectation to the solution of $\bm{A}\bm{x}=\bm{b}$ . The constant in the rate of convergence for the two-subspace Kaczmarz method is at most equal to that of the best known results for the randomized Kaczmarz method (1.1). When the matrix $\bm{A}$ has many correlated rows, the constant is significantly lower than that of the standard method, yielding substantially faster convergence. This has positive implications for many applications such as nonuniform sampling in Fourier analysis, as discussed in Section 1.

We emphasize that the bounds presented in our main theorems are weaker than what we actually prove, and that even when $\delta$ is small, if the rows of $\bm{A}$ still have many correlations, Lemmas 2.1 and 4.1 still guarantee improved convergence. For example, if the matrix $\bm{A}$ has correlated rows but contains a pair of identical rows and a pair of orthogonal rows, it will of course be that $\delta=0$ and $\Delta=1$ . However, we see from the proofs of our main theorems that the two-subspace method still guarantees substantial improvement over the standard method. Numerical experiments in cases like this produce results identical to those in Section 5.

It is clear both from the numerical experiments and Theorem 1.1 that the two-subspace Kaczmarz performs best when the correlations $\left\langle\bm{a_{r}},\bm{a_{s}}\right\rangle$ are bounded away from zero. In particular, the larger $\delta$ is the faster the convergence of the two-subspace method. The dependence on $\Delta$ , however, is not as straightforward. Theorems 1.1 and 4.3 suggest that when $\Delta$ is very close to $1$ the two-subspace method should provide similar convergence to the standard method. However, in the experiments of Section 5 we see this is not the case. This dependence on $\Delta$ appears to be only an artifact of the proof.

As is the case for many iterative algorithms, the presence of noise introduces complications both theoretically and empirically. Theorem 3.1 guarantees expected exponential convergence to the noise threshold. For pessimistic values of $D$ , the noise threshold provided by Theorem 3.1 is greater than that of the standard method, (1.2), by a factor of $\sqrt{R}$ . In addition, large values of $\Delta$ produce large error thresholds in this bound. As in the noiseless case, we believe this dependence on $\Delta$ to be an artifact of the proof.

A further and important complication that noise introduces is a semi-convergence effect, a well-known effect in Algebraic Reconstruction Technique (ART) methods (see e.g. ). For example, in Figure 5 (d), the estimation error for the two-subspace method decreases to a point and then begins to increase. It remains an open problem to determine an optimal stopping condition without knowledge of the solution $\bm{x}$ .

2. Future Work

The issue of detecting semiconvergence is a very deep problem. The simple solution would be to terminate the algorithm once the residual $\|\bm{A}\bm{x_{k}}-\bm{b}\|_{2}$ decreases below some threshold. However, the residual decreases in each iteration even when the estimation error begins to increase. Determining the residual threshold beyond which one should terminate is not an easy problem and work in this area continues to be done.

We also hope to improve the error threshold bound of Theorem 3.1 for the two-subspace method. We conjecture that the $(1-\Delta)$ term can be removed or improved, and that the dependence on $R$ can be reduced to $\sqrt{R}$ in the error term of Theorem 3.1.

Finally, a natural extension to our method would be to use more than two rows in each iteration. Indeed, extensions of the two-subspace algorithm to arbitrary subspaces can be analyzed .