Randomized Kaczmarz solver for noisy linear systems

Introduction

The Kaczmarz method is one of the most popular solvers of overdetermined linear systems and has numerous applications from computer tomography to image processing. It is an iterative method, and so therefore is practical in the realm of very large systems of equations. The algorithm consists of a series of alternating projections, and is often considered a type of Projection on Convex Sets (POCS) method. Given a consistent system of linear equations of the form

where $x_{k}$ is the $k^{th}$ iterate and $i=(k$ mod $m)+1$.

Although the Kaczmarz method is popular in practice, theoretical results on the convergence rate of the method have been difficult to obtain. Most known estimates depend on properties of the matrix $A$ which may be time consuming to compute, and are not easily comparable to those of other iterative methods (see e.g. , , ). Since the Kaczmarz method cycles through the rows of $A$ sequentially, its convergence rate depends on the order of the rows. Intuition tells us that the order of the rows of $A$ does not change the difficulty level of the system as a whole, so one would hope for results independent of the ordering. One natural way to overcome this is to use the rows of $A$ in a random order, rather than sequentially. Several observations were made on the improvements of this randomized version , but only recently have theoretical results been obtained .

In designing a random version of the Kaczmarz method, it is necessary to set the probability of each row being selected. Strohmer and Vershynin propose in to set the probability proportional to the Euclidean norm of the row. Their revised algorithm can then be described by the following:

In , Strohmer and Vershynin prove the following exponential bound on the expected rate of convergence for the randomized Kaczmarz method,

The first remarkable note about this result is that it is essentially independent of the number $m$ of equations in the system. Indeed, by the definition of $R$, $R$ is proportional to $n$ within a square factor of $\kappa(A)$, the condition number of $A$ ($\kappa(A)$ is defined as the ratio of the largest to smallest singular values of $A$). This bound also demonstrates, however, that the Kaczmarz method is an efficient alternative to other methods only when the condition number is very small. If this is not the case, then other alternative methods may offer improvements in practice.

Since the results of , there has been some further discussion about the benefits of this randomized version of the Kaczmarz method (see ). The Kaczmarz method has been studied for over seventy years, and is useful in many applications. The notion of selecting the rows randomly in the method has been proposed before (see ), and improvements over the standard method were observed. However, the work by Strohmer and Vershynin in provides the first proof on the rate of convergence. The rate is exponential in expectation and is in terms of standard matrix properties. We are not aware of any other Kaczmarz method that provably achieves exponential convergence.

It is important to note that the method of row selection proposed in this version of the randomized Kaczmarz method is not optimal, and an example that demonstrates this is given in . However, under this selection strategy, the convergence rates proven in are optimal, and there are matrices that satisfy the proven bounds exactly. The selection strategy in this method was chosen because it often yields very good results, allows a provable guarantee of exponential convergence, and is computationally efficient.

Since the algorithm selects rows based on their row norms, it is natural to ask whether one can simply scale the rows any way one wishes. Indeed, choosing the rows based on their norms is related to the notion of applying a diagonal preconditioner. However, since finding the optimal diagonal preconditioner for a system $Ax=b$ is itself a task that is often more costly than inverting the entire matrix, we select an easier, although not optimal, preconditioner that simply scales by the (square of the) row norms. This type of preconditioner yields a balance of computational cost and optimality (see ). The distinction between the effect of an alternative diagonal preconditioner on the Kaczmarz method versus the randomized method discussed here is important. If the system is multiplied by a diagonal matrix, the standard Kaczmarz method will not change, since the angles between all rows do not change. However, such a multiplication to the system in our randomized setting changes the probabilities of selecting the rows (by definition). It is then not a surprise that this will also affect the convergence rate proved for this method (since multiplication will affect the value of $R$ in (1.1)).

This randomized version of the Kaczmarz method provides clear advantages over the standard method in many cases. Using the selection strategy above, Strohmer and Vershynin were able to provide a proof for the expected rate of convergence that shows exponential convergence. No such convergence rate for any Kaczmarz method has been proven before. These benefits lead one to question whether the method works in the more realistic case where the system is corrupted by noise. In this paper we provide theoretical and empirical results to suggest that in this noisy case the method converges exponentially to the solution within a specified error bound. The error bound is proportional to $\sqrt{R}$, and we also provide a simple example showing this bound is sharp in the general setting.

Main Results

Theoretical and empirical studies have shown the randomized Kaczmarz algorithm to provide very promising results. Here we show that it also performs well in the case where the system is corrupted with noise. In this section we consider the consistent system $Ax=b$ after an error vector $r$ is added to the right side:

Note that we do not require the perturbed system to be consistent. First we present a simple example to gain intuition about how drastically the noise can affect the system. To that end, let $A$ be the $n\times n$ identity matrix, $b=0$, and suppose the error is the vector whose entries are all one, $r=(1,1,\ldots,1)$. Then the solution to the noisy system is clearly $x=r=(1,1,\ldots,1)$, and the solution to the unperturbed problem is $x=0$. By Jensen’s inequality, we have

Now considering the noisy problem, we may substitute $r$ for $x$ in (1.1). Combining this with Jensen’s inequality above, we obtain

Next, by taking expectation and using (2.1) above, we have

Finally by the definition of $r$ and $R$, this implies

This means that the limiting error between the iterates $x_{k}$ and the original solution $x$ is $\sqrt{R}$. In it is shown that the bound provided in (1.1) is optimal, so even this trivial example demonstrates that if we wish to maintain a general setting, the best error bound for the noisy case we can hope for is proportional to $\sqrt{R}$. Our main result proves this exact theoretical bound.

