An Asynchronous Parallel Randomized Kaczmarz Algorithm

Introduction

We consider the problem of finding a solution to a consistent linear system

Besides consistency of $Ax=b$, we assume that throughout that $A$ has no zero rows. In fact, we assume (to simplify the analysis) that that the rows of $A$ are normalized, that is,

although we define the algorithm as if normalization had not been applied.

We are interested in the casein which $A$ is extremely large and sparse. The randomized Kaczmarz (RK) is an algorithm for solving (1) that requires only $O(n)$ storage and has a linear (geometric) rate of convergence. In some situations, it is even more efficient than the conjugate gradient (CG) method (Strohmer and Vershynin 2009), which forms the basis of the most popular iterative algorithms for solving large linear systems. At iteration $j$, the RK algorithm randomly selects a row $i(j)\in\{1,2,\dotsc,m\}$ of the linear system (the probability of choosing row $i$ is ${\|a_{i}\|^{2}}/{\|A\|^{2}_{F}}$) and does an orthogonal projection of the current estimate vector onto the hyperplane $a_{i(j)}^{T}x=b_{i(j)}$:

This update formula can be derived also by applying the basic stochastic gradient algorithm to the objective $\frac{1}{2}\|Ax-b\|^{2}=\frac{1}{2}\sum_{i}(a_{i}^{T}x-b_{i})^{2}$ where $(a_{i}^{T}x-b_{i})a_{i}$ is the stochastic gradient corresponding to a random choice of index $i$ and $1/\|a_{i}\|^{2}$ is the steplength for that gradient estimate. The expected linear convergence rate of RK can be proved trivially as follows (Strohmer and Vershynin 2009; Needell 2010; Leventhal and Lewis 2010). Denoting by $x_{j}^{*}$ the projection of iterate $x_{j}$ onto the solution set of (1), we have from (2) that

Given the probability $\|a_{i}\|^{2}/\|A\|^{2}_{F}$ of choosing $i$, we have by taking expectations that

where $\lambda_{\min}$ is the smallest nonzero eigenvalue value of $A^{T}A$.

Recently, asynchronous parallel stochastic algorithms have received broad attention for solving large convex optimization problems. Niu et al. 2011 proposed a simple, effective asynchronous scheme to parallelize the stochastic gradient algorithm. In this approach, the unknown vector $x$ is stored in memory locations accessible to all cores of a multicore processor, and all cores are free to update $x$ in an asynchronous, uncoordinated fashion. It is assumed that there is a bound $\tau$ on the age of the updates, that is, no more that $\tau$ updates in total can be occur between the time at which any processor reads the current $x$ and the time at which it makes its update. Hogwild! (Niu et al. 2011) allows a lock-free implementation, since the update to a single element of $x$ is an atomic operation. Avron et al. 2014, Liu et al. 2013 and Sridhar et al. 2013 applied a similar asynchronous scheme to stochastic coordinate descent, and have proved attractive convergence properties.

We apply the same asynchronous parallel technique used in Hogwild! (Niu et al. 2011) to the standard RK algorithm. The unknown vector $x$ is stored in a shared location, and all cores simultaneously run a RK process, updating $x$ in an asynchronous fashion. Although our asynchronous parallel randomized Kaczmarz algorithm (AsyRK) can be viewed as an application of Hogwild! to the objective $\frac{1}{2}\|Ax-b\|^{2}$ (with a particular choice of step length), our analysis shows a linear convergence rate for AsyRK that outperforms the $1/t$ sublinear convergence rate for Hogwild!.

Our analysis also provides an indication of the maximum number of cores that can be involved in the computation while still yielding approximately linear speedup. This “bound” is expressed in terms of the number of equations $m$ and the maximal eigenvalue of $A^{T}A$.

An outline of the remainder of the paper is as follows. We review related work in Section 2. Section 3 illustrates details of the AsyRK algorithm. The convergence rate of AsyRK is described in Section 4, with proofs given in Appendix A. Some simple experiments illustrate linear speedup in Section 6. We discuss extensions to the inconsistent case in Section 7, and make some concluding observations in Section 8.

$\|X\|$ is the spectral norm of the matrix $X$, while $\|X\|_{F}$ is the Frobenius norm.

