Stochastic Coordinate Sampling

• Stochastic Coordinate Sampling is a family of methods that updates selected coordinates or blocks, reducing computational overhead while maintaining convergence properties.
• It incorporates diverse sampling designs—uniform, non-uniform, adaptive, and blockwise—to optimize convergence rates and manage bias in stochastic gradient estimations.
• SCS is applied in both optimization and lattice QCD, efficiently computing high-dimensional gradients and Wick contractions for isospin-breaking corrections.

Stochastic Coordinate Sampling (SCS) denotes a family of stochastic procedures in which only selected coordinates, coordinate blocks, or coordinate-space points are processed at a given iteration or estimator evaluation, rather than a full gradient, full variable vector, or full all-to-all sum. In optimization, the paradigm appears as randomized coordinate descent, stochastic coordinate gradient descent, and primal-dual coordinate methods embedded in broader proximal SGD or saddle-point frameworks; in lattice QCD, it is used to sample source, sink, and internal vertex positions in coordinate space when constructing Wick contractions for isospin-breaking corrections to hadronic vacuum polarization (Gorbunov et al., 2019, Leluc et al., 2021, Bruno et al., 27 Feb 2026).

1. Scope, terminology, and historical placement

The optimization literature does not present SCS as a single standardized algorithm. In a unified proximal-SGD theory, SCS/RCD is not treated as a separate optimization class; instead, it is one instance of proximal SGD in which the stochastic direction is built by sampling coordinates and, in some cases, maintaining a memory or control-variate state (Gorbunov et al., 2019). A related non-convex framework places stochastic coordinate gradient descent between classical SGD and coordinate descent: instead of updating the entire parameter vector with a full stochastic gradient, it updates only selected coordinates of a stochastic gradient estimate (Leluc et al., 2021).

Non-uniform and arbitrary sampling entered the coordinate-sampling literature early through explicit probability design. A 2013 paper introduced NSync, a nonuniform synchronous parallel coordinate descent method in which a random subset of coordinates is updated in parallel and the subsets may be chosen non-uniformly; a 2017 paper on SPDC then studied primal-dual type stochastic optimization algorithms with non-uniform sampling, presented a convergence analysis of SPDC with arbitrary sampling, and proposed Optimality Violation-based Sampling SPDC (ovsSPDC) together with heuristic variants ovsSPDC+ and ovsSPDC++ [(Richtárik et al., 2013); (Shibagaki et al., 2017)].

Terminology is not uniform across fields. In one 2025 paper on two-stage stochastic programming, the acronym SCS denotes Stochastic Conjugate Subgradient, not Stochastic Coordinate Sampling, and the paper explicitly is not about coordinate selection or coordinate-wise updates (Zhang et al., 26 Mar 2025). This acronym collision matters because the coordinate-sampling literature is otherwise closely tied to randomized coordinate descent, stochastic dual coordinate ascent, and coordinate-subsampled SGD.

2. Core mechanism in stochastic optimization

A convenient unifying abstraction is the proximal-SGD update

$$x^{k+1}=\operatorname{prox}_{\gamma R}\bigl(x^k-\gamma g^k\bigr), \qquad E[g^k\mid x^k]=\nabla f(x^k),$$

where SCS enters through the construction of the stochastic direction $g^k$. In the coordinate-descent representative SEGA, only one coordinate $i\in[d]$ is sampled at each iteration: $$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$ The subsequent step is again a proximal-SGD step, so the coordinate choice is encoded in the estimator rather than in a separate algorithmic class (Gorbunov et al., 2019).

In stochastic coordinate gradient descent, the update is written directly at the parameter level: $$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$ with conditional coordinate policy

$$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$

This separates data sampling and coordinate sampling, and the framework assumes

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$

The same framework allows biased zeroth-order gradient estimators through a parameter $h$ controlling the bias, with

$$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$

A central theoretical observation is that, under $c=0$, the only unbiased coordinate policy is uniform sampling,

$g^k$0

Hence adaptive coordinate policies generally introduce bias in the coordinate-selection mechanism (Leluc et al., 2021).

3. Sampling design: uniform, non-uniform, adaptive, and blockwise

In NSync, the sampling object is a random subset $g^k$1 with law

$g^k$2

and parallel coordinate update

$g^k$3

The induced marginal coordinate probabilities are

$g^k$4

and the convergence parameter is governed by

$g^k$5

Under the nonuniform Expected Separable Overapproximation assumption and weighted strong convexity, the choice of $g^k$6 directly controls the rate. In the serial case, the optimal probabilities are

