Richardson Debiasing in SGD

• Richardson debiasing is a method that uses Richardson extrapolation to cancel the leading-order gradient bias caused by imputation in stochastic gradient descent.
• It constructs a further-thinned mask to generate an extra gradient estimate, enabling cancellation of the bias from O(‖p‖) to O(‖p‖²) or higher.
• The approach applies broadly across models and missingness mechanisms, demonstrating improved test loss and parameter estimation in empirical studies.

Searching arXiv for the cited paper and closely related missing-data SGD debiasing work. Richardson debiasing, in the setting of learning with incomplete covariates, is a procedure that uses Richardson extrapolation to remove leading-order gradient bias introduced by imputation-based stochastic optimization. For supervised learning with differentiable loss $\ell(w;x,y)$, the central observation is that the population bias of an imputed stochastic gradient admits an exact first-order dependence on the missingness ratio vector $p$, with leading term $O(\|p\|)$. The method exploits that structure by constructing, from an already incomplete observation, a deliberately more incomplete version and combining the two corresponding gradients so that the linear bias term cancels. Under independent hMCAR or independent sMAR, one Richardson step reduces the gradient bias from $O(\|p\|)$ to $O(\|p\|^2)$, and under independent masking the construction extends to higher-order cancellation because the bias is a multilinear polynomial in $p$ (Genans et al., 19 May 2026).

1. Formal setting and missingness model

The framework considers covariates $X \in \mathcal X \subseteq \mathbb R^d$, response $Y \in \mathcal Y$, parameter $w \in \mathcal W \subseteq \mathbb R^q$, and differentiable loss $\ell(w;x,y)$. The population risk is

$p$0

with gradient

$p$1

SGD updates take the form

$p$2

where $p$3 is a stochastic gradient estimator.

Incomplete covariates are encoded by a mask $p$4, where $p$5 denotes that coordinate $p$6 is missing. Missing entries are filled by an imputation rule $p$7, applied per observation, producing an imputed covariate vector $p$8, where $p$9 denotes imputation randomness and satisfies $O(\|p\|)$0. The learner then uses

$O(\|p\|)$1

as a gradient estimator.

The missingness ratio vector is

$O(\|p\|)$2

Two classes of missingness mechanisms are treated. Under heterogeneous MCAR (hMCAR), $O(\|p\|)$3 and each feature has its own missingness probability $O(\|p\|)$4. Under scalable MAR (sMAR), there is a set $O(\|p\|)$5 of always-observed variables, denoted $O(\|p\|)$6, such that $O(\|p\|)$7 and

$O(\|p\|)$8

for known intensity functions $O(\|p\|)$9 satisfying $O(\|p\|)$0. Writing

$O(\|p\|)$1

yields the unified form $O(\|p\|)$2.

These definitions place Richardson debiasing within standard stochastic optimization, but with the missingness mechanism itself treated as a structured perturbation parameterized by $O(\|p\|)$3. That perspective is essential because the debiasing step acts on the missingness scale rather than directly on the loss or the model.

2. First-order bias structure

The population gradient bias induced by imputation-based SGD is

$O(\|p\|)$4

where the expectation is over $O(\|p\|)$5, masks $O(\|p\|)$6, and imputation randomness $O(\|p\|)$7. The core structural result is that, for any data-independent imputation rule $O(\|p\|)$8, under hMCAR or sMAR,

$O(\|p\|)$9

with $O(\|p\|^2)$0 independent of $O(\|p\|^2)$1. The $O(\|p\|^2)$2-th column of $O(\|p\|^2)$3 is

$O(\|p\|^2)$4

where $O(\|p\|^2)$5 is the gradient obtained after declaring only coordinate $O(\|p\|^2)$6 missing and imputing it (Genans et al., 19 May 2026).

This identifies the leading bias as the sum of single-coordinate missingness effects. The remainder term $O(\|p\|^2)$7 collects co-missingness effects, namely cases in which at least two coordinates are missing together. Under conditional independence of the mask components given $O(\|p\|^2)$8,

$O(\|p\|^2)$9

Accordingly, the leading bias is linear in $p$0.

The result is described as universal in the sense used in the source: it holds for any parametric model, any differentiable loss, and any data-independent imputation rule. Better imputation can reduce the magnitude of $p$1 but cannot remove the $p$2 scaling unless it perfectly recovers the missing features. A plausible implication is that improved imputers and Richardson debiasing address different components of the error: the former may shrink the coefficient of the linear term, whereas the latter is designed to cancel the term itself.

Under independent masking, the bias has a more explicit multilinear polynomial form,

$p$3

with coefficients

$p$4

