Skip to content

New ecorr sqrtsolve#8

Merged
vhaasteren merged 3 commits into
mainfrom
new-ecorr-sqrtsolve
Feb 10, 2026
Merged

New ecorr sqrtsolve#8
vhaasteren merged 3 commits into
mainfrom
new-ecorr-sqrtsolve

Conversation

@vhaasteren
Copy link
Copy Markdown
Member

Summary

  • Added new non-Cholesky sqrtsolve implementation that computes N^{-1/2} X using a closed-form diagonal-plus-rank-1 inverse square-root identity, replacing the previous Cholesky-based approach as the default
  • Preserved the original Cholesky-based whitening under cholsqrtsolve for all three classes (ShermanMorrison, FastShermanMorrison, and the test reference class)
  • The new sqrtsolve is ~2.4x faster than cholsqrtsolve for the whitening step, and computing Z^T N^{-1} Z via sqrtsolve + W^T W is competitive with (or slightly faster than) the direct BLAS-accelerated Sherman-Morrison 2D solve

Details

Mathematical approach: For each ECORR block with covariance D + j 11^T, the new sqrtsolve uses the rank-1 identity

(I + u u^T)^{-1/2} x = x + alpha u (u^T x), alpha = -1 / (sqrt(1+v) (sqrt(1+v) + 1))

where v = ||u||^2. This avoids forming or solving a triangular Cholesky factor entirely. The alpha formula is the numerically stable equivalent of (1/sqrt(1+v) - 1)/v, avoiding cancellation for all values of v.

What changed:

  • python_block_sqrtsolve_rank1 / python_idx_sqrtsolve_rank1: new pure-Python implementations using the closed-form formula
  • cython_block_sqrtsolve_rank1 / cython_idx_sqrtsolve_rank1: new Cython implementations of the same
  • The previous Cholesky-based implementations are renamed to *_cholsqrtsolve_* (all four variants: python block/idx, cython block/idx)
  • ShermanMorrison and FastShermanMorrison classes now expose both sqrtsolve() (non-Cholesky, default) and cholsqrtsolve() (Cholesky-based) public methods
  • _sqrtsolve_D2 on all classes now calls the non-Cholesky path; _cholsqrtsolve_D2 preserves the old behavior

Tests:

  • Existing test_sqrtsolve_D12 updated to test the new non-Cholesky sqrtsolve (element-wise cross-checks between Python reference and Cython, plus end-to-end Z^T N^{-1} Z = (WZ)^T (WZ) verification)
  • New test_cholsqrtsolve_D12 added for the Cholesky-based path, with full cross-validation across all three implementations (Python rank-1, scipy Cholesky, Cython rank-1) and the same end-to-end correctness check
  • Redundant element-wise comparisons between classes that now share the same algorithm are commented out

Benchmark (not included in codebase)

48,000 TOAs (1500 epochs x 32 obs/epoch), 500 parameters:

Operation Time
Direct Z^T N^{-1} Z (BLAS SM 2D) 319 ms
N^{-1/2} Z (sqrtsolve whitening) 193 ms
L^{-1} Z (cholsqrtsolve whitening) 513 ms
sqrtsolve + W^T W 277 ms

Numerical accuracy identical across all paths (~2e-15 relative error).

Test plan

  • All existing tests pass (test_solve_D1, test_solve_1D1, test_solve_2D2, test_sqrtsolve_D12, test_errors, test_logdet)
  • New test_cholsqrtsolve_D12 passes
  • End-to-end correctness: Z^T N^{-1} Z == (WZ)^T (WZ) verified for all code paths (block, indexed, shuffled)****

@codecov-commenter
Copy link
Copy Markdown

codecov-commenter commented Feb 10, 2026

Codecov Report

❌ Patch coverage is 93.10345% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 91.76%. Comparing base (ff28e9e) to head (5e62e78).

Files with missing lines Patch % Lines
fastshermanmorrison/fastshermanmorrison.py 93.10% 4 Missing ⚠️
Additional details and impacted files

Impacted file tree graph

@@            Coverage Diff             @@
##             main       #8      +/-   ##
==========================================
- Coverage   92.42%   91.76%   -0.66%     
==========================================
  Files           2        2              
  Lines         198      255      +57     
==========================================
+ Hits          183      234      +51     
- Misses         15       21       +6     
Flag Coverage Δ
unittests 91.76% <93.10%> (-0.66%) ⬇️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines Coverage Δ
fastshermanmorrison/fastshermanmorrison.py 92.40% <93.10%> (-0.87%) ⬇️

Continue to review full report in Codecov by Sentry.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update ff28e9e...5e62e78. Read the comment docs.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants