Struct InnerSolution 

pub struct InnerSolution<'dp> {

    pub log_likelihood: f64,
    pub penalty_quadratic: f64,
    pub hessian_op: Arc<dyn HessianOperator>,
    pub beta: Array1<f64>,
    pub penalty_coords: Vec<PenaltyCoordinate>,
    pub penalty_logdet: PenaltyLogdetDerivs,
    pub deriv_provider: Box<dyn HessianDerivativeProvider + 'dp>,
    pub firth: Option<ExactJeffreysTerm>,
    pub hessian_logdet_correction: f64,
    pub penalty_subspace_trace: Option<Arc<PenaltySubspaceTrace>>,
    pub rho_curvature_scale: f64,
    pub rho_prior: RhoPrior,
    pub n_observations: usize,
    pub nullspace_dim: f64,
    pub gaussian_weight_log_sum_half: f64,
    pub dp_floor_scale: f64,
    pub dispersion: DispersionHandling,
    pub ext_coords: Vec<HyperCoord>,
    pub ext_coord_pair_fn: Option<Box<dyn Fn(usize, usize) -> HyperCoordPair + Send + Sync>>,
    pub rho_ext_pair_fn: Option<Box<dyn Fn(usize, usize) -> HyperCoordPair + Send + Sync>>,
    pub fixed_drift_deriv: Option<FixedDriftDerivFn>,
    pub contracted_psi_second_order: Option<ContractedPsiSecondOrderFn>,
    pub barrier_config: Option<BarrierConfig>,
    pub kkt_residual: Option<ProjectedKktResidual>,
    pub active_constraints: Option<Arc<ActiveLinearConstraintBlock>>,
    pub stochastic_trace_state: Arc<Mutex<StochasticTraceState>>,
}

The unified inner solution produced by any inner solver.

Contains everything the outer REML/LAML evaluator needs. Produced by:

• Single-block PIRLS (via PirlsResult::into_inner_solution())
• Blockwise coupled Newton (via BlockwiseInnerResult::into_inner_solution())
• Sparse Cholesky (via SparsePenalizedSystem::into_inner_solution())

Fields

log_likelihood: f64

ℓ(β̂) — log-likelihood at the converged mode. For Gaussian: −0.5 × deviance (RSS). For GLMs: actual log-likelihood.

penalty_quadratic: f64

β̂ᵀS(ρ)β̂ — penalty quadratic form at the mode.

hessian_op: Arc<dyn HessianOperator>

The Hessian operator providing logdet, trace, and solve. Both cost and gradient use this same object.

IMPORTANT: This MUST encode the observed Hessian H_obs = X’W_obs X + S at the converged mode, where W_obs includes the residual-dependent correction for non-canonical links. Using expected Fisher H_Fisher = X’W_Fisher X + S would make this a PQL surrogate rather than the exact Laplace approximation. See response.md Section 3 for the mathematical justification.

beta: Array1<f64>

β̂ — coefficients at the converged mode (in the operator’s native basis).

penalty_coords: Vec<PenaltyCoordinate>

Penalty coordinates for the rho block.

Each coordinate represents one smoothing-parameter direction A_k = λ_k S_k through either a full-root or a block-local root.

penalty_logdet: PenaltyLogdetDerivs

Derivatives of log|S(ρ)|₊ — precomputed from penalty structure.

deriv_provider: Box<dyn HessianDerivativeProvider + 'dp>

Provider of third-derivative corrections for non-Gaussian families.

The c and d arrays (dW/deta, d^2W/deta^2) carried by this provider MUST be the observed derivatives, not the Fisher derivatives. For non-canonical links the observed c/d include residual-dependent corrections: c_obs = c_Fisher + h’*B - (y-mu)*B_eta d_obs = d_Fisher + h’‘B + 2h’*B_eta - (y-mu)*B_etaeta These corrections matter for the outer gradient (C[v] correction) and outer Hessian (Q[v_k, v_l] correction). See response.md Section 3.

firth: Option<ExactJeffreysTerm>

