Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Suppose we wish to solve a consistent linear system of equations; i.e., given ARn×d\vec{A}\in\R^{n\times d} and bRn\vec{b}\in\R^n, find x\vec{x}^* so that

Ax=b.\vec{A}\vec{x}^* = \vec{b}.

Each equation ajTx=bj\vec{a}_j^\T \vec{x} = b_j, where ajT\vec{a}_j^\T is the jj-th row of A\vec{A}, determines a hyperplane

Hj:={xRd:ajTx=bj},H_j := \{ \vec{x}\in\R^d : \vec{a}_j^\T \vec{x} = b_j \},

and the solution set of (7.5) is the intersection H1HnH_1 \cap \cdots \cap H_n.

The Kaczmarz method builds an approximate solution by repeatedly enforcing one equation at a time. Given a current iterate xi1\vec{x}_{i-1} and a row index ss, the Kaczmarz update orthogonally projects xi1\vec{x}_{i-1} onto the hyperplane HsH_s:

xi=xi1+bsasTxi1as2as.\vec{x}_{i} = \vec{x}_{i-1} + \frac{b_{s} - \vec{a}_{s}^\T \vec{x}_{i-1}}{\|\vec{a}_{s}\|^2} \, \vec{a}_{s}.

After the update, the ss-th equation is satisfied exactly (although it may be violated again by later updates). Note that the update (7.6) reads only a single row of A\vec{A} and a single entry of b\vec{b}, and costs just O(d)O(d) arithmetic operations. This makes Kaczmarz-type methods attractive when nn is very large or when we can only access A\vec{A} one row at a time (e.g. streaming or out-of-core settings).

The classical Kaczmarz method cycles through the rows in order s=1,2,,n,1,2,s = 1, 2, \ldots, n, 1, 2, \ldots, but its convergence behavior depends on the ordering of the rows and is difficult to analyze. The randomized Kaczmarz method of Strohmer & Vershynin, 2009 instead samples rows randomly, proportional to their squared row-norms. This is exactly row-norm sampling, which we encountered as a subsampling distribution in Chapter 2.

Analysis

The key result of Strohmer & Vershynin, 2009 is that randomized Kaczmarz converges linearly in expectation, at a rate governed by the ratio σmin(A)2/AF2\smin(\vec{A})^2 / \|\vec{A}\|_\F^2.

Proof

Write ei:=xix\vec{e}_{i} := \vec{x}_{i} - \vec{x}^* for the error at step ii, and let s=sis = s_i be the row sampled at step ii. Since the system is consistent, bs=asTxb_{s} = \vec{a}_{s}^\T \vec{x}^*, so the update (7.6) gives

ei=ei1asTei1as2as=(IasasTas2)ei1.\vec{e}_{i} = \vec{e}_{i-1} - \frac{\vec{a}_{s}^\T \vec{e}_{i-1}}{\|\vec{a}_{s}\|^2} \, \vec{a}_{s} = \left( \vec{I} - \frac{\vec{a}_{s}\vec{a}_{s}^\T}{\|\vec{a}_{s}\|^2} \right) \vec{e}_{i-1}.

The matrix in parentheses is an orthogonal projection, so by the Pythagorean theorem,

ei2=ei12(asTei1)2as2.\| \vec{e}_{i} \|^2 = \| \vec{e}_{i-1} \|^2 - \frac{(\vec{a}_{s}^\T \vec{e}_{i-1})^2}{\|\vec{a}_{s}\|^2}.

Now take the expectation over the choice of ss, conditioned on xi1\vec{x}_{i-1}. Since P[s=j]=aj2/AF2\PP[s = j] = \|\vec{a}_j\|^2 / \|\vec{A}\|_\F^2,

E[(asTei1)2as2|xi1]=j=1naj2AF2(ajTei1)2aj2=Aei12AF2σmin(A)2AF2ei12,\EE\left[ \frac{(\vec{a}_{s}^\T \vec{e}_{i-1})^2}{\|\vec{a}_{s}\|^2} \,\middle|\, \vec{x}_{i-1} \right] = \sum_{j=1}^{n} \frac{\|\vec{a}_j\|^2}{\|\vec{A}\|_\F^2} \cdot \frac{(\vec{a}_{j}^\T \vec{e}_{i-1})^2}{\|\vec{a}_{j}\|^2} = \frac{\|\vec{A}\vec{e}_{i-1}\|^2}{\|\vec{A}\|_\F^2} \geq \frac{\smin(\vec{A})^2}{\|\vec{A}\|_\F^2} \|\vec{e}_{i-1}\|^2,

where the last inequality uses that A\vec{A} has full column rank. Therefore,