$P_{t}$ is the square $n\times n$ matrix of all zeros, except for a $1$ in the $(t,t)$ position.

Several quantities characterize the rows and columns of $A$: $\theta_{i}:=\|a_{i}\|_{0}$, $\mu:=\max_{i}\|a_{i}\|_{0}$, $\nu:=\max_{j}\|\bar{a}_{j}\|_{0}$.

$\alpha:=\max_{i,t}\|A\theta_{i}P_{t}a_{i}\|$. One can verify that $\alpha\leq\sqrt{\nu}\mu$ and $\alpha\leq\|A\|\mu$.

The support index set of $x$ is defined as $\mbox{\rm supp}(x)$.

$\lambda_{\min}$ is defined as the minimal nonzero eigenvalue value of $A^{T}A$, while $\lambda_{\max}$ is defined as the maximal eigenvalue value of $A^{T}A$.

We make a few observations about $\lambda_{\max}$. If $A$ is a matrix whose elements are i.i.d Gaussian random variables from $\mathcal{N}(0,1)$, then row-normalized, fundamental results in random matrices (Vershynin 2011) yield that $\lambda_{\max}$ is bounded by $\left({\sqrt{m}+\sqrt{n}\over\sqrt{n}}\right)^{2}\leq O(1+m/n)$ with high probability. As long as $m/n$ is bounded by a constant, $\lambda_{\max}$ is bounded by a constant as well. If $A$ is a sparse matrix, then

$A$ is row-normalized, that is, $\|a_{i}\|=1~\forall i\in\{1,2,\cdots,m\}$.

Note that $\|A\|_{F}^{2}=m$ when the rows of $A$ are normalized.

Related Work

The original Kaczmarz algorithm (Kaczmarz 1937) used a cyclic projection procedure to solve consistent linear systems $Ax=b$. Kaczmarz proved convergence to the unique solution when $A$ is a square nonsingular matrix. The cyclic ordering of the iterates made it difficult to obtain iteration-based convergence results, but Galantai 2005 proved a linear convergence rate in terms of cycles. Since the 1980s, the Kaczmarz algorithm has found an important application area in Algebraic Reconstruction Techniques (ART) for image reconstruction; see for example Herman 1980; Herman 2009. It is sometimes referred to in this literature as the “sequential row-action ART algorithm.”

Strohmer and Vershynin 2009 studied the behavior of RK in the case of a consistent system $Ax=b$ in which $A$ has full column rank (making the solution unique). They proved linear convergence rate for RK in expectation. Needell 2010 also assumed full column rank, but dropped the assumption of consistency, showing that the RK algorithm converges linearly to a ball of fixed radius centered at the solution, where the radius is proportional to the distance of $b$ from the image space of $A$. Eldar and Needell 2011 presented a modified version of the randomized Kaczmarz method which selects the optimal projection from a randomly chosen set at each iteration. This technique improves the convergence rate, but requires more computation per iteration. Liu and Wright 2013 proposed an accelerated RK algorithm that uses a Nesterov-type accelerated scheme, improving the linear convergence rate constant from $1-\lambda_{\min}/m$ (corresponding to (3), after normalization of rows) to $1-\sqrt{\lambda_{\min}}/m$.

Leventhal and Lewis 2010 extended the RK algorithm for consistent linear equalities $Ax=b$ to the more general setting of consistent linear inequalities and equalities: $A_{I}x\geq b_{I}$, $A_{E}x=b_{E}$. The basic idea is quite similar to the RK algorithm: iteratively update $x_{k+1}$ by projecting $x_{k}$ onto the hyperplane or half space for a randomly selected equality or inequality constraint. The linear convergence rate was proven to be $1-1/{(L^{2}\|A\|^{2}_{F})}$, where $L$ is the Hoffman constant (Hoffman 1952) for the full system.

Zouzias and Freris 2012 considered the case of possibly inconsistent (1). They proposed a randomized extended Kaczmarz algorithm by first projecting $b$ orthogonally onto the image space of $A$ to obtain $b_{\bot}$, then orthogonally projecting the initial point $x_{0}$ onto the hyperplane $Ax=b_{\bot}$. Essentially, the RK algorithm is applied twice. The convergence rate is proven to be $1-{\lambda_{\min}/\|A\|^{2}_{F}}$, which is the same as the RK algorithm for consistent linear systems. This method can be considered as a randomized variant of the extended Kaczmarz method proposed by Popa 1999.

