Skip to content

perf: optimize every custom periodogram implementation across all backends - #15

Merged
tjayasinghe merged 9 commits into
release/v1.1.0from
feat/perf-optimizations
Jul 2, 2026
Merged

perf: optimize every custom periodogram implementation across all backends#15
tjayasinghe merged 9 commits into
release/v1.1.0from
feat/perf-optimizations

Conversation

@tjayasinghe

Copy link
Copy Markdown
Owner

Implements the full optimization plan across all custom (non-reference) periodogram implementations — numpy, numba, cupy RawKernel, and torch/array-API tiers — plus the extension of the multicore numba CPU tier from BLS to all six ported methods.

Algorithmic (all backends at once)

  • MHAOV: Gram matrix assembled from harmonic trig sums (Chebyshev recurrence, one cos/sin evaluation) instead of the (F,N,2H+1) design tensor — ~3x on numpy/cupy/torch, an order of magnitude less transient memory, gemm-free by construction (Blackwell cuBLAS workaround deleted).
  • GLS direct path: one trig evaluation instead of six (shared base-grid angles + double-angle identities) — ~1.8x on torch CPU/CUDA; frequency chunk auto-capped by N.
  • PDM: fold binned once into n_bins*n_covers fine bins; covers regrouped exactly — ~2.8x numpy, 3x fewer CUDA atomics.

New numba CPU tier (PDM, CE, String-Length, MHAOV, TLS)

prange kernels mirroring the CUDA designs; cpu/auto resolve to them via a shared fast_cpu_backend attribute (replaces the bespoke BLS override). Warm speedups vs the vectorized numpy paths: PDM 306x, CE 135x, MHAOV 57x, TLS 53x, SL 23x; parity <= 2.5e-11 (TLS bit-identical; SL stable-sort tie order preserved via mergesort).

GPU

  • float32-templated RawKernels (BLS/PDM/CE/TLS) honoring precision: ~8.6x for BLS (610->71 ms segmented) and TLS on GeForce-class FP64-limited cards; float64 remains the default; float32 BLS guards empty windows with a sum-scaled ivar floor.
  • String-Length CUDA kernel: stable in-block bitonic sort — 75x over numpy, ~2x the old vectorized cupy path; torch path fuses sort+gather via torch.sort(stable=True).
  • CE: exact int shared-memory count atomics (half the smem) + block-parallel entropy reduction.
  • BLS device cache: one upload per light curve across period segments (cupy+torch), host-side t_min (no per-segment sync), one stacked D2H for the 7 outputs.
  • GLS NUFFT: base-grid pair batched as one n_trans=2 transform (finufft + cufinufft plan cache); one-shot cufinufft assembles on device (194->290 lc/s).

Shims

scatter_add numpy via buffered bincount (vs unbuffered add.at), cupy via cupyx.scatter_add; new row-wise scatter helpers avoid materializing (P,N) broadcast copies on torch. ensure_shared_memory now opts in on dynamic+static totals (fixes a CUDA_ERROR_INVALID_VALUE edge at exactly 48 KB).

Not included (documented reasoning)

  • TLS Ofir grid: requires physical stellar-density durations, which breaks TLS's period-independent phase-bin widths — a scientific redesign, not an optimization.
  • TLS FFT correlation: break-even at n_bins=256 and rounding can flip the dip-admissibility test.

Verification

233 tests pass (12 new in tests/test_fast_backends.py, incl. GPU float32 and tie-stability checks) on a machine with numba + cupy + CUDA torch (Blackwell); ruff and mypy --strict clean; docs build with -W. All measured speedups are in the commit messages; checksums/parities quoted per commit.

🤖 Generated with Claude Code

tjayasinghe and others added 9 commits July 1, 2026 14:03
- scatter_add: numpy via bincount (buffered; ~10-30x over unbuffered
  add.at), cupy via cupyx.scatter_add
- new row-wise scatter_add_rows / scatter_counts_rows: torch scatters a
  stride-0 expanded view (no (P,N) broadcast copy, no flat index); numpy
  gets one fused bincount pass
- BLS/CE/TLS batch kernels use the row-wise shims; TLS also pads once to
  the widest template and accumulates the correlation in place; BLS
  hoists the per-duration column arange
- PDM bins each point once into n_bins*n_covers fine bins and regroups
  covers exactly (roll + group-sum): 3 scatters instead of 3*n_covers on
  the vectorized path, 3 shared-memory atomics per point instead of
  3*n_covers in the CUDA kernel. PDM numpy: 3652 -> 1289 ms on a
  3k-point/20k-period sweep, bit-identical theta.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…tensor

Every normal-equation entry is analytically a harmonic trig sum
(product-to-sum identities), so compute C_m/S_m for m<=2H via the
Chebyshev recurrence from ONE cos/sin evaluation and assemble the tiny
(F,d,d) Gram from them: O(F*N*d^2) -> O(F*N*H) work, O(d) less transient
memory (no more multi-GB design tensors at 1e5 points), and gemm-free on
every backend by construction, which retires the Blackwell cuBLAS
workaround branch. Same math, same ridge; checksums unchanged.

