Kaczmarz method

We have

${\displaystyle \forall z\in \mathbb {C} ^{n}:\quad \sum _{j=1}^{m}|\langle z,a_{j}\rangle |^{2}\geq {\frac {\|z\|^{2}}{\|A^{-1}\|^{2}}}}$

(2)

Using

${\displaystyle \|A\|^{2}=\sum _{j=1}^{m}\|a_{j}\|^{2}}$

we can write (2) as

${\displaystyle \forall z\in \mathbb {C} ^{n}:\quad \sum _{j=1}^{m}{\frac {\|a_{j}\|^{2}}{\|A\|^{2}}}\left|\left\langle z,{\frac {a_{j}}{\|a_{j}\|}}\right\rangle \right|^{2}\geq \kappa (A)^{-2}{\|z\|^{2}}}$

(3)

The main point of the proof is to view the left hand side in (3) as an expectation of some random variable. Namely, recall that the solution space of the ${\displaystyle j-th}$ equation of ${\displaystyle Ax=b}$ is the hyperplane

${\displaystyle \{y:\langle y,a_{j}\rangle =b_{j}\},}$

whose normal is ${\displaystyle {\tfrac {a_{j}}{\|a_{j}\|^{2}}}.}$ Define a random vector Z whose values are the normals to all the equations of ${\displaystyle Ax=b}$, with probabilities as in our algorithm:

${\displaystyle Z={\frac {a_{j}}{\|a_{j}\|}}}$ with probability ${\displaystyle {\frac {\|a_{j}\|^{2}}{\|A\|^{2}}}\qquad \qquad \qquad j=1,\ldots ,m}$

Then (3) says that

${\displaystyle \forall z\in \mathbb {C} ^{n}:\quad \mathbb {E} |\langle z,Z\rangle |^{2}\geq \kappa (A)^{-2}{\|z\|^{2}}}$

(4)

The orthogonal projection ${\displaystyle P}$ onto the solution space of a random equation of ${\displaystyle Ax=b}$ is given by ${\displaystyle Pz=z-\langle z-x,Z\rangle Z.}$

Now we are ready to analyze our algorithm. We want to show that the error ${\displaystyle {\|x_{k}-x\|^{2}}}$ reduces at each step in average (conditioned on the previous steps) by at least the factor of ${\displaystyle (1-\kappa (A)^{-2}).}$ The next approximation ${\displaystyle x_{k}}$ is computed from ${\displaystyle x_{k-1}}$ as ${\displaystyle x_{k}=P_{k}x_{k-1},}$ where ${\displaystyle P_{1},P_{2},\ldots }$ are independent realizations of the random projection ${\displaystyle P.}$ The vector ${\displaystyle x_{k-1}-x_{k}}$ is in the kernel of ${\displaystyle P_{k}.}$ It is orthogonal to the solution space of the equation onto which ${\displaystyle P_{k}}$ projects, which contains the vector ${\displaystyle x_{k}-x}$ (recall that ${\displaystyle x}$ is the solution to all equations). The orthogonality of these two vectors then yields

${\displaystyle \|x_{k}-x\|^{2}=\|x_{k-1}-x\|^{2}-\|x_{k-1}-x_{k}\|^{2}.}$

To complete the proof, we have to bound ${\displaystyle \|x_{k-1}-x_{k}\|^{2}}$ from below. By the definition of ${\displaystyle x_{k}}$, we have

${\displaystyle \|x_{k-1}-x_{k}\|=\langle x_{k-1}-x,Z_{k}\rangle }$

where ${\displaystyle Z_{1},Z_{2},\ldots }$ are independent realizations of the random vector ${\displaystyle Z.}$

Thus

${\displaystyle \|x_{k}-x\|^{2}\leq \left(1-\left|\left\langle {\frac {x_{k-1}-x}{\|x_{k-1}-x\|}},Z_{k}\right\rangle \right|^{2}\right){\|x_{k-1}-x\|^{2}}.}$

Now we take the expectation of both sides conditional upon the choice of the random vectors ${\displaystyle Z_{1},\ldots ,Z_{k-1}}$ (hence we fix the choice of the random projections ${\displaystyle P_{1},\ldots ,P_{k-1}}$ and thus the random vectors ${\displaystyle x_{1},\ldots ,x_{k-1}}$ and we average over the random vector ${\displaystyle Z_{k}}$). Then

