Struct ArrowFactorCache 

pub struct ArrowFactorCache {

    pub htt_factors: ArrowFactorSlab,
    pub htt_factors_undamped: ArrowUndampedFactors,
    pub schur_factor: Option<Array2<f64>>,
    pub schur_factor_is_undamped: bool,
    pub beta_schur_deflation: Option<BetaSchurDeflationSpectrum>,
    pub joint_hessian_log_det: Option<f64>,
    pub solver_mode: ArrowSolverMode,
    pub ridge_t: f64,
    pub ridge_beta: f64,
    pub htbeta: ArrowHtbetaCache,
    pub d: usize,
    pub row_dims: Arc<[usize]>,
    pub row_offsets: Arc<[usize]>,
    pub k: usize,
    pub manifold_mode_fingerprint: u64,
    pub row_hessian_fingerprint: u64,
    pub pcg_diagnostics: ArrowPcgDiagnostics,
    pub gauge_deflated_directions: usize,
    pub deflated_row_directions: Arc<[Vec<Array1<f64>>]>,
    pub deflation_row_spectra: Arc<[Option<RowDeflationSpectrum>]>,
    pub beta_gauge_quotient: Option<ArrowBetaGaugeQuotient>,
}

Fields

htt_factors: ArrowFactorSlab

Per-row lower-triangular Cholesky factors of H_tt^(i) + ridge_t·I.

These are the damped factors used inside the Newton solve. The IFT predictor must NOT use them — see Self::htt_factors_undamped.

htt_factors_undamped: ArrowUndampedFactors

Per-row lower-triangular Cholesky factors of the UNDAMPED H_tt^(i) (no ridge_t added).

The IFT predictor formula Δt_i = -(H_tt^(i))⁻¹ · (H_tβ^(i) Δβ + δg_t^(i)) is derived from ∂g_t/∂t = H_tt at the stationary point, with no LM damping term. Reusing the damped factors would bias the predicted shift toward zero in proportion to ridge_t. We pay one extra O(N d³) Cholesky per Newton solve — the same complexity class as the Newton solve itself — to make the IFT exact.

schur_factor: Option<Array2<f64>>

Lower-triangular Cholesky factor of the Schur complement when the selected BA mode formed/factored dense RCS. None for ArrowSolverMode::InexactPCG, where Agarwal-style inexact LM avoids the dense K × K factor.

schur_factor_is_undamped: bool

True iff schur_factor is the reduced Schur complement built from the undamped evidence row factors (H_tt, no LM ridge) and ridge_beta = 0.

A Newton solve may be damped while the cache still carries an undamped evidence Schur for logdet / selected-inverse consumers. When this is false, consumers must not combine schur_factor with Self::undamped_factor: that would mix two different bordered-arrow operators.

beta_schur_deflation: Option<BetaSchurDeflationSpectrum>

Authoritative original-coordinate spectrum and null mask used when the undamped evidence β-Schur was unit-deflated. The mask, rather than a threshold re-derived from L Lᵀ, defines which directions contribute log 1 = 0 to the value and zero to every inverse/trace contraction.

joint_hessian_log_det: Option<f64>

Exact undamped joint-Hessian log-determinant produced by the dense factorization path. REML evidence consumes this directly so the Laplace normalizer cannot miss the log-det even when later cache consumers only need solves/traces.

On the matrix-free large-k SAE evidence path this is set from the Stochastic Lanczos Quadrature reduced-Schur log-determinant (see Self::undamped_arrow_log_det_with_schur and crate::arrow_schur::slq_logdet) so no dense k × k Cholesky is ever formed; arrow_log_det_from_cache reads THIS field first, before any schur_factor diagonal fallback.

solver_mode: ArrowSolverMode

BA mode used to create this cache.

ridge_t: f64

Ridge values used to build the cached factors (recorded so the warm-start predictor knows whether the cache is still valid for a requested ridge level).

ridge_beta: f64

htbeta: ArrowHtbetaCache

Per-row cross-block access for H_tβ^(i) x.

Large caches retain a row matvec callback or disable β-coupled IFT prediction instead of cloning every dense d × K slab.

d: usize

Maximum per-row latent dim (upper bound; matches sys.d at creation).

row_dims: Arc<[usize]>