Among synchronous parallel methods, Censor et al. 2001 proposed a parallel component averaging method to solve (1). This approach parallel-projects the current $x$ onto all (or multiple) hyperplanes, then applies an averaging scheme to the projections to obtain the next iterate. This method is essentially a gradient descent method for solving ${1\over 2}\|Ax-b\|^{2}$, so is able to handle inconsistent problems. This paper also notes (Censor et al. 2001, Section 5.1) that for sparse problems, parallelism can be obtained by simultaneously projecting the current iterate onto a set of mutually orthogonal hyperplanes, obtained by considering equations whose nonzero components appear in disjoint locations. When obtained forom image reconstruction problems, such sets of equations can be obtained by considering parallel rays that are sufficiently far apart so as to pass through disjoint sets of pixels. This type of parallelism has small granularity, and the amount of communication required between processors may make it unattractive in practice.

A related approach is Block-Cimmino (or Block-AMS) algorithm (Aharoni and Censor 1989), which can be considered as the block version of Kaczmarz algorithm. Other variants of Block-Cimmino algorithm are described in (Elfving and Nikazad 2009; Nikazad 2008).

Another synchronous parallel approach (for general convex optimization) due to Ferris and Mangasarian 1994 distributes variables among multiple processors and optimizes concurrently over each subset. A synchronization step searches the affine hull formed by the current iterate and the partial optima found by each processor.

In discussing asynchronous parallel methods, we make a distinction according to whether it is assumed that the reading of $x$ by each processor is “consistent” or not. The term “consistent” in this context means that the $x$ used by each processor to evaluate its update is an iterate that actually existed at some point in time, whose components were not changed repeatedly by other processors during reading (yielding a hybrid of two or more iterates). Bertsekas and Tsitsiklis 1989 introduced an asynchronous parallel implementation for general fixed point problems $x=q(x)$ over a separable convex closed feasible region. The optimization problem of minimizing $f$ over a closed convex set $\Omega$ can be formulated as a fixed-point problem by defining $q(x):=\mathcal{P}_{\Omega}[(I-\alpha\nabla f)(x)]$, where $\mathcal{P}_{\Omega}$ denotes Euclidean projection onto $\Omega$. The vector $x$ is stored in memory accessible to all cores, and the cores update the value of $x$ without locking or coordination. Inconsistent reading of $x$ is allowed. Linear convergence is established — using admirably straightforward analysis — provided that $\nabla^{2}f(x)$ satisfies a diagonal dominance condition, guaranteeing that the iteration $x=q(x)$ is a maximum norm contraction mapping for sufficient small $\alpha$. We note, however, that this condition is even stronger than strong convexity.

Elsner et al. 1990 proposed an asynchronous parallel RK algorithm, again for a situation in which all processors have access to $x$ stored in commonly accessible memory. Each processor iteratively runs the following procedures, where $x$ denotes a globally shared version of the variable vector and $x^{\prime}$ and $x^{\prime\prime}$ denote copies stored locally on each processor: (a) read the current global $x$ into the local $x^{\prime}$; (b) write the convex combination of the local variables $x^{\prime}$ and $x^{\prime\prime}$ into $x^{\prime\prime}$, and also into the shared memory as a new $x$; (c) project the local $x^{\prime\prime}$ onto the selected hyperplane to get a new $x^{\prime\prime}$. Note that the algorithm requires locking the shared memory in step (b) because it does not allow two processors to access shared memory at the same time. Our computational experiences with related algorithms (e.g., AsySCD (Liu et al. 2013) and Hogwild! (Niu et al. 2011)) indicate that memory locking of this type degrades computational performance seriously. Moreover, the convergence analysis establishes convergence, but does not prove a linear convergence rate.

Hogwild! (Niu et al. 2011) is a lock-free, asynchronous parallel version of the stochastic gradient method. All processors share the same memory storing $x$ and update it simultaneously. Unlike Bertsekas and Tsitsiklis 1989, inconsistent reads of $x$ are not permitted by the analysis. When the updates satisfy a certain sparsity property, the convergence of Hogwild! approximately matches the $1/t$ rate of serial stochastic gradient, as described and analyzed by (Nemirovski et al. 2009). Recent work by Avron et al. 2014 concerned an asynchronous linear solver for $Ax=b$ (for symmetric positive definite $A$) using the same asynchronous scheme as Hogwild!, proving a linear convergence rate.