${\displaystyle \mathbb {E} _{Z_{1},\ldots ,Z_{k-1}}{\|x_{k}-x\|^{2}}=\left(1-\mathbb {E} _{Z_{1},\ldots ,Z_{k-1},Z_{k}}\left|\left\langle {\frac {x_{k-1}-x}{\|x_{k-1}-x\|}},Z_{k}\right\rangle \right|^{2}\right){\|x_{k-1}-x\|^{2}}.}$

By (4) and the independence,

${\displaystyle \mathbb {E} _{Z_{1},\ldots ,Z_{k-1}}{\|x_{k}-x\|^{2}}\leq (1-\kappa (A)^{-2}){\|x_{k-1}-x\|^{2}}.}$

Taking the full expectation of both sides, we conclude that

${\displaystyle \mathbb {E} \|x_{k}-x\|^{2}\leq (1-\kappa (A)^{-2})\mathbb {E} {\|x_{k-1}-x\|^{2}}.\blacksquare }$

The superiority of this selection was illustrated with the reconstruction of a bandlimited function from its nonuniformly spaced sampling values. However, it has been pointed out[10] that the reported success by Strohmer and Vershynin depends on the specific choices that were made there in translating the underlying problem, whose geometrical nature is to find a common point of a set of hyperplanes, into a system of algebraic equations. There will always be legitimate algebraic representations of the underlying problem for which the selection method in[8] will perform in an inferior manner.[8][10][11]

The Kaczmarz iteration (1) has a purely geometric interpretation: the algorithm successively projects the current iterate onto the hyperplane defined by the next equation. Hence, any scaling of the equations is irrelevant; it can also be seen from (1) that any (nonzero) scaling of the equations cancels out. Thus, in RK, one can use ${\displaystyle \|a_{i}\|}$ or any other weights that may be relevant. Specifically, in the above-mentioned reconstruction example, the equations were chosen with probability proportional to the average distance of each sample point from its two nearest neighbors — a concept introduced by Feichtinger and Gröchenig. For additional progress on this topic, see,[12][13] and the references therein.

Interesting new insights about the randomized Kaczmarz method that can be gained from the analysis of the method include:

• The general rate of the Gower-Richtarik algorithm precisely recovers the rate of the randomized Kaczmarz method in the special case when it reduced to it.
• The choice of probabilities for which the randomized Kaczmarz algorithm was originally formulated and analyzed (probabilities proportional to the squares of the row norms) is not optimal. Optimal probabilities are the solution of a certain semidefinite program. The theoretical complexity of randomized Kaczmarz with the optimal probabilities can be arbitrarily better than the complexity for the standard probabilities. However, the amount by which it is better depends on the matrix ${\displaystyle A}$. There are problems for which the standard probabilities are optimal.
• When applied to a system with matrix ${\displaystyle A}$ which is positive definite, Randomized Kaczmarz method is equivalent to the Stochastic Gradient Descent (SGD) method (with a very special stepsize) for minimizing the strongly convex quadratic function ${\displaystyle f(x)={\tfrac {1}{2}}x^{T}Ax-b^{T}x.}$ Note that since ${\displaystyle f}$ is convex, the minimizers of ${\displaystyle f}$ must satisfy ${\displaystyle \nabla f(x)=0}$, which is equivalent to ${\displaystyle Ax=b.}$ The "special stepsize" is the stepsize which leads to a point which in the one-dimensional line spanned by the stochastic gradient minimizes the Euclidean distance from the unknown(!) minimizer of ${\displaystyle f}$, namely, from ${\displaystyle x^{*}=A^{-1}b.}$ This insight is gained from a dual view of the iterative process (below described as "Optimization Viewpoint: Constrain and Approximate").

The Gower-Richtarik method enjoys six seemingly different but equivalent formulations, shedding additional light on how to interpret it (and, as a consequence, how to interpret its many variants, including randomized Kaczmarz):

• 1. Sketching viewpoint: Sketch & Project
• 2. Optimization viewpoint: Constrain and Approximate
• 3. Geometric viewpoint: Random Intersect
• 4. Algebraic viewpoint 1: Random Linear Solve
• 5. Algebraic viewpoint 2: Random Update
• 6. Analytic viewpoint: Random Fixed Point

