Struct MultinomialFamily 

pub struct MultinomialFamily {
    pub y_one_hot: Array2<f64>,
    pub weights: Array1<f64>,
    pub total_classes: usize,
    pub design: Arc<Array2<f64>>,
    pub penalties: Arc<Vec<PenaltyMatrix>>,
    /* private fields */
}

Joint-coupled multinomial-logit family with shared design and shared smoothing penalty across active classes.

Block layout

K − 1 parameter blocks, indexed a = 0..K-1, each carrying coefficient vector β_a ∈ ℝ^P. Class K − 1 is the reference (β_{K-1} ≡ 0) and does not appear in the block list.

Invariants

• y_one_hot.dim() == (N, K), with K = total_classes ≥ 2.
• weights.len() == N, finite and non-negative.
• design.nrows() == N, design.ncols() == P.
• every penalty in penalties has shape (P, P) (symmetric, PSD).

All are validated by MultinomialFamily::new.

Fields

y_one_hot: Array2<f64>

Categorical response matrix Y ∈ ℝ^{N × K}. Each row must be a point on the probability simplex (y_c ≥ 0, Σ_c y_c = 1): a one-hot indicator or a label-smoothed probability vector. Rows whose mass departs from 1 are rejected by MultinomialFamily::new — the softmax residual and Fisher block are the derivatives of Σ_c y_c log p_c only under the simplex constraint. Column K − 1 is the reference class.

weights: Array1<f64>

Per-row weights w ∈ ℝ^N, finite and non-negative.

total_classes: usize

Total class count K ≥ 2. Active classes are 0..K-1; class K − 1 is the reference.

design: Arc<Array2<f64>>

Shared design matrix X ∈ ℝ^{N × P}, identical across all active classes. Carried as Arc<Array2<f64>> so the per-block specs and the family share storage with zero copies.

penalties: Arc<Vec<PenaltyMatrix>>