$g^k$7

while for uniform sampling

$g^k$8

The same paper notes that the optimal serial method can be faster than the fully parallel method in iteration complexity when the partial separability degree $g^k$9 is large enough (Richtárik et al., 2013).

Adaptive sampling replaces fixed probabilities by state-dependent probabilities. In AdaSDCA, the sampling rule is built from the dual residue

$i\in[d]$0

which vanishes exactly when coordinate $i\in[d]$1 satisfies the dual optimality relation. The relaxed optimal adaptive distribution is

$i\in[d]$2

The practical variant AdaSDCA+ computes $i\in[d]$3 only at epoch boundaries and then decreases the probability of a sampled coordinate by a factor $i\in[d]$4 within the epoch (Csiba et al., 2015).

Arbitrary and non-uniform sampling also arise in primal-dual coordinate methods. The SPDC paper of 2017 studies arbitrary sampling and proposes Optimality Violation-based Sampling SPDC, explicitly tying the sampling distribution to optimality violation rather than to fixed importance weights (Shibagaki et al., 2017). This suggests that in SCS, probability design is itself a first-class algorithmic component rather than a minor implementation choice.

4. Convergence theory and provable comparisons

A general convergence template for coordinate-sampled proximal methods is given by the two-recursion assumption

$i\in[d]$5

$i\in[d]$6

where

$i\in[d]$7

and $i\in[d]$8 tracks the quality of the coordinate estimator. Under $i\in[d]$9-strong quasi-convexity, the Lyapunov function

$$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$0

contracts linearly. For exact coordinate methods with $$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$1, the method converges linearly to the exact optimum; if $$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$2 or $$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$3 is positive, it converges only to a neighborhood. For SEGA, the framework yields

$$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$4

and the corollary

$$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$5

for $$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$6 and $$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$7 (Gorbunov et al., 2019).

NSync gives a complementary high-probability linear guarantee. If

$$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$8

then

$$h^{k+1}=h^k+e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr), \qquad g^k=d\,e_i\bigl(\nabla_i f(x^k)-h_i^k\bigr)+h^k.$$9

The same analysis proves the lower bound

$$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$0

showing that the expected subset size alone does not determine complexity; the interaction between sampling probabilities, coordinate curvature, and the ESO constants remains decisive (Richtárik et al., 2013).

Sampling without replacement can alter the asymptotic theory. For positive-definite quadratics with unit diagonal, a 2025 paper proves that random-permutation coordinate descent (RPCD) has a strictly better asymptotic contraction factor than random coordinate descent (RCD) on a structured class of permutation-invariant Hessians. The general RCD lower bound is

$$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$1

the stronger lower bound on the structured class is

$$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$2

and the RPCD upper bound is

$$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$3

For $$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$4, the RPCD bound is strictly smaller than the RCD lower bound on the same class. The paper also emphasizes the limitations: the exact RPCD upper bound is proved only for a structured quadratic class, and the extension to all positive definite quadratics is conjectural (Kim et al., 29 May 2025).

5. Major algorithmic families and specialized extensions

Primal-dual SCS is exemplified by the doubly stochastic primal-dual coordinate method for ERM. Starting from the bilinear saddle-point formulation

$$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$5

the method samples a subset $$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$6 of dual coordinates with $$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$7 and a subset $$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$8 of primal coordinates with $$\theta_{t+1}=\theta_t-\gamma_{t+1}D(\zeta_{t+1})g_t, \qquad D(k)=e_k e_k^\top,$$9, both uniformly. The coupling strength is summarized by

$$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$0

The iteration complexity is summarized as

$$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$1

The method is reported to have lower overall complexity than existing coordinate methods when the data matrix has factorized structure or when the proximal mapping on each block is computationally expensive, such as an eigenvalue decomposition (Yu et al., 2015).

Adaptive SCS has also been developed for non-convex and zeroth-order regimes. In stochastic coordinate gradient descent, MUSKETEER implements the “all for one / one for all” mechanism: past noisy gradients are aggregated into gains $$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$2, normalized either by $$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$3 normalization or softmax, and then mixed with uniform exploration,

$$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$4

The framework proves almost sure convergence under an extended Robbins–Monro condition and gives non-asymptotic bounds under the Polyak–Łojasiewicz condition. The experiments reported gains on synthetic block-structured problems, zeroth-order optimization, and neural-network training; on MNIST, the reported test accuracy improves from $$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$5 for SGD to $$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$6 with $$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$7 (Leluc et al., 2021).