Optional exact Jeffreys/Firth term in the active coefficient basis.

hessian_logdet_correction: f64

Additive correction for the Hessian logdet when hessian_op encodes a uniformly rescaled exact curvature matrix.

penalty_subspace_trace: Option<Arc<PenaltySubspaceTrace>>

When the cost uses log|U_Sᵀ H U_S|_+ (rank-deficient LAML fix), this carries the matching projected kernel so the gradient trace tr(K · Ḣ) agrees with the cost’s derivative. See PenaltySubspaceTrace for the full derivation.

rho_curvature_scale: f64

Uniform scale s applied to rho-coordinate penalty derivatives in the H-dependent trace / solve parts of the outer calculus.

Contract (CRITICAL — gradient/cost consistency)

rho_curvature_scale is NOT a free knob. It encodes the convention that the supplied hessian_op represents the rescaled curvature H_op = s · (∇²(-ℓ) + Σ_k e^{ρ_k} S_k), i.e. every contribution to the curvature (likelihood Hessian AND penalty λ_k S_k) has been uniformly multiplied by s before reaching the evaluator. Under this convention:

• ∂H_op/∂ρ_k = s · λ_k S_k (matches the curvature_lambdas = s · λ drift used inside the gradient’s trace term),
• K = H_op⁻¹ = (1/s) · (∇²(-ℓ) + λS)⁻¹,
• tr(K · ∂H_op/∂ρ_k) = tr((∇²(-ℓ) + λS)⁻¹ · λ_k S_k) (the analytic gradient of the unscaled log|H|),
• log|H_op| = log|∇²(-ℓ) + λS| + p · log(s), which the caller MUST un-scale by supplying hessian_logdet_correction += −p · log(s) so that hop.logdet() + hessian_logdet_correction evaluates the same unscaled log|H| whose derivative the gradient trace computes.

Callers that set rho_curvature_scale ≠ 1 without ALSO pre-scaling hessian_op AND adding the matching −p·log(s) term to hessian_logdet_correction will get a gradient that is off by the factor s from dV/dρ_k. The unified evaluator does not scale hop for the caller — that would defeat the purpose of the curvature-conditioning trick survival families use to keep the outer eigendecomposition numerically stable.

See survival::location_scale::exact_newton_outer_curvature for the canonical example: rho_curvature_scale = exp(-log_scale) paired with hessian_logdet_correction = p · log_scale = −p · log(scale).

The evaluator enforces rho_curvature_scale > 0 and finite; pass 1.0 (the documented default) when no curvature conditioning is in play.

rho_prior: RhoPrior

Configured prior over rho coordinates. The evaluator receives the realized cost/gradient tuple separately; this copy lets EFS use the conjugate Gamma rate in its multiplicative denominator.

n_observations: usize

Number of observations.

nullspace_dim: f64

M_p: dimension of the penalty null space (unpenalized coefficients).

gaussian_weight_log_sum_half: f64

½·Σᵢ log(wᵢ) — half the sum of log prior weights.

This is the per-observation Gaussian normalization constant that the log_likelihood (computed by [calculate_loglikelihood_omitting_constants]) deliberately drops. The full weighted-Gaussian negative log-likelihood normalization is ½·Σᵢ log(2π·φ/wᵢ) = (n/2)·log(2πφ) − ½·Σᵢ log(wᵢ), because Var(yᵢ) = φ/wᵢ under inverse-variance prior weights.

