Struct ClosedFormPenaltyOperator 

pub struct ClosedFormPenaltyOperator { /* private fields */ }

Matrix-free closed-form anisotropic Duchon penalty operator.

Stores the parameters of the closed-form pair-block (q, m, s, κ, η, knot centers, and optional constraint factors). The hot matvec path stays matrix-free; cached_dense is populated only by dense_form().

Implementations

impl ClosedFormPenaltyOperator

pub fn new( centers: ArrayView2<'_, f64>, q: usize, m: usize, s: usize, kappa: f64, aniso_log_scales: Option<&[f64]>, kernel_nullspace: Option<&Array2<f64>>, polynomial_block_cols: usize, outer_identifiability: Option<&Array2<f64>>, ) -> Self

Build an operator with the same closed-form parameters that basis::closed_form_operator_penalty_in_total_basis consumes.

pub fn raw_dim(&self) -> usize

Number of raw kernel rows (K).

pub fn dim(&self) -> usize

Number of columns after applying constraint composition: the dimension that callers see when invoking matvec/dense_form.

pub fn matvec(&self, w: ArrayView1<'_, f64>, out: ArrayViewMut1<'_, f64>)

Evaluate (S w) writing the result into out.

With constraints composed, S' = T^T diag(Z, I_poly)^T S_raw diag(Z, I_poly) T. We apply the chain right-to-left:

1. u = T w (dim → total_pre)
2. u_kernel = u[..kernel_cols]; u_poly = u[kernel_cols..]
3. v = Z u_kernel (kernel_cols → K)
4. y = S_raw v (K → K), via on-the-fly pair-block evaluation
5. y_kernel = Z^T y (K → kernel_cols)
6. compose with zero polynomial block, then out = T^T [y_kernel; 0]

pub fn diag(&self) -> Array1<f64>

Diagonal S[i,i] for i in 0..dim. In the raw layout this is the analytic self-pair repeated K times. With constraint composition the diagonal is not the K-space diagonal; we extract it via e_i^T S' e_i = matvec(e_i)[i].

pub fn trace(&self) -> f64

Trace tr(S'). In raw layout this is K times the analytic self-pair; otherwise it uses the composed-basis diagonal.

pub fn log_det_plus_lambda_i(&self, lambda: f64) -> Result<f64, String>

Exact log det(S' + λI). S' is rank-deficient under typical constraints (kernel/polynomial nullspace), so the regularization λ > 0 is mandatory.

pub fn dense_form(&self) -> Array2<f64>

Materialize the full constrained operator as a dense Array2 for callers that explicitly request a matrix or for validation against closed_form_operator_penalty_in_total_basis. Uses the internal cache: the first call builds, subsequent calls clone from the cache.

Trait Implementations

impl Clone for ClosedFormPenaltyOperator

fn clone(&self) -> Self

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl PenaltyOp for ClosedFormPenaltyOperator

fn dim(&self) -> usize

Side length of the (square) operator.

fn matvec(&self, w: ArrayView1<'_, f64>, out: ArrayViewMut1<'_, f64>)

Apply the operator: out = S w.

fn diag(&self) -> Array1<f64>

Diagonal entries S[i,i].

fn trace(&self) -> f64

Trace tr(S) = Σ_i S[i,i].

fn log_det_plus_lambda_i(&self, lambda: f64) -> Result<f64, String>

Exact log det(S + λI) for λ > 0. S is allowed to be rank-deficient; the regularization makes the regularized operator strictly positive definite.

fn as_dense(&self) -> Array2<f64>

Materialize the operator as a dense matrix. Required for the existing analyze_penalty_block path and for callers that need a &Array2 view (Cholesky factorization, etc.). Implementations that already hold a dense form should return it cheaply.

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

Symmetric eigendecomposition S = V diag(λ) V^T. The default implementation materializes via as_dense and runs the same FaerEigh path the existing dense pipeline uses, which preserves numerical agreement with analyze_penalty_block. Implementations that have a faster path (Lanczos top-k for very large K) may override.