Liu et al. 2013 followed the model of Hogwild! to propose an asynchronous parallel stochastic coordinate descent (AsySCD) algorithm and proved sublinear ($1/t$) convergence on general convex functions and a linear convergence rate on functions that satisfy an essential strong convexity property. Richtárik and Takáč 2012 proposed a parallel coordinate descent method for minimization of a composite convex objective with separable nonsmooth part. Their method is a synchronous parallel approach (in contrast to AsySCD, which is asynchronous), but it is implemented in an asynchronous fashion. Another distinction between the two approaches is found in the convexity assumptions, which are slightly weaker in Liu et al. 2013.

Algorithm

Each thread in our AsyRK algorithm performs the following simple steps: (1) Choose an index $i$ randomly from $\{1,2,\dotsc,m\}$; (2) read the components of $x$ that correspond to the nonzeros in $a_{i}$ from shared memory; (3) calculate $a_{i}^{T}x-b_{i}$; (4) select $t\in\text{supp}(a_{i})$; (5) update component $t$ of $x$ in the shared memory by a multiple of $(a_{i})_{t}(a_{i}^{T}x-b_{i})$. In principle, no memory locking takes place during either read or write, but we assume that the reads are “consistent,” according to the discussion above. (We note that inconsistent reading is expected to be rather rare in the case of sparse $A$, because only those elements of $x$ that correspond to nonzero locations in $a_{i}$ need to be read, and inconsistency possibly occurs only when this subset of elements is updated at least twice by other processors while it is being read.) The update to component $t$ of $x$ can be implemented as a unitary operation, requiring no memory locking.

Algorithm 1 gives a global, aggregated view of this multithreaded process. An iteration counter $j$ is incremented each time $x$ is updated by a thread. We use $k(j)$ to denote the iterate at which $x$ was read by the thread that updated $x_{j}$ to $x_{j+1}$. (We always have $k(j)\leq j$, and strict inequality holds when other threads have updated $x$ between the time it is read and the time the update is performed by this thread.) The index $i(j)\in\{1,2,\dotsc,m\}$ denotes the row that was selected by the thread that updated $x_{j}$ to $x_{j+1}$. The index $t(j)\in\text{supp}(a_{i(j)})$ denotes the component of $x$ that is chosen (randomly) to be updated at iteration $j$. We assume that the delay between reading and update for each thread is not too long, that is,

for some integer $\tau\geq 1$. $\tau$ can be assumed to be similar to the number of processors that are involved in the computation. Note that the step depends on $\theta_{i(j)}$ (the cardinality of the chosen row) and a parameter $\gamma$ which is critical to the analysis of the following sections.

Main Results

This section presents the convergence analysis for AsyRK. The key issue for AsyRK is to choose an appropriate steplength parameter $\gamma$. At an intuitive level, we would like $\gamma$ to be large enough to make significant progress in the approximate gradient direction. On the other hand, we want to keep it small enough that the approximate gradient information computed at the earlier iterate $k(j)$ is still relevant when the time comes to do the update at iteration $j$. That is, the difference between $x_{k(j)}$ and $x_{j}$ should not be too large. Along these lines, we require the ratios of expected residuals at any two successive iterations to be bounded above and below, as follows:

where $\rho$ is a user defined parameter, usually set to be slightly larger than $1$. The steplength $\gamma$ depends strongly on $\rho$.

We state a result about convergence of AsyRK in Algorithm 1.

Assume that Assumption 1 is satisfied. Let $\rho$ be any number greater than 1 and define the quantity $\psi$ as follows:

Suppose the steplength parameter $\gamma>0$ in Algorithm 1 satisfies the following three bounds:

This theorem indicates a linear rate of convergence, outperforming the sublinear “$1/j$” convergence rate for the asynchronous stochastic gradient method Hogwild!. The key reason for this improvement is that because $a_{i}^{T}x^{*}-b_{i}=0$ for all $i$, the stochastic gradient estimates all approach zero as $x$ approaches $x^{*}$, a property that does not hold for general stochastic gradient algorithms.