Per-smooth-term penalty components, each a P × P operator expressed in block-local form (PenaltyMatrix::Blockwise embedding the term’s local S_t at its col_range within the shared P-column coefficient space). Every active class block receives this entire list, so the outer REML/LAML loop selects an independent smoothing parameter per (class, term) — matching mgcv/VGAM per-term smoothing. The full block-replicated penalty is I_{K-1} ⊗ (Σ_t λ_{a,t} S_t); pre-summing the terms (one fused λ per class) is exactly the multi-term fusion that over-smooths one term while under-smoothing another (#561). Carried as Arc<Vec<…>> so per-block specs share storage with zero copies.

Implementations

impl MultinomialFamily

pub const fn active_classes(&self) -> usize

Total number of active blocks, M = K − 1.

pub fn new( y_one_hot: Array2<f64>, weights: Array1<f64>, total_classes: usize, design: Arc<Array2<f64>>, penalties: Arc<Vec<PenaltyMatrix>>, ) -> Result<Self, String>

Validate inputs and construct the family.

All shape and finiteness invariants are checked here so the CustomFamily methods can rely on pre-validated geometry.

pub fn with_joint_jeffreys_span(self, span: Option<Arc<Array2<f64>>>) -> Self

Install the MEASURED Jeffreys/Firth span (gam#2612).

The caller is promising the basis is orthonormal, is expressed in the raw joint coefficient order, and does not change for the lifetime of the fit. See Self::joint_jeffreys_span and CustomFamily::jeffreys_span_basis.

pub fn with_joint_jeffreys_term(self, enabled: bool) -> Self

Select whether this multinomial adapter instance contributes the full-span Jeffreys/Firth correction.

pub fn with_initial_log_lambda(self, log_lambda: f64) -> Self

Seed the warm-start log λ carried into the reference-symmetric joint smoothing penalties (gam#1587). The formula REML driver sets this from its init_lambda so the joint-penalty outer ρ starts at the same seed the per-block path used historically; the outer loop then selects the optimum.

pub fn with_joint_initial_log_lambdas(self, seeds: Vec<f64>) -> Self

Seed PER-SPEC warm-start log λ values for the joint smoothing penalties, in the builders’ term-major spec order (equivariant carrier: s = t·K + c; shared centered carrier: s = t). Overrides the shared Self::with_initial_log_lambda seed entry-by-entry; the spec builders reject a wrong length. This is the resume path for a joint-penalty rho_checkpoint and the fixed-ρ pin for criterion diagnostics (#2349).

pub fn build_block_specs(&self) -> Vec<ParameterBlockSpec>

Build the canonical block specs for this family.

One ParameterBlockSpec per active class, all sharing the same design (zero-copy through Arc<Array2<f64>>) and an independent PenaltyMatrix::Dense copy of S. The gauge_priority is set so that the active class closest to the reference owns shared affine / null-space directions: class a gets priority 100 + (M − a). Class 0 (farthest from the reference) is the most likely to retain a shared direction in canonicalisation; class M − 1 is the least likely. This matches the task’s “descending priorities” gauge convention.

initial_log_lambdas is initialised to zeros (one entry per penalty term per block: each block carries one λ_{a,t} per smooth term t). Callers that want a custom warm start override per-block before passing to fit_custom_family_with_rho_prior.

pub fn beta_flat_dim(&self) -> usize

Total stacked-coefficient dimension (K − 1) · P.

pub fn centered_joint_penalty_specs( &self, ) -> Result<Vec<JointPenaltySpec>, String>

Build the reference-symmetric (“centered”) full-width smoothing penalties λ_t · (M ⊗ S_t), one per smooth term t, in raw stacked (class-major) coordinates [β_0; …; β_{K-2}] (gam#1587).

M = I_{K-1} − J_{K-1}/K is the closed-form CLR whitening metric of the softmax class gauge (the multinomial analogue of the resolved ALR sibling #1549). The quadratic form βᵀ (M ⊗ S_t) β equals the symmetric CLR penalty Σ_{k=0}^{K-1} β̃_{k}ᵀ S_t β̃_{k} over centered coefficients β̃_k = β_k − (1/K)Σ_b β_b (β_{K-1} ≡ 0), a symmetric function of all K classes — so the penalized fit no longer depends on which class is the arbitrary softmax reference. Block (a, b) of the returned (M·P)×(M·P) matrix is M[a,b]·S_t; M is SPD (eigenvalues 1 with multiplicity K−2 and 1/K once), so each M ⊗ S_t is PSD with nullspace_dim = (K−1)·nullspace_dim(S_t).

Every spec carries the per-term precision label multinomial_term_{t} so the outer loop ties one shared λ_t across all classes (the gauge the centered metric requires; an untied per-(class,term) λ is itself a second source of reference dependence).

pub fn joint_smoothing_dimension(&self) -> usize

Build the permutation-EQUIVARIANT heterogeneous smoothing penalties: for each smooth term t, K per-class penalties λ_{t,c} · γ_cᵀ S_t γ_c on the CENTERED class functions γ_c = β_c − (1/K)Σ_b β_b (with β_ref ≡ 0), one λ per class — including the softmax reference class.

This is the resolution of the #1587 (reference invariance) vs #1855 (heterogeneous per-class smoothness) tension. The reverted per-block carrier penalized the ALR contrasts β_a = γ_a − γ_ref, whose “per-class” smoothness is an artifact of which class is the baseline (the family of diagonal ALR precisions is not closed under reference changes). Penalizing the centered functions is reference-free by construction: relabeling classes permutes the (γ_c, λ_{t,c}) pairs together, so the fitted probabilities after label alignment are identical, while REML still selects genuinely heterogeneous per-class smoothness (a wiggly class takes a small λ_c, an easy class shrinks its centered deviation toward the mean function).

In stacked ALR coordinates [β_0; …; β_{m−1}] (m = K−1), class c’s centering row is C_a = e_aᵀ − 𝟙ᵀ/K for an active class and C_ref = −𝟙ᵀ/K for the reference, so spec (t, c) carries the PSD rank-rank(S_t) matrix (C_cᵀC_c) ⊗ S_t. With all λ_{t,c} equal the sum collapses exactly to the shared centered metric: Σ_c C_cᵀC_c = I − J/K = M, so this family strictly generalizes Self::centered_joint_penalty_specs.

K = 2 is the degenerate case: γ_ref = −γ_0, both centered functions have identical wiggliness, and the two per-class metrics are proportional (only λ_0 + λ_1 would be identified). The shared centered spec is the correct model there, so this builder returns it. The number of smoothing coordinates the OUTER search actually has.

This is the length of the joint penalty spec list, and it is NOT (K − 1) · n_penalties. Under the equivariant carrier (#1587) each penalty component emits ONE spec PER CLASS (s = t·K + c), so a K = 3 model carries 3·n_penalties coordinates, not 2·n_penalties; the K ≤ 2 arm emits one shared centered spec per component. Any policy keyed on “how many ρ are there” — the exact-outer-curvature dimension gate, a cost estimate, a box — must read THIS and not the per-block count the pre-#1587 layout had, which is a different number for every K > 2 model.

Computed from the shapes alone, so a caller deciding a policy does not have to materialize n_penalties · K dense (m·p)² matrices to find out how many there will be.

pub fn equivariant_class_penalty_specs( &self, ) -> Result<Vec<JointPenaltySpec>, String>

Trait Implementations

impl Clone for MultinomialFamily

fn clone(&self) -> MultinomialFamily

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl CustomFamily for MultinomialFamily

fn jeffreys_span_aggregate_penalty(&self) -> Result<Option<Array2<f64>>, String>

The directions the multinomial’s smoothing reaches, as one λ-free aggregate (#2612).

Every joint spec is (r_c r_cᵀ) ⊗ S_t for a term t and a class contrast r_c, and Σ_c r_c r_cᵀ is exactly the centered class metric M = I − J/K, which is positive definite. So

    Σ_{t,c} (r_c r_cᵀ) ⊗ S_t  =  M ⊗ Σ_t S_t

and, because M is PD, ker(M ⊗ Σ_t S_t) = ℝ^{K−1} ⊗ ker(Σ_t S_t) — the unpenalized columns of the shared design, replicated across the active classes. Every positive combination Σ λ_{t,c} (r_c r_cᵀ) ⊗ S_t has the SAME kernel, so this aggregate answers “which directions does the smoothing reach” without knowing a single λ; it is therefore constant in both β and ρ, which is what keeps the Jeffreys derivative tower valid unchanged.

Built directly rather than by summing equivariant_class_penalty_specs, which materialises K · T dense (K−1)P-square matrices and is called on paths that run per inner-Newton cycle: this is one P-square sum and one (K−1)-square metric.

fn joint_jeffreys_term_required(&self) -> bool

Whether the family’s inner/outer solves need the full-span Jeffreys curvature H_Φ and score ∇Φ.

fn joint_jeffreys_term_strength(&self) -> f64

Multiplicative strength of the Jeffreys/Firth contribution.

fn jeffreys_span_basis(&self) -> Result<Option<Array2<f64>>, String>

The directions the Jeffreys/Firth term is allowed to act on, stated DIRECTLY as an orthonormal total_p × m basis, when the family has measured them rather than derived them from a penalty’s kernel (gam#2612).

fn coefficient_mode_homotopy_member( &self, progress: f64, ) -> Result<Option<Self>, String>

Return the family member at progress on a canonical homotopy from a uniquely selected coefficient objective (progress = 0) to self (progress = 1).

fn joint_penalty_specs(&self) -> Result<Vec<JointPenaltySpec>, String>

Full-span cross-block smoothing penalties, in raw (pre-canonicalisation) coordinates over the entire stacked parameter vector Σ_b p_b_raw.

fn exact_newton_joint_hessian_beta_dependent(&self) -> bool

Whether the joint likelihood Hessian H_L depends on β.

fn inner_coefficient_objective_is_globally_convex(&self) -> bool

Whether the complete inner coefficient objective is globally convex on its feasible set for every admissible smoothing parameter.

fn pseudo_logdet_mode(&self) -> PseudoLogdetMode

How the penalized Hessian’s log-determinant and its derivatives should handle eigenvalues below the numerical-stability floor.

fn has_explicit_joint_hessian(&self) -> bool

Whether the family has an explicit override of exact_newton_joint_hessian (or its _with_specs variant) that returns the true coupled joint Hessian rather than the trait’s block-diagonal default.

fn requires_joint_outer_hyper_path(&self) -> bool

Whether outer hyper-derivative evaluation must use a joint exact path.

fn levenberg_on_ill_conditioning(&self) -> bool

Whether the coupled-joint inner Newton should engage its self-vanishing Levenberg–Marquardt damping μ on a FULL-RANK-but-ILL-CONDITIONED penalized Hessian (cond > COND_NEWTON_SAFETY), not only on a rank-deficient one (nullity > 0). Default false (binary / AFT / others byte-identical). Survival marginal-slope overrides to true (#808: full-rank but cond ≈ 5.8e6; the self-vanishing μ shapes only the trajectory, so the converged β is unbiased and the log-slope target is preserved). Survival-local by trait override so the shared spectral-range solver stays byte-identical for every other family — in particular AFT (survival_location_scale), whose intercept-only-scale fits can be high-cond and which a shared (unconditional) gate would regress (#735/#736).

fn inner_coefficient_hessian_hvp_available( &self, specs: &[ParameterBlockSpec], ) -> bool

Explicit name for the inner coefficient-space Hessian HVP capability.

fn inner_joint_workspace_gradient_available( &self, specs: &[ParameterBlockSpec], ) -> bool

fn inner_joint_workspace_log_likelihood_available( &self, specs: &[ParameterBlockSpec], ) -> bool

fn coefficient_hessian_cost(&self, specs: &[ParameterBlockSpec]) -> u64

Per-evaluation arithmetic cost of forming or applying the inner coefficient-space Hessian once, in flop-equivalent units. This is used for diagnostics, seed-budget policy, and first-order iteration caps when a family genuinely lacks analytic second-order support. It is not allowed to hide an analytic Hessian from the outer optimizer.

fn evaluate( &self, block_states: &[ParameterBlockState], ) -> Result<FamilyEvaluation, String>

Evaluate log-likelihood and per-block working quantities at current block predictors.

fn log_likelihood_only( &self, block_states: &[ParameterBlockState], ) -> Result<f64, String>

Compute only the log-likelihood without building working sets.

fn exact_newton_joint_hessian( &self, block_states: &[ParameterBlockState], ) -> Result<Option<Array2<f64>>, String>

Optional exact joint coefficient-space Hessian across all blocks.

fn exact_newton_joint_gradient_evaluation( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<ExactNewtonJointGradientEvaluation>, String>

Optional exact joint log-likelihood / score evaluation in flattened coefficient space without building per-block Hessian working sets.

fn exact_newton_joint_hessian_workspace( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<Arc<dyn ExactNewtonJointHessianWorkspace>>, String>

Optional per-evaluation workspace for exact joint Hessian operators and directional derivatives.

fn exact_newton_joint_hessian_directional_derivative( &self, block_states: &[ParameterBlockState], d_beta_flat: &Array1<f64>, ) -> Result<Option<Array2<f64>>, String>

Optional exact directional derivative of the joint coefficient-space Hessian.

fn joint_jeffreys_information_directional_derivative_all_axes_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<Vec<Array2<f64>>>, String>

BATCHED all-axes FIRST beta-directional derivative of Self::joint_jeffreys_information_with_specs: with the direction sweeping every canonical axis e_a, return the p dense matrices {Hdot[e_a]}_{a=0..p}.

fn joint_jeffreys_information_second_directional_all_axes_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_u_flat: &Array1<f64>, ) -> Result<Option<Vec<Array2<f64>>>, String>

BATCHED all-axes second beta-directional derivative of Self::joint_jeffreys_information_with_specs: with d_beta_u fixed and the second direction sweeping every canonical axis e_a, return the p dense matrices {H²dot[d_beta_u, e_a]}_{a=0..p}.

fn joint_jeffreys_information_contracted_trace_hessian_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], weight: &Array2<f64>, ) -> Result<Option<Array2<f64>>, String>

Contracted second beta-derivative matching Self::joint_jeffreys_information_with_specs:

fn joint_jeffreys_information_contracted_trace_hessian_available(&self) -> bool

Whether Self::joint_jeffreys_information_contracted_trace_hessian_with_specs can supply the wide-p Jeffreys completion without the pairwise H'' fallback. Default false preserves the historical width cap exactly.

fn exact_newton_joint_hessiansecond_directional_derivative( &self, block_states: &[ParameterBlockState], d_beta_u_flat: &Array1<f64>, d_beta_v_flat: &Array1<f64>, ) -> Result<Option<Array2<f64>>, String>

Optional exact second directional derivative of the joint Hessian.

fn outer_derivative_pilot_schedule( &self, ) -> Option<OuterDerivativePilotSchedule>

Optional sampled-derivative pilot owned by this family.

fn persistent_warm_start_fingerprint( &self, specs: &[ParameterBlockSpec], options: &BlockwiseFitOptions, ) -> Option<String>

Family-owned fingerprint for persistent coefficient warm-starts.

fn current_identifiability_family_scalars( &self, block_states: &[ParameterBlockState], ) -> Result<Option<Arc<dyn Any + Send + Sync>>, String>

Build family-owned current-state scalars for a post-fit identifiability audit.

fn identifiability_probit_frailty_scale(&self) -> f64

Scale stored on the converged identifiability linearization state. Families returning dynamic scalars override this when their callbacks use a non-unit probit/frailty scale.

fn log_likelihood_only_with_options( &self, block_states: &[ParameterBlockState], options: &BlockwiseFitOptions, ) -> Result<f64, String>

Options-aware log-likelihood evaluation for line search.

fn supports_log_likelihood_early_exit(&self) -> bool

Whether log_likelihood_only_with_options can use BlockwiseFitOptions::early_exit_threshold to reject line-search trials without computing the full log-likelihood.

fn exact_newton_outerobjective(&self) -> ExactNewtonOuterObjective

Selects the outer objective semantics for exact-Newton families.

fn use_projected_penalty_logdet(&self) -> bool

Whether the outer REML/LAML logdet term ½ log|H + Sλ| and its analytic trace gradient ½ tr((H+Sλ)⁺ ∂Sλ) are evaluated over the FULL identifiable subspace range(H + Sλ) (mgcv’s generalized determinant, gam#752) rather than the penalty-range subspace range(Sλ).

fn coefficient_gradient_cost(&self, specs: &[ParameterBlockSpec]) -> u64

Per-evaluation arithmetic cost of one analytic-gradient outer evaluation, in flop-equivalent units. Used only when the family genuinely has no analytic outer Hessian and the planner must use a first-order optimizer.

fn exact_outer_derivative_order( &self, specs: &[ParameterBlockSpec], options: &BlockwiseFitOptions, ) -> ExactOuterDerivativeOrder

Declares how much exact outer calculus this family wants to expose for the current realized problem size.

fn outer_derivative_policy( &self, specs: &[ParameterBlockSpec], psi_dim: usize, options: &BlockwiseFitOptions, ) -> OuterDerivativePolicy

Realized outer-derivative policy at the current problem size.

fn outer_derivative_subsample_capable(&self) -> bool

Whether this family’s outer-only paths honour HT-weighted partial sums over options.outer_score_subsample.

fn outer_seed_config(&self, n_params: usize) -> SeedConfig

Family-specific outer seeding policy.

fn output_channel_assignment( &self, specs: &[ParameterBlockSpec], ) -> Option<Vec<usize>>

Per-block output-channel assignment for the identifiability audit.

fn block_geometry( &self, block_states: &[ParameterBlockState], spec: &ParameterBlockSpec, ) -> Result<(DesignMatrix, ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>), String>

Optional dynamic geometry hook for blocks whose design/offset depend on current values of other blocks.

fn block_geometry_is_dynamic(&self) -> bool

Whether block_geometry(...) can change with the current block state.

fn block_geometry_directional_derivative( &self, block_states: &[ParameterBlockState], block_index: usize, block_spec: &ParameterBlockSpec, d_beta: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<BlockGeometryDirectionalDerivative>, String>

Optional directional derivative of the effective block geometry wrt the current block coefficients.

fn post_update_block_beta( &self, block_states: &[ParameterBlockState], block_index: usize, block_spec: &ParameterBlockSpec, beta: ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, String>

Optional per-block coefficient projection applied after each block update.

fn max_feasible_step_size( &self, block_states: &[ParameterBlockState], block_index: usize, delta: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<f64>, String>

Optional barrier-aware maximum feasible step size for a block update.

fn max_feasible_joint_step_size( &self, block_states: &[ParameterBlockState], delta: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<f64>, String>

Optional barrier-aware maximum feasible step size for the JOINT update.

fn joint_trust_metric_block_floor( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>>, String>

Optional scale-aware floor for the joint trust-region metric D.

fn block_linear_constraints( &self, block_states: &[ParameterBlockState], block_index: usize, block_spec: &ParameterBlockSpec, ) -> Result<Option<ConstraintSet>, String>

Optional inequality constraints for a block update: `A * beta_block

fn block_coefficient_coordinate( &self, block_states: &[ParameterBlockState], block_index: usize, block_spec: &ParameterBlockSpec, ) -> CoefficientCoordinate

Is block block_index’s COEFFICIENT COORDINATE model content, or is only its column space? (#2748)

fn exact_newton_hessian_directional_derivative( &self, block_states: &[ParameterBlockState], block_index: usize, d_beta: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional exact directional derivative of a block’s ExactNewton Hessian.

fn exact_newton_hessian_second_directional_derivative( &self, block_states: &[ParameterBlockState], block_index: usize, d_beta_u: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_beta_v: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional exact second directional derivative of a block’s ExactNewton Hessian.

fn exact_newton_joint_loglik_gradient( &self, block_states: &[ParameterBlockState], ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>>, String>

Optional block-concatenated log-likelihood gradient g = nabla l(theta) assembled from the SAME single source of truth as Self::exact_newton_joint_hessian (e.g. a per-row jet-tower kernel), so the damped Newton H delta = g is solved on a consistent (objective, gradient, Hessian) triple. The default returns None, leaving the caller on its legacy hand-assembled gradient.

fn exact_newton_joint_hessian_workspace_with_options( &self, states: &[ParameterBlockState], specs: &[ParameterBlockSpec], options: &BlockwiseFitOptions, ) -> Result<Option<Arc<dyn ExactNewtonJointHessianWorkspace>>, String>

Outer-aware variant of exact_newton_joint_hessian_workspace.

fn batched_outer_gradient_terms( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, rho: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, options: &BlockwiseFitOptions, hessian_workspace: Option<Arc<dyn ExactNewtonJointHessianWorkspace>>, ) -> Result<Option<BatchedOuterGradientTerms>, String>

Optional batched analytic-gradient hook.

fn batched_outer_hessian_terms( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, rho: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, hessian_workspace: Option<Arc<dyn ExactNewtonJointHessianWorkspace>>, ) -> Result<Option<BatchedOuterHessianTerms>, String>

Optional batched analytic-Hessian / HVP hook.

fn prefers_matrix_free_inner_joint( &self, specs: &[ParameterBlockSpec], block_states: &[ParameterBlockState], ) -> bool

Opt families in to the matrix-free inner-Newton/PCG path on top of the generic use_joint_matrix_free_path heuristic.

fn outer_hyper_hessian_hvp_available( &self, specs: &[ParameterBlockSpec], ) -> bool

True only when the family has a real profiled outer Hessian-vector product over θ = (ρ, ψ), without enumerating all θ_i θ_j pairs.

fn outer_hyper_hessian_dense_available( &self, specs: &[ParameterBlockSpec], ) -> bool

True when the family can expose the dense profiled outer Hessian. Generic custom-family pairwise derivative paths default to dense availability; families with only inner HVP support should override this if dense θθ assembly is not a valid capability for their path.

fn outer_hyper_hessian_operator( &self, specs: &[ParameterBlockSpec], ) -> Option<Arc<dyn HessianOperator>>

Family-supplied exact outer Hessian operator over θ = (ρ, ψ).

fn exact_newton_joint_hessian_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional spec-aware exact joint Hessian.

fn joint_hessian_is_structurally_coupled( &self, block_states: &[ParameterBlockState], ) -> Result<bool, String>

Structural-coupling probe shared by the _with_specs joint dispatch gates: is the family’s exact_newton_joint_hessian a genuinely coupled matrix (nonzero off-diagonal blocks), as opposed to the trait’s block-diagonal default? This is the marker-free signal that lets the engine trust a coupled multi-block family that overrode the joint Hessian without hand-setting has_explicit_joint_hessian(). Returns false when no joint Hessian is available or it is block-diagonal.

fn likelihood_blocks_uncoupled(&self) -> bool

Whether the family’s log-likelihood Hessian is block-diagonal in the joint coefficient vector — i.e. ∂²L/∂β_a∂β_b = 0 for every pair of distinct blocks a ≠ b. Default false (assume coupling, the safe answer); families whose blocks share no η/W coupling override to true to opt into the default working-set joint-Hessian assembly for multi-block specs.

fn joint_jeffreys_information_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional Tier-B Jeffreys information matrix.

fn joint_jeffreys_information_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

First beta-directional derivative of Self::joint_jeffreys_information_with_specs.

fn joint_jeffreys_information_second_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_u_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_betav_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Second beta-directional derivative of Self::joint_jeffreys_information_with_specs.

fn exact_newton_joint_contracted_trace_hessian( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], weight: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional contracted second beta-derivative of the observed joint Newton information:

fn joint_jeffreys_information_matches_observed_hessian(&self) -> bool

Whether Self::joint_jeffreys_information_with_specs is the SAME object as the observed joint Newton Hessian (exact_newton_joint_hessian_with_specs).

fn joint_jeffreys_information_depends_on_psi(&self) -> bool

Whether Self::joint_jeffreys_information_with_specs depends EXPLICITLY on the ψ hyperparameters — i.e. whether ∂_ψ H_info|_β ≠ 0.

fn inner_objective_is_self_concordant(&self) -> bool

Whether this family’s penalized inner objective −ℓ(β) + ½βᵀS(λ)β is SELF-CONCORDANT in the coefficients, so the coupled-joint inner Newton may globalize with a self-concordant DAMPED step (α = 1/(1+λ_N), λ_N the Newton decrement) instead of the trust-region ratio search.

fn outer_default_trustworthy_for_joint_hessian( &self, specs: &[ParameterBlockSpec], ) -> bool

Internal helper: do the outer-REML _with_specs defaults trust the inner-fit’s block-diagonal-from-blocks output for this family?

fn exact_newton_outer_curvature( &self, block_states: &[ParameterBlockState], ) -> Result<Option<ExactNewtonOuterCurvature>, String>

Optional scale-aware exact joint curvature for the outer REML calculus.

fn exact_newton_outer_curvature_directional_derivative( &self, block_states: &[ParameterBlockState], d_beta_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional first directional derivative matching exact_newton_outer_curvature.

fn exact_newton_outer_curvature_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Spec-aware variant of exact_newton_outer_curvature_directional_derivative.

fn exact_newton_outer_curvature_second_directional_derivative( &self, block_states: &[ParameterBlockState], d_beta_u_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_beta_v_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional second directional derivative matching exact_newton_outer_curvature.

fn exact_newton_outer_curvature_second_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_u_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_beta_v_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Spec-aware variant of exact_newton_outer_curvature_second_directional_derivative.

fn exact_newton_joint_hessian_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional spec-aware exact first directional derivative of the joint Hessian.

fn exact_newton_joint_hessian_second_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_u_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_betav_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional spec-aware exact second directional derivative of the joint Hessian.

fn joint_outer_hyper_surrogate_hessian_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional joint multi-block outer-hyper surrogate Hessian over the flattened coefficient vector.

fn joint_outer_hyper_surrogate_hessian_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional first beta-directional derivative of the joint surrogate outer-hyper Hessian.

fn joint_outer_hyper_surrogate_hessian_second_directional_derivative_with_specs( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], d_beta_u_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_betav_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional second beta-directional derivative of the joint surrogate outer-hyper Hessian.

fn diagonalworking_weights_directional_derivative( &self, block_states: &[ParameterBlockState], block_index: usize, d_eta: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>>, String>

Optional exact directional derivative of diagonal working weights along a predictor-space direction d_eta for BlockWorkingSet::Diagonal.

fn diagonalworking_weights_second_directional_derivative( &self, block_states: &[ParameterBlockState], block_index: usize, d_eta_u: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, d_eta_v: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>>, String>

Optional exact second directional derivative of diagonal working weights.

fn exact_newton_joint_psi_terms( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, psi_index: usize, ) -> Result<Option<ExactNewtonJointPsiTerms>, String>

Optional exact first-order joint psi terms over the flattened coefficient vector.

fn exact_newton_joint_psisecond_order_terms( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, psi_index_i: usize, psi_index_j: usize, ) -> Result<Option<ExactNewtonJointPsiSecondOrderTerms>, String>

Optional exact second-order joint psi terms over the flattened coefficient vector.

fn exact_newton_joint_psi_workspace( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, ) -> Result<Option<Arc<dyn ExactNewtonJointPsiWorkspace>>, String>

Optional per-evaluation workspace for exact joint ψ derivatives.

fn exact_newton_joint_psi_workspace_with_options( &self, states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, options: &BlockwiseFitOptions, ) -> Result<Option<Arc<dyn ExactNewtonJointPsiWorkspace>>, String>

Outer-aware variant of exact_newton_joint_psi_workspace.

fn exact_newton_joint_psi_workspace_for_first_order_terms(&self) -> bool

Whether the family’s exact joint ψ workspace should also be built for first-order ψ terms during outer gradient evaluation.

fn exact_newton_joint_psihessian_directional_derivative( &self, block_states: &[ParameterBlockState], specs: &[ParameterBlockSpec], hyper_layout: &CustomFamilyHyperLayout, psi_index: usize, d_beta_flat: &ArrayBase<OwnedRepr<f64>, Dim<[usize; 1]>>, ) -> Result<Option<ArrayBase<OwnedRepr<f64>, Dim<[usize; 2]>>>, String>

Optional mixed beta/psi Hessian drift D_beta H_psi[u].

impl Debug for MultinomialFamily

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

Formats the value using the given formatter.