Sampling rules

FiniteSum(n, M, Fi, Gi!, Hi)

f(x) = (1/M) Σ_{i=1}^{M} F_i(x), sampled by subsetting the index set.

• Fi(i, x) -> Real
• Gi!(i, x, g) -> nothing, accumulating ∇F_i(x) into g
• Hi(i, x) -> AbstractMatrix

The empirical-risk case. Subsampling the same index set at x_k and x_k + s_k gives common random numbers for free, and it is worth having: the noise in f̂(x_k) − f̂(x_k + s_k) then cancels to the extent that the F_i are smooth, which is what keeps ρ̂ usable at small radii.

PerturbedSum(base, M; σg = 1.0, σH = 0.0, seed = 0)

A finite sum whose mean is exactly a known model: with mean-zero b_i and symmetric mean-zero C_i,

\[F_i(x) = f_{\text{base}}(x) + b_i^\top x + \tfrac12 x^\top C_i x , \qquad \tfrac1M \textstyle\sum_i F_i \equiv f_{\text{base}} .\]

The perturbations are smooth in x, so common random numbers behave as they do on a real empirical risk rather than cancelling artificially, while the exact f, ∇f and ∇²f remain available through base.

That combination is what makes a controlled experiment possible: the estimation error at every iterate is known, not merely bounded.

Contrast PerturbedExpectation, the same construction without the finite population — there the perturbations are drawn fresh and their mean is zero only in expectation, so no sample size makes the estimate exact.

PerturbedExpectation(base; σg = 1.0, σH = 0.0)

f(x) = E_ξ[ f_base(x) + b(ξ)ᵀx + ½ xᵀC(ξ)x ] with b ~ N(0, σg²I) and C symmetric mean-zero: an expectation whose mean is exactly base.

The expectation twin of PerturbedSum, and the difference between them is the point of having both:

• PerturbedSum centres its M perturbations, so at N = M the estimate is exact and the accuracy hypotheses of the stochastic theory are discharged.
• PerturbedExpectation draws fresh, so the sample mean of b is O(σ_g/√N) at every N. There is no sample size at which the estimate is exact, and no full-batch iteration to fall back on.

has_truth is true here — base is available — which makes it the right expectation problem for a controlled experiment. A general expectation has no such handle, and TRParams(true_stop = true) is rejected against it.

GaussianDraw

One expectation batch: N realisations (b_i, C_i) drawn i.i.d., held as an n × N matrix and an optional vector of matrices.

Unlike a finite-sum batch this is data, not indices — which is the whole difference between the two classes at the evaluation level, and the reason batch_obj and friends take an opaque batch argument rather than a Vector{Int}.

batch_obj(prob, x, batch) -> Float64

Estimate f(x) from one batch.

Must be a deterministic function of (x, batch): the same batch is reused for f̂(x_k) and f̂(x_k + s_k) within an iteration, and that is what makes the noise in the difference shrink with ‖s_k‖. An implementation that draws randomness inside this function defeats it.

batch_grad!(prob, x, batch, g) -> nothing

Estimate ∇f(x) from one batch, writing it into g. Deterministic in (x, batch), for the reason given under batch_obj.

batch_hess(prob, x, batch) -> Matrix

Estimate ∇²f(x) from one batch. Cached per (iterate, batch) by the oracles, since it is a pure function of both and on the finite-difference problems it costs O(n) gradient evaluations.

grad_variance(prob, x, batch) -> Float64

tr Var(∇F_i) over the batch — the σ_g² the NormTest and RadiusProportional consume.

obj_variance(prob, x, batch) -> Float64

Var(F_i) over the batch — the σ_f² that SequentialEstimation and the objective half of RadiusProportional consume.

SamplingRule

How N_k^grad and N_k^obj are chosen. A concrete subtype implements

grad_sample_size(rule, st::SamplingState) -> Int
obj_sample_size(rule, st::SamplingState)  -> Int

and may implement couples_to_radius, needs_scores, user_cap, requires_finite_population and reset_sampling_rule!.

Implementations must be pure: the solver calls both queries in the same iteration, so a rule that mutates itself in one is invoked twice per step. SequentialEstimation caches on st.k for exactly this reason.

FullBatch()

N_k = M at every iteration: the whole finite sum.

