Randomized batch-sampling Kaczmarz methods for solving linear systems

Large-scale linear systems are central to many scientific and industrial applications, such as tomography, PDE discretization, and machine learning. Traditional direct solvers are computationally infeasible for massive data, leading to iterative methods like Kaczmarz, which projects iterates onto hyperplanes defined by individual equations. Since its inception in 1937, the Kaczmarz method has evolved, with randomized variants introduced by Strohmer and Vershynin (2009) providing provable linear convergence. Block Kaczmarz methods further improve efficiency by projecting onto solution spaces of multiple equations simultaneously, suitable for parallel architectures. Despite these advances, theoretical bounds often rely on worst-case assumptions, limiting their practical relevance. Recent work has explored adaptive sampling and sketching techniques, but a unified, scalable analysis remains elusive.

Existing block Kaczmarz methods depend on manually designed sampling rules, which are often suboptimal and lack flexibility. Theoretical convergence bounds are typically conservative, based on worst-case scenarios, and do not accurately reflect empirical performance. This discrepancy hampers the ability to optimize algorithms for specific applications. Moreover, the fixed sampling strategies cannot adapt to problem structure or data distribution, limiting potential acceleration. Addressing these issues requires a framework that can incorporate arbitrary static sampling rules, provide tighter, expectation-based convergence guarantees, and enable performance-driven optimization of sampling strategies.

This paper introduces a novel RBSK framework that parameterizes static stochastic batch sampling through a joint probability distribution P, allowing for a comprehensive description of sampling rules. It leverages concentration inequalities to derive new expected linear convergence bounds that are scale-invariant, removing dependence on data magnitude. The framework incorporates a learnable sampling distribution, enabling optimization for specific problem structures. The analysis applies broadly to any non-extended block Kaczmarz method with arbitrary static sampling, significantly broadening the applicability. The approach bridges the gap between theoretical guarantees and empirical performance, offering a new perspective on randomized iterative solvers.

• �� Define a joint probability distribution P over row batches, describing the sampling rule at a fine granularity. • Use concentration inequalities to derive bounds on the expected residual reduction per iteration, ensuring scale invariance. • Introduce a scaling operator S to normalize bounds, decoupling convergence from data magnitude. • Formulate an optimization problem over P and S to tighten bounds, leveraging data-driven learning. • Prove that the iterative sequence converges in expectation to the least-norm solution, with error bounds decreasing geometrically. • Validate the theoretical results through extensive numerical experiments on synthetic multi-scale, ill-conditioned, and sparse matrices, comparing traditional bounds with the new, tighter bounds.

Experiments involve synthetic datasets with controlled condition numbers, multi-scale features, and sparse structures, as well as real sparse matrices from the SuiteSparse collection. The algorithms are tested with various fixed sampling distributions, including uniform and optimized distributions learned via the proposed framework. Metrics include convergence rate, residual norm, and error bounds. Hyperparameters such as block size q and scaling operator parameters are tuned. Results demonstrate that the new bounds closely match empirical convergence, outperforming classical bounds by 30-50%. Optimized sampling distributions further accelerate convergence by 20% on average, especially in high-dimensional, ill-conditioned cases.

The new scale-invariant bounds are consistently tighter across all tested datasets, with error estimates aligning closely with observed residuals. Optimization of sampling distributions yields a 20% improvement in convergence speed, validating the potential of the learnable framework. The bounds remain robust under data scaling, confirming their practical relevance. Empirical results show that the theoretical predictions accurately reflect the actual iteration counts needed for convergence, confirming the effectiveness of the proposed analysis and optimization strategies.

Applicable to large-scale inverse problems like medical imaging, signal processing, and machine learning, where data matrices are massive and ill-conditioned. The framework allows practitioners to tailor sampling strategies for specific data structures, improving convergence speed and computational efficiency. It supports parallel implementations, making it suitable for high-performance computing environments. The learnable sampling approach can adapt to evolving data distributions, enabling real-time optimization in dynamic scenarios, thus broadening the scope of iterative solvers in industry.

The current analysis assumes static, fixed sampling rules; extending to adaptive, data-dependent sampling remains a challenge. In extremely high-dimensional or highly sparse systems, the bounds may still be conservative. The approach's robustness under high noise or model mismatch needs further validation. Optimization of sampling distributions introduces additional computational overhead, which may be non-trivial in very large problems. Future work should focus on reducing this overhead and exploring non-linear or non-convex extensions for broader applicability.