Kaczmarz method

The Kaczmarz method or Kaczmarz's algorithm is an iterative algorithm for solving linear equation systems formula_1. It was first discovered by the Polish mathematician Stefan Kaczmarz, and was rediscovered in the field of image reconstruction from projections by Richard Gordon, Robert Bender, and Gabor Herman in 1970, where it is called the Algebraic Reconstruction Technique (ART). ART includes the positivity constraint, making it nonlinear.
The Kaczmarz method is applicable to any linear system of equations, but its computational advantage relative to other methods depends on the system being sparse. It has been demonstrated to be superior, in some biomedical imaging applications, to other methods such as the filtered backprojection method.
It has many applications ranging from computed tomography (CT) to signal processing. It can be obtained also by applying to the hyperplanes, described by the linear system, the method of successive projections onto convex sets (POCS).
Algorithm 1: Kaczmarz algorithm.
Let formula_2 be a system of linear equations, let formula_3 be the number of rows of "A", formula_4 be the formula_5th row of complex-valued matrix formula_6, and let formula_7 be arbitrary complex-valued initial approximation to the solution of formula_8. For formula_9 compute:
where formula_10 and formula_11 denotes complex conjugation of formula_12.
If the system is consistent, formula_13 converges to the minimum-norm solution, provided that the iterations start with the zero vector.
A more general algorithm can be defined using a relaxation parameter formula_14
There are versions of the method that converge to a regularized weighted least squares solution when applied to a system of inconsistent equations and, at least as far as initial behavior is concerned, at a lesser cost than other iterative methods, such as the conjugate gradient method.
Algorithm 2: Randomized Kaczmarz algorithm.
In 2009, a randomized version of the Kaczmarz method for overdetermined linear systems was introduced by Thomas Strohmer and Roman Vershynin in which the "i"-th equation is selected randomly with probability proportional to formula_16
This method can be seen as a particular case of stochastic gradient descent.
Under such circumstances formula_17 converges exponentially fast to the solution of formula_18 and the rate of convergence depends only on the scaled condition number formula_19.
Proof.
We have
Using
we can write () as
The main point of the proof is to view the left hand side in () as an expectation of some random variable. Namely, recall that the solution space of the formula_25 equation of formula_8 is the hyperplane
whose normal is formula_28 Define a random vector "Z" whose values are the normals to all the equations of formula_29, with probabilities as in our algorithm:
Then () says that
The orthogonal projection formula_32 onto the solution space of a random equation of formula_8 is given by formula_34
Now we are ready to analyze our algorithm. We want to show that the error formula_35 reduces at each step in average (conditioned on the previous steps) by at least the factor of formula_36 The next approximation formula_37 is computed from formula_38 as formula_39 where formula_40 are independent realizations of the random projection formula_41 The vector formula_42 is in the kernel of formula_43 It is orthogonal to the solution space of the equation onto which formula_44 projects, which contains the vector formula_45 (recall that formula_46 is the solution to all equations). The orthogonality of these two vectors then yields
To complete the proof, we have to bound formula_48 from below. By the definition of formula_37, we have
where formula_51 are independent realizations of the random vector formula_52
Thus
Now we take the expectation of both sides conditional upon the choice of the random vectors formula_54 (hence we fix the choice of the random projections formula_55 and thus the random vectors formula_56 and we average over the random vector formula_57). Then
By () and the independence,
Taking the full expectation of both sides, we conclude that
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 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 will perform in an inferior manner.
The Kaczmarz iteration () 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 () that any (nonzero) scaling of the equations cancels out. Thus, in RK, one can use formula_61 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 . For additional progress on this topic, see, and the references therein.
Algorithm 3: Gower-Richtarik algorithm.
In 2015, Robert M. Gower and Peter Richtarik developed a versatile randomized iterative method for solving a consistent system of linear equations formula_62 which includes the randomized Kaczmarz algorithm as a special case. Other special cases include randomized coordinate descent, randomized Gaussian descent and randomized Newton method. Block versions and versions with importance sampling of all these methods also arise as special cases. The method is shown to enjoy exponential rate decay (in expectation) - also known as linear convergence, under very mild conditions on the way randomness enters the algorithm. The Gower-Richtarik method is the first algorithm uncovering a "sibling" relationship between these methods, some of which were independently proposed before, while many of which were new.
Insights about Randomized Kaczmarz.
Interesting new insights about the randomized Kaczmarz method that can be gained from the analysis of the method include:
Six Equivalent Formulations.
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):
We now describe some of these viewpoints. The method depends on 2 parameters:
1. Sketch and Project.
Given previous iterate formula_77 the new point formula_78 is computed by drawing a random matrix formula_75 (in an iid fashion from some fixed distribution), and setting
That is, formula_78 is obtained as the projection of formula_82 onto the randomly sketched system formula_83. The idea behind this method is to pick formula_75 in such a way that a projection onto the sketched system is substantially simpler than the solution of the original system formula_85. Randomized Kaczmarz method is obtained by picking formula_72 to be the identity matrix, and formula_75 to be the formula_88 unit coordinate vector with probability formula_89 Different choices of formula_72 and formula_75 lead to different variants of the method.
2. Constrain and Approximate.
A seemingly different but entirely equivalent formulation of the method (obtained via Lagrangian duality) is
where formula_93 is also allowed to vary, and where formula_94 is any solution of the system formula_21 Hence, formula_78 is obtained by first constraining the update to the linear subspace spanned by the columns of the random matrix formula_97, i.e., to
and then choosing the point formula_20 from this subspace which best approximates formula_100. This formulation may look surprising as it seems impossible to perform the approximation step due to the fact that formula_94 is not known (after all, this is what we are trying the compute!). However, it is still possible to do this, simply because formula_78 computed this way is the same as formula_78 computed via the sketch and project formulation and since formula_94 does not appear there.
5. Random Update.
The update can also be written explicitly as
where by formula_106 we denote the Moore-Penrose pseudo-inverse of matrix formula_107. Hence, the method can be written in the form formula_108, where formula_109 is a random update vector.
Letting formula_110 it can be shown that the system formula_111 always has a solution formula_112, and that for all such solutions the vector formula_113 is the same. Hence, it does not matter which of these solutions is chosen, and the method can be also written as formula_114. The pseudo-inverse leads just to one particular solution. The role of the pseudo-inverse is twofold:
6. Random Fixed Point.
If we subtract formula_94 from both sides of the random update formula, denote
and use the fact that formula_117 we arrive at the last formulation:
where formula_119 is the identity matrix. The iteration matrix, formula_120 is random, whence the name of this formulation.
Convergence.
By taking conditional expectations in the 6th formulation (conditional on formula_121), we obtain
By taking expectation again, and using the tower property of expectations, we obtain
Gower and Richtarik show that
where the matrix norm is defined by
Moreover, without any assumptions on formula_75 one has formula_127 By taking norms and unrolling the recurrence, we obtain
Theorem [Gower & Richtarik 2015].
"Remark". A sufficient condition for the expected residuals to converge to 0 is formula_129 This can be achieved if formula_6 has a full column rank and under very mild conditions on formula_131 Convergence of the method can be established also without the full column rank assumption in a different way.
It is also possible to show a stronger result:
Theorem [Gower & Richtarik 2015].
The expected squared norms (rather than norms of expectations) converge at the same rate:
"Remark". This second type of convergence is stronger due to the following identity which holds for any random vector formula_20 and any fixed vector formula_94:
Convergence of Randomized Kaczmarz.
We have seen that the randomized Kaczmarz method appears as a special case of the Gower-Richtarik method for formula_136 and formula_75 being the formula_88 unit coordinate vector with probability formula_139 where formula_12 is the formula_88 row of formula_142 It can be checked by direct calculation that