Note that the upper bound on steplength parameter $\gamma$ decreases as the bound $\tau$ on the age of the iterates increases. This dependency allows us to figure out how many threads can be executed in parallel without significantly degrading the convergence behavior.

This following corollary proposes an interesting particular choice for the parameters for which the convergence expressions become more comprehensible. The result requires a condition on the delay bound $\tau$ in terms of $m$ and $\lambda_{\max}$.

Suppose that Assumption 1 is satisfied and that

and set $\gamma=1/\psi$, where $\psi$ is defined as in (5), we have that

Over a span of $m$ iterations, (11) implies a decrease factor of approximately $1-{\lambda_{\min}}/{(\mu+1)}$. This rate estimate indicates that for a delay $\tau$ (and hence a number of processors) in the range implied by (9), the number of iterations required for convergence is not affected much by the delay, so we can expect an almost linear speedup from the multicore implementation in this regime.

We conclude this section with a high-probability estimate for convergence of $\{\|x_{j}-x_{j}^{*}\|^{2}\}_{j=1,2,\dots}$.

Suppose that the assumptions of Corollary 2 hold, and that $\rho$ and $\psi$ are defined as there. For $\epsilon>0$ and $\eta\in(0,1)$, if

The proofs of all results in this section appear in Section A.

Comparison

This section compares the theoretical performance of RK, AsySCD (applied to minimization of $\frac{1}{2}\|Ax-b\|^{2}$) and AsyRK. In Table 1, we show the complexities (per iteration) and convergence rates (with respect to number of iterations) of three algorithms in the first and second rows. The third row gives the maximal possibly number of processors to parallelize three algorithms respectively. The last row computes the convergence rate in term of the operation using the possibly maximal number of processors, which can be roughly understood as the running time comparison. In reporting statistics for AsySCD, we consider two alternative implementations: (1) randomly choose a coordinate $i$ and compute it by $(a_{i}A)x$, which needs $O(\delta^{2}mn)$ operations per iteration; and (2) compute $A^{T}A:=Q$ offline and randomly choose an coordinate $i$ to compute it by $Q_{i.}x$, which needs $O(n)$ operations per iteration. We report the complexity per iteration of AsySCD as the minimum of these two estimates.

We perform a comparison of convergence behavior of these three algorithms on a Gaussian random matrix $A$ with i.i.d. elements generated from $\mathcal{N}(0,1/n)$. All rows of $A$ have norm approximately $1$, and $\lambda_{\max}$ is approximately $1+m/n$. For these values, the convergence rates per iteration of AsySCD and RK are similar. By comparison, the convergence rate of AsyRK seems worse than RK and AsyRK by a factor $(\mu+1)$. This is because AsyRK only updates a single coordinate rather than all coordinates corresponding to the nonzero elements in the stochastic gradient. If we modify Algorithm 1 to updated all components in $\mbox{\rm supp}(a_{i(j)})$ (rather than just the component $t(j)$), the convergence rate for AsyRK becomes quite similar to the other two methods, without an appreciable increase in cost.

Next we compare the parallel implementations. From the last row of Table 1, we see that AsyRK improves the rate of RK if $\delta mn\gg\delta^{2}n^{2}\lambda_{\max}$, or equivalently $m\gg\delta n\lambda_{\max}$. Assuming the Gaussian ensemble for $A$ and that $m$ and $n$ are comparable (so that $\lambda_{\max}=O(1)$), there is a potential factor of improvement in runtime of $O(1/\delta)$ for AsyRK over RK. To compare AsyRK and AsySCD, we note that $\lambda_{\max}=O(L_{\text{res}})$ under the same scenario for $m$, $n$, and $A$. Comparing the rates (running time) in the last row of Table 1, we find that when the mild condition $\delta<O(n^{-1/4})$ holds, AsyRK converges much faster than AsySCD. Overall, AsyRK has a clear advantage in complexity when applied to sparse problems.

Experiments

We illustrate the behavior of AsyRK on sparse synthetic data. Our chief interest is the efficiency of multicore implementations (one thread per core), compared to a single-thread implementation.

