Block Kaczmarz Method with Inequalities

Jonathan Briskman, Deanna Needell

Introduction

The Kaczmarz method is an iterative algorithm for solving linear systems of equations. It is usually applied to large-scale overdetermined systems because of its simplicity and speed (but also converges in the underdetermined case to the least-norm solution under appropriate initial conditions). Each iteration projects onto the solution space corresponding to one row in the system, in a sequential fashion. Strohmer and Vershynin prove that when the rows are selected from a certain random distribution rather than sequentially, that the randomized method converges to the solution at a linear rate . The method has been applied to fields including image reconstruction, digital signal processing, and computer tomography . Leventhal and Lewis modify the randomized Kaczmarz method to apply to systems of linear equalities and inequalities , thereby extending results on the standard method in this setting (see e.g. and references therein). Unlike the traditional randomized algorithm which enforces a single constraint at each iteration, the block Kaczmarz approach recently analyzed by Needell and Tropp enforces multiple constraints simultaneously and thus offers computational advantages. Here we demonstrate convergence for a system of linear equalities and inequalities by combining a randomized block Kaczmarz method for the equalities with a randomized Kaczmarz algorithm for the inequalities. These results indicate that the block Kaczmarz method can be used for a system of equalities and inequalities, and in some cases may quicken convergence. We also consider the case of utilizing blocking in both the equalities and inequalities, although this can be detrimental unless the geometry of the system meets certain conditions.

where A\bm{A} is a real (or complex) n×dn\times d matrix, typically with n≫dn\gg d.

We define the eigenvalues λmin⁡(A),…,λmax⁡(A)\lambda_{\min}(\bm{A}),\ldots,\lambda_{\max}(\bm{A}) of a matrix analogously. For convenience we will assume that each row ai\bm{a_{i}} of A\bm{A} has unit norm, ∥ai∥2=1\lVert{\bm{a_{i}}}\rVert_{2}=1, and we call such matrices standardized.

Now we consider a system of linear equalities and inequalities and denote by SS its non-empty set of feasible solutions. We thus consider the matrix A\bm{A} whose rows can be arranged such that

and we will write I=I_{=} and I≤I_{\leq} to denote the row indices of A=\bm{A}_{=} and A≤\bm{A}_{\leq}, respectively. Therefore, we ask that

We will assume that the set of rows {1,2,...,n}\{1,2,...,n\} is partitioned such that the first nen_{e} rows correspond to equalities, and the remaining ni=n−nen_{i}=n-n_{e} rows to inequalities. Thus A=\bm{A}_{=} is an ne×dn_{e}\times d matrix and A≤\bm{A}_{\leq} is ni×dn_{i}\times d.

The error bound for this system of linear inequalities uses the function e:Rn→Rne:\textbf{R}^{n}\rightarrow\textbf{R}^{n} defined as in by

2. Details of Kaczmarz

The simple Kaczmarz method is an iterative algorithm that approximates a least-squares minimizer x⋆\bm{x}_{\star} to the problem in (1.1). It takes an arbitrary initial approximation x0\bm{x}_{0}, and at each iteration jj the current iterate is projected orthogonally onto the solution hyperplane {⟨ai,x⟩=bi}\{\left\langle\bm{a}_{i},\bm{x}\right\rangle=\bm{b}_{i}\}, using the update rule

where i=ji=j mod n+1n+1 . With an unfortunate ordering of the rows, this method as-is can produce very slow convergence. However, it has been well known that using randomized selection often eliminates this effect . The randomized Kaczmarz method put forth by Strohmer and Vershynin uses a random selection method for the selection of row ii such that each row ii is selected with probability proportional to ∥ai∥22\lVert{\bm{a}_{i}}\rVert_{2}^{2}. This randomization provides an algorithm that is both simple to analyze and enforce in many cases. In this paper we assume each row has unit norm, so each row is selected uniformly at random from {\{1,…,n}\} in the simple randomized Kaczmarz approach. This assumption is both for notational convenience, and because the use of matrix pavings discussed below only hold for standardized matrices. In practice, one can employ pre-conditioning on non-standardized systems, or extend the construction of matrix pavings to non-standardized systems . Strohmer and Vershynin prove a linear rate of convergence for consistent systems that depends on the scaled condition number of A\bm{A}, and not on the number of equations nn ,