We now describe some of these viewpoints. The method depends on 2 parameters:

• a positive definite matrix ${\displaystyle B}$ giving rise to a weighted Euclidean inner product ${\displaystyle \langle x,y\rangle _{B}:=x^{T}By}$ and the induced norm

${\displaystyle \|x\|_{B}=\left(\langle x,x\rangle _{B}\right)^{\frac {1}{2}},}$

• and a random matrix ${\displaystyle S}$ with as many rows as ${\displaystyle A}$ (and possibly random number of columns).

Given previous iterate ${\displaystyle x^{k},}$ the new point ${\displaystyle x^{k+1}}$ is computed by drawing a random matrix ${\displaystyle S}$ (in an iid fashion from some fixed distribution), and setting

${\displaystyle x^{k+1}={\underset {x}{\operatorname {arg\ min} }}\|x-x^{k}\|_{B}{\text{ subject to }}S^{T}Ax=S^{T}b.}$

That is, ${\displaystyle x^{k+1}}$ is obtained as the projection of ${\displaystyle x^{k}}$ onto the randomly sketched system ${\displaystyle S^{T}Ax=S^{T}b}$. The idea behind this method is to pick ${\displaystyle S}$ in such a way that a projection onto the sketched system is substantially simpler than the solution of the original system ${\displaystyle Ax=b}$. Randomized Kaczmarz method is obtained by picking ${\displaystyle B}$ to be the identity matrix, and ${\displaystyle S}$ to be the ${\displaystyle i^{th}}$ unit coordinate vector with probability ${\displaystyle p_{i}=\|a_{i}\|_{2}^{2}/\|A\|_{F}^{2}.}$ Different choices of ${\displaystyle B}$ and ${\displaystyle S}$ lead to different variants of the method.

A seemingly different but entirely equivalent formulation of the method (obtained via Lagrangian duality) is

${\displaystyle x^{k+1}={\underset {x}{\operatorname {arg\ min} }}\left\|x-x^{*}\right\|_{B}{\text{ subject to }}x=x^{k}+B^{-1}A^{T}Sy,}$

where ${\displaystyle y}$ is also allowed to vary, and where ${\displaystyle x^{*}}$ is any solution of the system ${\displaystyle Ax=b.}$ Hence, ${\displaystyle x^{k+1}}$ is obtained by first constraining the update to the linear subspace spanned by the columns of the random matrix ${\displaystyle B^{-1}A^{T}S}$, i.e., to

${\displaystyle \left\{h:h=B^{-1}A^{T}Sy,\quad y{\text{ can vary }}\right\},}$

and then choosing the point ${\displaystyle x}$ from this subspace which best approximates ${\displaystyle x^{*}}$. This formulation may look surprising as it seems impossible to perform the approximation step due to the fact that ${\displaystyle x^{*}}$ is not known (after all, this is what we are trying the compute!). However, it is still possible to do this, simply because ${\displaystyle x^{k+1}}$ computed this way is the same as ${\displaystyle x^{k+1}}$ computed via the sketch and project formulation and since ${\displaystyle x^{*}}$ does not appear there.

The update can also be written explicitly as

${\displaystyle x^{k+1}=x^{k}-B^{-1}A^{T}S\left(S^{T}AB^{-1}A^{T}S\right)^{\dagger }S^{T}\left(Ax^{k}-b\right),}$

where by ${\displaystyle M^{\dagger }}$ we denote the Moore-Penrose pseudo-inverse of matrix ${\displaystyle M}$. Hence, the method can be written in the form ${\displaystyle x^{k+1}=x^{k}+h^{k}}$, where ${\displaystyle h^{k}}$ is a random update vector.

Letting ${\displaystyle M=S^{T}AB^{-1}A^{T}S,}$ it can be shown that the system ${\displaystyle My=S^{T}(Ax^{k}-b)}$ always has a solution ${\displaystyle y^{k}}$, and that for all such solutions the vector ${\displaystyle x^{k+1}-B^{-1}A^{T}Sy^{k}}$ is the same. Hence, it does not matter which of these solutions is chosen, and the method can be also written as ${\displaystyle x^{k+1}=x^{k}-B^{-1}A^{T}Sy^{k}}$. The pseudo-inverse leads just to one particular solution. The role of the pseudo-inverse is twofold:

• It allows the method to be written in the explicit "random update" form as above,
• It makes the analysis simple through the final, sixth, formulation.

If we subtract ${\displaystyle x^{*}}$ from both sides of the random update formula, denote

${\displaystyle Z:=A^{T}S\left(S^{T}AB^{-1}A^{T}S\right)^{\dagger }S^{T}A,}$

and use the fact that ${\displaystyle Ax^{*}=b,}$ we arrive at the last formulation:

${\displaystyle x^{k+1}-x^{*}=\left(I-B^{-1}Z\right)\left(x^{k}-x^{*}\right),}$

where ${\displaystyle I}$ is the identity matrix. The iteration matrix, ${\displaystyle I-B^{-1}Z,}$ is random, whence the name of this formulation.

By taking conditional expectations in the 6th formulation (conditional on ${\displaystyle x^{k}}$), we obtain

${\displaystyle \mathbb {E} \left.\left[x^{k+1}-x^{*}\right|x^{k}\right]=\left(I-B^{-1}\mathbb {E} [Z]\right)\left[x^{k}-x^{*}\right].}$

By taking expectation again, and using the tower property of expectations, we obtain

${\displaystyle \mathbb {E} \left[x^{k+1}-x^{*}\right]=(I-B^{-1}\mathbb {E} [Z])\mathbb {E} \left[x^{k}-x^{*}\right].}$

Gower and Richtarik[14] show that

${\displaystyle \rho$ :=\left\|I-B^{-{\frac {1}{2}}}\mathbb {E} [Z]B^{-{\frac {1}{2}}}\right\|_{B}=\lambda _{\max }\left(I-B^{-1}\mathbb {E} [Z]\right),}

where the matrix norm is defined by

${\displaystyle \|M\|_{B}:=\max _{x\neq 0}{\frac {\|Mx\|_{B}}{\|x\|_{B}}}.}$

Moreover, without any assumptions on ${\displaystyle S}$ one has ${\displaystyle 0\leq \rho \leq 1.}$ By taking norms and unrolling the recurrence, we obtain

${\displaystyle \left\|\mathbb {E} \left[x^{k}-x^{*}\right]\right\|_{B}\leq \rho ^{k}\|x^{0}-x^{*}\|_{B}.}$

Remark. A sufficient condition for the expected residuals to converge to 0 is ${\displaystyle \rho <1.}$ This can be achieved if ${\displaystyle A}$ has a full column rank and under very mild conditions on ${\displaystyle S.}$ Convergence of the method can be established also without the full column rank assumption in a different way.[15]

It is also possible to show a stronger result:

The expected squared norms (rather than norms of expectations) converge at the same rate:

${\displaystyle \mathbb {E} \left\|\left[x^{k}-x^{*}\right]\right\|_{B}^{2}\leq \rho ^{k}\left\|x^{0}-x^{*}\right\|_{B}^{2}.}$

Remark. This second type of convergence is stronger due to the following identity[14] which holds for any random vector ${\displaystyle x}$ and any fixed vector ${\displaystyle x^{*}}$:

${\displaystyle \left\|\mathbb {E} \left[x-x^{*}\right]\right\|^{2}=\mathbb {E} \left[\left\|x-x^{*}\right\|^{2}\right]-\mathbb {E} \left[\|x-\mathbb {E} [x]\|^{2}\right].}$

We have seen that the randomized Kaczmarz method appears as a special case of the Gower-Richtarik method for ${\displaystyle B=I}$ and ${\displaystyle S}$ being the ${\displaystyle i^{th}}$ unit coordinate vector with probability ${\displaystyle p_{i}=\|a_{i}\|_{2}^{2}/\|A\|_{F}^{2},}$ where ${\displaystyle a_{i}}$ is the ${\displaystyle i^{th}}$ row of ${\displaystyle A.}$ It can be checked by direct calculation that

${\displaystyle \rho =\|I-B^{-1}\mathbb {E} [Z]\|_{B}=1-{\frac {\lambda _{\min }(A^{T}A)}{\|A\|_{F}^{2}}}.}$