Dropping −½·Σ log(wᵢ) does not move the ρ-argmin in exact arithmetic (it is constant in ρ), but it makes the ProfiledGaussian objective VALUE scale-dependent: under a global rescale w → c·w the invariance- preserving smoothing λ → c·λ leaves the cost SHAPE fixed but inflates its absolute value by (n/2)·log c. That inflation breaks the exact weight-scale invariance of the selected λ̂ / EDF / fit (issue #877). Restoring this term makes the ProfiledGaussian cost value exactly invariant to w → c·w (with σ̂² absorbing the c factor), matching mgcv.

Only consumed by the ProfiledGaussian arm; the Fixed-dispersion arm already omits the Gaussian normalization constant by design and is not affected.

dp_floor_scale: f64

Deviance scale D₀ used as the relative reference for the smooth penalized-deviance floor (see [crate::estimate::smooth_floor_dp]).

Set to the weighted null deviance of the Gaussian response, D₀ = Σ wᵢ(yᵢ − ȳ_w)², which is the natural upper reference for D_p and — crucially — transforms as D₀ → a²·D₀ under a response rescale y → a·y, exactly as D_p does. Flooring D_p at a fixed fraction of D₀ therefore keeps the profiled Gaussian REML criterion exactly scale-equivariant (issue #1127); an absolute floor does not.

Only consumed by the ProfiledGaussian arm. Defaults to 1.0, which reproduces the historical absolute floor byte-for-byte for every caller that does not supply a response scale.

dispersion: DispersionHandling

How the dispersion parameter is handled.

ext_coords: Vec<HyperCoord>

External (non-ρ) hyperparameter coordinates with their fixed-β objects. These are appended after the ρ coordinates in the gradient/Hessian output.

ext_coord_pair_fn: Option<Box<dyn Fn(usize, usize) -> HyperCoordPair + Send + Sync>>

Callback to compute second-order fixed-β objects for a pair (i, j) of external coordinates (or external × ρ cross pairs). Arguments: (ext_index_i, ext_index_j) → HyperCoordPair. When None, the outer Hessian is not computed for extended coordinates.

rho_ext_pair_fn: Option<Box<dyn Fn(usize, usize) -> HyperCoordPair + Send + Sync>>

Callback for ρ × ext cross pairs: (rho_index, ext_index) → HyperCoordPair.

fixed_drift_deriv: Option<FixedDriftDerivFn>

M_i[u] = D_β B_i[u] callback for extended coordinates. Arguments: (ext_index, direction) → correction matrix.

contracted_psi_second_order: Option<ContractedPsiSecondOrderFn>

Direction-contracted second-order ψ hook for the profiled θ-HVP (#740). When present, the outer-Hessian operator builder skips the K² per-pair base_h2 ψψ assembly and instead applies this once per matvec to obtain every output row’s tr(K · D²_ψ H_L[ψ_i, ψ(α)]) in a single family row pass. None keeps the exact per-pair assembly. See ContractedPsiSecondOrderFn.

barrier_config: Option<BarrierConfig>

Optional log-barrier configuration for monotonicity-constrained coefficients. When present, the barrier cost and Hessian corrections are added to the outer REML/LAML objective.

kkt_residual: Option<ProjectedKktResidual>

Optional inner KKT residual r = ∇_β L_pen(β̂) at the converged β̂, already projected onto the free subspace (see ProjectedKktResidual for the invariant and why the type wraps this). Some activates the implicit-function-theorem corrections in reml_laml_evaluate (cost gets −½ rᵀ H⁻¹ r, ρ-gradient and ρρ Hessian get the matching first and second derivatives of that same scalar correction). None keeps the envelope-only behaviour for callers that genuinely guarantee exact KKT.

active_constraints: Option<Arc<ActiveLinearConstraintBlock>>

Optional active linear-inequality constraints at the converged inner iterate. Some(rows) means the joint constraint matrix’s row indices in rows.active_indices are pinned (treated as equality constraints at the cert point). The unified evaluator combines this with the penalty_subspace_trace to form the constraint-aware kernel K_T = K_S − K_S Aᵀ (A K_S Aᵀ)⁻¹ A K_S for per-coordinate IFT mode responses v_k = ∂β/∂ρ_k. See ConstrainedSubspaceKernel for the full derivation and consistency with log|U_Tᵀ H U_T|.

None is the legacy/unconstrained path (no active inequality constraints to project against).

stochastic_trace_state: Arc<Mutex<StochasticTraceState>>

Fit-level stochastic trace state. Shared by stochastic trace batches so CRN probe prefixes stay fixed and matrix-free trace CG can warm-start from the previous solve of the same probe id.