Let $A$ have full column rank and assume the system $Ax=b$ is consistent. Let $x_{k}^{*}$ be the $k^{th}$ iterate of the noisy randomized Kaczmarz method run with $Ax\approx b+r$, and let $a_{1},\ldots a_{m}$ denote the rows of $A$. Then we have

where $R=\|A^{-1}\|^{2}\|A\|_{F}^{2}$, $\gamma=\max_{i}\frac{|r_{i}|}{\|a_{i}\|_{2}}$, and the expectation is taken over the choice of the rows in the algorithm.

In the case discussed above, note that we have $\gamma=1$, so the example indeed shows the bound is sharp.

By applying the bound $\sqrt{R}\leq\kappa(A)\sqrt{n}$ to Theorem 2.1 above, we obtain the bound

These bounds look similar in spirit, providing some more reassurance to the sharpness of the error bound. It is important to note though that the first is obtained by applying the left inverse rather than an iterative method, which explains why the bounds are not exactly equal. Of course for problems of large sizes, applying the inverse may not even be computationally feasible.

Before proving the theorem, it is important to first analyze what happens to the solution spaces of the original equations $Ax=b$ when the error vector is added. Letting $a_{1},\ldots a_{m}$ denote the rows of $A$, we have that each solution space $\langle{a_{i}},{x}\rangle=b_{i}$ of the original system is a hyperplane whose normal is $\frac{a_{i}}{\|a_{i}\|_{2}}$. When noise is added, each hyperplane is translated in the direction of $a_{i}$. Thus the new geometry consists of hyperplanes parallel to those in the noiseless case. A simple computation provides the following lemma which specifies exactly how far each hyperplane is shifted.

Note that this lemma does not imply that the noisy system is consistent. By definition of $H_{i}^{*}$ it is clear that each subspace is non-empty, but we are not requiring that the intersection of all $H_{i}^{*}$ be non-empty.

First, if $w\in H_{i}$ then $\langle{a_{i}},{w+\alpha a_{i}}\rangle=\langle{a_{i}},{w}\rangle+\alpha\|a_{i}\|_{2}^{2}=b_{i}+r_{i}$, so $w+\alpha a_{i}\in H_{i}^{*}$. Next let $u\in H_{i}^{*}$. Set $w=u-\alpha a_{i}$. Then $\langle{a_{i}},{w}\rangle=\langle{a_{i}},{u}\rangle-r_{i}=b_{i}+r_{i}-r_{i}=b_{i}$, so $w\in H_{i}^{*}$. This completes the proof. ∎

We will also utilize the following lemma which is proved in the proof of Theorem 2 in .

where $R=\|A^{-1}\|^{2}\|A\|_{F}^{2}$, and the expectation is taken over the choice of the rows in the algorithm.

We are now prepared to prove Theorem 2.1.

Let $x_{k-1}^{*}$ denote the $(k-1)^{th}$ iterate of noisy randomized Kaczmarz. Using notation as in Lemma 2.2, let $H_{i}^{*}$ be the solution space chosen in the $k^{th}$ iteration. Then $x_{k}^{*}$ is the orthogonal projection of $x_{k-1}^{*}$ onto $H_{i}^{*}$. Let $x_{k}$ denote the orthogonal projection of $x_{k-1}^{*}$ onto $H_{i}$ (see Figure 1).

By Lemma 2.2 and the fact that $a_{i}$ is orthogonal to $H_{i}$ and $H_{i}^{*}$, we have that $x_{k}^{*}-x=x_{k}-x+\alpha_{i}a_{i}$. Again by orthogonality, we have $\|x_{k}^{*}-x\|_{2}^{2}=\|x_{k}-x\|_{2}^{2}+\|\alpha_{i}a_{i}\|_{2}^{2}$. Then by Lemma 2.3 and the definition of $\gamma$, we have

where the expectation is conditioned upon the choice of the random selections in the first $k-1$ iterations. Then applying this recursive relation iteratively and taking full expectation, we have

Numerical Examples

In this section we describe some of our numerical results for the randomized Kaczmarz method in the case of noisy systems. Figure 2 depicts the error between the estimate by randomized Kaczmarz and the actual signal, in comparison with the predicted threshold value for several types of matrices. The first study was conducted for 100 trials using $2000\times 100$ Gaussian matrices (matrices who entries are i.i.d. Gaussian with mean and variance $1$) and independent Gaussian noise of norm $0.02$. The systems were homogeneous, meaning $x=0$ and $b=0$. The thick line is a plot of the threshold value, $\gamma\sqrt{R}$ for each trial. The thin line is a plot of the error in the estimate after the given amount of iterations for the corresponding trial. The scatter plot displays the convergence of the method over several randomly chosen trials from this study, and clearly shows exponential convergence. The second study is a similar study but for the experiments in which we used partial Fourier matrices. In this case we use $m=700$ and $n=101$. For $j=1\ldots 700$ and $k=-50\ldots 50$, we set $A_{j,k}=\exp(2\pi ikt_{j})$, where $t_{j}$ are generated uniformly at random on $$. This type of generation is used to create nonuniformly spaced sampling values, and is used in many applications in signal processing, such as in the reconstruction of bandlimited signals. The third study is similar but used matrices whose entries are Bernoulli ($0/1$each with probability$0.5$). All of these experiments were conducted to demonstrate that the error found in practice is close to that predicted by the theoretical results. As is evident by the plots, the error is quite close to the threshold in all cases.

I would like to thank Roman Vershynin for suggestions that simplified the proofs and for many thoughtful discussions. I would also like to thank Thomas Strohmer for his very appreciated guidance.

References