Native sparse FPCA/PACE validation¶

This page records validation strategy for the 0.10 development branch. It does not claim that the native sparse estimator has been integrated into the 0.9.1 maintenance line.

Evidence is deliberately separated¶

The sparse estimator uses four distinct evidence classes:

Evidence	Question answered	Current use
Analytical truth	Is a numerical/statistical primitive implemented correctly?	weighted covariance eigensystem
Simulation recovery	Can known latent structure be recovered under finite sparse sampling?	full native sparse FPCA/PACE
Independent implementation equivalence	Does the same mathematical specification agree with an independent implementation?	one-dimensional local-linear pooled mean smoother
Cross-implementation sensitivity	Do related implementations remain scientifically coherent when their contracts cannot be made identical?	whole-estimator FDApy comparison

The last category is intentionally not called equivalence.

Analytical weighted-operator truth¶

The covariance eigensolver is tested using a known rank-two covariance operator on a quadrature grid. Qualification checks generating eigenvalues, sign-invariant weighted component alignment, and weighted eigenfunction orthonormality.

The operator being diagonalized is

\[ W^{1/2}\widehat G W^{1/2}, \]

not the unweighted grid covariance matrix.

Full-estimator simulation recovery¶

The private 0.10 truth generator produces

\[ X_i(t) = \mu(t) + \sum_{k=1}^2 \xi_{ik}\phi_k(t) + \epsilon_i(t) \]

with known mean, eigenfunctions, eigenvalues, latent scores, measurement noise, curve-specific observation times, and sample counts.

Recovery is evaluated with mean and covariance integrated squared error, relative eigenvalue error, absolute weighted component similarity, principal cosines for the retained functional subspace, sign-aligned PACE score correlation/RMSE, and score-system failure rate.

Near-tied eigenvalues are evaluated primarily as a subspace problem rather than by requiring arbitrary individual eigenvectors to remain fixed.

Difficult-regime stress matrix¶

The dedicated sparse-native-validation workflow runs deterministic repeated scenarios covering:

2–4 samples per curve;
variable 5–10 and richer sample regimes;
low and high measurement noise;
well-separated and near-tied eigenvalues;
uniform and strongly center-clustered observation times;
boundary-poor sampling with an explicitly restricted analysis interval;
diagonal-difference noise estimation on an explicit interior support;
deliberately narrow bandwidths used as a fail-closed probe.

Each replicate records whether fitting succeeded or failed closed, the structured failure code, sample-count scale, covariance-pair count, grid size, and recovery metrics when available. Failure frequency is treated as part of the estimator's scientific behavior rather than discarded by conditioning only on successful fits.

The JSON artifact is generated by:

python scripts/run_sparse_fpca_stress_validation.py \
  --replicates 3 \
  --output sparse-fpca-stress-validation.json

Predeclared thresholds in that script are qualification guards for the seeded CI scenarios; they are not universal statistical guarantees.

Measurement-noise diagonal qualification¶

The diagonal-difference measurement-noise estimator is qualified separately from population-subspace and PACE-score recovery.

For the PACE-style measurement-error decomposition, the latent covariance diagonal is not taken directly from the diagonal of the generic local-linear covariance surface. Off-diagonal covariance pairs are rotated by 45 degrees and a diagonal-specific local polynomial is fitted: locally linear along the diagonal direction and quadratic perpendicular to it. The adjusted latent diagonal uses the declared covariance bandwidth. The observed/raw diagonal is smoothed separately using the declared noise_bandwidth.

The resulting signed estimate is averaged only over the analyst-declared noise_support. A non-positive estimate remains a structured noise_variance_invalid failure; it is not silently clipped to zero or a small positive constant.

A dedicated recovery matrix holds the generating process fixed where possible and varies mean, covariance and raw-diagonal bandwidths, curve count, per-curve sample count, noise magnitude, and averaging support. It records both the historical generic-surface diagonal estimate and the adjusted PACE-diagonal estimate so that numerical changes remain auditable.

The evidence is generated by:

python scripts/run_sparse_noise_variance_validation.py \
  --replicates 5 \
  --output-json sparse-noise-variance-validation.json \
  --output-csv sparse-noise-variance-validation.csv