where x⋆\bm{x}_{\star} is the solution to the consistent system (1.1) and K=∥A∥F2/σmin⁡2(A)K=\|\bm{A}\|_{F}^{2}/\sigma^{2}_{\min}(\bm{A}) denotes the scaled condition number. Needell extended this work to the inconsistent case and proves linear convergence to the least-squares solution within some fixed radius ,

where e=Ax⋆−b\bm{e}=\bm{Ax_{\star}}-\bm{b} denotes the residual vector. Because the Kaczmarz method projects directly onto each solution hyperplane, such a convergence radius is unavoidable without adding a relaxation parameter.

The randomized Kaczmarz method can be adapted to the case of a linear system of equalities and inequalities described in (1.3). Leventhal and Lewis apply the Kaczmarz method to a consistent system of linear equalities and inequalities (here consistent simply means the feasible set SS is non-empty). At each iteration jj, the previous iterate only projects onto the solution hyperplane if the inequality is not already satisfied. If the inequality is satisfied for row ii selected at iteration jj (aiTx≤bi)(\bm{a}_{i}^{T}\bm{x}\leq\bm{b}_{i}), the approximation xj\bm{x}_{j} is set as xj−1\bm{x}_{j-1} . The update rule for this algorithm is thus

This algorithm converges linearly in expectation , with

In order to bound the right hand side of this expression, the authors rely on a lemma due to Hoffman . This result states that for any system (1.3) with non-empty solution set SS, there exists a constant LL independent of b\bm{b} such that for all x\bm{x},

When A==A\bm{A}_{=}=\bm{A} is full column rank, the Hoffman constant is the inverse of the smallest singular value, L=σmin⁡−1(A)L=\sigma_{\min}^{-1}(\bm{A}).

which coincides with (1.5) for consistent systems of equalities.

3. Block Kaczmarz

A block variant of the randomized Kaczmarz method due to Elfving has been recently analyzed by Needell and Tropp and can improve the convergence rate in certain cases. The block Kaczmarz method first partitions the rows {1,...,n}\{1,...,n\} into mm blocks, denoted τ1,…τm\tau_{1},\ldots\tau_{m}. Instead of selecting one row per iteration as done with the simple Kaczmarz method, the block Kaczmarz algorithm chooses a block uniformly at random at each iteration. Thus the block Kaczmarz method enforces multiple constraints simultaneously. At each iteration, the previous iterate xj−1\bm{x}_{j-1} is projected onto the solution space to Aτx=bτ\bm{A}_{\tau}\bm{x}=\bm{b}_{\tau}, which enforces the set of equations in block τ\tau . Aτ\bm{A}_{\tau} and bτ\bm{b}_{\tau} are written as the row submatrix of A\bm{A} and the subvector of b\bm{b} indexed by τ\tau respectively, yielding an iterative rule of

The pseudoinverse used in (1.9) returns the solution to the underdetermined least squares problem for a wide or square row submatrix Aτ\bm{A}_{\tau}.

Depending on the characteristics of the submatrix Aτ\bm{A}_{\tau}, the block method can provide better convergence than the simple method. If we assume that the submatrices Aτ\bm{A}_{\tau} are well conditioned, the additional cost of computing their pseudo-inverse can be overcome by the gain in utilizing block multiplications (see our experiments in Section 4). In fact, if the blocks admit a fast multiply (for example if the matrix is built of DFT or circulant blocks), then the computational cost of the block iteration (1.9) is similar to the cost of the simple update rule in (1.4). Since the convergence depends heavily on the conditioning of each submatrix, one seeks partitions of the rows into blocks for which each block is well-conditioned. The notion of a row-paving allows one to do precisely that.

We define an (mm, β\beta) row paving The standard definition of a row paving also includes a constant α\alpha which serves as a lower bound to the smallest singular value. We ignore that parameter here since it will not be utilized. of matrix A\bm{A} as a partition T={τ1,...τm}T=\{\tau_{1},...\tau_{m}\} of the row indices such that

The size of the paving, or number of blocks, is mm. The value of β\beta is the upper paving bound, which controls the spectral norms of the submatrices. Needell and Tropp show that these parameters determine the performance of the algorithm, with convergence for a consistent system admitting an (m,β)(m,\beta) paving given by

Therefore the convergence rate depends on the size mm and upper bound β\beta; the algorithm’s performance improves with low values of mm and β\beta, and large σmin⁡2(A)\sigma_{\min}^{2}(\bm{A}). The authors also prove convergence for inconsistent systems, with the same convergence rate and convergence radius which depends also on the minimum of all λmin⁡(AτAτ∗)\lambda_{\min}(\bm{A}_{\tau}\bm{A}_{\tau}^{*}), see for details.