A finite-sum run under this rule is exactly the deterministic run — ρ̂ = ρ, the stopping test is the real one, no accuracy hypothesis is needed — which makes it the reference every other sampling rule should be scored against, and a sharp regression test: FiniteSumTRSolver with FullBatch() must reproduce DeterministicTRSolver on FullBatchNLP iterate for iterate.

Rejected on an ExpectationProblem: there is no full batch. That refusal is the type system carrying the substantive difference between the two classes, not a technicality — an expectation's accuracy hypotheses can never be discharged by sampling harder, and this rule is what discharging them looks like.

SamplingState

What a sampling rule may condition on: the iteration index k, the radius Δ, the gradient-estimate norm g_norm, the running variances σg², σf², their parallel/orthogonal split ip², orth², the previous predicted reduction pred, the previous sample size N_prev, and the population N_pop.

N_pop is M for a finite sum and typemax(Int) for an expectation. Only FullBatch reads it, and it raises on an expectation — where "all of it" does not exist.

SampleStats(σg², σf², ip², orth²)

The four second moments the adaptive rules consume:

• σg² = tr Var(∇Fᵢ) — the NormTest;
• σf² = Var(Fᵢ) — SequentialEstimation and the objective half of RadiusProportional;
• ip² = Var(∇Fᵢᵀĝ) — the variance along the estimated gradient;
• orth² = E‖∇Fᵢ − (∇Fᵢᵀĝ/‖ĝ‖²)ĝ‖² — the variance orthogonal to it.

Splitting σg² = ip²/‖ĝ‖² + orth² is the whole content of the inner-product tests: only the component along ĝ decides whether the estimated direction is still a descent direction, so bounding the sum is stricter than descent requires.

batch_stats(prob, x, batch, g) -> SampleStats

All four moments in one pass. Generic for any problem exposing per-observation scores; falls back to the gradient and objective variances alone otherwise, in which case the inner-product family is rejected at oracle construction rather than at the first iteration.

batch_stats(p::ScoredProblem, x, batch, g) -> SampleStats

The parallel/orthogonal split, exact rather than estimated separately: with S the n × N score matrix and ĝ its column mean,

ip²   = Var_i(Sᵢᵀ ĝ),
orth² = (1/(N−1)) Σᵢ ‖Sᵢ − ĝ‖² − ip²/‖ĝ‖² ,

so σ_g² = ip²/‖ĝ‖² + orth² holds to rounding and the tests are comparable on the same batch. σf² comes back on the same pass, so the caller need not traverse twice.

FixedSample(N; N_obj = N)

Sample-average approximation: the same N at every iteration.

The baseline every adaptive rule should be scored against, and the only rule for which total work is proportional to the iteration count — which is what iteration-count comparisons quietly assume.

RadiusProportional(; κ_g = 1.0, κ_f = 1.0, N_min = 2, N_max = nothing,
                     σ_floor = 1e-12)

The STORM requirement, made operational:

\[N_k^{\mathrm{grad}} = \lceil (\sigma_g/(\kappa_g \Delta_k))^2 \rceil, \qquad N_k^{\mathrm{obj}} = \lceil (\sigma_f/(\kappa_f \Delta_k^2))^2 \rceil ,\]

so the gradient estimate is accurate to O(Δ_k) and the objective estimates to O(Δ_k²) with fixed probability.

This is the rule that makes the mechanism pay for its own radius: N_k grows like Δ_k^{-2}, so a criticality-anchored mechanism drives the sample size up without bound near a solution while a criticality-blind one does not.

N_max is an expectation-only budget. On a finite sum the cap is M and the solver's budget keyword is where a sub-population limit belongs.

NormTest(; θ = 0.5, N_min = 2, N_max = nothing, N_obj = nothing)

Byrd, Chin, Nocedal & Wu (2012): σ_g²/N_k ≤ θ²‖ĝ_k‖².

Never consults Δ_k, so the sampling rule and the radius mechanism stay independent — which makes it the right control for isolating what the coupling in RadiusProportional actually does. It still drives N_k → ∞, but through ‖ĝ_k‖ → 0 rather than through the radius, so the growth is a property of the problem rather than of the mechanism.

InnerProductTest(; θ = 0.9, N_min = 2, N_max = nothing, N_obj = nothing)

Bollapragada, Byrd & Nocedal (2018): control the variance of the sampled gradient along the estimated direction,

