Block-streaming Schur pseudo-likelihood: avoid materializing the full p×p covariance

[microprediction]

Summary

The Schur pseudo-likelihood factorizes into a sum of block-conditional terms, so it never needs the dense p×p covariance. The current SchurCovariance / SchurLikelihood path materializes and damps the full matrix, leaving a memory (and compute) saving on the table — exactly in the high-dimensional regime the method targets.

Current behavior

• precise/schurcov.py: maintains a full p×p EWA covariance, then damps cross-block entries (out[cross] *= gamma, or the geodesic variant) before use.
• SchurLikelihood scores a full cov passed in.

So storage and the relevant inverses are O(p²) / O(p³) regardless of gamma.

Why the full matrix is unnecessary

ℓ_γ = Σ_k log N(x_k; μ_{k|<k}(γ), S_k(γ)), and each block term only touches its own block plus its conditioning set:

• γ = 0 (composite / block-diagonal): only the diagonal blocks are needed — no cross-block entry is ever formed. Memory O(p·b) (block size b), K small independent inverses instead of one p³.
• 0 < γ < 1 with bounded (Vecchia-style) conditioning: if each block conditions on only a fixed neighbourhood of m earlier variables, the full inverse is never formed; memory/compute scale with b+m, not p. This is the same structure as vecchia1988 / katzfuss2021 (already cited in the paper) that makes million-point GPs tractable.
• γ = 1 with full conditioning: the last block conditions on everything → no saving (expected).

So the saving grows smoothly as γ → 0 and as the conditioning bandwidth is capped.

Proposal

A block-streaming estimator/assessor that accumulates only within-block (and optionally bounded-conditioning) sufficient statistics and never allocates the dense p×p matrix:

• new estimator variant (or a flag on SchurCovariance) that keeps per-block EWA scatter + the cross-statistics needed for its conditioning set only;
• optional n_neighbors / bandwidth parameter for Vecchia-style ordered conditioning;
• SchurLikelihood computed block-by-block from a streamed test set.

Acceptance

• Produces ℓ_γ (and covariance_/partial_fit behaviour where applicable) matching the dense path within numerical tolerance at equal settings.
• Peak memory O(p·(b+m)) rather than O(p²); verified on a high-p case (e.g. a few hundred to a few thousand variables).

Motivation / context

Surfaced while building large high-dimensional experiments (spatial GP fields, ~hundreds-of-points weather forecast-residual panels, large asset universes) where p×p storage/inversion dominates. The math already supports streaming; only the implementation materializes the matrix.