Surprisingly, every standardized matrix admits a good row paving. The following result is due to which builds off the foundational work of .

For any δ∈(0,1)\delta\in(0,1) and standardized n×dn\times d matrix A\bm{A}, there is a row paving satisfying

Although this is an existential result, there are constructive methods to obtain such pavings, and for certain classes of matrices, they can even be obtained by a random partitioning of the rows .

With such a paving in tow, the convergence of (1.10) becomes

Although often comparable to the convergence rate for the simple method (1.5), numerical results confirm that the block method offers significant reduction in computation time due to the speed of matrix–vector multiplication (see e.g. ).

4. Contribution

This paper analyzes the system with matrix described in (1.2) using an algorithm with the block Kaczmarz approach for the equalities given by A=\bm{A}_{=} and the simple method for the inequalities given by A≤\bm{A}_{\leq}. A paving is created for A=\bm{A}_{=}, with the inequalities excluded. At each iteration, we select from A=\bm{A}_{=} with a fixed probability pp and from A≤\bm{A}_{\leq} with probability 1−p1-p. In the former case, we select a block τ\tau from paving TT uniformly at random, and in the latter case we select a row ii of A≤\bm{A}_{\leq} uniformly at random. In the case of a block of equalities being selected, the algorithm proceeds by updating xj\bm{x}_{j} using (1.9). When an inequality row is selected, xj\bm{x}_{j} is updated using the rule (1.6). We prove that this method yields linear convergence to the solution set SS. We also include a discussion about paving both A=\bm{A}_{=} and A≤\bm{A}_{\leq}, which identifies a geometric property of the system which allows for improved convergence by utilizing two pavings. We show that when this property is not satisfied, utilizing both pavings can be detrimental to convergence.

5. Organization

Section 2 lays out our main result, Theorem 2.1, and provides a proof. We discuss blocking the full matrix in Section 3 and Section 4 explains numerical experiments and results. We conclude with discussion and related work in Section 5.

Analysis of the Block Kaczmarz Algorithm for a System of Inequalities

In this section we analyze the convergence of the described method, which is detailed in Algorithm 2.1.

Notice that the probability of selecting a block of A=\bm{A}_{=} is βmni+βm\frac{\beta m}{n_{i}+\beta m}. This quantity corresponds to the relative size of A=A_{=} in the system, where the size is measured in terms of the paving quantities βm\beta m. This value may be difficult to compute precisely, and the simpler threshold of ne/nn_{e}/n appears to also work well in practice. We provide no evidence that our selection of this threshold is most efficient, nor any more efficient than using one proportional to the number of equality rows nen_{e}. We find that this algorithm yields linear convergence in expectation with a rate that only depends on the number of inequalities nin_{i}, paving size mm, and upper bound β\beta.

Our main result is described in Theorem 2.1.

Remarks. 1. Note that when there are no block projections, no inequalities, or neither, Theorem 2.1 recovers the results of the standard randomized Kaczmarz for inequalities , the standard randomized block Kaczmarz method or the standard randomized Kaczmarz method , respectively. We thus view this result as a completely generalized convergence bound.

2. If we let ρs\rho_{s} and ρb\rho_{b} be the convergence rates of the simple and block methods for mixed systems, respectively, then by (1.8) and Theorem 2.1,

It is evident that our expected convergence rate will be faster per iteration than the simple method when ni+βm<nn_{i}+\beta m<n. Since β\beta can be chosen close to 11 and m<nem<n_{e} is then number of rows in A=\bm{A}_{=}, this holds quite easily.

3. Since a single iteration using a block Aτ\bm{A}_{\tau} in general may cost more than an iteration utilizing a single row, it is more fair to compare per epoch, rather than per iteration. An epoch is typically the minimum number of iterations needed to visit each row of the matrix. When there are inequalities present that are already satisfied in a given iteration, that iteration may make no contribution and cost very little computationally. Thus the notion of epoch may be slightly skewed here, but if we ignore this subtlety the simple method will have approximately nn iterations per epoch, compared to ni+mn_{i}+m iterations per epoch with the block method. The approximate per epoch convergence rates can thus be compared as