Our experiments run on $1$ to $10$ threads on an Intel Xeon machine, with all threads sharing a single memory socket. Our implementations deviate modestly from the version of AsyRK analyzed here. First, $A$ is partitioned into slices (row submatrices) of equal size, and each thread is assigned one slice. Each thread then selects the rows in its slice to update in order, with the order being reshuffled after each scan. This scheme essentially changes from sampling with replacement (as analyzed) to sampling without replacement, which has empirically better performance. (The same advantage is noted in implementations of Hogwild!.) The second deviation from the analyzed version is that all coordinates corresponding to nonzeros in the selected row $a_{i(j)}$ are updated, not just the $t(j)$ component. This scheme makes a single thread behave like $|a_{i(j)}|$ threads, thus implicitly increasing the number of cores involved in the computation. Note that this variant represents the obvious extension of randomized RK. In fact, when implemented on a single thread, it is precisely the usual randomized RK scheme.

For the plots in Figures 1 and 2, we choose $m=80000$ and $n=100000$, with $\delta=0.001$, and set the steplength $\gamma$ as $1$ in Figure 1 and $\delta=0.003$ in Figure 2. The left-hand graph in each figure indicates the number of threads / cores and plots residual (defined as $\|Ax-b\|^{2}$) vs epoch count, where one epoch is equivalent to $n$ iterations. Note that the curves tend to merge, indicating that the workload required for AsyRK is almost independent of the number of cores. This observation validates our result in Corollary 2, which indicates that provided it is below a certain threshold, the value of $\tau$ does not affect convergence rate. The right-hand graph in each figure shows speedup over different numbers of cores. Near-linear speedup is observed for $\delta=0.001$ (Figure 1), while for $\delta=0.003$ there is a dropoff for larger numbers of cores (Figure 2). This can perhaps be explained by the difference between our implementation from the analyzed version, in that the nonzeros in the full row $a_{i(j)}$ are updated rather than just a single element. The effect of this policy can be incorporated into the analysis roughly by increasing the value of the maximum delay parameter $\tau$. In this case, a matrix that is three times more dense could be modeled by a value of $\tau$ that is three times larger. The effect may be to raise $\tau$ above the threshold for which linear speedup can be expected, thus explaining the (graceful) degradation in speedup for larger numbers of cores in Figure 2.

Next, we compare AsyRK to AsySCD (Liu et al. 2013) on sparse synthetic data sets, on $10$ cores (single socket) of the Intel Xeon. Various values of $m$, $n$, and $\delta$ are chosen for comparison in Table 2. A similar number of epochs is required by both algorithms, reflecting the similarity of their theoretical convergence rates; see Section 5. However, AsyRK is one order of magnitude faster than AsySCD to achieve the same accuracy. The main reason is that, as we showed in Table 1, the per-iteration complexity of AsySCD is much higher than AsyRK, for these values of the parameters.

Extension to Inconsistent Systems

Although this paper assumes that the linear system is consistent, we can extend the algorithm described above to find the least-squares solution of inconsistent linear systems.

The minimizer of the least-squares objective $\|Ax-b\|^{2}$ is equivalent to the linear system $A^{T}Ax=A^{T}b$, which can be stated as the following square, consistent system of linear equations:

for any positive values of $\zeta$ and $\phi$. Similar reformulations have appeared previously in the literature; see Eggermont 1981, for example. Here we are mainly interested in the optimal values for $\zeta$ and $\phi$. We can choose $\zeta$ and $\phi$ to maximize the critical quantity in the analysis of Algorithm 1, which is the ratio of the minimum nonzero eigenvalue of $A^{T}A$ to its squared Frobenius norm. (In Theorem 1, this ratio appears as $\lambda_{\min}/m$, because of the normalization of the rows of $A$.) To show how this quantity depends on $\zeta$ and $\phi$, we note first that the coefficient matrix in (14) is

For fixed $\phi$, the first term is monotonically increasing with respect to $\zeta$ while the second term is monotonically decreasing with respect to $\zeta$. We can thus express the optimal $\zeta$ value as $\zeta^{*}=\sigma_{r}\sqrt{\phi/2}$, which is the value for which these two terms are equal. By substituting this value into (15), we obtain

It is clear from the last expression that the optimal value for $\phi$ is $\phi^{*}=1$, giving the following maximal value for (16):