The qualification result is deliberately interpreted as a sensitivity surface, not as a universal bandwidth recommendation. In particular, measurement-noise variance can be difficult to recover when it is small relative to latent functional variance because it is obtained as the difference between two separately smoothed diagonal quantities. Failure frequency and signed estimates are therefore retained rather than conditioning only on positive estimates.

Evaluation-grid convergence¶

A separate estimator-level test refits the same sparse dataset at \(G=21,41,81\) grid points and compares relative eigenvalue change, functional- subspace principal cosines, and PACE score correlations. This targets stabilization under grid refinement; it does not designate one grid size as scientifically correct.

Ridge sensitivity¶

score_ridge is part of the estimator:

\[ \widehat\Sigma_i(\gamma) = \widehat G(T_i,T_i) + \widehat\sigma_\epsilon^2 I + \gamma I. \]

Validation explicitly demonstrates that increasing \(\gamma\) can reduce the condition number while changing and shrinking the returned conditional scores. It is therefore never described as invisible numerical housekeeping.

FDApy: what can be matched¶

FDApy remains an external reference/compatibility implementation. With FDApy 1.0.x, its LocalPolynomial implementation can be configured with the same one-dimensional observations, evaluation points, Epanechnikov kernel, local-linear degree, and numerical bandwidth.

That pooled-mean primitive is therefore tested with tight pointwise agreement and is recorded as independent implementation equivalence for that primitive only.

FDApy: why the full estimator is sensitivity, not equivalence¶

The whole sparse pipeline cannot currently be made identical without changing one of the implementations. In FDApy 1.0.x:

multidimensional local-polynomial kernel weights use Euclidean distance in the \((s,t)\) plane; native eyetrajectoriespy covariance smoothing uses a product Epanechnikov kernel;
FDApy has its own covariance-diagonal and curvewise difference-based noise machinery, and its PACE transform applies a tolerance floor to the noise variance;
the irregular PACE path smooths/interpolates each sparse curve to a dense representation before applying the dense PACE calculation;
after UFPCA fitting, the covariance exposed to PACE is reconstructed from the retained eigenvalues/eigenfunctions;
native eyetrajectoriespy evaluates fitted population objects at native observation times and uses the full fitted covariance surface in the conditional score system.

The full comparison therefore uses broad invariant guards: finite relative eigenvalue differences, absolute weighted component similarity, principal cosines, and sign-aligned score correlation/RMSE. Agreement is interpreted as cross-implementation sensitivity evidence, not proof of numerical equivalence.

Cross-language comparison¶

The current CRAN fdapace implementation was inspected as the most relevant R reference candidate. Its FPCA() interface can accept user-supplied mean, covariance, and measurement-error variance, which is useful for controlled cross-language sensitivity work. However, its current eigenanalysis diagonalizes the regular-grid covariance matrix and applies its own grid-step normalization, whereas eyetrajectoriespy solves the trapezoidal quadrature-weighted operator problem

\[ W^{1/2}\widehat G W^{1/2}. \]

The R implementation also constructs a fitted covariance from its retained positive eigensystem before conditional-expectation scoring. Those differences mean there is currently no defensible exact R equivalence case for the full native estimator.

A future R comparison may still be added as explicitly labeled cross-implementation sensitivity, especially by supplying the same fitted population mean/covariance/noise inputs, but it must not be promoted to equivalence unless the discretized operator and PACE covariance contracts can also be matched.

Sparse-specific performance characterization¶

Sparse covariance work depends on within-curve pair count as well as curve count. The dedicated performance artifact therefore records

\[ \left( n,\; \operatorname{median}N_i,\; P=\sum_i N_i(N_i-1),\; G \right) \longrightarrow \left( \text{runtime}, \text{peak RSS} \right). \]

The current qualification cases separately cover very sparse, moderately sparse, and irregular-rich sampling. The script does not compare package speed against another library:

python scripts/run_sparse_fpca_performance.py \
  --repeats 3 \
  --output sparse-fpca-performance.json

If optimization becomes necessary, acceptable routes include indexing compact kernel neighborhoods, chunking, caching, and vectorization. Randomly discarding covariance pairs is not an acceptable hidden optimization.

Bandwidth selection remains out of scope¶

The 0.10 estimator requires declared numeric bandwidths. CV/GCV is not needed for 0.10 completion because its resampling unit is itself a methodological choice for repeated sparse data. A future selector must explicitly declare and audit its resampling unit, candidate grid, criterion, selected value, and failure behavior.