Per-row latent dims: row_dims[i] is the active dim for row i.

row_offsets: Arc<[usize]>

Flat-buffer row offsets for delta_t / IFT output vectors. row_offsets[i] is the start of row i; row_offsets[n] is the total length.

k: usize

β dimensionality K.

manifold_mode_fingerprint: u64

Geometry tag for the row-local factors and cross-blocks.

row_hessian_fingerprint: u64

Row-system tag for the cached per-row factors, cross-blocks, and shared-block diagonal used to build the Schur factor.

pcg_diagnostics: ArrowPcgDiagnostics

PCG instrumentation from the solve that produced this cache.

Zero-valued (default) when the selected mode did not use PCG (i.e. Direct or SqrtBA).

gauge_deflated_directions: usize

Number of row-local gauge directions stiffened in an undamped evidence factorization.

Each direction is stiffened at UNIT stiffness kappa = 1.0, so it contributes log(1) = 0 to the row-block logdet through the returned Cholesky factor: the gauge orbit is a criterion null direction and adds nothing to the Laplace normalizer (the quotient pseudo-determinant convention, cf. PenaltyPseudologdet). Zero theta/rho dependence.

deflated_row_directions: Arc<[Vec<Array1<f64>>]>

Per-row unit-norm directions vᵢ (in each row’s d-dim latent block coordinates) that an undamped evidence factorization stiffened to UNIT stiffness λ̃ = 1 (gauge or spectral deflation). Indexed by row; empty for every PD row factored without deflation, and empty overall on the non-deflating solver paths (streaming / cross-row-penalty CG / device).

A deflated direction contributes log(1) = 0 to the row-block log-det and is ρ/θ-INDEPENDENT, so its true contribution to ∂log|H|/∂ρ is 0. The analytic outer-gradient traces (assignment_log_strength_hessian_trace and logdet_theta_adjoint) contract ∂H_raw/∂ρ (the RAW, pre-deflation block derivative) against the DEFLATED inverse, which assigns 1/λ̃ = 1 to each vᵢ and therefore spuriously adds ½ vᵢᵀ (∂H_raw/∂ρ) vᵢ. Those traces subtract this per-row term (kept-subspace restriction) using these directions; without them the REML outer ρ-gradient is biased by +Σ_deflated ½ vᵢᵀ ∂H_raw/∂ρ vᵢ.

deflation_row_spectra: Arc<[Option<RowDeflationSpectrum>]>

Per-row RAW spectral decomposition of an undamped evidence H_tt block that underwent SPECTRAL deflation, surfaced so the outer ρ/θ-gradient traces can apply the EXACT deflation-map (Daleckii–Krein) derivative correction, not just the within-row kept-subspace term.

The criterion VALUE re-deflates H_tt at every ρ, so its gradient is tr(H_deflated⁻¹ DΦ[∂H_raw/∂ρ]), where Φ is the spectral pin-to-unit map. By Daleckii–Krein DΦ[Ȧ] = U (F ∘ UᵀȦU) Uᵀ with the divided- difference matrix F_{ml} = (λ̃ₘ − λ̃ₗ)/(λₘ − λₗ) (raw λ in the denominator, conditioned λ̃ in the numerator). The kept×kept block of F is 1 (the kept subspace contracts the raw derivative unchanged), the deflated×deflated block is 0, and the kept(m)×deflated(i) block is (λₘ − 1)/(λₘ − λᵢ) — this last, ROTATION, term is what the per-row kept-subspace correction alone misses; it couples to the β-block through the Schur back-substitution carried in (H⁻¹)_tt.

Some(spectrum) only for spectrally-deflated rows; None for PD rows, gauge-only deflation (ρ-independent structural null — within-row term suffices), and every non-SAE-evidence solver path (streaming / device / cross-row CG). Empty overall when no row deflated spectrally.

beta_gauge_quotient: Option<ArrowBetaGaugeQuotient>

Shared-border scale gauge used by the evidence factor.

When present, schur_factor factors P S P + Q Q^T, and every public inverse primitive projects both its border RHS and result with P. The unit-pinned orbit contributes zero to arrow_log_det and zero to every analytic trace, so value and gradient live on the same quotient.

Implementations

impl ArrowFactorCache

pub fn n_rows(&self) -> usize

pub fn htbeta_available(&self) -> bool