\[\frac{1}{N_k}\operatorname{Var}(\nabla F_i^\top \hat g) \le \theta^2 \|\hat g\|^4 .\]

This is the condition that keeps −ĝ a descent direction with fixed probability. It is weaker than the NormTest, which also pays for error orthogonal to ĝ — error that rotates the direction without threatening descent — so it reaches the same guarantee at a materially smaller sample.

θ is not on the same scale as the norm test's: values near 1 are usual here, against 0.1–0.5 there. Requires a ScoredProblem.

OrthogonalityTest(; ν = 2.0, N_min = 2, N_max = nothing, N_obj = nothing)

The companion test: bound the variance orthogonal to the estimated direction, E‖proj⊥∇Fᵢ‖²/N_k ≤ ν²‖ĝ‖².

The inner-product test keeps the sampled direction pointing downhill; this one keeps it from rotating too far while doing so. Neither implies the other, and the norm test implies both at a larger sample — the sense in which it is conservative rather than wrong. Requires a ScoredProblem.

AugmentedInnerProduct(; θ = 0.9, ν = 2.0, N_min = 2, N_max = nothing, N_obj = nothing)

Both conditions at once — the form Bollapragada, Byrd & Nocedal actually recommend. Taking the maximum controls the sign and the rotation of the sampled direction simultaneously while still costing less than bounding the total variance. Requires a ScoredProblem.

GeometricSample(; N₀ = 8, rate = 1.1, N_max = nothing)

N_k = ⌈N₀·rate^k⌉: a schedule fixed in advance, independent of everything the run observes.

The other control. It reaches any accuracy eventually and cannot respond to the radius at all, so comparing it against RadiusProportional at matched total sample budget separates "spending more" from "spending it at the right iterations".

rate^k overflows to Inf around k ≈ 7000 for rate = 1.1, inside the default iteration budget, which is why the conversion goes through _sample_size.

SequentialEstimation(; κ = 0.25, α = 0.05, N_min = 8, N_max = nothing,
                       monotone = true, growth = 2.0, N_start = 8)

Bastin, Cirillo & Toint (2006): make the sample large enough that the noise in the estimated decrease is small beside the decrease the model predicts,

\[N_k \ge \frac{2 z_{\alpha/2}^2 \hat\sigma_f^2}{\kappa^2 \mathrm{pred}_{k-1}^2}.\]

Two things distinguish it from the gradient-based tests.

It tracks progress, not criticality. The accuracy demanded falls as progress becomes harder to make, rather than as the gradient becomes small.

It is monotone by default. A batch that shrinks makes f̂ jump for reasons unrelated to the step, and ρ̂ then measures the change of estimator rather than of objective. growth caps how fast it may rise in one iteration.

CertifiedDecrease(; p = 0.9, N_min = 2, N_max = nothing, N_start = 8,
                    monotone = false, growth = 4.0, variance = :empirical)

Sample until the achieved decrease is certified against the noise in the estimate of that decrease, using paired differences under common random numbers.

With the same realisations at both ends of the step,

\[D_i = F(x_k, \xi_i) - F(x_k + s_k, \xi_i), \qquad N_{k+1} = \left\lceil z_p^2 \hat\sigma_N^2 / \hat\delta_N^2 \right\rceil ,\]

where δ̂_N and σ̂_N² are the sample mean and variance of the D_i.

Why pairing is the whole point

Holding ξ_i fixed at both points makes D_i → 0 pathwise as s_k → 0, so σ_D² = O(‖s_k‖²). Two independent batches give O(1) instead, and the rule is then unusable near a solution. On a heavy-tailed sample problem the measured variance ratio between the paired and unpaired estimators runs to several orders of magnitude.

What the rule costs, and why it does not see the radius

Both the numerator and the denominator are linear in the step: δ̂ ∼ ‖g_k‖‖s_k‖ and σ̂ ∼ C‖s_k‖, so

\[N_k \sim \frac{z_p^2 C^2 \|s_k\|^2}{\|g_k\|^2 \|s_k\|^2} = \frac{z_p^2 C^2}{\|g_k\|^2} ,\]

and the step length cancels. Pairing buys a large constant, not a better rate, and the per-iteration sampling cost is insensitive to the radius rule: mechanisms differ in how many iterations they take, not in what each costs to certify. This is the opposite of RadiusProportional, whose N_k ∼ Δ_k^{-2} makes the mechanism pay for its own radius, and the two together bracket the question of whether a sampling rule can see the radius at all.

