Struct MatrixFreeSpdOperator 

pub struct MatrixFreeSpdOperator { /* private fields */ }

Operator-backed SPD Hessian with exact spectral REML algebra.

The operator closure is still useful for construction paths that naturally expose HVPs, but REML cost/gradient/Hessian terms must all come from one exact decomposition so ∂ log|H| = tr(H⁻¹ ∂H) holds. We therefore materialize the coefficient Hessian by canonical-basis HVPs under an explicit memory cap and delegate logdet, traces, and solves to DenseSpectralOperator.

Implementations

impl MatrixFreeSpdOperator

pub fn new_with_mode<F>(dim: usize, apply: F, mode: PseudoLogdetMode) -> Self

where F: Fn(&Array1<f64>) -> Array1<f64> + Send + Sync + 'static,

pub fn new_with_mode_and_dense_assemble<F>( dim: usize, apply: F, mode: PseudoLogdetMode, dense_assemble: Option<Arc<dyn Fn() -> Option<Array2<f64>> + Send + Sync>>, ) -> Self

where F: Fn(&Array1<f64>) -> Array1<f64> + Send + Sync + 'static,

Like [new_with_mode], but additionally accepts an optional single-pass dense assembly of the same penalized operator. When present and it yields a matrix, materialize_dense_operator uses it instead of the dim canonical-basis matvecs. See the field doc on dense_assemble.

Trait Implementations

impl HessianOperator for MatrixFreeSpdOperator

fn prefers_stochastic_trace_estimation(&self) -> bool

The operator delegates logdet, trace_hinv_*, trace_logdet_*, solve, and solve_multi to a lazily-built DenseSpectralOperator whenever the exact-dense materialization fits the configured byte cap (see exact_dense_spectral_budget_ok / EXACT_DENSE_SPECTRAL_MAX_BYTES). In that regime the algebra is exact spectral — there is no stochastic preference to advertise, and forcing the caller to take the Hutchinson path would replace an O(p²) exact reduction with O(k·apply) noisy probes.

When the budget is exceeded the dense factor cannot be built and the CG trace-solve path added in 2bd6af68 is the only feasible route; the flag flips to true so stochastic_trace_solve* callers route through cg_trace_solve instead of crashing in exact_dense_spectral().expect.

fn logdet_traces_match_hinv_kernel(&self) -> bool

Mirror the prefers_stochastic_trace_estimation gate: when the dense factor is reachable the operator’s logdet / trace_hinv reductions all resolve through DenseSpectralOperator, whose logdet_traces_match_hinv_kernel is false for the smooth-spectral regularization variants we run. Reporting true here would let the outer evaluator route logdet-gradient/Hessian traces through the Hutchinson H⁻¹ kernel which does not satisfy ∂ log|H| = tr(H⁻¹ ∂H) under smooth-spectral. The CG-only regime (budget exceeded) lacks a dense reference so falling back to the stochastic kernel is acceptable as a best-effort estimate.

fn logdet(&self) -> f64

log|H|₊ — pseudo-logdet using only active eigenvalues/pivots.

fn as_exact_dense_spectral(&self) -> Option<&DenseSpectralOperator>

Exact dense spectral representation, when this backend has one.

fn trace_hinv_product(&self, a: &Array2<f64>) -> f64

tr(H₊⁻¹ A) — trace of pseudo-inverse times a symmetric matrix. Uses the SAME decomposition as logdet.

fn trace_hinv_operator(&self, op: &dyn HyperOperator) -> f64

tr(H₊⁻¹ B) for an operator-backed Hessian drift.

fn trace_hinv_product_cross(&self, a: &Array2<f64>, b: &Array2<f64>) -> f64

tr(H⁻¹ A H⁻¹ B) for dense symmetric Hessian drifts.

fn trace_hinv_matrix_operator_cross( &self, matrix: &Array2<f64>, op: &dyn HyperOperator, ) -> f64

tr(H⁻¹ A H⁻¹ B) for a dense drift A and an operator-backed drift B.

fn trace_hinv_operator_cross( &self, left: &dyn HyperOperator, right: &dyn HyperOperator, ) -> f64

tr(H⁻¹ A H⁻¹ B) for operator-backed Hessian drifts.

fn trace_logdet_operator(&self, op: &dyn HyperOperator) -> f64

tr(G_ε(H) B) for an operator-backed Hessian drift.

fn solve(&self, rhs: &Array1<f64>) -> Array1<f64>

H⁻¹ v — linear solve using the active decomposition.

fn solve_multi(&self, rhs: &Array2<f64>) -> Array2<f64>

H⁻¹ M — multi-column solve.

fn stochastic_trace_solve(&self, rhs: &Array1<f64>, rel_tol: f64) -> Array1<f64>

H⁻¹ v for stochastic trace probes.

fn stochastic_trace_solve_for_probe( &self, rhs: &Array1<f64>, rel_tol: f64, probe_id: u64, trace_state: Option<&Arc<Mutex<StochasticTraceState>>>, ) -> Array1<f64>

H⁻¹ v for a deterministic stochastic trace probe id.

fn stochastic_trace_solve_multi( &self, rhs: &Array2<f64>, rel_tol: f64, ) -> Array2<f64>

H⁻¹ M for stochastic trace probes.

fn trace_logdet_hessian_cross( &self, h_i: &Array2<f64>, h_j: &Array2<f64>, ) -> f64

Cross-trace for the logdet Hessian: ∂²_{ij} log|R_ε(H)| = tr(G_ε Ḧ_{ij}) + spectral_cross(Ḣ_i, Ḣ_j).

fn trace_logdet_hessian_cross_matrix_operator( &self, h_i: &Array2<f64>, h_j: &dyn HyperOperator, ) -> f64

Operator-backed mixed form of [trace_logdet_hessian_cross].

fn trace_logdet_hessian_cross_operator( &self, h_i: &dyn HyperOperator, h_j: &dyn HyperOperator, ) -> f64

Operator-backed form of [trace_logdet_hessian_cross].

fn active_rank(&self) -> usize

Number of active dimensions (rank of pseudo-inverse).

fn dim(&self) -> usize

Full dimension of H.

fn is_dense(&self) -> bool

Whether this operator is backed by a dense factorization.

fn as_dense_spectral(&self) -> Option<&DenseSpectralOperator>

Access the dense spectral backend when this operator is powered by a single eigendecomposition.

fn has_matrix_free_trace_cg_operator(&self) -> bool

Whether this backend exposes a matrix-free operator usable by trace CG.

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

Assemble the raw dense Hessian represented by this backend for active-constraint tangent projection.

fn trace_logdet_gradient(&self, a: &Array2<f64>) -> f64

tr(G_ε(H) A) — trace for the logdet gradient ∂_i log|R_ε(H)|.

fn xt_logdet_kernel_x_diagonal(&self, x: &DesignMatrix) -> Array1<f64>

diag(X · G_ε(H) · Xᵀ) — the leverage corresponding to trace_logdet_gradient. trace_logdet_gradient(Xᵀ diag(w) X) = Σᵢ wᵢ · h^G[i].

fn trace_logdet_h_k( &self, a_k: &Array2<f64>, third_deriv_correction: Option<&Array2<f64>>, ) -> f64

Efficient computation of tr(G_ε(H) Hₖ) for the logdet gradient.

fn trace_logdet_block_local( &self, block: &Array2<f64>, scale: f64, start: usize, end: usize, ) -> f64

tr(G_ε(H) · A_block) where A_block is a p_block × p_block matrix embedded at rows/columns [start..end].