This result is similar to that found by Needell and Tropp , with the block convergence rate at best equal to that of the simple convergence rate when β=1\beta=1. However, as already noted, the block method is quite advantageous computationally.

Combining the paving result of Prop. 1.2 with Theorem 2.1 yields the following corollary.

Instate the assumptions and notation of Theorem 2.1 and let A=\bm{A}_{=} be equipped with an (m,β)(m,\beta) row-paving as in Proposition 1.2. Then the iterates of Algorithm 2.1 satisfy

where γ=(1−1L2(ni+C∥A=∥2log⁡(1+n)))\gamma=\left(1-\frac{1}{L^{2}(n_{i}+C\|\bm{A}_{=}\|^{2}\log(1+n))}\right) and CC is some absolute constant.

Fix an iteration jj of Algorithm 2.1. We proceed as in and . First, we suppose that q≤βmni+βmq\leq\frac{\beta m}{n_{i}+\beta m}, so that a block τ\tau of equalities is selected this iteration. Then writing PSP_{S} as the orthogonal projection onto SS, we have bτ=AτPSxj−1\bm{b}_{\tau}=\bm{A}_{\tau}P_{S}\bm{x}_{j-1} since PSxj−1∈SP_{S}\bm{x}_{j-1}\in S. We then have

Taking expectation (over the choice of the block τ\tau, conditioned on previous choices), and using the fact that Aτ†Aτ\bm{A}_{\tau}^{\dagger}\bm{A}_{\tau} is an orthogonal projector, along with the properties of the paving yields

Since d(xj−1,S)=∥xj−1−PSxj−1∥2d(\bm{x}_{j-1},S)=\lVert{\bm{x}_{j-1}-P_{S}\bm{x}_{j-1}}\rVert_{2} and d(xj,S)≤∥xj−PSxj−1∥2d(\bm{x}_{j},S)\leq\lVert{\bm{x}_{j}-P_{S}\bm{x}_{j-1}}\rVert_{2}, this means that

Next suppose that instead i∈I≤i\in I_{\leq} is selected. Then since each row ai\bm{a_{i}} has unit norm,

where the last line follows from the fact that ⟨ai,PSxj−1⟩≤bi\langle\bm{a_{i}},P_{S}\bm{x}_{j-1}\rangle\leq b_{i} and e(Axj−1−b)i≥0e(\bm{A}\bm{x}_{j-1}-\bm{b})_{i}\geq 0. Now taking expectation again we have

Combining these results and letting E=E_{=} and E≤E_{\leq} denote the events that a block from TT and a row from I≤I_{\leq} is selected, respectively, we have

Since p=βm ni+βmp=\frac{\beta m}{\ n_{i}+\beta m}, we have 1−pni=1ni+βm\frac{1-p}{n_{i}}=\frac{1}{n_{i}+\beta m} and we can simplify

where we have utilized the Hoffman bound (1.7) in the second inequality.

Utilizing independence of the random selections and recursing on this relation yields the desired result. ∎

A Discussion about Blocking Inequalities

It is natural to ask whether one can benefit by blocking both the equalities as above and also the inequalities, as described by Algorithm 3.2. Indeed, Section 4 will show dramatic improvements in computational time when the rows of A=\bm{A}_{=} are paved and block projections as in Algorithm 2.1 are used. So can one benefit even more by paving also the rows of A≤\bm{A}_{\leq}? The answer to this question heavily depends on the structure of the matrix A\bm{A}.

If we only consider A=\bm{A}_{=}, a block projection as in (1.9) enforces all the equations indexed by τ\tau to be satisfied. This is of course desirable when the rows indexed by τ\tau correspond to equalities. Also, if a single inequality corresponding to row ii in A≤\bm{A}_{\leq} is not satisfied and we perform a single projection as in (1.4), we are again enforcing that inequality to hold with equality. However, this improves the estimation since in this case we know the solution set SS lies on the opposite side of the hyperplane {x:⟨x,ai⟩=bi}\{\bm{x}:\langle\bm{x},\bm{a_{i}}\rangle=b_{i}\} as the current estimation (see Figure 1 (a)).

On the other hand, if we employ a block projection as in (1.9) to a set of inequalities indexed by τ\tau which are not satisfied by the current estimation xj−1\bm{x}_{j-1} then we enforce all of them to hold with equality simultaneously. Depending on the geometry of the involved rows, this may result in an improved estimation or actually one much farther from the solution set. Of course, one might alternatively want to solve the convex program to project onto the intersection of the corresponding half-spaces, but we would like to maintain the efficiency and simplicity of the block Kaczmarz method.