couples_to_radius is true because the statistic is read off the step; the cancellation above is a statement about the asymptotic rate, not about the rule ignoring s_k.

monotone = false by default, unlike SequentialEstimation: the certified size is meaningful on its own at every iteration, and forcing it upward hides the ‖g_k‖^{-2} growth this rule exists to exhibit.

monotone = true is the special case a = 1 of SmoothedSize, which is where a two-sided bound a N_k ≤ N_{k+1} ≤ b N_k with a < 1 belongs; wrap this rule rather than adding a floor to it.

Which variance

variance = :empirical (default) sizes the batch from the sample variance of the D_i. variance = :outer_product uses s_kᵀB̃_ks_k − (g̃_kᵀs_k)² instead, the first-order Taylor approximation formed from the score matrix; see paired_variance_kind for what separates them and paired_op_variance for the computation. The outer-product form needs per-observation scores, so needs_scores becomes true and the rule is then rejected on a problem that has none.

Both estimates are recorded on the trace whenever the problem supplies scores, under :paired_variance_trajectory and :paired_op_variance_trajectory, so the choice can be audited against the one not taken.

SmoothedSize(rule; a = 0.75, b = 2.0, N_floor = 1, N_incr = 0, N_start = 8)

Wrap a SamplingRule and bound how fast the size it asks for may move.

The inner rule proposes a candidate N⁺_{k+1}; this reports

\[N_{k+1} = \min\bigl\{ b N_k , \; \max\{ a N_k , \; N^{\min}_k , \; N^{+}_{k+1} \} \bigr\} , \qquad N^{\min}_{k+1} = N^{\min}_k + N_{incr} ,\]

with 0 < a ≤ 1 and b > 1. The two-sided bound is the naive smoothing of the PreDoc, written there as Naive(a, b); N_floor and N_incr are the growing floor N^{min}_k of its sample-size update step.

Why it exists

CertifiedDecrease already offers monotone = true, which is a = 1: the size may never fall. That is a strong assumption, and it is not the one the PreDoc makes. It also hides what the certified formula does, since a size that can only rise tells you nothing about the iterations at which the formula asked for less.

monotone = false is the other extreme, a = 0 and b = ∞, and there the raw requirement swings over orders of magnitude between consecutive iterations. a and b interpolate, and the interesting settings are strictly between the two.

What it does not do

It does not change the inner rule's statistic, only the size read off it. A run under SmoothedSize and a run under the bare rule see the same estimates of the same quantities; they differ in what they do with them. In particular the smoothing is not a variance reduction and should never be described as one.

When the floor and the ceiling conflict

N_incr > 0 makes the floor grow linearly while b N_k grows geometrically from whatever N_k happens to be, so the floor can exceed the ceiling — most easily right after a sharp fall. The floor wins, and N_{k+1} = N^{min}_k. It is a requirement the caller stated outright, whereas the ceiling is a smoothing preference; and the alternative would let a run sit below its own stated floor for as long as it took b to catch up.

Composition

Every trait is forwarded, so the wrapper is transparent to check_rule_problem, to the oracle's cap logic and to the solver: needs_paired, needs_scores, couples_to_radius, requires_finite_population, user_cap and paired_variance_kind all read through to rule, and record_paired! is passed straight down.

The gradient and objective channels are smoothed separately, each against its own previous size, so a rule that distinguishes them — RadiusProportional asks for Δ^{-2} and Δ^{-4} — keeps the distinction. A rule that returns one size for both, which is every rule that pairs, is unaffected: the two channels then carry identical state.

Example

The PreDoc's setting, a certified size smoothed by Naive(0.75, 2) with a floor starting at 10 and rising by 10 each iteration:

rule = SmoothedSize(CertifiedDecrease(p = 0.95, N_min = 10);
                    a = 0.75, b = 2.0, N_floor = 10, N_incr = 10, N_start = 100)

record_prediction!(m::SampledNLP, pred) -> pred

Store the model's predicted reduction. SequentialEstimation sizes the next batch against it; no other rule reads it.

reset_sampling_rule!(rule) -> nothing

Restore mutable state. Default: no-op.

