Struct ImplicitDesignPsiDerivative 

pub struct ImplicitDesignPsiDerivative { /* private fields */ }

Implicit representation of ∂X/∂ψ_d that supports matrix-vector products without materializing the full (n x p) derivative matrices.

For anisotropic Matern / Duchon terms with D axes, the dense path creates D matrices of size (n x p_smooth) for dX/dpsi_d. At n=400K, p=2000, D=16, that is ~100 GB.

Two storage modes:

Materialized (small-to-medium problems): stores pre-computed arrays

• phi_values[i*n_knots + j] = phi(r_{ij})
• q_values[i*n_knots + j] = phi’(r_{ij}) / r_{ij}
• t_values[i*n_knots + j] = (phi’’(r_{ij}) - q_{ij}) / r_{ij}^2
• axis_components[i*n_knots + j, d] = exp(2 eta_d) * (x_{id} - c_{jd})^2 Memory: O(n * k * (D + 2)).

Streaming (large scale): stores only data/centers/eta/kernel params and recomputes (q, t, s_a) on the fly during each matvec. Memory: O(nd + kd) – no per-(data,knot) storage.

The raw-psi chain rule: shape_a = q * s_a shape_ab = t * s_a * s_b + 2 q s_a 1[a=b] dphi/dpsi_a = shape_a + c * phi d2phi/(dpsi_a dpsi_b) = shape_ab + c (shape_a + shape_b) + c^2 phi where c = 0 for Matérn and c = delta / d for hybrid Duchon.

Implementations

impl ImplicitDesignPsiDerivative

pub fn new( phi_values: Array1<f64>, q_values: Array1<f64>, t_values: Array1<f64>, axis_components: Array2<f64>, ident_transform: Option<Array2<f64>>, full_ident_transform: Option<Array2<f64>>, n: usize, n_knots: usize, n_poly: usize, n_axes: usize, ) -> Self

Construct from pre-computed radial jet scalars.

Arguments

• q_values: (n * n_knots,) — φ’(r)/r for each (data, knot) pair.
• t_values: (n * n_knots,) — (φ’’(r) - q) / r² for each pair.
• axis_components: (n * n_knots, D) — s_{d,ij} = exp(2η_d) · h_d² for each pair/axis.
• ident_transform: optional (n_knots × p_constrained) constraint projection.
• full_ident_transform: optional further projection after padding.
• n, n_knots, n_poly, n_axes: dimensions. Construct from pre-computed (materialized) radial jet scalars. This is the original path for small-to-medium problems where O(nk(d+2)) storage is acceptable.

pub fn new_streaming( data: Arc<Array2<f64>>, centers: Arc<Array2<f64>>, eta: Vec<f64>, radial_kind: RadialScalarKind, ident_transform: Option<Array2<f64>>, full_ident_transform: Option<Array2<f64>>, n_poly: usize, ) -> Self

Construct a streaming operator that recomputes (q, t, s_a) on the fly from raw data/centers/eta during each matvec. No O(n*k) arrays are stored. This is the large-scale path.