20k points x 30k frequencies, H=3: numpy 81.6s -> 30.1s, cupy 2.04s ->
0.71s, torch:cuda 2.12s -> 0.66s.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The base-grid sums for w and w*y share one angle matrix, and the
doubled-frequency sums follow exactly from the double-angle identities
(cos2t = 2cos^2 t - 1, sin2t = 2 sin t cos t), so the portable direct
path now evaluates cos/sin once per frequency instead of six times.
The frequency chunk is also auto-capped by N so the (chunk, N)
transients stay ~64 MB on any device.

5k points x 100k frequencies: torch:cpu 3.06s -> 1.66s, torch:cuda
264ms -> 150ms; max deviation vs finufft unchanged at 2.2e-11.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…t on device

The w and w*y transforms share their nonuniform points, so they now run
as one n_trans=2 NUFFT (finufft and cufinufft, plan cache keyed by
(nf, n_trans)). The one-shot cufinufft path keeps all three sums on the
GPU and assembles the power there, so only the final spectrum crosses
back: 194 -> 290 lc/s on 800-point curves; engine path 869 -> 888 lc/s.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
One-iteration-per-trial prange kernels mirroring the CUDA designs (PDM
uses the fine-bin regroup; SL uses a stable mergesort argsort so tied
phases pair identically to the array-API backends; MHAOV runs the
Chebyshev harmonic recurrence in scalar registers with a per-frequency
d x d solve). cpu/auto now resolve to numba when the [fast] extra is
installed via a shared fast_cpu_backend attribute, which also replaces
the bespoke BLS resolve_backend override.

3k points, warm kernels, vs the vectorized numpy paths: PDM 306x, CE
135x, MHAOV 57x, TLS 53x, SL 23x; parity <= 2.5e-11 (TLS bit-identical).

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…e cache

- BLS/PDM/CE/TLS kernel sources are templated on the working precision
  (REAL) and honor settings.precision on the cupy backends: float64 stays
  the default; float32 is an opt-in that runs the FLOP-bound box/matched-
  filter scans at full consumer-GPU rate (Blackwell GeForce: BLS 610 ->
  71 ms, TLS 34 -> 4 ms, ~8.6x). PDM keeps double accumulators (the
  ssd = sq - sum^2/cnt form is cancellation-prone) and CE counts are now
  exact int shared atomics (half the shared memory) with the entropy
  reduction spread across the block instead of one thread.
- TLS caches its concatenated templates in shared memory.
- BLS gains a caller-owned per-light-curve device cache so the segmented
  search uploads tau/yw/ivar once per run (cupy + torch), computes t_min
  on the host (no per-segment device sync), and returns all seven outputs
  in one stacked D2H copy. Multiband keeps one cache per band.
- numba BLS zeroes only the n_bins+1 bin prefix actually used per period.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- New one-block-per-period RawKernel: folds, bitonic-sorts the (phase,
  original-index) pairs in shared memory (index tie-break = the same
  stable order as the array-API paths), and reduces the string length
  with no (P, N) intermediates. Curves beyond the shared-memory
  capacity (~8k points on Blackwell) fall back to the vectorized path.
  3k points x 50k periods: 123 ms vs 9.3 s numpy, parity 1e-13 with
  forced phase ties.
- torch path uses torch.sort(stable=True) to fuse the argsort +
  row-gather pair (57 ms on torch:cuda).
- ensure_shared_memory now opts in whenever dynamic+static shared
  exceeds the 48 KB default; dynamic == 48 KB alongside static arrays
  previously failed with CUDA_ERROR_INVALID_VALUE.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
- tests/test_fast_backends.py: numba-vs-numpy parity for PDM/CE/SL/
  MHAOV/TLS (SL with forced phase ties), float32 CUDA kernel checks
  (BLS tolerates isolated near-tie box flips but requires the same
  detected period), scatter-shim unit tests (numpy + torch), BLS
  device-cache reuse equivalence, and cpu->numba resolution for all six
  methods; new requires_numba marker.
- CHANGELOG Performance section; backends/installation/methods/index
  docs updated for the numba tier and float32 kernels; mypy override
  for cupyx.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
The ubuntu-latest image stopped shipping libEGL.so.1, so pytest-qt now
fails at configure importing PySide6 QtGui (INTERNALERROR before any
test runs). Unrelated to the code changes on this branch.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
@tjayasinghe
tjayasinghe merged commit 75eee04 into release/v1.1.0 Jul 2, 2026
6 checks passed
tjayasinghe added a commit that referenced this pull request Jul 2, 2026
…ons merge

Re-ran the full validation + benchmark suite against merged PR #15 (numba CPU
tier for PDM/CE/String-Length/MHAOV/TLS, algorithmic rewrites, precision-
templated CUDA kernels). Correctness holds at the same tolerances as before,
but the new CPU tier is now fast enough to match or beat the GPU on a single
light curve for PDM/CE/MHAOV, overturning the old "50-190x GPU speedup"
figures. Updated REPORT.md, docs/benchmarks.md, docs/guide/backends.md,
docs/guide/batch.md, and README.md to reflect the current numbers, and fixed
the stale hardcoded prose in make_report.py so future re-runs stay accurate.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
@tjayasinghe
tjayasinghe deleted the feat/perf-optimizations branch August 14, 2026 23:37
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.

1 participant