The degree-$p$5 terms are the leading linear effects; terms with $p$6 encode higher-order interactions among missing coordinates. The abstract states that, when missing indicators are independent, the population gradient bias is a multilinear polynomial in $p$7 and depends only on population gradient errors induced by declaring a single coordinate missing (Genans et al., 19 May 2026).

3. Richardson extrapolation by controlled thinning

Classical Richardson extrapolation starts from an expansion

$p$8

If both $p$9 and $X \in \mathcal X \subseteq \mathbb R^d$0 are available for some $X \in \mathcal X \subseteq \mathbb R^d$1, then

$X \in \mathcal X \subseteq \mathbb R^d$2

cancels the linear term and leaves bias of order $X \in \mathcal X \subseteq \mathbb R^d$3. In the present setting, the role of $X \in \mathcal X \subseteq \mathbb R^d$4 is played by the population gradient $X \in \mathcal X \subseteq \mathbb R^d$5, with $X \in \mathcal X \subseteq \mathbb R^d$6.

The practical difficulty is that, for any given observation, one ordinarily sees only one mask realization at scale $X \in \mathcal X \subseteq \mathbb R^d$7. The method resolves this by deliberately adding missingness. Given an observed mask $X \in \mathcal X \subseteq \mathbb R^d$8 at scale $X \in \mathcal X \subseteq \mathbb R^d$9, choose $Y \in \mathcal Y$0 such that $Y \in \mathcal Y$1 for all $Y \in \mathcal Y$2. For each coordinate subject to missingness, draw

$Y \in \mathcal Y$3

conditionally independently, and define

$Y \in \mathcal Y$4

If a coordinate is already missing under $Y \in \mathcal Y$5, it remains missing. If it is observed, it is further hidden with a controlled probability. The resulting mask satisfies

$Y \in \mathcal Y$6

and the components remain conditionally independent.

From the two masks $Y \in \mathcal Y$7 and $Y \in \mathcal Y$8, the construction imputes once at the more-thinned level: $Y \in \mathcal Y$9 The imputed vector at the original scale is then recovered by restoring the entries that were only hidden by the second thinning: $w \in \mathcal W \subseteq \mathbb R^q$0 This yields

$w \in \mathcal W \subseteq \mathbb R^q$1

and the Richardson-corrected gradient

$w \in \mathcal W \subseteq \mathbb R^q$2

A technically important feature is that genuinely missing entries shared by both masks use the same imputed values in $w \in \mathcal W \subseteq \mathbb R^q$3 and $w \in \mathcal W \subseteq \mathbb R^q$4. The source states that this sharing is necessary for first-order cancellation of the bias operator $w \in \mathcal W \subseteq \mathbb R^q$5, and that re-imputing independently at the two levels would break the effect.

4. One-step and multi-step cancellation

Under independent hMCAR or independent sMAR, one-step Richardson debiasing satisfies

$w \in \mathcal W \subseteq \mathbb R^q$6

The mechanism is direct. Since

$w \in \mathcal W \subseteq \mathbb R^q$7

the bias of the Richardson gradient is

$w \in \mathcal W \subseteq \mathbb R^q$8

because the linear terms satisfy $w \in \mathcal W \subseteq \mathbb R^q$9. Under conditional independence, both remainder terms are $\ell(w;x,y)$0, so the corrected bias is also $\ell(w;x,y)$1 (Genans et al., 19 May 2026).

The same logic extends to higher order because, under independent hMCAR or independent sMAR, $\ell(w;x,y)$2 is a multilinear polynomial in $\ell(w;x,y)$3. Grouping by degree,

$\ell(w;x,y)$4

Choose scales

$\ell(w;x,y)$5

construct iteratively thinned gradients $\ell(w;x,y)$6, and define

$\ell(w;x,y)$7

where the coefficients solve the Vandermonde system

$\ell(w;x,y)$8

Then

$\ell(w;x,y)$9

If $p$00, then the $p$01-th order Richardson estimator cancels the bias exactly: $p$02 The reason is that the polynomial contains no monomials of degree larger than $p$03, so annihilating degrees $p$04 through $p$05 removes the full bias.

A special case emphasized in the source is linear regression with squared loss, for which the bias polynomial degree is at most $p$06. First-order Richardson therefore reduces bias from $p$07 to $p$08, while second-order Richardson cancels the bias exactly. This sharp finite-degree behavior clarifies why higher-order Richardson is especially attractive when the number of missing coordinates is small or the induced bias polynomial is low degree.

5. Algorithmic realization and empirical behavior

At iteration $p$09 of Richardson-SGD, each sampled incomplete observation is processed by using the original mask $p$10, generating a further-thinned mask $p$11, imputing once on the more-thinned sample, restoring artificially hidden entries to recover $p$12, computing $p$13 and $p$14, forming

$p$15

and updating

$p$16