We conclude that by normalizing the rows of $A$, estimating its minimum singular value $\sigma_{r}$, and setting $\phi=1$ and $\zeta=\sigma_{r}/\sqrt{2}$, we obtain an optimally conditioned system (14), to which the approach of this section can be applied.

Conclusion

We have proposed a simple asynchronous parallel randomized Kaczmarz algorithm, and proved linear convergence. Our analysis also indicates the proposed method can be expected to yield near-linear speedup if the number of processors is bounded by a multiple of the number of equations in the system. Computational results, including comparison with an asynchronous stochastic coordinate descent method, confirm the effectiveness of the approach.

Acknowledgements

The authors acknowledge support of National Science Foundation Grants DMS-0914524 and DMS-1216318, ONR Award N00014-13-1-0129, AFOSR Award FA9550-13-1-0138, and a Wisconsin Alumni Research Foundation 2011-12 Fall Competition Award. The authors would like to sincerely thank Yijun Huang for her implementation of algorithms AsyRK and AsySCD used in this paper.

Appendix A Proofs

This section provides proofs for our main results in Section 4.

We start with the following useful results, noting that the random variable $t(j)$ is distributed uniformly over the set $\mbox{\rm supp}(a_{i(j)})$:

We prove each of the two inequalities in (7) by induction. We start from the right-hand inequality. First we consider the expansion of $\|Ax_{j+1}-b\|^{2}$ for any values of $j$:

Next we consider the expectation of three terms $T_{1}$, $T_{2}$, and $T_{3}$ in (18). For $T_{1}$, we have

where the third line is from the observation that $t(d)$ and $t(j)$ are conditionally independent given $i(d)$ and $i(j)$; the fifth line uses the result that $i(d)$ only affects $x_{d+1}$ and subsequent iterates and $k(j)$ is less than $d+1$ (so $x_{k(j)}$ and $i(d)$ are independent to each other). Combining (19), (20), (21), and (18), we obtain

We can use this bound to show that the right-hand inequality in (7) holds for $j=0$. By setting $j=0$ in (22) and noting that $k(0)=0$ and that the last summation is vacuous, we obtain

For the inductive step, we use (22) again, assuming that the right-hand inequality in (7) holds up to stage $j$, and thus that

provided that $0\leq j-k(j)\leq\tau$ and $0\leq j-k(d)\leq 2\tau$, as assumed. By substituting into the right-hand side of (22) again, we obtain

where the last inequality uses the third bound on $\gamma$ from (6). We conclude that the right-hand side inequality in (7) holds for all $j$.

We now work on the left-hand inequality in (7). For all $j$, we have the following:

We can use this bound to show that the left-hand inequality in (7) holds for $j=0$. By setting $j=0$ in (22) and noting that $k(0)=0$, we obtain

provided that $0\leq j-k(j)\leq\tau$, as assumed (4). By substituting into the left-hand side of (24) again, we obtain

We conclude that the left-hand side inequality in (7) holds for all $j$.

At this point, we have shown that both inequalities in (7) are satisfied for all $j$.

We next prove (8). Consider the expansion of $\|x_{j+1}-x_{j+1}^{*}\|^{2}$:

Next, we estimate the expectations of $T_{4}$, $T_{5}$, and $T_{6}$. For $T_{4}$, we have

By following a derivation similar to (21) for $T_{6}$, we obtain

Since for $d=k(j),k(j)+1,\dotsc,j-1$, we have

By substituting (27), (28), and (30) into (26), we obtain

Note first that for $\rho$ defined by (10), and using (9), we have

Thus from the definition of $\psi$ (5), and using (9) again, we have

We show now that the steplength parameter choice $\gamma={1/\psi}$ satisfies all the bounds in (6), by showing that the second and third bounds are implied by the first. For the second bound, we have

where the second inequality follows from (10) and the final inequality follows from the definition of $\psi$ in (5) and the fact that $\psi>\mu\geq 1$.

For the third bound in (6), we have (by taking squares) that

We can thus set $\gamma=1/\psi$, and by substituting this choice into (8) and using (31), we obtain (11). ∎

where the second inequality applies (11), the third inequality uses the definition of $j$ (12), and the second last inequality uses the inequality $(1-c)^{1/c}\leq e^{-1}~\forall c\in(0,1)$, which completes the proof. ∎

References