1
Adaptive Row Selection Meets Asynchrony in Randomized Kaczmarz Evan Coleman Department of Computer Science University of Mary Washington Fredericksburg, Virginia, USA
arXiv:2607.08313v1 [cs.DC] 9 Jul 2026
Abstract Randomized Kaczmarz is a natural fit for large sparse least-squares and tomographic reconstruction, and adaptive row selection can reduce iteration counts. However, deploying adaptive selection on a shared-memory machine means sampling from a residual that lock-free workers are concurrently modifying, often using stale data. We present the first systematic study of this regime: residual-weighted and greedy Kaczmarz under asynchronous execution, measured across 339 runs on a 96-core node with realized (not injected) delays. Four findings carry directly to practice. (i) Stability is governed by a boundary ℓ∗ (T ) between sampling aggressiveness and thread count; below it, more aggressive sampling is strictly better, so one should tune to just inside the cliff. (ii) Threshold-greedy selection (the standard accelerated rule) is unstable at high thread counts, diverging almost immediately. (iii) Under-relaxation buys back the cliff at a predictable cost, giving a usable safety knob. (iv) Consistent-snapshot reads admit a rare, scheduling-dependent divergence that live (inconsistent) reads never exhibited and that is also cheaper, making inconsistent reads the right default. We validate the implementation against published sequential results and outline the distributed two-level sampler these measurements motivate.
Index Terms asynchronous iterative methods, randomized Kaczmarz, greedy selection, residual-weighted sampling, shared-memory parallelism
2
I. I NTRODUCTION Randomized Kaczmarz [1] is a row-action solver: each step projects the current iterate onto one equation’s hyperplane, costing O(nnz) of a single row with no factorization. That structure makes it a natural choice for large sparse least-squares and for applications such as tomographic reconstruction, where it is the algebraic backbone of the classical Algebraic Reconstruction Technique (ART) method. Adaptive selection (e.g., greedy [2] or residual-power-weighted [3], [4], [5]) accelerates convergence by preferring high-residual rows; asynchronous execution [6], [7] extracts parallelism by letting lock-free workers update shared state without barriers. Each of these mechanisms is well understood in isolation, but deploying them together creates a situation neither analysis covers. Adaptive selection must rank rows by residual magnitude; but on an asynchronous machine the residual is a shared vector that other workers are concurrently modifying (possibly with stale views). The selector is therefore choosing from an imperfect distribution, where the discrepancy from ideal grows with the number of in-flight updates. Existing adaptive-Kaczmarz theory assumes the residual is current; existing asynchronous theory assumes the sampling distribution is fixed. What actually happens between those assumptions has not been measured, and existing parallel implementations [8], [9], [10] do not analyze it. Using the power-weighted family P (i) ∝ |ri |ℓ (a single knob from uniform, ℓ → 0, through greedy, ℓ → ∞) we instrument residual-weighted and threshold-greedy Kaczmarz across 339 runs on a 96-core node with delays that are realized by the hardware rather than synthetically injected, on: (1) a benign dense problem, (2) a coherent tomographic problem, and (3) the standard sparse suite. We find that the read-consistency effect known from asynchronous Jacobi does not govern here [11]. What governs is a stability boundary between how aggressively selection chases residual peaks and how many workers are running, a cliff that aggressive sampling races toward as concurrency grows. The contributions of this paper are: •
A stability boundary ℓ∗ (T ): We map a boundary in the (sampling exponent ℓ, thread count T ) plane that decreases in T : more concurrency tolerates less aggressive sampling. Below the boundary greedier sampling is strictly faster so tuning amounts to a race to a concurrency-dependent cliff.
•
Greedy selection is asynchrony-incompatible at scale: Threshold-greedy selection [2], the standard accelerated rule, sits outside the boundary at high thread count: it diverges within a fraction of a sweep at T = 96 on every test problem, under either read semantics.
3
•
Step-size control recovers the cliff: Under-relaxation restores convergence past the boundary in our experiments: every configuration that diverges at full step (β = 1) converges at β ≤ 0.5 (21 of 21 runs), at a rate cost of β(2 − β) times a measured interference premium of roughly two.
•
Read-mode semantics as tail risk: Consistent-snapshot reads admit a rare, schedulingdependent divergence (one of 36 runs at a coherent cell; the same-seed survivor differs only in thread interleaving) that inconsistent (live) reads never exhibited (zero of 36) and that is also cheaper per step, which makes inconsistent reads the right default.
•
A validated, exact implementation: Sequential greedy/uniform speedups reproduce published results [2] where the row-normalization protocols coincide, and the maintained residual agrees with a recomputed b − Ax at every exit (to within 2.4 × 10−12 on converged runs and to within 3.8 × 10−15 in relative terms at the divergence flag), so every reported divergence is numerical behavior, not a data race. II. R ELATED W ORK
The Kaczmarz method projects the iterate onto one equation’s hyperplane per step. Its tomography connection is classical: Algebraic Reconstruction Technique (ART) applies this row-action projection view to tomographic projection equations [12]. Strohmer and Vershynin [1] showed that sampling row i with probability proportional to ∥ai ∥2 yields expected linear convergence at a rate set by the scaled condition number; Leventhal and Lewis [13] cast this as randomized coordinate descent on the normal equations, which is the view we use here. Adaptive rules sharpen the constant by steering selection toward large residuals. Bai and Wu [2] introduced greedy randomized Kaczmarz (GRK), which thresholds the residual and samples among the rows above the cutoff; Steinerberger [3] analyzed the continuous family P (i) ∝ |ri |p , recovering uniform sampling at p = 0 and the maximal-residual Motzkin rule as p → ∞, with the p = 1 case appearing in adaptive sketch-and-project [5] and the same power family extended to sparse recovery in a Bregman setting [14]. The Motzkin and sampling-Kaczmarz–Motzkin line [15], [16], [4] quantifies the acceleration through the residual’s dynamic range. Across all of this work, selection is computed from the current residual (a synchronous assumption). Parallelism and hardware optimization, where considered, typically remain synchronous: averaging independent projections [17] changes the iteration, while randomized block Kaczmarz [18] processes multiple
4
rows simultaneously to increase arithmetic intensity and mitigate the memory-bound nature of single-row updates. Asynchronous linear solvers trace to chaotic relaxation [19] and were given a modern asynchronous and randomized analysis by Avron, Druinsky, and Gupta [6], who proved convergence of shared-memory asynchronous updates under bounded staleness. More broadly, lockfree shared-memory stochastic updates were popularized by HOGWILD! [20], where sparsity makes overwrites tolerable. For Kaczmarz specifically, Liu, Wright, and Sridhar [7] analyzed an asynchronous parallel variant, but under two assumptions we revisit: a fixed (non-adaptive) sampling distribution, and (effectively) consistent reads, argued to be benign because each update touches only one row’s sparse support. Component-averaged row projection (CARP) [21] parallelizes row actions by a different route. Most relevant, recent parallel and distributed implementations [8], [9], [10] operate in exactly this setting, and RGDBEK [10] explicitly parallelizes greedy selection, but reports performance engineering without a staleness-aware analysis of when the two cooperate. The closest analytical precedent is our companion study of asynchronous residual-weighted Jacobi [11], where read consistency is the binding axis. III. A LGORITHM AND I MPLEMENTATION A. Power-weighted sampling We solve a consistent m × n system Ax = b, with rows normalized so ∥ai ∥ = 1 and b is scaled accordingly, following the Steinerberger/AsyRK convention [3], [7]. A relaxed Kaczmarz step on row i is x ← x + βri ai ,
r = b − Ax,
so β = 1 exactly satisfies equation i. Rows are sampled as P (i) ∝ |ri |ℓ . The exponent ℓ interpolates from uniform sampling (ℓ = 0, equivalent to Strohmer–Vershynin for unit rows) to the maximal-residual Motzkin rule (ℓ → ∞). We also test threshold-greedy GRK [2], and a deterministic cyclic baseline that visits rows in fixed order. To summarize how concentrated the residual is, we track its inverse participation ratio (IPR) ν 2 (r) = m∥r∥44 /∥r∥42 ∈ [1, m], which gives the effective number of large residual entries.
5
Algorithm 1 Per-worker loop (power-weighted, exponent ℓ, step β) 1: while not converged do 2:
d ← atomic load of global commit counter
3:
obtain residual view of r
4:
sample row i with P (i) ∝ |ri |ℓ
5:
riused ← if inconsistent then re-read ri else snapshot ri
6:
λ ← β riused
7:
atomic x += λ ai
▷ support of row i
8:
atomic r −= λ Qi,:
▷ Gram row = stencil
9:
c ← fetch-and-add commit counter; record delay c − d
▷ dispatch ▷ snapshot or live reads ▷ inverse CDF
10: end while
B. Residual maintenance Rather than recompute r = b − Ax, workers maintain it exactly via the Gram-row update of [3]. Committing λ = βriused on row i updates x ← x + λai ,
r ← r − λQi,: ,
where Q = AA⊤ and Qii = 1. Thus row i of Q is the set of equations coupled to row i through overlapping support, analogous to a stencil row in Jacobi. For this measurement study Q is precomputed once in sparse form. Each update costs O(m) for inverse-CDF sampling plus O(degQ (i)) for the commit; the sampling pass dominates. As a check, every run joins its workers at termination, recomputes b − Ax, and compares it to the maintained residual. C. Read semantics Workers are lock-free, so a worker may sample from a stale residual. We compare two policies. With consistent reads, the worker snapshots r once and uses that snapshot both to sample i and to form riused . With inconsistent reads, the sampling CDF is built from live atomic reads and ri is read again after i is selected. Thus the value used for the update may differ from the value that caused the row to be selected. Consistent reads pay an O(m) copy per step; GRK adds another source of staleness because its threshold and sampling passes may see different live residuals.
6
TABLE I T EST PROBLEMS AND MEASURED STATIC STRUCTURE ( EMPTY ROWS DROPPED , ROWS NORMALIZED , b = Ax⋆ ). degQ IS off
THE MEAN G RAM ROW DEGREE ( COUPLING DENSITY ); ρ̄2
THE MEAN OFF - DIAGONAL G RAM ROW 2- NORM ; χ̄ = E|Qij |
2
THE MEAN PAIRWISE COUPLING ; ν (b) THE INITIAL RESIDUAL CONCENTRATION . T HE S UITE S PARSE VALIDATION 3
MATRICES SPAN m = 10–2063, n = 10–24310, degQ = 4.76–400, AND 10 χ̄ = 1.8–768.
degQ ρ̄2off 103 χ̄
ν 2 (b)
Problem
m×n
gaussian
4,000 × 1,000
4,000 2.00
25.2
2.9
gaussian-4x
16,000 × 1,000 16,000 4.00
25.2
2.9
tomo
7,546 × 4,096
3,016 2.58
12.0 4.5–6.7
tomo (point)
7,546 × 4,096
3,016 2.58
12.0
47.0
29,980 × 16,384 11,866 2.82
6.0
5.3
tomo-4x
At T =96: τ χ̄ = 2.40 (gaussian), 1.14 (tomo), 0.57 (tomo-4x).
IV. E XPERIMENTAL S ETUP All experiments ran on a 96-core ACES node (Texas A&M), with threads pinned using OMP_PROC_BIND=spread and OMP_PLACES=cores. Staleness was measured, not injected: each worker records the global commit counter when it begins a step and again when it commits, and we report the realized dispatch-to-commit gap. Thus the delay distribution is a property of the hardware and schedule. All test systems are consistent by construction: we fix x⋆ and set b = Ax⋆ , so both the relative residual ∥r∥/∥b∥ and the relative solution error ∥x − x⋆ ∥2 /∥x⋆ ∥2 are available. Runs stop at ∥r∥/∥b∥ < 10−6 , are flagged divergent when ∥r∥/∥b∥ > 108 , and are capped by a problem-dependent sweep limit; one sweep denotes m commits, a full pass over the rows in expectation. Paired read-mode comparisons use matched seeds, so each consistent/inconsistent pair shares the matrix, right-hand side, and per-thread random streams and differs only in read policy. Table I summarizes the test problems: a benign dense system, coherent tomographic systems at two scales, and 16 SuiteSparse matrices from the Bai–Wu test suite [2]. The tomographic systems use a parallel-beam geometry with 90 angles × 93 detector bins on a 64 × 64 grid (tomo-4x: 180 × 185 on 128 × 128); rays that miss the grid are dropped. In total we report 339 instrumented runs. At termination, every run joins its workers, recomputes b − Ax, and compares it to the maintained residual; the maximum discrepancy across converged runs is 2.4 × 10−12 .
7
V. R ESULTS A. Sampling aggressiveness under asynchrony We begin with the observation that motivates the rest of the paper. Figure 1 compares cyclic, uniform, power-weighted sampling with ℓ ∈ {0.5, 1, 2, 4}, and threshold-greedy GRK at T = 96 under inconsistent reads on the two built-in problems. Synchronous theory suggests that stronger residual bias should reduce iteration counts. Under asynchrony, that ordering holds only up to a stability cliff. On the gaussian problem, cyclic, uniform, and ℓ ≤ 1 converge; ℓ = 1 reaches tolerance in 5.8 sweeps, roughly twice as fast as uniform’s 11.4 sweeps. In contrast, ℓ ≥ 2 diverges. GRK is the extreme case: on both problems and under either read policy, it fails within a fraction of a sweep, reaching residuals of 1012 –1021 . Its threshold rule concentrates selection on the largest residual entries, so many workers repeatedly update the same coupled neighborhoods. The ℓ = 2 gaussian run shows the mechanism more clearly because it first behaves like a successful solve. The residual falls to about 3 × 10−2 within one sweep while the concentration statistic ν 2 remains in the range 3–6; then the run diverges within the next half-sweep as ν 2 rises to 130–170. We interpret this as an interference floor. While the true residual is large, it dominates the perturbations caused by the Θ(T ) in-flight updates. Once the residual reaches the same scale as those perturbations, a residual-weighted sampler begins ranking interference rather than signal. The selected neighborhoods are then over-updated, the perturbation grows, and the feedback loop closes. Uniform and cyclic sampling do not rank residual entries, so they do not amplify this floor in the same way. This first experiment therefore points to sampling aggressiveness, rather than read consistency alone, as the axis that must be mapped next. B. The stability boundary ℓ∗ (T ) Figure 2 maps the outcome over the (ℓ, T ) plane, with three seeds per cell. A concurrencydependent boundary ℓ∗ (T ) separates convergence from divergence, and it moves downward as T grows: more workers tolerate less aggressive residual bias. On gaussian, ℓ = 3 is stable at T = 8, mixed at T = 32 (two seeds converge and one diverges), and unstable at T = 96; the high-concurrency cliff lies between ℓ = 1.5 and ℓ = 2. On tomo, the boundary is shifted upward: ℓ = 3 remains stable through T = 32 and fails only at T = 96. In the tomo panel of Fig. 2, ℓ ≤ 2 cells report the residual at the 200-sweep cap, essentially at tolerance; the ℓ = 3, T ≤ 32 cells converge outright in ∼177 sweeps.
8
Fig. 1. Relative-residual trajectories by sampling rule at T = 96 under inconsistent reads. On gaussian, ℓ = 1 is the fastest survivor while ℓ ≥ 2 and GRK diverge. On tomo, the stability cliff occurs at a larger exponent, but GRK again fails rapidly. The optimum is therefore interior: residual bias helps until asynchrony drives the sampler past a problem-dependent cliff.
Below the boundary, the synchronous ordering largely survives. On gaussian, median sweep counts decrease as ℓ increases within the stable region, so the best setting is not uniform sampling but the largest exponent still inside the cliff. Thus tuning becomes a race to a problem- and concurrency-dependent boundary. The static statistics in Table I indicate what moves that boundary. The normalized off-diagonal √ Gram-row norm ρ̄off 2 / m is nearly the same for gaussian and tomo and is m-independent for gaussian, yet their cliffs differ by about one unit of ℓ. The better discriminator is the mean pairwise coupling χ̄ = E|Qij |: 103 χ̄ = 25.2 for gaussian, 12.0 for tomo, and 6.0 for tomo4x, exactly matching the observed stability ordering. With mean realized delay τ ≈ T − 1, the corresponding interference scale τ χ̄ at T = 96 is 2.40, 1.14, and 0.57; Fig. 6 tests this discriminator suite-wide. This also explains why gaussian-4x rescues ℓ = 2 at T = 96, changing from 0/6 converged runs to 6/6. Its pairwise coupling is the same as gaussian, so the interference scale has not changed. What changes is the convergence clock: with m/n = 16, the system reaches tolerance in roughly one sweep, before the interference floor that destroys the smaller gaussian instance can dominate. The stability cliff is therefore a race between convergence and interference growth.
9
Fig. 2. Stability boundary ℓ∗ (T ) over the (sampling exponent, thread count) plane; three seeds per cell, green annotated with median sweeps to tolerance (or relative residual at the 200-sweep cap). The boundary decreases with concurrency, and below it the synchronous ordering survives (greedier is strictly faster) so tuning is a race to a concurrency-dependent cliff. The mixed cell (two converge, one diverges) shows the boundary is probabilistic.
C. Step-size control at the cliff If instability occurs because workers over-commit to an interference-contaminated residual, then reducing the projection step should improve stability. Figure 3 reruns the cliff cells from Figure 2 with smaller relaxation parameter β. Every configuration that diverges at full step, β = 1, converges for some β ≤ 0.5; across the 21 tested rescue runs, none diverged. The metastable gaussian cell at (ℓ = 3, T = 32) also improves from two-of-three to three-of-three converged seeds. The recovery is not free. In the synchronous Kaczmarz iteration, under-relaxation predicts a sweep-count penalty proportional to 1/[β(2 − β)]. The observed penalty is larger by about 1.7–2.3×, with the excess growing as β decreases. Thus β damps the instability, but it also slows the correction of the rows that the stale sampler continues to over-select. This extra cost is the asynchronous interference premium. The engineering conclusion is that β is a reliable safety knob: when concurrency pushes
10
Fig. 3. Sweeps to tolerance versus relaxation parameter β at cliff cells that diverge for β = 1. Dotted curves show the synchronous 1/[β(2 − β)] prediction. Under-relaxation restores convergence in all tested rescue runs, but with a 1.7–2.3× asynchronous premium beyond the synchronous prediction.
Fig. 4. Scaling on the headline tomo instance. Wall-clock speedup saturates at 4.4× because the O(m) global sampling pass is memory-bound. In contrast, sweeps to tolerance and sustained ν 2 remain nearly flat through T = 96, so asynchrony is cheap in iteration terms but not in wall-clock time.
the desired sampling exponent beyond the stability boundary, under-relaxation can buy back convergence at a measurable rate cost. The tomo β = 0.25 runs are censored at the 600-sweep cap but continue descending through 1.3–1.5 × 10−5 , and are marked accordingly in Figure 3. D. Read-mode semantics as tail risk Read consistency does not set the stability boundary, but it does affect the tail behavior. In the paired tomo runs at T = 96, ℓ = 2, de-censored to 400 sweeps, one consistent-read run diverges catastrophically at sweep 2.3 (∥r∥/∥b∥ ≈ 5 × 1027 ), while its matched inconsistent-read twin converges. The pair uses the same matrix, right-hand side, and per-thread random streams;
11
TABLE II C ATASTROPHIC DIVERGENCES AT THE COHERENT CELL ( TOMO , T =96, ℓ=2), WITH EXACT 95% C LOPPER –P EARSON INTERVALS . T HE INTERVALS OVERLAP : THE CLAIM IS EXISTENCE AND ASYMMETRY, NOT A RATE DIFFERENCE .
read mode
catastrophic / runs
rate
95% CI
consistent
1 / 36
2.8%
[0.1%, 14.5%]
inconsistent
0 / 36
0%
[0%, 9.7%]
it differs only in read policy and thread interleaving. Across all matched tomo runs, consistent reads had 1 catastrophic divergence in 36 trials, while inconsistent reads had 0 in 36 (Table II). The confidence intervals overlap, so we do not estimate a failure rate. The result is instead an existence proof: snapshot reads admit a scheduling-dependent early-divergence mode that live reads did not exhibit in our runs. This observation also favors live reads on cost. Consistent reads copy the full residual once per step, whereas live reads build the sampling CDF directly from atomic loads and re-read the selected entry before committing. Since live reads were cheaper and showed no tail failures in these matched tests, they are the better default for the broad-coupling regime studied here. Finally, we tested whether the collision-dominated Jacobi regime appears in tomography by using a point phantom with much higher initial concentration, ν 2 (b) ≈ 47. Both read policies converged in all six runs, with similar sweep counts. The protection is structural: a point object’s sinogram is supported on only a few rays per angle but across all angles, so the Radon transform caps achievable concentration at roughly m/nangles (≈ 84 here, against the measured 47) — orders of magnitude below the ν 2 ∼ m that same-row pile-up requires. Thus, for these Kaczmarz problems, the main failure mode is not read inconsistency itself but residual-adaptive sampling of an interference floor. E. Implementation soundness and scaling The divergences above are properties of the asynchronous iteration, not residual-maintenance errors. On every converged run, the maintained residual agrees with a recomputed b − Ax to within 2.4 × 10−12 ; and in validation re-runs of every divergent (ℓ, T ) cell and every greedy configuration (33 runs spanning both read modes and all five test problems), the maintained and recomputed residuals at the divergence flag (where ∥r∥/∥b∥ reaches 108 –1082 ) agree to within
12
Fig. 5. Six same-configuration consistent-reads runs, tomo, T =96, ℓ=2. Five converge; one diverges catastrophically at sweep 2.3 (residual 5 × 1027 ), while its same-seed twin survives; the runs differ only in thread interleaving, so the failure is schedulingdependent, not deterministic.
3.8 × 10−15 in relative terms: the explosive trajectories are computed exactly. The sequential implementation also reproduces published GRK-over-uniform iteration-count ratios [2] on the near-constant-row-norm matrices where the normalization protocols coincide: bibd gives ratios 2.9/3.0/3.0 versus their 3.2/3.1/3.3, df2177 gives 4.9 versus 5.0, crew1 gives 5.0 versus 4.5, WorldCities gives 6.1 versus 5.7, and cari gives 4.2 versus 4.7. The main disagreements occur on heterogeneous-norm or ill-conditioned matrices, where row normalization changes the effective problem and residual stopping need not track solution error. Figure 4 separates iteration behavior from wall-clock behavior on the headline tomo instance. Sweeps to tolerance are nearly flat from T = 1 to T = 96, and the sustained concentration statistic remains ν 2 ≈ 3.7–4.0, so asynchrony is essentially free in iteration terms. Wall-clock scaling is much weaker: speedup reaches only 4.4× on 96 threads, and T = 2 is slower than T = 1. The bottleneck is the O(m) atomic-read pass used to build the global sampling CDF, which is memory-bound and dominates the sparse Gram-row commit. Rows per thread and coupling are therefore joint stability resources, and both are computable from A before any run. Figure 6 plots the 16 SuiteSparse matrices in the (m/T, τ χ̄) plane
13
Fig. 6. Asynchronous survival predicted from two static quantities (uniform sampling, T =96): rows-per-thread m/T vs. in-flight coupling τ χ̄; green converged, red diverged. crew1 and bibd_17_8 share m/T =1.4 but differ in coupling and fate.
against their T = 96 uniform-sampling outcomes. Below roughly one row per thread nothing survives, regardless of coupling; above that floor, survival tracks τ χ̄. The sharpest evidence is a natural controlled pair: crew1 and bibd_17_8 sit at identical m/T = 1.4 with opposite outcomes, separated only by coupling (τ χ̄ of 4.4 versus 19.0); WorldCities, the lone high-m/T casualty, carries the largest coupling in the suite (τ χ̄ = 32). VI. C ONCLUSION Asynchronous Kaczmarz with residual-adaptive selection is governed by coupling structure. Unlike the narrow-coupling, collision-bound regime of asynchronous Jacobi [11], the broadcoupling regime of Kaczmarz (e.g., dense least-squares and tomography) is bound by the product of sampling aggressiveness and coupling density, which creates a neighborhood interference floor. Measurements identify the discriminator concretely: the mean pairwise coupling χ̄ =
14
√ E|Qij | strictly orders the stability boundary, whereas spectral proxies (ρ̄off 2 / m) and worst-case coherence fail to predict it. The results here suggest it would be optimal to default to inconsistent reads, which showed no tail divergences in our runs and cost less per step; tune the sampling exponent just inside the boundary ℓ∗ (T ) where aggressiveness is strictly beneficial; recover the cliff with under-relaxation if concurrency forces a smaller ℓ than desired; and budget at least O(10) rows per thread. Two directions seem ideal for further exploration: (1) developing the staleness-aware theory behind the ℓ∗ (T ) boundary and super-quadratic under-relaxation premium, and (2), the observed O(m) global sampling bottleneck motivates work on a distributed two-level sampler, where block-level masses refreshed by asynchronous reductions explicitly decouple sampling staleness from iterate staleness. ACKNOWLEDGMENT This work leveraged the ACES cluster at Texas A&M University under allocation CIS250436 from the ACCESS program (NSF grants #2138259, #2138286, #2138307, #2137603, #2138296). R EFERENCES [1] T. Strohmer and R. Vershynin, “A randomized Kaczmarz algorithm with exponential convergence,” Journal of Fourier Analysis and Applications, vol. 15, no. 2, pp. 262–278, 2009. [2] Z.-Z. Bai and W.-T. Wu, “On greedy randomized Kaczmarz method for solving large sparse linear systems,” SIAM Journal on Scientific Computing, vol. 40, no. 1, pp. A592–A606, 2018. [3] S. Steinerberger, “A weighted randomized Kaczmarz method for solving linear systems,” Mathematics of Computation, vol. 90, no. 332, pp. 2815–2826, 2021. [4] J. Haddock and A. Ma, “Greed works: An improved analysis of sampling Kaczmarz–Motzkin,” SIAM Journal on Mathematics of Data Science, vol. 3, no. 1, pp. 342–368, 2021. [5] R. M. Gower, D. Molitor, J. Moorman, and D. Needell, “On adaptive sketch-and-project for solving linear systems,” SIAM Journal on Matrix Analysis and Applications, vol. 42, no. 2, pp. 954–989, 2021. [6] H. Avron, A. Druinsky, and A. Gupta, “Revisiting asynchronous linear solvers: Provable convergence rate through randomization,” Journal of the ACM, vol. 62, no. 6, pp. 51:1–51:27, 2015, preprint: arXiv:1304.6475. [7] J. Liu, S. J. Wright, and S. Sridhar, “An asynchronous parallel randomized Kaczmarz algorithm,” arXiv preprint arXiv:1401.4780, 2014. [8] J. Wang, M. Tian, Y. Wang, G. He, and T. Liu, “PRKP: A parallel randomized iterative algorithm for solving linear systems,” in 2022 IEEE 24th Int Conf on High Performance Computing & Communications; 8th Int Conf on Data Science & Systems; 20th Int Conf on Smart City; 8th Int Conf on Dependability in Sensor, Cloud & Big Data Systems & Application (HPCC/DSS/SmartCity/DependSys).
IEEE, 2022, pp. 244–249.
[9] E. S. Bölükbaşı, F. Ş. Torun, and M. Manguoğlu, “A distributed memory parallel randomized Kaczmarz for sparse system of equations,” Concurrency and Computation: Practice and Experience, vol. 36, no. 25, p. e8274, 2024.
15
[10] A. Panchal and R. Behera, “RGDBEK: Randomized greedy double block extended Kaczmarz algorithm with hybrid parallel implementation and applications,” arXiv preprint arXiv:2509.19267, 2025. [11] E. Coleman, “Residual-weighted randomized Jacobi: Sharpened bounds via residual concentration and asynchronous extension,” arXiv preprint arXiv:2606.01232, 2026. [12] R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography,” Journal of Theoretical Biology, vol. 29, no. 3, pp. 471–481, 1970. [13] D. Leventhal and A. S. Lewis, “Randomized methods for linear constraints: convergence rates and conditioning,” Mathematics of Operations Research, vol. 35, no. 3, pp. 641–654, 2010. [14] L. Zhang, Z. Yuan, H. Wang, and H. Zhang, “A weighted randomized sparse Kaczmarz method for solving linear systems,” Computational and Applied Mathematics, vol. 41, p. 383, 2022. [15] J. Nutini, B. Sepehry, I. Laradji, M. Schmidt, H. Koepke, and A. Virani, “Convergence rates for greedy Kaczmarz algorithms, and faster randomized Kaczmarz rules using the orthogonality graph,” in Proceedings of the 32nd Conference on Uncertainty in Artificial Intelligence (UAI), 2016. [16] J. Haddock and D. Needell, “On Motzkin’s method for inconsistent linear systems,” BIT Numerical Mathematics, vol. 59, no. 2, pp. 387–401, 2019. [17] J. D. Moorman, T. K. Tu, D. Molitor, and D. Needell, “Randomized Kaczmarz with averaging,” BIT Numerical Mathematics, vol. 61, no. 1, pp. 337–359, 2021. [18] D. Needell, R. Zhao, and A. Zouzias, “Randomized block Kaczmarz method with projection for solving least squares,” Linear Algebra and its Applications, vol. 484, pp. 322–343, 2015. [19] D. Chazan and W. Miranker, “Chaotic relaxation,” Linear Algebra and its Applications, vol. 2, no. 2, pp. 199–222, 1969. [20] B. Recht, C. Re, S. Wright, and F. Niu, “HOGWILD!: A lock-free approach to parallelizing stochastic gradient descent,” Advances in Neural Information Processing Systems, vol. 24, 2011. [21] D. Gordon and R. Gordon, “Component-averaged row projections: a robust, block-parallel scheme for sparse linear systems,” SIAM Journal on Scientific Computing, vol. 27, no. 3, pp. 1092–1117, 2005.