For minibatches, the corrected gradients are averaged before the parameter update. The method requires one more mask via thinning, one imputation at level $p$17, and one extra gradient evaluation per sample.

The construction is model-agnostic in the sense stated in the source: it applies whenever the loss is differentiable, the per-sample gradient $p$18 can be computed from imputed covariates, and the imputation rule is data-independent with $p$19. The source explicitly lists compatibility with linear and generalized linear models, logistic regression, Poisson regression, and any parametric model where gradients can be computed on imputed covariates. It also lists zero imputation, mean imputation, MICE, Random-Forest MICE, k-NN, and neural imputers, provided they are used per-sample (Genans et al., 19 May 2026).

The empirical study covers generalized linear models—linear (Gaussian), logistic, and Poisson—on synthetic Gaussian covariates with various correlation structures and on real datasets including Diabetes, California housing, Covertype, Breast Cancer, and Bike Sharing. Missingness mechanisms include hMCAR, heterogeneous MCAR, and sMAR, with average missingness levels typically $p$20 and additional sweeps over $p$21. Across models and datasets, the reported pattern is that plain imputation plus SGD suffers from bias due to missing data, while Richardson-SGD always improves or matches the corresponding imputation baseline, often substantially. The source specifically reports improved test loss and parameter MSE for logistic regression on Covertype under hMCAR across a wide range of missingness levels, significant reductions in final parameter MSE for multiple logistic-regression datasets, consistent gains in linear and Poisson regression under MCAR, heterogeneous MCAR, and sMAR, and synthetic examples in which second-order Richardson in linear regression essentially overlaps with the complete-data trajectory while third-order Richardson gives additional improvement for logistic regression.

The source also gives asymptotic optimization consequences. For plain imputation-based SGD,

$p$22

and one-pass SGD satisfies

$p$23

Using Richardson-SGD and a decreasing stepsize $p$24 under standard smoothness and strong convexity assumptions yields

$p$25

so the missingness-induced error floor improves by two powers of $p$26.

Practical guidance in the source favors modest values such as $p$27, noting that larger $p$28 reduce variance inflation because $p$29 decreases with $p$30, but are constrained by $p$31. First-order Richardson is recommended as the default compromise between bias reduction and variance, while higher-order versions are suggested when missingness is small or localized, or when the bias polynomial degree is small. When $p$32 and, in sMAR, $p$33 are estimated rather than known, the resulting bias is

$p$34

which remains dominated by $p$35 if the estimates are sufficiently accurate.

6. Relation to prior corrections, limitations, and interpretive issues

Richardson debiasing differs from plain imputation-based SGD by acting directly on the missingness scale. Most missing-data methods treat the mask as fixed and attempt to reconstruct missing entries more accurately or to model the joint distribution of $p$36. Here the central operation is different: the algorithm deliberately changes the missingness level by further thinning, then uses extrapolation to cancel bias terms. The source identifies this as conceptually novel in the missing-data context (Genans et al., 19 May 2026).

The method also differs from earlier debiasing procedures for linear regression. The source states that earlier work, exemplified by Sportisse et al. 2020, derives closed-form bias corrections for linear regression under MCAR with zero imputation and independent masks, and that those corrections are model-specific and rely heavily on linearity. Richardson-SGD is presented as a generalization: it works for any parametric loss and any data-independent imputation, extends beyond MCAR to sMAR, and in linear regression reduces to the previously known exact-correction regime because the bias polynomial has degree at most $p$37.

Several limitations are explicit. The clean $p$38 result and the multilinear polynomial structure rely on conditional independence of the mask components given $p$39. With strong dependencies, including anticorrelation, co-missingness probabilities can be first-order in $p$40, so the remainder $p$41 may also be $p$42 and the effectiveness of Richardson may be reduced. The analysis also assumes data-independent imputation, although the source notes that iterative imputers such as MICE can be adapted by training them on an auxiliary dataset and then treating them as fixed per-sample imputers. The stated convergence theory is for one-pass SGD; multi-epoch behavior, with correlated gradient noise across epochs, is not fully analyzed. In sMAR, the construction requires known intensity functions $p$43, and estimation error or misspecification can degrade performance, even though empirical results indicate that the method often remains beneficial.

A common misconception would be to interpret the procedure as arbitrary corruption of already incomplete data. The source supports a narrower and more technical interpretation. The added missingness is controlled so that the further-thinned mask has the same distribution as a mask drawn directly at scale $p$44, and the two gradients are coupled through shared imputations on the genuinely missing entries. The debiasing effect therefore comes not from noise injection per se, but from a structured cancellation argument tied to the exact dependence of $p$45 on $p$46. This suggests a broader methodological principle: when a nuisance parameter induces a tractable expansion of the estimation bias, deliberately changing that parameter and extrapolating across scales may produce a higher-order estimator.