Posterior weighting

posterior_scale(sws, sq_residuals, k[, prior::IGPrior = JeffreysPrior()];
                model_dof, codimension=1) -> σ̂

The LINEAR scale σ̂ = √ŝ² of a sweep! partition, where ŝ² is fit_scale over the partition's RSS and counted dof ν = (k − model_dof) · codimension:

σ̂² = (RSS + 2β₀) / (ν + 2α₀)

Pseudo-observation form: 2α₀ extra DOF and 2β₀ extra RSS contributed by the prior; the default Jeffreys prior (α₀ = β₀ = 0) gives √(RSS/ν). Matches the closed-form ŝ² in the IG-marginal scoring paper (Sec. 4, eq. (614) of main.tex). Square the return value to report σ̂².

A thin sqrt over fit_scale, so the two cannot drift: this was the last hand-inlined copy of that formula. Its Jeffreys overload floored differently — max(RSS, eps)/max(ν, 1) floors the RSS, fit_scale floors the QUOTIENT, a floor on σ̂² consistent with every other fit_scale caller — so for a real partition (ν ≥ 1) the two agree wherever RSS ≥ ν·eps. Below that edge (for ν ≥ 2; at ν = 1 both floors land on eps) the old form gave √(max(RSS, eps)/ν) and this one gives √eps: at RSS = 0, ν = 100, 1.49e-9 → 1.49e-8; at RSS = 1e-15, ν = 10, 1.0e-8 → 1.49e-8. That is a noise-free or near-noise-free partition, and no caller in the stack reaches it. (At ν ≤ 0, k ≤ model_dof, there is no partition to scale; the result is fit_scale's there, pinned by the tests.)

fit_scale(prior::IGPrior, rss, dof) -> s2

Fit the squared scale ŝ² as the inverse-gamma MAP, (rss + 2β₀)/(dof + 2α₀).

Returns ŝ², the SQUARED scale (the PosteriorWeightState.s_hat_sq field) — never a linear σ̂, which is what posterior_scale returns.

ŝ² is PER COMPONENT: rss sums squared NORMS of dg-vector residuals, so dof counts components — `νdata = (k − modeldof)·dgfor a hard partition (the seed constructors,posteriorscale, RobustVisualGeometry'sfinalize), ordg·Σπᵢin the IRLS M-step (updatescale!). The2α₀/2β₀pseudo-counts make a refinement inherit the scoring prior (β₀floors ŝ² at low inlier counts);IGPrior(0, 0)(Jeffreys) givesrss/dof`.

The first argument is THE ESTIMATOR: IGPrior selects the inverse-gamma MAP, its prior being that estimator's parameters. This is the fit contract: the residuals are TAKEN as inliers, no initial guess is needed, and no inlier decision is made here. fit_prior in this package is the precedent.

(A 4-argument form, fit_scale(prior, sq_residuals, weights, n) = Σwᵢrᵢ²/Σwᵢ, is deleted: its one caller was update_scale!, and that per-point moment WAS the missing ÷dg — for dg-vector residuals it estimates d_g·σ².)

The result is floored at eps(), so a zero rss under Jeffreys returns eps(), not a NaN. A caller that would rather HOLD its previous scale in that case must test the dof itself at the call site.

PosteriorWeightState(sq_residuals, n, model_dof, k_prev::Int;
               sws, w, inlier_posterior, sq_r, log_two_outlier_halfwidth, codimension=1)

Pre-sorted init: reuses an existing partition in sws.perm instead of running a fresh sweep. k_prev is the optimal inlier count from a preceding sweep on these same squared residuals (typically the RANSAC scoring sweep in _score_and_sweep!).

Skips the O(n log n) sort entirely — the sweep already ran once this hypothesis cycle. Prior is not a parameter here: the partition already reflects whatever prior the caller's sweep used, so (ν₀, ŝ²₀) inherits it automatically.

Precondition: sws.perm[1:k_prev] is the sort order of sq_residuals[1:n] from the preceding sweep, and k_prev > model_dof.

PosteriorWeightState(n; ν, s_hat_sq, log_two_outlier_halfwidth, codimension = 1, inlier_prior = 1/2,
                     prior = JeffreysPrior())

A state with FRESH length-n buffers at the given mixture, for a caller that fills the residuals itself and runs the E-step (e.g. mixture_ecm!) before reading any weight. The buffers are uninitialised; ŝ² is floored at eps and the inlier prior clamped to [1e-3, 1 − 1e-3] (the one clamp, update_inlier_prior!'s).

median_seeded_posterior(sq, n; outlier_halfwidth, codimension = 2, ν = 1e4,
                        prior = IGPrior()) -> PosteriorWeightState

A PosteriorWeightState seeded from the MEDIAN residual, for a caller with no sweep! partition to seed from. Every point starts an inlier candidate at the scale ŝ² = median(r²)/median(χ²_d) (floored at eps), with ν large (a near-Gaussian start) and the E-step run once. Allocates its own buffers (sq[1:n] is copied).

codimension is 1 or 2. median(χ²_d) is VisualGeometryCore's chisq_median(d), the stack's one definition of it.

Moved here from RobustVisualGeometry (scoring_glue.jl, used by fit_division_camera), verbatim: it builds this package's state from this package's E-step, so it belongs beside the other two constructors.

partition_seeded_posterior(residuals, n, scoring::MarginalScoring; model_dof,
                           sws = alloc_scoring_workspace(n, codimension(scoring), scoring),
                           w, inlier_posterior, sq_r)
    -> Union{PosteriorWeightState, Nothing}

A PosteriorWeightState seeded from the SCORER's own partition of the residuals residuals[1:n] (unsquared; norms for a codimension-2 residual). score_residuals! squares them into sq_r and sweeps them under scoring — the scale-free partition, at the scoring's own model_dof, outlier half-width, codimension and scale prior — and the pre-sorted PosteriorWeightState(sq_r, n, model_dof, k*; …) seeds the Student-t of the selected k* points (partition_student_t: ŝ² by fit_scale at the residual dof (k* − model_dof)·d_g) and runs the E-step once.

The seed's ν is the partition's posterior-predictive dof, 2α₀ + (k* − model_dof)·d_g — the dof a Gaussian with an inverse-gamma scale posterior marginalizes to. An EM that holds it there (PartitionDof) keeps ν a statement of how well the partition determines σ rather than of the inlier tail's shape.

model_dof is the SEED's: the scoring's minimal sample size plus one for a seed that also estimates the scale (RobustVisualGeometry's LO seed), or the residual-dof count of a fixed fit (a bundle's free parameters per correspondence). Returns nothing when the partition is too thin to carry a scale (k* ≤ model_dof: no residual dof).

THE residuals-only partition seeding — the same sweep and the same seed rule as the RANSAC loop's (which reuses the sweep it already ran, through the pre-sorted constructor directly). The buffers default to fresh length-n vectors; sws must hold at least n points. There is no default outlier half-width: it is the scoring's, which requires one.

soft_posterior_weights!(w, inlier_posterior, sq_residuals, ν, s_hat_sq, log_two_outlier_halfwidth, n, d_g=1)

E-step: compute posterior inlier probability and combined IRLS weight.

Stores inlier_posterior[i] = P(inlier | r_i) separately for the ν M-step, and w[i] = inlier_posterior[i] · u_i (combined weight for the θ solve).

d_g is the constraint codimension (1 for scalar residuals, 2 for 2D). The Student-t density uses (ν + d_g)/2 instead of (ν + 1)/2 for the multivariate case, and the normalizer becomes (ν π s²)^{d_g/2}.

The separation is critical: the ν Newton update must use fixed E-step weights (not recomputed from the changing ν), otherwise ν and inlier_posterior couple in a positive feedback loop that drives ν → ∞.

update_weights!(pws, n)

E-step: recompute posterior inlier probabilities and combined IRLS weights from current ν, ŝ², and squared residuals.

update_scale!(pws, n)

M-step, scale block: overwrite ŝ² from the E-step's weights and the current squared residuals, and nothing else. Call after update_residuals! fills pws.sq_residuals; the ν block is update_nu!, called after this.

ŝ² is the ECM step of the d_g-variate Student-t with scale ŝ²·I:

ŝ²_new = Σ πᵢ uᵢ rᵢ² / (d_g · Σ πᵢ),     uᵢ = (ν + d_g)/(ν + rᵢ²/ŝ²),

with πᵢ the inlier posteriors and rᵢ² squared NORMS. The E-step's weight is precision-scaled, wᵢ = πᵢ·(ν + d_g)/(νŝ² + rᵢ²) = πᵢuᵢ/ŝ², so the numerator is ŝ²·Σwᵢrᵢ² at the ŝ² that E-step used, and the inverse-gamma prior enters as fit_scale(prior, ŝ²·Σwᵢrᵢ², d_g·Σπᵢ). That is the per-component convention of every counted-dof fit_scale.

Until 2026-09-22 this was Σwᵢrᵢ²/Σwᵢ, whose fixed point is dg·σ² (ŝ ≈ 1.40σ at dg = 2 on Gaussian data), and whose prior pseudo-counts entered in the units of the weights. At d_g = 1 under Jeffreys the fixed point is unchanged (both are stationary where Σπu = Σπ); the path to it moves.

RSDR's own function, NOT VGC's policy verb — VGC's 1-arg sweep verb on AbstractWeightPolicy is update_nuisance!, the composite block whose parts here are update_scale! then update_nu!. RVG's PosteriorLikelihood supplies update_nuisance! by calling the two in that order. Returns ŝ².

update_scale!(pws, n, policy) -> ŝ²

The scale block under its update policy. Estimate() is update_scale!(pws, n); VisualGeometryCore's FixedScale(s2) holds ŝ² at s2 (a positive finite SQUARED scale; FixedScale() carries none and is refused).

update_nu!(pws, n)

M-step, ν block: overwrite the Student-t dof ν with next_nu (the safeguarded Newton MLE) from the E-step's inlier posteriors, the current squared residuals and the CURRENT ŝ², and nothing else. Call after update_scale!, which it reads. Returns ν.

update_nu!(pws, n, policy) -> ν

The ν block under its update policy: Estimate() is update_nu!(pws, n); FixedNu holds ν.

next_nu(ν, sq_residuals, s_hat_sq, w_fixed, n; d_g=1, kwargs...) -> ν_new

M-step: Newton update for ν maximizing the EM Q-function with fixed E-step weights.

Maximizes Q(ν) = Σ wi · log t{dg}(ri; ν, ŝ²) w.r.t. ν, where wi are the fixed posterior inlier probabilities from the E-step (not recomputed), ri² the squared NORM of a d_g-vector residual and ŝ² the per-component scale:

log t_d(r; ν, ŝ²) = logΓ((ν+d)/2) − logΓ(ν/2) − (d/2)·log(νπŝ²) − ((ν+d)/2)·log(1 + z), z = r²/(νŝ²)

∂Q/∂ν = Σ wi · [½ψ((ν+d)/2) − ½ψ(ν/2) − d/(2ν) − ½log(1 + zi) + (ν+d)zi/(2ν(1+zi))]

∂²Q/∂ν² = Σ wi · [¼ψ₁((ν+d)/2) − ¼ψ₁(ν/2) + d/(2ν²) + zi/(2ν(1+zi)) − zi(ν + 2d + d·zi)/(2ν²(1+zi)²)]

(_nu_derivatives), with z_i the SCALE-NORMALISED squared residual — it was written with the symbol this package reserves for a decision threshold. The last term is the derivative of the gradient's (ν+d)z/(2ν(1+z)) through dz/dν = −z/ν. Until 2026-09-22 it read −z(ν+d+1+z)/(2ν²(1+z)²), which equals it only at d = 1: at d = 2 on t₄ data at ν = 4 the true curvature is ≈ −206 and that form gave ≈ +3.3, so Newton stepped AWAY from the maximum.

Up to max_newton Newton steps from the current ν, each clamped to [ν_min, ν_max]. Newton is taken only where Q is concave (∂²Q < 0), where its step is uphill; elsewhere the step is geometric, uphill along the gradient — ν ← 2ν if ∂Q/∂ν > 0, ν ← ν/2 otherwise (Q is not concave in ν everywhere: at d = 2 on t₄ data it is convex near ν = 10).

update_inlier_prior!(pws, n)

M-step, inlier-prior block: overwrite the mixing proportion π with the mean inlier posterior of the E-step, and nothing else. The mean Σπᵢ/n is clamped to [1e-3, 1 − 1e-3] (the E-step reads it through log((1 − π)/π)). Call after update_weights! filled the posteriors it reads. The ECM of a mixture whose π is estimated rather than fixed at ½ runs update_weights!, update_scale!, then this. Returns π.

update_inlier_prior!(pws, n, policy) -> π

The inlier-prior block under its update policy: Estimate() is update_inlier_prior!(pws, n); FixedPrior holds π.

mixture_ecm!(pws, r, n)

One ECM cycle of the Student-t + uniform mixture with the mixing proportion ESTIMATED (ν held fixed), from the residuals r[1:n] (unsquared). It squares them into sq_residuals(pws), then runs the E-step (update_weights!), the scale block (update_scale!) and the inlier-prior block (update_inlier_prior!), in that order. Each block is an exact conditional maximisation, so the mixture NLL (mixture_nll) never rises. Returns nothing.

posterior_weights!(w, sws, sq_residuals, k; model_dof, log_two_outlier_halfwidth, codimension=1)

Compute Student-t posterior weights w[i] = inlier_posterior[i] · u_i from the partition produced by sweep!. inlier_posterior is P(inlier | r) under the Student-t vs uniform mixture, u_i is the Student-t IRLS weight (ν + dg) / (RSS + r²_i).

sws.perm must still hold the sort from the preceding sweep! call.

The scale prior is sws.scale_prior — the prior that sweep scored the partition under, to which the workspace is committed (see ScoringWorkspace). The Student-t is that prior's posterior predictive over the partition, partition_student_t: ν = (k − model_dof)·dg + 2α₀ degrees of freedom (at least 2) and squared scale ŝ² = (RSS + 2β₀)/((k − model_dof)·dg + 2α₀), so u_i = (ν + dg)/(ν·ŝ² + r²_i) — the predictive PosteriorWeightState seeds from. At the Jeffreys prior (α₀ = β₀ = 0) the pseudo-counts vanish. Per point this is posterior_weight of that StudentTMixture.

The prior used to be IGNORED here: an Inverse-Gamma workspace swept its partition under IG and then weighted it under Jeffreys.

posterior_inlier_probs!(w, sws, sq_residuals, k; model_dof, log_two_outlier_halfwidth, codimension=1)

Compute posterior inlier probabilities w[i] ∈ [0,1] without the IRLS weight factor. Same preconditions — and the same sws.scale_prior — as posterior_weights!, so posterior_weights! is this times u_i.

Estimate()

The update policy of a nuisance part that is ESTIMATED: its M-step block runs (update_scale!(pws, n), update_inlier_prior!(pws, n), update_nu!(pws, n)). One value for every part.

FixedPrior(inlier_prior)

The update policy of an inlier prior HELD at inlier_prior ∈ (0, 1): the π block writes it (through the E-step's one clamp, [1e-3, 1 − 1e-3]).

FixedNu(ν)

The update policy of a Student-t dof HELD at ν (positive, finite): the ν block writes it.

PartitionDof()

The update policy of a Student-t dof HELD at the value its partition seed gave it. That value is the posterior-predictive dof of the marginal scorer's partition, ν = 2α₀ + (k* − model_dof)·d_g (partition_student_t, the one statement of that count, with its ν ≥ 2 floor), where α₀ is the scale prior's shape (0 under Jeffreys). partition_seeded_posterior and the pre-sorted constructor set it once; the ν block never moves it afterwards.

Why hold it there: a Gaussian whose scale carries an inverse-gamma posterior marginalizes to a Student-t with exactly that dof — ν measures how well the partition DETERMINES σ (it grows with the residual dof), not the shape of the inlier tail. Estimating ν instead lets the fit trade one for the other: on a contaminated set a heavy inlier tail (ν → 2) absorbs the outliers, the mixing proportion runs to its clamp and almost nothing is classified an outlier — the identifiability collapse this policy prevents.

inlier_prior(pws::PosteriorWeightState) -> Float64

The mixing proportion π the E-step weighs with: the prior probability that a residual is an inlier (1/2 unless estimated or held elsewhere).

StudentTMixture(ν, s_hat_sq, log_two_outlier_halfwidth, d_g[, inlier_prior = 1/2])
StudentTMixture{T}(m::StudentTMixture)

The Student-t inlier + uniform outlier mixture of one E-step, with every per-residual constant precomputed. Its parameters are ν degrees of freedom, squared scale ŝ² = s_hat_sq, outlier support of half-width outlier_halfwidth per axis (log_two_outlier_halfwidth = log(2·outlier_halfwidth), density p_out = (2·outlier_halfwidth)^-d_g), inlier mixing proportion π = inlier_prior. The loggammas are evaluated HERE, once.

An immutable isbits value, so it can be an argument of a KernelAbstractions kernel, and the three scalar functions of it — inlier_posterior, posterior_weight, mixture_nll — are allocation-free plain functions of a squared residual, callable inside a kernel as VGC's ρ is. StudentTMixture{Float32}(m) converts the constants for a Float32 kernel.

These are THE formulas of the mixture: soft_posterior_weights! (the E-step over a buffer) and posterior_weights!/posterior_inlier_probs! (over a sweep! partition) build one of these and call the same scalar kernel.

inlier_posterior(m::StudentTMixture, sq_r) -> P(inlier | r)

The posterior probability that a residual with squared norm sq_r is an inlier, π·t / (π·t + (1−π)·p_out). Scalar and kernel-callable.

posterior_weight(m::StudentTMixture, sq_r) -> w

The IRLS weight of the mixture at sq_r: P(inlier | r) · (ν + d_g)/(ν·ŝ² + r²), the responsibility times the Student-t's own weight. Scalar and kernel-callable.

posterior_weight_derivative(m::StudentTMixture, sq_r) -> dw/d(sq_r)

The slope of posterior_weight in the squared residual, −w/(ν·ŝ² + r²) · (1 + (ν + d_g)/2 · (1 − P(inlier | r))) — the Student-t weight's own decay plus the responsibility's. w is strictly decreasing in sq_r, so this is negative wherever w > 0; it is what a Newton inversion of the weight (a tabulated half-quadratic conjugate, VisualGeometryCore's TabulatedLift) needs. Scalar and kernel-callable.

mixture_nll(m::StudentTMixture, sq_r) -> −log(π·t_ν(r; ŝ²) + (1−π)·p_out)

The negative log-density of the mixture at one residual — the per-point summand of the mixture NLL. Computed against the outlier term in log space, −(log((1−π)p_out) + softplus(log π·t − log((1−π)p_out))), so it is finite for an infinite residual (it is the outlier term's) and for an exact inlier. Scalar and kernel-callable.

mixture_nll(pws::PosteriorWeightState, r, n) -> Float64

The mixture negative log-likelihood Σ_{i ≤ n} −log[π·t_ν(rᵢ; ŝ²) + (1−π)·p_out] of the residuals r[1:n] (UNSQUARED) at the state's current (ν, ŝ², π). It sums the per-point mixture_nll(m, rᵢ²) of StudentTMixture(pws) in Float64. Reads the state; writes nothing.

partition_student_t(prior::IGPrior, rss, k; model_dof, codimension = 1) -> (ν, ŝ²)

The Student-t of a hard partition of k points with residual sum of squares rss. It is the scoring marginal's posterior predictive over it under prior, and THE seed of every Student-t in the stack: counted dof ν_data = (k − model_dof)·codimension, then

ν  = max(ν_data + 2α₀, 2),     ŝ² = fit_scale(prior, rss, ν_data).

The ν ≥ 2 clamp keeps the seed's variance finite (a Student-t has one only for ν > 2) at a thin partition, and fit_scale floors the QUOTIENT at eps. One rule for every caller: both PosteriorWeightState seeds, the partition mixture of posterior_weights!/posterior_inlier_probs!, and RobustVisualGeometry's LO seed (SeedPartition) and result scale (_finalize). model_dof is the caller's: the sweep's minimal_sample_size for the scoring marginal, one more for an LO seed that also estimates the scale.

square_residuals!(sq, r, n) -> sq

sq[i] = r[i]² for i ≤ n, squared in sq's OWN precision (x = eltype(sq)(r[i]), then x·x): a Float64 buffer squares a Float32 residual exactly. THE one residual-squaring loop of the stack — the refill of a PosteriorWeightState's sq_residuals (the ECM below, RobustVisualGeometry's posterior seeding and IRLS residual block), the scores score_residuals! sweeps, and the pool staging (precompute_pair).

weighted_sq(w, sq_r) -> w·sq_r, exactly 0 at w = 0

One term of weighted_rss, for a caller summing weights it computes on the fly rather than holding in a buffer (an ECM scale over a loss's own IRLS weights). A zero weight contributes exactly 0 even against an infinite squared residual, where 0·Inf = NaN.

weighted_rss(w, sq_residuals, n) -> Float64

The weighted residual sum of squares Σ_{i ≤ n} wᵢ·rᵢ² of SQUARED residuals, accumulated in Float64 left to right from 0.0. A zero-weight term contributes exactly 0: an excluded point may carry an infinite squared residual (the Sampson degeneracy rule — a vanishing projected variance), and 0·Inf = NaN would poison the sum and everything read from it (the ŝ² M-step of update_scale!, an IRLS cost and its stop rule).

s_hat_sq(pws::PosteriorWeightState) -> Float64

The current squared scale estimate ŝ² of the Student-t inlier density (a SQUARE; see fit_scale).

student_t_dof(pws::PosteriorWeightState) -> Float64

The Student-t degrees of freedom ν of the inlier density — not the model's DOF (minimal_sample_size), and not Distributions' dof.

partition_rss(sws::ScoringWorkspace, residuals, k) -> Float64

The residual sum of squares of the partition the last sweep! selected. It is Σ_{j ≤ k} residuals[sweep_order(sws)[j]]², from UNSQUARED residuals (signed scalars), accumulated in Float64 whatever their storage eltype.

For a caller whose sweep ran on a DIFFERENT per-point key than these residuals — RobustVisualGeometry's seed partition under predictive scoring sweeps the Mahalanobis scores but reads its scale off the plain residuals — and which so cannot use the squared buffer the sweep sorted. Pair it with fit_scale(scale_prior(sws), partition_rss(sws, r, k), ν_data).

partition_rss(perm, residuals, k) -> Float64

The same sum over an explicit sweep permutation — one a caller copied out of the workspace (RobustVisualGeometry's incumbent best_perm).

sq_residuals(pws::PosteriorWeightState) -> Vector{Float64}

The squared residuals the E-step weighs — the buffer update_scale! and update_weights! read, refilled by the caller between iterations.