pub fn used_device(&self) -> bool

Whether the Newton solve that produced this cache actually executed on the device: the device-resident Direct dense solve or the device-resident matrix-free SAE PCG (whose matvec runs in CUDA kernels). This does NOT include the injected host-procedural reduced-Schur matvec, whose arithmetic runs on the CPU even when a CUDA context was opened to build per-row factors (#1209) — that path sets ArrowPcgDiagnostics::injected_host_procedural_matvec instead. Read-only routing provenance: lets a fit result record device-vs-CPU as ground truth instead of inferring it from the runtime probe. Mirrors ArrowPcgDiagnostics::used_device_arrow.

pub fn undamped_factor(&self, row: usize) -> ArrayView2<'_, f64>

pub fn undamped_factor_count(&self) -> usize

pub fn undamped_factors_iter( &self, ) -> impl Iterator<Item = ArrayView2<'_, f64>> + '_

pub fn compute_undamped_arrow_log_det(&self) -> Option<f64>

pub fn undamped_arrow_log_det_with_schur( &self, schur_log_det: f64, ) -> Option<f64>

Undamped joint log-determinant log|H| = Σ_i log|H_tt^(i)| + log|S| using an EXTERNALLY-supplied reduced-Schur term schur_log_det = log|S| instead of a dense schur_factor diagonal sum.

This is the matrix-free large-k SAE evidence path: the reduced Schur is never Cholesky-factored, so schur_factor is None and log|S| comes from Stochastic Lanczos Quadrature (crate::arrow_schur::slq_logdet). The per-row latent-block term is computed exactly as in Self::compute_undamped_arrow_log_det, with the same ridge, positivity, and finiteness guards.

pub fn delta_t_len(&self) -> usize

The total length of delta_t / IFT output vectors for this cache.

pub fn apply_htbeta_row( &self, row: usize, delta_beta: ArrayView1<'_, f64>, out: &mut Array1<f64>, ) -> bool

pub fn apply_htbeta_row_transpose( &self, row: usize, v: ArrayView1<'_, f64>, out: &mut Array1<f64>, fallback_op: Option<&RowHtbetaMatvec>, ) -> bool

Accumulate out[a] += H_βt^(row)[a, :] · v for all a in 0..k.

v has length row_dims[row]; out has length k. The caller must zero out before the first call if it needs a fresh result. Returns false when the cache is Disabled and no fallback_op is provided; callers must treat the accumulator as invalid in that case.

pub fn arrow_log_det(&self) -> Option<f64>

Authoritative evidence joint log-determinant for the exact operator this cache exposes to selected-inverse and adjoint consumers.

The factorization path computes this once into Self::joint_hessian_log_det. This accessor deliberately does not reconstruct Σ_i log|H_tt^(i)| + log|Schur_β| from loose pieces: a damped Newton cache can otherwise pair undamped row factors with a damped Schur solve and silently describe no live operator. Returning only the stored joint value keeps REML evidence, fixed-state tests, selected inverse, and logdet_theta_adjoint on the same factorization branch.

pub fn latent_block_inverse_diagonal( &self, ) -> Result<Array1<f64>, ArrowSchurError>

Diagonal of the latent (t-block) of the full bordered-arrow inverse (H⁻¹)_tt, in delta_t layout (length Self::delta_t_len).

For the bordered arrow Hessian H = [[A, B], [Bᵀ, H_ββ]] with A = H_tt (block-diagonal per row, A_i = H_tt^(i)) and B = H_tβ, the standard block-inverse identity gives the t-block (H⁻¹)_tt = A⁻¹ + A⁻¹ B S⁻¹ Bᵀ A⁻¹, where S = H_ββ − Bᵀ A⁻¹ B is the Schur complement on β. Because A is block-diagonal, the (i, j) diagonal entry of (H⁻¹)_tt is computed purely from row i’s factor and cross-block:

a    = A_i⁻¹ e_j                       (chol_solve on the per-row factor)
[A_i⁻¹]_{jj} = a[j]
w    = B_iᵀ a = H_βt^(i) a             (a K-vector)
z    = S⁻¹ w                           (chol_solve on the Schur factor)
diag = a[j] + w · z

The UNDAMPED per-row factors (Self::undamped_factor) are used so the result is the inverse of the true H_tt, not the LM-damped H_tt + ridge_t·I — same rationale the IFT predictor docstring gives at the top of this struct.

Consuming the diagonal as a per-(atom, axis) trace

(H⁻¹)_tt is the latent covariance block. The selected-inverse trace for a contiguous group of latent coordinates (e.g. one atom’s rows, or one axis across rows) is simply the sum of the returned diagonal entries over those row_offsets[i] + j indices — no off-diagonal terms are needed for the trace tr[(H⁻¹)_tt · D] against a diagonal selector D.

Errors

Returns ArrowSchurError::SchurFactorFailed when this cache has no dense Schur factor or no usable H_βt coupling — i.e. it was produced by an ArrowSolverMode::InexactPCG solve (no dense K × K factor) or by a Disabled htbeta cache. The selected-inverse block-trace is not yet supported for the matrix-free PCG mode; that branch needs a separate Lanczos/Hutchinson estimator.

pub fn full_inverse_apply( &self, w_t: ArrayView1<'_, f64>, w_beta: ArrayView1<'_, f64>, ) -> Result<(Array1<f64>, Array1<f64>), ArrowSchurError>

Solve the full bordered-arrow system H·u = w on the cached factor (#1006): w arrives in arrow layout — w_t flat per Self::delta_t_len / row_offsets, w_beta of length K — and the solution comes back in the same layout. Standard block elimination on the SAME factors whose log-determinant the evidence reports:

  y_i      = H_tt^(i)⁻¹ · w_t^(i)
  r_β      = w_β − Σ_i H_βt^(i) · y_i
  u_β      = Schur⁻¹ · r_β
  u_t^(i)  = y_i − H_tt^(i)⁻¹ · (H_tβ^(i) · u_β)

This is the IFT / adjoint back-solve the analytic outer ρ-gradient consumes: u_j = H⁻¹ (∂g/∂ρ_j) per outer coordinate and the H⁻¹-side of the third-order correction −½·Γᵀ·H⁻¹·(∂g/∂ρ_j). Contract: the cache must be the ridge-0 Direct evidence factor (undamped per-row factors + dense Schur), so the solve is against the criterion’s own H — never a damped surrogate (that would desync the gradient from the reported evidence).

pub fn schur_inverse_apply( &self, rhs: ArrayView1<'_, f64>, ) -> Result<Array1<f64>, ArrowSchurError>

Apply the β-block of the full inverse, (H⁻¹)_ββ · rhs = S_β⁻¹ · rhs, where S_β is the Schur complement on β whose Cholesky factor this cache holds in Self::schur_factor.

For the bordered arrow Hessian H = [[A, B], [Bᵀ, H_ββ]], the β-block of H⁻¹ is exactly the inverse of the Schur complement S_β = H_ββ − Bᵀ A⁻¹ B. One Cholesky back-substitution per call, reusing the cached factor; rhs and the returned vector both have length K.

This is the general single-solve primitive for the β border. Callers that need a Schur-inverse trace tr(S_β⁻¹ M) against a structured penalty M (e.g. the SAE λ_smooth Fellner-Schall step, where M = blockdiag_k(λ_k S_k ⊗ I_p)) build it as Σ_col e_colᵀ S_β⁻¹ M e_col — apply this to each column of M (exploiting whatever sparsity M has) and read off result[col]. Keeping M’s layout on the caller side avoids coupling this solver to penalty-op types.

Errors

Returns ArrowSchurError::SchurFactorFailed when this cache has no dense Schur factor (an ArrowSolverMode::InexactPCG solve) — the same not-yet-supported branch as Self::latent_block_inverse_diagonal — or when rhs.len() != k.

pub fn schur_inverse_block( &self, block: Range<usize>, ) -> Result<Array2<f64>, ArrowSchurError>

Dense principal sub-block of the β-block of the full inverse, (H⁻¹)_ββ[block, block] = S_β⁻¹[block, block], shape (W, W) with W = block.len().

For the bordered arrow Hessian H = [[A, B], [Bᵀ, H_ββ]], the β-block of H⁻¹ is exactly S_β⁻¹ (the inverse of the Schur complement whose Cholesky factor this cache holds). This returns the contiguous block × block sub-block — e.g. one SAE atom’s decoder coefficients via gam_terms::sae::manifold::SaeManifoldTerm::beta_block_offsets — by solving S_β x = e_j for each j ∈ block (reusing the cached factor) and gathering the block rows of each solution column. W back-substitutions of size K; the result is symmetrized to clear back-substitution rounding asymmetry. Up to a dispersion scale φ, this block is the joint posterior covariance Cov(β_block) of those coefficients with the latent coordinates already marginalized out (that is precisely what Schur-eliminating the per-row t-blocks does).

Same dense-Schur requirement / error contract as Self::schur_inverse_apply; additionally errors when block runs past K.

pub fn schur_inverse_apply_deflated( &self, rhs: ArrayView1<'_, f64>, ) -> Result<Array1<f64>, ArrowSchurError>

Deflation-aware selected inverse of the cached β-Schur complement — a drop-in for Self::schur_inverse_apply that pseudo-inverts across the numerically-null curvature directions instead of dividing by them.

Why this exists (the λ→0 EDF divergence)

The REML EDF/log-det-trace term contracts (H⁻¹)_ββ against λS. At the ρ lower face a decoder direction can be null in BOTH the data (J_ββ ≈ 0) AND the penalty (s ≈ 0), making S_β = J + λS singular along it. The plain Self::schur_inverse_apply then divides by a ~zero pivot and returns Inf/NaN (the value stays finite — only this H⁻¹-contraction blows up). This method instead forms the spectral pseudo-inverse M⁺ of the SAME operator M = L Lᵀ the plain path inverts, dropping every eigen-direction at or below the solver’s canonical rank floor SPECTRAL_DEFLATION_REL_FLOOR · max|λ| (the exact threshold factor_spectral_deflated_criterion_row and the per-row gauge deflation already use — NOT a new epsilon and NOT a λ-smoothing floor). A doubly-null direction (j ≈ 0 ∧ s ≈ 0) deflates to 0 (it is unidentifiable, not a real DOF); a penalty-only direction survives. The result is finite by construction.

Interior equivalence

Away from the boundary every eigenvalue of M sits orders of magnitude above the floor, so NO direction deflates and M⁺ = M⁻¹ to round-off — this returns the plain selected inverse with no silent bias. Only the λ→0 face deflates. The exact-Newton path keeps calling the plain Self::schur_inverse_apply and is byte-for-byte unchanged.

Errors

Same dense-Schur / undamped-factor / rhs.len() != K contract as Self::schur_inverse_apply, plus a failed symmetric eigendecomposition of the reconstructed M.

pub fn schur_deflated_applier( &self, ) -> Result<impl Fn(ArrayView1<'_, f64>) -> Array1<f64> + '_, ArrowSchurError>

Precompute the deflated spectral pseudo-inverse ONCE and return a reusable applier — the many-RHS form of Self::schur_inverse_apply_deflated. The EDF trace contracts (H⁻¹)_ββ against one λS⊗I column per basis coefficient (Σ_k M_k·r_k columns total); recomputing the O(K³) eigendecomposition per column would multiply that cost by the border width for no reason. Each apply through the returned closure is O(K²) (two dense mat-vecs through the eigenbasis), identical in complexity to the plain Self::schur_inverse_apply back-substitution it replaces.

Same deflation semantics, contract, and errors as Self::schur_inverse_apply_deflated; the closure itself is infallible (rhs length is the caller’s loop invariant — a wrong length panics in the underlying gemv shape check rather than dividing by a null pivot).

pub fn schur_inverse_block_deflated( &self, block: Range<usize>, ) -> Result<Array2<f64>, ArrowSchurError>

Deflation-aware dense principal sub-block of (H⁻¹)_ββ — the drop-in for Self::schur_inverse_block used by the per-atom EDF trace. Identical contract, but each column is solved through the spectral pseudo-inverse (see Self::schur_inverse_apply_deflated) so a boundary atom with a doubly-null decoder direction yields a finite block instead of NaN.

The eigendecomposition of M = L Lᵀ is computed ONCE and reused across all W = block.len() columns.

Trait Implementations

impl Clone for ArrowFactorCache

fn clone(&self) -> ArrowFactorCache

Returns a duplicate of the value.

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

Performs copy-assignment from source.

impl Debug for ArrowFactorCache

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

Formats the value using the given formatter.