couples_to_radius(rule) -> Bool

Whether N_k depends on Δ_k. When it does, the radius mechanism and the sampling rule stop being independent axes: for a criticality-anchored rule, ĝ_k sets Δ_k sets N_k sets the accuracy of ĝ_{k+1}. RDelta has no such loop.

needs_scores(rule) -> Bool

Whether the rule requires per-observation scores. Checked once, at oracle construction, against has_scores on the problem.

requires_finite_population(rule) -> Bool

Whether the rule is meaningful only on a FiniteSumProblem.

true only for FullBatch: "all of it" is not a sample size an expectation has. Checked at oracle construction, so pairing it with an expectation is a constructor error rather than an N_k = typemax(Int) allocation.

sample_cap(rule) -> Int

The rule's own cap resolved for arithmetic: the user's N_max if they set one, typemax(Int) otherwise. The solver combines it with its own budget and with population(prob).

See user_cap, which preserves the distinction between "unset" and "set to a large number" — the one check_population_cap needs.

grad_sample_size(rule, st) -> Int

Terms to draw for the gradient estimate. The quantity the stochastic experiments vary. st carries Δ_k, so a rule may couple the sample size to the radius — RadiusProportional does, which is what makes the sampling cost a property of the radius mechanism rather than of the problem alone.

obj_sample_size(rule, st) -> Int

Terms to draw for the objective estimates. Usually larger: the accuracy requirement on f is O(Δ_k²) against O(Δ_k) on the gradient, so the sample scales as Δ_k^{-4} rather than Δ_k^{-2}. On a run whose radius becomes small it is the objective estimates, not the gradients, that dominate the bill.

resample!(m::SampledNLP, k, Δ, g_norm) -> (Ng, Nf)

Choose N_k from the sampling rule and draw the batches for iteration k.

Called once per iteration by the solver, before anything is evaluated. The variance estimates fed to the rule are those of the previous batch, since the new one does not exist yet — the usual plug-in, and the reason N_min matters: a rule cannot estimate a variance from one sample.

The rule states a requirement; the oracle applies the cap. For an expectation that is the budget, for a finite sum min(budget, M).

update_variances!(m::SampledNLP, x, g) -> (σg², σf²)

Refresh the variance estimates from the current batch, at the current iterate.

Takes the batch gradient g the solver has already computed instead of recomputing it: this previously evaluated the score matrix twice — once through batch_grad! and once inside batch_stats, which then ignored the g it was handed — on top of the solver's own gradient call.

samples_used(m) -> (grad = …, obj = …, confirm = …, total = …)

Cumulative term evaluations — the cost measure for a sampled comparison. Iteration counts are proportional to work only under FixedSample.

confirm is what confirm_gradient_norm! spent verifying candidate stops, and is zero unless TRParams(true_stop = …) asked for verification. It is separated out because it is not optimisation work — comparing a mechanism that confirms against one that does not, on total, would charge the first for being honest — but it is included in total, because it is work that was done.

reset_sampling!(m::SampledNLP) -> m

Restore the counters, the histories, the batches and the random stream to their initial state, so the same oracle can be reused across configurations and each run sees the same realisations — which is what makes a comparison of mechanisms under noise a comparison of mechanisms rather than of random seeds.

draw_batch(prob, rng, N; replace = true)

Produce a batch of size N.

For a FiniteSumProblem this subsets 1:M and replace decides whether with or without replacement. For an ExpectationProblem it draws N fresh realisations and replace is meaningless (an unbounded population is always sampled with replacement in effect).

sample_schemes() -> NTuple{3, Symbol}

The three ways a batch may be carried from one iteration to the next.

Within an iteration the batch is fixed and every evaluation shares it; that is the common random numbers the whole design rests on and no scheme touches it. What the scheme decides is what happens between iterations, and the three answers are the ones the PreDoc states.

:independent

𝒩_k and 𝒩_{k+1} are independent. The estimate f̃_k then moves between iterations for two reasons at once — the iterate moved, and the sample changed — and early on, when N_k is small and the variance is large, the second dominates.

:nested

𝒩_k ⊆ 𝒩_{k+1} on a growth: the previous batch is kept and N_{k+1} − N_k fresh realisations are appended. On a shrink the survivors are drawn uniformly from the previous batch, so 𝒩_{k+1} ⊂ 𝒩_k and nesting holds in that direction too. Consecutive estimates are then positively correlated and f̃ moves mostly because the iterate moved.

