Fast Compute via MC Boosting

Modern training and inference pipelines in statistical learning and deep learning repeatedly invoke linear-system solves as inner loops, yet high-accuracy deterministic solvers can be prohibitively expensive when solves must be repeated many times or when only partial information (selected components or linear functionals) is required. We position \emph{Monte Carlo boosting} as a practical alternative in this regime, surveying random-walk estimators and sequential residual correction in a unified notation (Neumann-series representation, forward/adjoint estimators, and Halton-style sequential correction), with extensions to overdetermined/least-squares problems and connections to IRLS-style updates in data augmentation and EM/ECM algorithms. Empirically, we compare Jacobi and Gauss–Seidel iterations with plain Monte Carlo, exact sequential Monte Carlo, and a subsampled sequential variant, illustrating scaling regimes that motivate when Monte Carlo boosting can be an enabling compute primitive for modern statistical learning workflows.

💡 Research Summary

The paper “Fast Compute via MC Boosting” addresses a pervasive bottleneck in modern statistical learning and deep‑learning pipelines: the repeated solution of large linear systems inside inner loops of training, inference, and model compression. Traditional deterministic solvers (direct factorization, high‑order iterative methods) require O(m³) or O(m²·n) work and compute the full solution even when only a few components or linear functionals are needed. The authors propose Monte Carlo (MC) boosting as a practical alternative that directly estimates selected components or functionals, is embarrassingly parallel, and offers a tunable accuracy‑compute trade‑off through sampling and sequential residual correction.

The core technical development begins by rewriting a linear system AX = B (or an over‑determined least‑squares problem) in a fixed‑point Neumann‑series form X = L + H X, where H = I − GA and G is a nonsingular preconditioner that can be interpreted as a Markov transition kernel. Convergence is guaranteed when the spectral radius ρ(H) < 1. This representation enables two families of random‑walk estimators:

1. Forward estimator – Random walks start from a prescribed initial distribution p and follow transition probabilities P derived from H. The weight of a walk of length ℓ is the product of H entries along the path multiplied by the product of transition probabilities. The Monte Carlo estimator X̂ = ∑ℓ Wℓ b_{kℓ} is unbiased for the linear functional ⟨h, x⟩. It is most efficient when the goal is a single scalar quantity (e.g., a prediction at a test point).

2. Adjoint estimator – By swapping the roles of the right‑hand side b and the functional h, and using Hᵀ as the transition matrix, a single set of walks can simultaneously estimate all components of the solution vector. The initial distribution is aligned with b, and the walk termination state determines which component receives the contribution. This dramatically reduces the number of required simulations when the full solution is needed.

Both estimators suffer from high variance when used alone. To accelerate convergence, the authors revisit Halton’s sequential residual‑correction scheme, which they term Monte Carlo boosting. Given a current approximation Y, the residual D = L + HY − Y defines a new fixed‑point problem Z = D + HZ. By applying a Monte Carlo estimator to Z and updating Y←Y + Z, the residual magnitude shrinks geometrically, leading to far fewer samples for a target error. Two concrete algorithms are described:

• Exact Sequential MC – At each boosting iteration the full residual vector D is used, yielding the theoretically fastest error decay.
• Subsampled Sequential MC – Only a random subset of rows/columns is used to form an approximate residual, dramatically lowering per‑iteration cost while retaining most of the convergence benefit.

The paper extends the framework to over‑determined systems. The weighted least‑squares solution (\hat{x} = (L^{\top}\Omega L)^{-1} L^{\top}\Omega f) is shown to be equivalent to an iteratively re‑weighted least‑squares (IRLS) step. By embedding Monte Carlo boosting within the IRLS loop, the authors obtain a stochastic EM/ECM‑style M‑step that can handle non‑Gaussian likelihoods and data‑augmentation schemes. They also discuss a numerically stable block‑matrix formulation that removes rows/columns whose weights approach zero, effectively reducing problem dimension on the fly.

Empirical evaluation covers three representative settings:

• 2‑D Poisson equation – Comparing Jacobi, Gauss‑Seidel, plain MC, exact sequential MC, and subsampled sequential MC. Results show that for modest accuracy (≈10⁻²) plain MC is fastest, but as the error tolerance tightens (≈10⁻⁴) sequential MC outperforms deterministic iterations by a factor of 5–10.
• Kernel ridge regression – Demonstrates that the adjoint estimator can recover the entire solution vector with 2–3× fewer walks than the forward estimator.
• Post‑training quantization (GPTQ) – Shows that curvature‑aware updates requiring (H + λI)⁻¹ v or selected inverse columns are efficiently handled by MC boosting, enabling scalable quantization of large language models.

Table 1 summarizes theoretical complexities: direct factorization O(m³), stationary iterations O(m²·s), plain MC O(c·N) where c is the number of queried components, exact sequential MC O(c·N·s) with rapid error decay, and subsampled sequential MC O(c·N·s·α) with α < 1 reflecting subsampling.

The authors conclude that Monte Carlo boosting offers a new paradigm: rather than solving the full linear system, one iteratively refines stochastic estimates of only the information needed. This paradigm is especially attractive for large‑scale training, model compression, Bayesian posterior approximations, and real‑time inference where memory and latency constraints preclude full factorization. Future directions include extensions to nonlinear operators, adaptive sampling policies, and highly optimized GPU/TPU implementations.