As an illustrative example, Figure 1 (b) and (c) demonstrate two possible scenarios in two dimensions. Here, the solution space is a single point marked SS, and we draw two hyperplanes Hi1H_{i_{1}} and Hi2H_{i_{2}} where Hi={x:⟨ai,x⟩=bi}H_{i}=\{\bm{x}:\langle\bm{a_{i}},\bm{x}\rangle=b_{i}\} . The yellow shaded regions denote areas where both inequalities hold true: {x:⟨ai1,x⟩≤bi1 and ⟨ai2,x⟩≤bi2}\{\bm{x}:\langle\bm{a}_{i_{1}},\bm{x}\rangle\leq b_{i_{1}}\text{ and }\langle\bm{a}_{i_{2}},\bm{x}\rangle\leq b_{i_{2}}\}. Notice that in (b), when the angle between xj−1−xj\bm{x}_{j-1}-\bm{x}_{j} and s−xj\bm{s}-\bm{x}_{j} is obtuse, the orthogonal projection of estimation xj−1\bm{x}_{j-1} onto their intersection is guaranteed to be closer to the solution set S. On the other hand, when that angle is acute we see exactly the opposite, as in (c). We can quantify this notion by the following definition.

In other words, the angle between w−z\bm{w}-\bm{z} and s\bm{s} (and thus s−z\bm{s}-\bm{z}) is obtuse.

We will see that performing block projections on the inequalities in the system only makes sense when one can obtain an obtuse row paving. We will use w=xj−1\bm{w}=\bm{x}_{j-1}, z=xj\bm{z}=\bm{x}_{j}, and s∈S\bm{s}\in S. Notice that if i1,i2∈τ∈Ti_{1},i_{2}\in\tau\in T, then the partition used in the system depicted in Figure 1 (c) does not constitute an obtuse row paving.

We conduct two simple experiments to demonstrate the different behavior of the algorithm. In all cases the matrix A\bm{A} is a 300×100300\times 100 matrix with standard normal entries, 100100 rows correspond to inequalities, and b\bm{b} is generated so that the solution set SS is non-empty. We measure the residual error which we define as ∥e(Axj−b)∥2\lVert{e(\bm{A}\bm{x}_{j}-\bm{b})}\rVert_{2}. Figure 2 (a) shows the behavior of the block method with this matrix and a row paving obtained via a random row partition of 3030 blocks (1010 rows per block). This generation will create a matrix with paving that with very high probability is not an obtuse row paving. As Figure 2 demonstrates, the block method does not converge to a solution in this case. However, as Figure 2 (c) shows, the simple Kaczmarz method succeeds in identifying a point in the solution space. Next, we create a matrix in the exact same way, and create the same random row paving. Then, however, we iterate through every block in the paving corresponding to inequalities and if two rows ii and kk in a block satisfy ⟨ai,ak⟩>0\langle\bm{a_{i}},\bm{a}_{k}\rangle>0, we replace row ai\bm{a_{i}} with −ai-\bm{a_{i}} and entry bib_{i} with −bi-b_{i}. This guarantees every block in the paving yields a geometry like that shown in Figure 1 (b), and gives an obtuse row paving. Note that of course this changes the solution space as well so one cannot employ this strategy in general. We then add positive values to the entries in b\bm{b} corresponding to inequalities to ensure the solution set SS is non-empty. With this new system and paving, we again run the block method and see that the method now converges to a point in the solution set, as seen in Figure 2 (b).

With this definition we obtain the following result, whose proof can be found in the appendix.

Let A\bm{A} satisfy the assumptions of Theorem 2.1 and in addition have an obtuse (m′,β′)(m^{\prime},\beta^{\prime}) row paving of A≤\bm{A}_{\leq}. Let x1,…x_{1},\ldots denote the iterates of Algorithm 3.2. Then using the notation of Theorem 2.1,

Note that row pavings of standardized matrices can be obtained readily, often by random partitions , whereas obtuse row pavings may be much more challenging to obtain in general. Of course, by default the trivial paving which assigns each set τ\tau to a single row always admits an obtuse row paving. We focus on Algorithm 2.1 which paves only A=\bm{A}_{=}, and leave further analysis of Algorithm 3.2 and constructions of obtuse row pavings for future work.

Experiments