E[ei2xi1](1σmin(A)2AF2)ei12.\EE\big[ \|\vec{e}_{i}\|^2 \,\big|\, \vec{x}_{i-1} \big] \leq \left( 1 - \frac{\smin(\vec{A})^2}{\|\vec{A}\|_\F^2} \right) \|\vec{e}_{i-1}\|^2.

The result follows by iterating this bound and using the law of total expectation.

Several remarks are in order:

  • The quantity AF2/σmin(A)2\|\vec{A}\|_\F^2 / \smin(\vec{A})^2 is sometimes called the Demmel condition number (squared) of A\vec{A}. Since AF2dσmax(A)2\|\vec{A}\|_\F^2 \leq d \, \smax(\vec{A})^2, Theorem 7.3 implies that

t=O(dcond(A)2log(1/ε))t = O\big( d \cond(\vec{A})^2 \log(1/\varepsilon) \big)

iterations suffice to guarantee E[xtx2]εx0x2\EE[\|\vec{x}_t - \vec{x}^*\|^2] \leq \varepsilon \|\vec{x}_0 - \vec{x}^*\|^2. Each iteration costs O(d)O(d), so the total cost is O(d2cond(A)2log(1/ε))O(d^2 \cond(\vec{A})^2 \log(1/\varepsilon)) operations — potentially sublinear in the size ndnd of A\vec{A} when ndn \gg d.

  • The rate is independent of nn; the number of rows enters only through how the sampling probabilities are defined.

  • The proof reveals why row-norm sampling is the natural choice: the as2\|\vec{a}_s\|^2 appearing in the sampling probability exactly cancels the as2\|\vec{a}_s\|^{-2} in the projection, so the expected per-step progress becomes Ae2/AF2\|\vec{A}\vec{e}\|^2 / \|\vec{A}\|_\F^2, which can be bounded using the singular values of A\vec{A}. This is the same cancellation that made row-norm sampling effective for approximate matrix multiplication.

Connection to stochastic gradient descent

Randomized Kaczmarz is closely related to stochastic gradient descent (SGD), the workhorse algorithm of modern machine learning. Consider the least-squares objective, written as a sum over the rows of A\vec{A}:

f(x)=12bAx2=j=1nfj(x),fj(x):=12(bjajTx)2.f(\vec{x}) = \frac{1}{2} \| \vec{b} - \vec{A}\vec{x} \|^2 = \sum_{j=1}^{n} f_j(\vec{x}), \qquad f_j(\vec{x}) := \frac{1}{2} \big( b_j - \vec{a}_j^\T \vec{x} \big)^2.

Computing the full gradient f(x)\nabla f(\vec{x}) requires reading all of A\vec{A}. Instead, SGD samples a random index ss from some distribution (p1,,pn)(p_1, \ldots, p_n) and takes a step in the direction of an unbiased estimate of the gradient:

xi=xi1ηpsfs(xi1)=xi1+ηps(bsasTxi1)as,\vec{x}_{i} = \vec{x}_{i-1} - \frac{\eta}{p_{s}} \nabla f_{s}(\vec{x}_{i-1}) = \vec{x}_{i-1} + \frac{\eta}{p_{s}} \big( b_{s} - \vec{a}_{s}^\T \vec{x}_{i-1} \big) \vec{a}_{s},

where η>0\eta > 0 is a step size. Indeed, E[ps1fs(x)]=jfj(x)=f(x)\EE[ p_s^{-1} \nabla f_s(\vec{x}) ] = \sum_j \nabla f_j(\vec{x}) = \nabla f(\vec{x}), regardless of the choice of sampling distribution.

Comparing (7.7) with (7.6), we see that randomized Kaczmarz is exactly SGD with row-norm importance sampling pj=aj2/AF2p_j = \|\vec{a}_j\|^2 / \|\vec{A}\|_\F^2 and the fixed step size η=1/AF2\eta = 1 / \|\vec{A}\|_\F^2. This perspective, developed in Needell et al., 2016, explains several features of the method. For consistent systems, each fjf_j is minimized (with value zero) at the solution x\vec{x}^*, so the stochastic gradients vanish at x\vec{x}^* and a constant step size yields linear convergence. For inconsistent systems (i.e. genuine least-squares problems), the stochastic gradients do not vanish at the minimizer, and randomized Kaczmarz with a fixed step size only converges to within a ball around the least-squares solution whose radius depends on the residual; decaying step sizes or tail averaging are needed to converge to the minimizer itself.

References
  1. Strohmer, T., & Vershynin, R. (2009). A Randomized Kaczmarz Algorithm with Exponential Convergence. Journal of Fourier Analysis and Applications, 15(2), 262–278. 10.1007/s00041-008-9030-4
  2. Needell, D., Srebro, N., & Ward, R. (2016). Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm. Mathematical Programming, 155(1–2), 549–573. 10.1007/s10107-015-0864-7