Struct ConstrainedPosteriorCorrection 

pub struct ConstrainedPosteriorCorrection {
    pub lift: Array2<f64>,
    pub removed_normal_variance: Array2<f64>,
    pub normal_mean_shift: Array1<f64>,
    pub rows: Vec<usize>,
    pub normal_upper_limits: Vec<f64>,
}

The low-rank correction that turns the unconstrained Laplace covariance into the truncated-posterior covariance.

Carrying the factored form rather than a dense p × p matrix lets a consumer that never materializes Σ (the factorized inference path, the prediction backends) apply the same correction with q extra solves: xᵀΣ_π x = xᵀΣx − ‖Δ^{1/2} Gᵀ x‖².

Fields

lift: Array2<f64>

G = Σ Aᵀ W⁻¹, p × q.

removed_normal_variance: Array2<f64>

Δ = W − Cov[u] ⪰ 0, q × q: the variance the truncation removes from the constraint-normal coordinates.

normal_mean_shift: Array1<f64>

E[u] − E_untrunc[u], q: how far truncation moves the posterior mean in constraint-normal coordinates. Positive componentwise for a half-line coordinate — the posterior mean is interior even when the mode is not — and of either sign for a coordinate that also carries an upper limit, where the far face pulls the mean back down.

rows: Vec<usize>

Indices, into the caller’s constraint system, of the rows retained.

normal_upper_limits: Vec<f64>

Upper limit on each retained coordinate: u_k ≤ normal_upper_limits[k], with f64::INFINITY where the coordinate is a half-line.

A two-sided coefficient bound l ≤ β_j ≤ u arrives as two rows whose normals are exactly anti-parallel. The second carries no constraint-normal DIRECTION the first does not already carry, so the rank filter drops it — correctly, as a direction. It is still a constraint, and this is where it is kept (#2523).

#[serde(default)] with an empty vector reading as “every retained coordinate is a half-line” is the encoding of a model saved before upper limits existed, which is exactly what those models meant. Live constructions always carry one entry per retained row.

+∞ is the COMMON case here — every one-sided shape or box constraint contributes one — and JSON has no literal for it. Without an explicit codec serde_json writes the entry as null and Vec<f64> then refuses its own output on the way back in, which made every shape-constrained fit that retains a face unloadable (#2601). serde_extended_real::vec_f64 makes null mean +∞ in both directions; that is byte-identical to what was already being written, so models saved before the codec existed read back with the meaning they always had.

Implementations

impl ConstrainedPosteriorCorrection

pub fn apply_to_covariance_in_place(&self, covariance: &mut Array2<f64>)

Σ ← Σ − G Δ Gᵀ, in place. The correction is rank q, so this never allocates a second p × p matrix next to the one being corrected.

pub fn apply_to_covariance(&self, covariance: &Array2<f64>) -> Array2<f64>

Σ_π = Σ − G Δ Gᵀ.

pub fn truncated_covariance_psd( &self, covariance: &Array2<f64>, constraints: &LinearInequalityConstraints, ) -> Result<Array2<f64>, String>

The same Σ_π, assembled as a SUM OF TWO GRAMS instead of as a subtraction — so its diagonal cannot be a cancellation and cannot come out negative (#2705 group A).

Σ − GΔGᵀ is the difference of two nearly equal numbers exactly where the answer matters most. A coordinate the constraint PINS has essentially all of its variance removed: on y ~ s(x, shape=convex) the measured entry went from Σ_ii = 2.30e-2 to 6.23e-13 — eleven digits gone — and on the neighbouring sqrt fixture the same subtraction lands at −3.09e-15, which is not a small variance but a rounding residue with a sign. Everything downstream (se_from_covariance, the dense SE loop) then has to argue about whether that sign is real.

Split the correction at Δ = W − C, C = Cov[u] ⪰ 0 the truncated constraint-normal covariance, and the same quantity is two Grams:

    Σ − GΔGᵀ = (Σ − G W Gᵀ) + G C Gᵀ = P Σ Pᵀ + G C Gᵀ,   P = I − G A.

With Σ = L Lᵀ and C = L_C L_Cᵀ that is (P L)(P L)ᵀ + (G L_C)(G L_C)ᵀ, and every diagonal entry is a sum of squares. The cancellation does not disappear — it moves INSIDE P L, where each entry carries an absolute error O(ε‖L‖) and is then SQUARED, so a pinned coordinate’s variance picks up O(p ε² Σ_ii) instead of O(ε Σ_ii): sixteen orders smaller, and non-negative by construction rather than by luck.

covariance must be the SPD matrix this correction was built for, and constraints the system its rows index. Both are exactly what the caller already holds where a dense Σ exists at all.

C is a cubature result, so it can carry a small negative eigenvalue — certify_removed_variance admits Δ_ii up to slack·W_ii past W_ii. Eigenvalues inside that certified band are read as the zero they are approximating; anything below it is refused, because a C that is materially indefinite is a broken moment computation and not a rounding question.

pub fn removed_variance_diagonal(&self) -> Array1<f64>

diag(G Δ Gᵀ) — the per-coefficient variance the truncation removes, for consumers that only ever build the covariance diagonal.

pub fn diagonal_uncertainty(&self) -> Array1<f64>

The absolute per-coefficient uncertainty this correction contributes to diag(Σ − GΔGᵀ).

Δ is a cubature result, not an exact quantity: it is certified to ORTHANT_MOMENT_RELATIVE_TOLERANCE relative (certify_removed_variance), and (GΔGᵀ)_ii = g_iᵀ Δ g_i is monotone in Δ in the PSD order, so a relative error ε in Δ moves the removed variance by at most ε · (GΔGᵀ)_ii. That product — and NOT the floating-point backward error of the subtraction — is the resolution at which diag(Σ − GΔGᵀ) can be read.

This accessor exists so that ORTHANT_MOMENT_RELATIVE_TOLERANCE is declared once and converted into a consumer-facing allowance once. #2705 group A is fits refused because the consumer (se_from_covariance, budget 16·n·eps ≈ 1e-14 relative) and the producer (this cubature, budget 1e-3 relative) hold two independent budgets for one number, ~11 orders apart, with nothing carrying the producer’s across the boundary.

pub fn posterior_mean(&self, unconstrained_center: &Array1<f64>) -> Array1<f64>

E_π[β] = β_unc + G·(E[u] − E_untrunc[u]).

pub fn upper_limits(&self) -> Vec<f64>

Upper limit per retained coordinate, materializing the legacy encoding of Self::normal_upper_limits.

Trait Implementations

impl Clone for ConstrainedPosteriorCorrection

fn clone(&self) -> ConstrainedPosteriorCorrection

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl Debug for ConstrainedPosteriorCorrection

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<'de> Deserialize<'de> for ConstrainedPosteriorCorrection

fn deserialize<__D>(__deserializer: __D) -> Result<Self, __D::Error>

where __D: Deserializer<'de>,

Deserialize this value from the given Serde deserializer.

impl Serialize for ConstrainedPosteriorCorrection

fn serialize<__S>(&self, __serializer: __S) -> Result<__S::Ok, __S::Error>

where __S: Serializer,

Serialize this value into the given Serde serializer.