We use Matlab to run some experiments using random matrices to test the convergence of the block Kaczmarz method applied to a system of equalities and inequalities. In each experiment, we create a random 500500 by 100100 matrix A\bm{A} where each element is an independent standard normal random variable. Each entry is then divided by the norm of its row so that the matrix is standardized. The first 400 rows of matrix A\bm{A} compose A=\bm{A}_{=}, and the remaining 100 rows are set as inequalities of A≤\bm{A}_{\leq} in the method described by (1.3). The experiments are run using the following procedure. For each of 100 trials,

Create matrix A\bm{A} in the manner described above.

Create x⋆\bm{x}_{\star} where each entry is selected independently from a standard normal distribution. Set b=Ax⋆\bm{b}=\bm{A}\bm{x}_{\star}.

Pave submatrix A=\bm{A}_{=} into 16 blocks with 25 equalities per block by a random partitioning of the rows.

Set initial approximations x0block=x0simp=A∗b\bm{x}_{0}^{\text{block}}=\bm{x}_{0}^{\text{simp}}=\bm{A}^{*}\bm{b}.

If q≤ne nq\leq\frac{n_{e}}{\ n}, choose block {1,...,m}\{1,...,m\} uniformly at random and update iterate xjblock\bm{x}_{j}^{\text{block}} using (1.9). (Note that the threshold ne n\frac{n_{e}}{\ n} is different than that given in the main algorithm and theorem, but it is easier to calculate and seems to work fine in practice.)

Else, choose a row uniformly at random from {401,...,500}\{401,...,500\} and update iterate xjblock\bm{x}_{j}^{\text{block}} using (1.6).

Update iterate xjsimp\bm{x}_{j}^{\text{simp}} using (1.6).

For both the simple and block algorithms, the median, minimum, and maximum values of the residual ∥e(Axj−b)∥22\lVert{e(\bm{A}\bm{x}_{j}-\bm{b})}\rVert_{2}^{2} of the 100 trials are recorded for each iteration jj.

Figure 3 compares the performance of the block Kaczmarz method used in this paper and the standard Kaczmarz method described by Leventhal and Lewis . The plot in Figure 3 (a) compares convergence per iteration. As the block Kaczmarz method enforces multiple equalities per iteration, it is unsurprising that it performs better in this experiment. Figure 3 (b) displays the convergence of the two methods per epoch. The block Kaczmarz algorithm has an epoch of m+nim+n_{i} iterations, and the standard Kaczmarz method has an epoch of size nn. Here, to be fair we only count an iteration towards an epoch if the estimated solution xj≠xj−1\bm{x}_{j}\neq\bm{x}_{j-1}. Thus in the case where a chosen inequality is already satisfied for iteration jj, this iteration does not count towards an epoch since no computation is being performed. We noticed, however, that whether or not we modified the count in this way, the behavior still produces results very similar to Figure 3. Once again the experiments yielded faster convergence with the block Kaczmarz approach. It is interesting to compare the results of Figure 3 (b) and those of Figure 2 (b) and (c). The per-epoch convergence of the methods and whether the block or standard appears faster varies slightly and depends on both the number of rows and columns. In general, the per-epoch convergence rates are reasonably comparable, as the analysis suggests. However, Figure 3 (c) compares the rate of convergence of the two algorithms by plotting the residual against the CPU time expended in the simulation. We believe that the ability to utilize efficient matrix–vector multiplication gives the method significantly improved convergence per second relative to the standard Kaczmarz algorithm, although other mechanisms may certainly be at work as well.

Conclusion and Related Work

The Kaczmarz algorithm was first proposed in . Kaczmarz demonstrated that the method converged to the solution of linear system Ax=b\bm{A}\bm{x}=\bm{b} for square, non-singular matrix A\bm{A}. Since then, the method has been utilized in the context of computer tomography as the Algebraic Reconstruction Technique (ART) . Empirical results suggested that randomized selection offered improved convergence over the cyclic scheme . Strohmer and Vershynin were the first to prove an expected linear convergence rate using a randomized Kaczmarz algorithm with specific random control. This result was extended by Needell to apply to inconsistent systems, which shows a linear convergence rate to within a fixed radius around the least-squares solution. Almost-sure convergence guarantees were recently proved by Chen and Powell . Zouzias and Freris analyze a modified version of the method in the inconsistent case, using a variant motivated by Popa to reduce the residual and thereby converge to the least squares solution. Relaxation parameters can also be introduced to obtain convergence to the least squares solution, see e.g. , and partially weighted sampling can lead to a tradeoff between convergence rate and radius . Liu, Wright, and Sridhar discuss applying a parallelized variant of the randomized Kaczmarz method, demonstrating that the convergence rate can be increased almost linearly by bounding the number of processors by a multiple of the number of rows of A\bm{A}.