Privacy introduces an additional constraint specific to coordinate methods: auxiliary memory. In DP-SCD for generalized linear models,

$$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$8

must be privatized together with the model variable because it is explicitly data-dependent. The key observation is that with independent zero-mean Gaussian noise added to both sequences, consistency holds in expectation,

$$d_t^{(k)}=\mathbb P(\zeta_{t+1}=k\mid \mathcal F_t), \qquad D_t=\operatorname{Diag}(d_t^{(1)},\ldots,d_t^{(p)}).$$9

The independent mini-batch coordinate construction yields sensitivity

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$0

whereas the sequential correlated baseline has sensitivity

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$1

The privacy guarantee is stated as

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$2

with $$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$3. The paper argues that decoupling and parallelizing coordinate updates is essential for utility (Damaskinos et al., 2020).

For linear systems, the SSD framework recovers Sampling Coordinate Descent and Capped Coordinate Descent through sketch choices such as $$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$4, introduces Greedy Sketching Rule and Greedy Capped Sketching Rule, incorporates Polyak momentum, and proves both global convergence and $$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$5 convergence for Cesàro averages (Morshed et al., 2020).

6. Coordinate-space SCS in lattice QCD

Outside optimization, stochastic coordinate sampling has a distinct coordinate-space meaning in lattice QCD. In the RBC/UKQCD calculation of isospin-breaking corrections to the hadronic vacuum polarization contribution to $$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$6, SCS is the practical Monte Carlo strategy used to evaluate the many coordinate-space Wick contractions needed for the electromagnetic and strong isospin-breaking terms without computing a prohibitively expensive full all-to-all propagator on every configuration (Bruno et al., 27 Feb 2026).

The starting observable is the time-momentum representation

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$7

with

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$8

and

$$\zeta_{t+1}\ \perp\!\!\!\perp\ \xi_{t+1}\mid \mathcal F_t.$$9

The RM123 expansion decomposes $h$0 into isoQCD, QED, and strong isospin-breaking pieces plus higher-order terms. The electromagnetic correction requires a four-point function with two internal electromagnetic vertices $h$1, while the strong isospin-breaking correction is implemented by inserting scalar-density operators and expanding around the isoQCD point. In both sectors, the relevant objects are correlators with extra operator insertions, and SCS samples subsets of source, sink, and internal coordinates together with correction factors that account for the finite number of samples.

The implementation uses a dataset of precomputed point-source propagators for randomly sampled tuples of source and sink positions. For the connected diagram $h$2, the sampled coordinates $h$3 are summed with a normalization factor $h$4 correcting for the number of stochastic samples. The paper emphasizes that the sample density can be tuned in different subdomains, for example to impose distance cuts between internal vertices or to oversample regions that contribute more strongly. For the disconnected diagram $h$5, the method uses a stochastic source position $h$6, constructs a spatially summed two-point correlator $h$7, and subtracts the gauge average to form $h$8; the SCS estimator then combines these subtracted quantities with a normalization $h$9. The paper notes that even in QED$$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$0, subtracting the vacuum expectation value of each summand in the SCS estimator significantly reduces the statistical error, although the exact subtraction is theoretically zero there.

The same coordinate-sampling philosophy is applied to tadpole-type disconnected diagrams. Tadpole fields

$$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$1

are split into low and high modes; the low-mode part is reconstructed exactly on over 2000–5000 low Dirac modes using multi-grid Lanczos, while the high-mode difference is estimated with a random sparse $$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$2 grid. The paper states that a $$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$3 sparse grid is often sufficient, corresponding to $$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$4 propagator computations. The whole SCS Wick-contraction machinery is validated by computing the full all-to-all propagator on a single $$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$5 gauge configuration and checking convergence of the stochastic GPT implementation as the number of stochastic propagators increases.

Finite-volume QED systematics are handled by comparing three photon propagator prescriptions: QED$$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$6, QED$$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$7, and QED$$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$8. The paper reports that the internal-vertex-distance dependence of diagram $$\|\mathbb E_\xi[g_h(\theta,\xi)]-\nabla f(\theta)\|_2\le ch.$$9 differs between QED$c=0$0 and QED$c=0$1, showing that the choice of QED prescription affects the long-distance structure and hence the systematic error budget. In this setting, SCS is not a descent method but the computational engine that turns formally exact all-to-all coordinate-space expressions for $c=0$2 and strong-isospin-breaking corrections into estimators based on sampled coordinates.