pub like the sibling new_* constructors: after the engine crate carve (#1521) the REML planner tests live in gam-solve and build streaming operators as fixtures, so this constructor is part of the cross-crate surface, not a crate-private helper.

pub fn n_data(&self) -> usize

Number of data points.

pub fn n_axes(&self) -> usize

Number of axes (D).

pub fn is_duchon_family(&self) -> bool

pub fn p_out(&self) -> usize

Output dimension: total basis columns in the final space.

pub fn append_full_transform( self, transform: &Array2<f64>, ) -> Result<Self, BasisError>

pub fn unproject_matrix(&self, u: &ArrayView2<'_, f64>) -> Array2<f64>

Batched unproject for a (p_out × rank) coefficient matrix. Returns (n_knots × rank) via two BLAS3 matmuls — the same algebra as unproject, but amortized across all rank columns of u. Used by forward_mul_matrix so per-axis trace evaluations can be a single chunked GEMM rather than rank-many forward_mul calls.

pub fn transpose_mul( &self, axis: usize, v: &ArrayView1<'_, f64>, ) -> Result<Array1<f64>, BasisError>

Compute (∂X/∂ψ_d)^T v for a given axis d and vector v of length n.

Returns a vector of length p_out (total basis dimension after all transforms).

Formula in raw knot space: [raw]j = Σ_i v_i · q{ij} · s_{d,ij} then project through Z and pad.

Note: q = φ_r/r and s_d = exp(2ψ_d)·h_d² are UNNORMALIZED axis components. With this convention, q·s_d = (φ_r/r)·(exp(2ψ_d)·h_d²) = φ_r·(s_d/r), which equals the correct ∂φ/∂ψ_d = φ_r·∂r/∂ψ_d = φ_r·s_d/r. No r² correction is needed — that would be required only if s_d were the fractional quantity s_d/r².

pub fn forward_mul( &self, axis: usize, u: &ArrayView1<'_, f64>, ) -> Result<Array1<f64>, BasisError>

Compute (∂X/∂ψ_d) u for a given axis d and vector u of length p_out.

Returns a vector of length n.

Formula: for each data point i, result_i = Σ_j q_{ij} · s_{d,ij} · u_knot_j where u_knot = Z · u_smooth (unprojected back to knot space).

pub fn transpose_mul_second_diag( &self, axis: usize, v: &ArrayView1<'_, f64>, ) -> Result<Array1<f64>, BasisError>

Compute (∂²X/∂ψ_d²)^T v — diagonal second derivative, same axis.

Matrix-free variant of materialize_second_diag: avoids forming the full (n × p_out) matrix when only a single adjoint matvec is needed.

pub fn transpose_mul_second_cross( &self, axis_d: usize, axis_e: usize, v: &ArrayView1<'_, f64>, ) -> Result<Array1<f64>, BasisError>

Compute (∂²X/∂ψ_d∂ψ_e)^T v — cross second derivative (d ≠ e).

pub fn forward_mul_second_diag( &self, axis: usize, u: &ArrayView1<'_, f64>, ) -> Result<Array1<f64>, BasisError>

Compute (∂²X/∂ψ_d²) u — forward diagonal second derivative.

pub fn forward_mul_second_cross( &self, axis_d: usize, axis_e: usize, u: &ArrayView1<'_, f64>, ) -> Result<Array1<f64>, BasisError>

Compute (∂²X/∂ψ_d∂ψ_e) u — forward cross second derivative.

pub fn materialize_first(&self, axis: usize) -> Result<Array2<f64>, BasisError>

Materialize the full (n × p_out) first-derivative matrix for axis d.

Efficient O(n * k) construction: builds the raw (n × k) kernel derivative matrix directly, then projects through identifiability transforms. This is used when the dense matrix is needed temporarily (e.g., for HyperCoord construction) while avoiding simultaneous storage of all D axes.

pub fn materialize_second_diag( &self, axis: usize, ) -> Result<Array2<f64>, BasisError>

Materialize the full (n × p_out) second diagonal derivative matrix for axis d.

pub fn materialize_second_cross( &self, axis_d: usize, axis_e: usize, ) -> Result<Array2<f64>, BasisError>

Materialize the full (n × p_out) cross second derivative matrix for axes (d, e).

Dense materialization of the t · s_d · s_e cross coupling.

pub fn row_chunk_first( &self, axis: usize, rows: Range<usize>, ) -> Result<Array2<f64>, BasisError>

pub fn row_chunk_first_raw( &self, axis: usize, rows: Range<usize>, ) -> Result<Array2<f64>, BasisError>

Raw (chunk × n_knots) first-order kernel scalars for axis d, without the identifiability/padding projection. Pairs with unproject_matrix in forward_mul_matrix: the kernel scalars stay in raw knot space while the rank side (F) is unprojected to knot space, so the per-chunk GEMM is (chunk × n_knots) · (n_knots × rank) rather than (chunk × p_out) · (p_out × rank). Saves both flops and a (chunk × p_out) intermediate.

pub fn row_chunk_second_diag( &self, axis: usize, rows: Range<usize>, ) -> Result<Array2<f64>, BasisError>

pub fn row_chunk_second_cross( &self, axis_d: usize, axis_e: usize, rows: Range<usize>, ) -> Result<Array2<f64>, BasisError>

pub fn row_vector_first_into( &self, axis: usize, row: usize, out: ArrayViewMut1<'_, f64>, ) -> Result<(), BasisError>

Single-row specialization of row_chunk_first(axis, row..row+1) that writes the length-p_out row directly into the caller-provided buffer.

This is the row-local API used by CustomFamilyPsiLinearMapRef::row_vector for survival rowwise exact-Hessian paths, which previously applied a unit-vector transpose_mul trick (O(n·K) per row) to recover a single row. Avoids allocating a temporary (1 × p_out) matrix per row call.

Trait Implementations

impl Clone for ImplicitDesignPsiDerivative

fn clone(&self) -> ImplicitDesignPsiDerivative

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl Debug for ImplicitDesignPsiDerivative

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

Formats the value using the given formatter.