Growth without replacement draws from the complement of the current batch, which costs O(M) per growth; with replacement it appends rand(1:M, ·) and duplicates are allowed, as they are in the initial draw.

:prefix

𝒩_k = {1, …, N_k}: literally the first N_k individuals, in population order. Every batch is a prefix of every larger one, so this is the most correlated scheme there is, and it is the only one that needs no random numbers at all after construction.

It requires a finite population. On an ExpectationProblem there is no first individual to take, and the constructor rejects it.

sample_scheme(m::SampledNLP) -> Symbol

Which of sample_schemes this oracle is carrying its batch under.

true_objective(prob, x)

The exact f(x), by evaluating every term rather than a batch.

Defined for every FiniteSumProblem — one pass over all M terms — and for those expectations whose construction supplies the mean. Guarded by has_truth.

true_gradient(prob, x)

The exact ∇f(x).

The counterpart of true_objective, and the one that matters for scoring: a stopping test on ‖ĝ_k‖ is a statement about one batch, and a mechanism that shrinks the radius fast enough will meet it on noise alone. Every reported number should come from this, and TRParams(true_stop = true) makes the solver's own stopping test do so.

Costs a full pass, so it belongs in the reporting and in trace, not in the inner loop. A run on a problem without truth cannot use true_stop and gets no :true_grad_trajectory.

needs_paired(rule) -> Bool

Whether rule consumes the paired decrease statistic, and so whether the solver must evaluate the trial point on the same batch and difference it observation by observation. false by default, so no rule pays for a statistic it does not read.

record_paired!(rule, δ, σ², N) -> nothing

Hand the rule the paired decrease statistic of the iteration just taken: the sample mean δ and sample variance σ² of

D_i = F(x_k, ξ_i) − F(x_k + s_k, ξ_i)

over the N observations of the current batch. A no-op for every rule that does not declare needs_paired.

The statistic belongs to the rule rather than to the oracle because the rule is its only consumer, and because it is state carried across iterations in exactly the way SequentialEstimation carries the predicted reduction.

paired_decrease_stats(m::SampledNLP, x, x_cand) -> (δ, σ², N)

The sample mean and variance of the paired differences

D_i = F(x, ξ_i) − F(x_cand, ξ_i),   i ∈ batch_f,

over the objective batch of the current iteration.

Both points are evaluated on the same realisations, which is what makes σ_D² = O(‖x_cand − x‖²); see obs_objective. Consumed by CertifiedDecrease through record_paired!.

Cost

None beyond what the iteration already pays. The solver has already evaluated f(x) and f(x_cand) as batch means over this batch; this recomputes the same two quantities at per-observation granularity. It is not charged to samples_used, because no new realisation is drawn.

Returns (NaN, NaN, N) when the batch is too small for a sample variance.

paired_variance_kind(rule) -> Symbol

Which estimate of Var(D_i) the rule wants handed to record_paired!, where D_i = F(x_k, ξ_i) − F(x_k + s_k, ξ_i).

• :empirical — the sample variance of the D_i themselves, computed by paired_decrease_stats. Assumption-free, and it costs one pass over the per-observation objective values at both ends of the step.
• :outer_product — s_kᵀ B̃_k s_k − (g̃_kᵀ s_k)², computed by paired_op_variance from the score matrix. This is the first-order Taylor approximation: D_i ≈ −∇F(x_k, ξ_i)ᵀ s_k, whose variance is s_kᵀ Var[∇F] s_k with Var[∇F] replaced by the outer-product estimate B̃_k − g̃_k g̃_kᵀ.

The two are different quantities and a rule may prefer either. The approximation is exact only in the limit of a vanishing step, and it inherits whatever error the outer product carries as an estimate of Var[∇F] — which is a different requirement from the information identity, since it does not ask B̃ to approximate the Hessian, only the score covariance.

Default :empirical. Only rules that declare needs_paired are asked.

paired_op_variance(m::SampledNLP, x, x_cand) -> Float64

The outer-product approximation to Var(D_i), where D_i = F(x, ξ_i) − F(x_cand, ξ_i), over the objective batch of the current iteration.

With s = x_cand − x, a first-order expansion of F(·, ξ) about x gives D_i ≈ −∇F(x, ξ_i)ᵀ s, so