The block Kaczmarz updating method was introduced by Elfving as a special case of the more general framework by Eggermont et.al. . The notion of using blocking in projection methods is certainly not new, and there is a large amount of literature on these types of methods, see e.g. and references therein. Needell and Tropp provide the first analysis showing an expected linear convergence rate which depends on the properties of the matrix A\bm{A} and of the submatrices Aτ\bm{A_{\tau}} resulting from the paving, connecting pavings and the block Kaczmarz scheme. The use of specialized blocks appears elsewhere, in particular, the works of Popa use blocks with orthogonal rows that are beneficial for the block Kaczmarz method . Needell, Zhao, and Zouzias expand on the results from and to demonstrate convergence to the least-squares solution for an inconsistent system using the block Kaczmarz method. Again the block approach can yield faster convergence than the simple method.

The Kaczmarz method was first applied to a system of equalities and inequalities by Leventhal and Lewis , who also consider polynomial constraints with the method. They give a linear convergence rate to the feasible solution space SS, using ∥A∥F2\lVert{\bm{A}}\rVert_{\rm F}^{2} and the Hoffman constant . We apply the block Kaczmarz scheme to the system described in , combining their result with that of Needell and Tropp to acquire a completely generalized result. We highlight several important complications which arise when attempting to apply the block scheme to inequalities. Nonetheless, whether a paving is used only partially or for the complete system, significant reduction in computational time can be achieved.

There are many interesting open problems related to the block Kazcmarz method and linear systems with inequalities. It has been well observed in the literature that selecting rows (or blocks) without replacement rather than with replacement as in the theoretical results leads to faster a convergence rate empirically . When selecting without replacement, independence between iterations vanishes, making a theoretical analysis more challenging. Secondly, it would be interesting to further investigate the use of obtuse row pavings. In systems with a large number of inequalities, the ability to pave the submatrix A≤\bm{A}_{\leq} with an obtuse row paving would lead to significantly faster convergence. In that case, one may like to identify a more general geometric property about the system that permits such pavings or an alternative formulation that offers convergence of the full block method.

Appendix A Proof of Theorem 3.2

Then since σ\sigma is part of an obtuse paving, the angle between xj−xj−1\bm{x}_{j}-\bm{x}_{j-1} and s−xj−1\bm{s}-\bm{x}_{j-1} must be obtuse. There thus exists a point t\bm{t} on the line segment L={γxj−1+(1−γ)s:0≤γ≤1}L=\{\gamma\bm{x}_{j-1}+(1-\gamma)\bm{s}:0\leq\gamma\leq 1\} such that xj−1−xj\bm{x}_{j-1}-\bm{x}_{j} and t−xj\bm{t}-\bm{x}_{j} are orthogonal (see Figure 4).

Now since t∈L\bm{t}\in L, we have ∥t−xj−1∥2≤∥xj−1−s∥2\lVert{\bm{t}-\bm{x}_{j-1}}\rVert_{2}\leq\lVert{\bm{x}_{j-1}-\bm{s}}\rVert_{2}, and thus letting θ\theta denote the angle between xj−xj−1\bm{x}_{j}-\bm{x}_{j-1} and t−xj−1\bm{t}-\bm{x}_{j-1} (see Figure 4), we have

By the definition of xj\bm{x}_{j}, this means that

Using this along with the paving properties we see that

Thus, taking expectation (over the choice of τ′\tau^{\prime}, conditioned on previous choices), yields

Combining this with (2) and letting E=E_{=} and E≤E_{\leq} denote the events that a block from TT and a block from T′T^{\prime} is selected, respectively, we have

Since p=βmβ′m′+βmp=\frac{\beta m}{\beta^{\prime}m^{\prime}+\beta m}, we have 1−pβ′m′=1β′m′+βm\frac{1-p}{\beta^{\prime}m^{\prime}}=\frac{1}{\beta^{\prime}m^{\prime}+\beta m} and we can simplify

where we have utilized the Hoffman bound (1.7) in the second inequality.

Iterating this relation along with independence of the random control completes the proof.

References