Enum PenaltyCoordinate 

pub enum PenaltyCoordinate {
    DenseRoot(Array2<f64>),
    DenseRootCentered {
        root: Array2<f64>,
        prior_mean: Array1<f64>,
    },
    BlockRoot {
        root: Array2<f64>,
        start: usize,
        end: usize,
        total_dim: usize,
    },
    BlockRootCentered {
        root: Array2<f64>,
        start: usize,
        end: usize,
        total_dim: usize,
        prior_mean: Array1<f64>,
    },
    KroneckerMarginal {
        eigenvalues: Vec<Array1<f64>>,
        dim_index: usize,
        marginal_dims: Vec<usize>,
        total_dim: usize,
    },
}

A rho-coordinate always contributes

A_k = λ_k S_k, S_k = R_k^T R_k.

For single-block/small problems it is fine to store the full-root R_k in the joint basis. For exact-joint multi-block paths that scaling is wasteful: the root is naturally block-local. This enum lets the unified evaluator consume both forms through one interface.

Variants

DenseRoot(Array2<f64>)

DenseRootCentered

Fields

root: Array2<f64>

prior_mean: Array1<f64>

BlockRoot

Fields

root: Array2<f64>

start: usize

end: usize

total_dim: usize

BlockRootCentered

Fields

root: Array2<f64>

start: usize

end: usize

total_dim: usize

prior_mean: Array1<f64>

KroneckerMarginal

Kronecker-factored penalty coordinate for tensor-product smooths.

In the reparameterized (eigenbasis) representation, the penalty I ⊗ ... ⊗ S_k ⊗ ... ⊗ I becomes I ⊗ ... ⊗ Λ_k ⊗ ... ⊗ I where Λ_k = diag(μ_{k,0}, ..., μ_{k,q_k-1}). This is diagonal in each mode, so apply/quadratic/trace operations avoid O(p²).

Fields

eigenvalues: Vec<Array1<f64>>

Marginal eigenvalues for ALL dimensions: eigenvalues[j] has length q_j.

dim_index: usize

Which marginal dimension this penalty coordinate corresponds to.

marginal_dims: Vec<usize>

Marginal basis dimensions: [q_0, ..., q_{d-1}].

total_dim: usize

Total joint dimension: ∏ q_j.

Implementations

impl PenaltyCoordinate

pub fn from_dense_root(root: Array2<f64>) -> Self

pub fn from_dense_root_with_mean( root: Array2<f64>, prior_mean: Array1<f64>, ) -> Self

pub fn from_block_root( root: Array2<f64>, start: usize, end: usize, total_dim: usize, ) -> Self

pub fn from_block_root_with_mean( root: Array2<f64>, start: usize, end: usize, total_dim: usize, prior_mean: Array1<f64>, ) -> Self

pub fn rank(&self) -> usize

pub fn dim(&self) -> usize

pub fn uses_operator_fast_path(&self) -> bool

pub fn project_into_subspace(&self, z: &Array2<f64>) -> Self

Restrict this penalty coordinate onto the free subspace spanned by the orthonormal columns of z (shape p × m, m ≤ p, zᵀz = I).

When a linear-inequality active set is non-empty, the inner solve and the penalized Hessian are reduced to the free subspace β = z β_f of dimension m = p − active_set_size. The penalty must move in lockstep: the quadratic βᵀ S_k β = β_fᵀ (zᵀ S_k z) β_f, and since S_k = R_kᵀ R_k the reduced root is R_k z (shape rank_k × m). For a block-local root R_k acting on β[start..end] the same identity gives reduced dense root R_k · z[start..end, :], so the reduced coordinate is always a (dimension-m) DenseRoot / DenseRootCentered — the block structure does not survive an arbitrary subspace rotation. A centered mean μ_k maps to zᵀ μ_k, the representation of μ_k in the free subspace.

This keeps dim() equal to the reduced beta.len(), which InnerSolutionBuilder::build asserts.

pub fn apply_penalty(&self, beta: &Array1<f64>, scale: f64) -> Array1<f64>

pub fn apply_penalty_view_into( &self, beta: ArrayView1<'_, f64>, scale: f64, out: ArrayViewMut1<'_, f64>, )

pub fn scaled_add_penalty_view( &self, beta: ArrayView1<'_, f64>, scale: f64, out: ArrayViewMut1<'_, f64>, )

pub fn quadratic(&self, beta: &Array1<f64>, scale: f64) -> f64

pub fn apply_shifted_penalty( &self, beta: &Array1<f64>, scale: f64, ) -> Array1<f64>

pub fn shifted_quadratic(&self, beta: &Array1<f64>, scale: f64) -> f64

pub fn scaled_dense_matrix(&self, scale: f64) -> Array2<f64>

pub fn scaled_block_local(&self, scale: f64) -> (Array2<f64>, usize, usize)

Returns the block-local scaled penalty matrix (p_block × p_block) along with the embedding range, WITHOUT materializing into total_dim × total_dim. For DenseRoot (full-rank, no block structure), returns (matrix, 0, p).

pub fn is_block_local(&self) -> bool

Whether this coordinate has block structure (not full-rank dense).

pub fn scaled_matvec(&self, v: &Array1<f64>, scale: f64) -> Array1<f64>

Apply λ_k S_k to a vector v without materializing the full matrix. For BlockRoot: extracts v[start..end], multiplies by local S_k, embeds result.

pub fn canonical_structural_key(&self) -> u64

A stable, formula-order-independent signature of this penalty coordinate’s STRUCTURAL CONTENT.

Two penalty coordinates that represent the same smoothing structure — the same wiggliness root, the same null-space ridge, the same tensor margin — produce the same key regardless of which block of the joint coefficient vector they happen to occupy or which order the user typed the terms in. It is derived ENTIRELY from rotation/placement-invariant content (rank, block width, the spectrum of the block-local penalty Sₖ = RₖᵀRₖ, or the marginal eigenvalue spectrum for a Kronecker margin), and NEVER from a coordinate’s position (start/dim_index) in the joint layout. Swapping s(x)+s(z) ↔ s(z)+s(x) or te(x,z) ↔ te(z,x) permutes the coordinates but leaves each coordinate’s key fixed.

This is the key the outer REML driver sorts on to present an identical canonical coordinate layout to the smoothing-parameter optimizer regardless of term/margin order, so the flat double-penalty REML valley is resolved order-invariantly (#1538/#1539). Values are quantized to a coarse relative grid so that floating-point round-off in the roots does not split an otherwise-identical key.

Trait Implementations

impl Clone for PenaltyCoordinate

fn clone(&self) -> PenaltyCoordinate

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl Debug for PenaltyCoordinate

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

Formats the value using the given formatter.