\[\mathrm{Var}(D_i) \approx s^\top \mathrm{Var}[\nabla F(x, \xi)] s \approx s^\top \tilde B s - (\tilde g^\top s)^2 ,\]

with B̃ = (1/N) Σ ∇F_i ∇F_iᵀ the outer-product matrix and g̃ the batch gradient, both over batch_f at x. The right-hand side is the sample variance of the scalars ∇F_iᵀs, and that is how it is computed here: one pass over the score matrix, forming Var({sᵀ∇F_i}) directly rather than materialising B̃. So the result is never negative, which the difference of the two quadratic forms can be in floating point when the two nearly cancel.

What it costs, and what it assumes

One score matrix over batch_f at x: O(n·N) storage and whatever the problem charges for scores. No new realisation is drawn, so like paired_decrease_stats it is not charged to samples_used — the comparison it enables is between two estimates of one quantity on one batch.

Two approximations, and they fail in different places. The Taylor step needs ‖s‖ small relative to the curvature of F(·, ξ) along s, so it degrades on a large step, which is where the radius is large and the batch is small. The outer product is used here as an estimate of Var[∇F] alone: it does not need the information identity, and it is legitimate on a misspecified model where BHHH as a Hessian is not. Those are separate questions and conflating them is the standard error.

Returns NaN when the batch is too small for a sample variance, matching paired_decrease_stats, and when the problem supplies no scores.

supports_paired(prob) -> Bool

Whether prob implements obs_objective, and so whether the paired statistic of CertifiedDecrease is available on it.

obs_objective(prob, x, batch) -> AbstractVector{Float64}

The per-observation objective values F(x, ξ_i) on batch, whose mean is the batch estimate of f(x).

obj_variance already needs these values; exposing them is what makes a paired statistic possible. Evaluating the same batch at two points and differencing gives

D_i = F(x, ξ_i) − F(x + s, ξ_i),

whose variance is O(‖s‖²) because the realisation is held fixed, against O(1) for two independent batches. That is the whole basis of certified decrease: see CertifiedDecrease.

A problem that cannot supply per-observation values cannot support the paired rules; supports_paired reports which can.

confirm_gradient_norm!(m, x) -> Float64

‖∇f(x)‖ as accurately as this oracle can produce it, for confirming a candidate stopping test. Mutates m: the cost is added to m.samples_confirm.

Two paths, and which one applies is a property of the problem, not a choice:

• has_truth — the exact gradient, via true_gradient. Every finite sum qualifies (one pass over M), as does an expectation supplying a closed-form mean, which costs nothing per term and is charged nothing.
• otherwise — the gradient over a freshly drawn batch of population_cap terms, the largest the budget allows. This is still an estimate; it is simply the best one available, and on an expectation there is nothing better.

The finite-sum truth path draws no random numbers, so a confirmed run and an unconfirmed one see the same realisations and remain comparable iterate for iterate. The cap-batch path does consume the stream, so an expectation without truth cannot be compared that way — confirmation there changes the run it is confirming.

grad_standard_error(m) -> Float64

The plug-in standard error of the batch gradient, from the variance estimate the oracle already maintains and the batch size that produced it.

This is what makes the statistical stopping modes free — no extra evaluation, only arithmetic on numbers update_variances! has already computed. It is a plug-in estimate and inherits the usual caveat: σ_g² comes from the previous batch, and on a batch of two it is itself a one-degree-of-freedom quantity. It is used to decide whether a stop is believable, never to decide the step.

The finite population correction is not optional here

For an expectation the population is unbounded and the answer is σ_g/√N_k. For a finite sum it is not, and using σ_g/√N_k there makes the statistical test useless rather than merely conservative: it would report a positive sampling error even at N_k = M, where the batch is the population and the gradient is exact. The practical effect is that no tolerance below z·σ_g/√M could ever be certified — on a population of 20 000 with σ_g ≈ 1 that is about 1.4e-2, so every realistic tol is unreachable and :both degenerates into always paying for :full.

So: zero at N_k ≥ M, since draw_batch returns the whole population there whatever replace says; the textbook √((M−N)/(M−1)) correction below that when drawing without replacement; and no correction when drawing with replacement, where the draws are genuinely i.i.d. and σ_g/√N_k is right.