Skip to content

ppvm-lindblad: deterministic row-parallel matvec for the cached expm generator - #252

Open
AlexSchuckert wants to merge 2 commits into
mainfrom
perf/deterministic-csr-matvec
Open

AlexSchuckert wants to merge 2 commits into
mainfrom
perf/deterministic-csr-matvec

Conversation

@AlexSchuckert

@AlexSchuckert AlexSchuckert commented Oct 16, 2026 •

Copy link
Copy Markdown
Collaborator

Summary

This makes exp(dt·M)·b in ppvm-lindblad bit-identical for every rayon thread count. It does so by replacing the column-chunked CscOp::dot with a row-major (CSR) cache and a row-parallel matvec. I'm framing it as a determinism and scaling fix, not a speed-up.

Problem. CscOp::dot split the columns into rayon::current_num_threads() chunks, accumulated each chunk into a dense length-n partial vector, and summed the partials on one core. The summation order, and hence the round-off, depended on the thread count. A rank-capped evolution amplifies that round-off to the size of the truncation error. On the 82-qubit XX ladder (open two-leg ladder, centre-rung spin density, A = 3B, drop_tol = 0), changing RAYON_NUM_THREADS moved the MSD by up to 1.6e-2 by Jt = 2. A stored run could therefore only be reproduced at the same core count. The design also allocated threads × n scalars per matvec and did an O(threads × n) serial reduction.

Change.

  • The action pass is still column-wise. Each column block stages its entries in exactly-sized buckets, one per row range (16 B/entry).
  • A second parallel pass counting-sorts every range into CSR (12 B/entry), dropping the buckets as it goes. The cache is never held in both layouts at once.
  • The range size is the next power of two ≥ max(64, ⌈n/256⌉) rows, so it depends on n only.
  • dot computes each output entry on one thread, over its row in ascending column order.
  • Together with a canonical term order in compute_action_terms (sorted by word hash), every sum in the cache build, the matvec and the leakage merge runs in an order fixed by the basis alone. Results differ from main only at round-off level (≤ 3.6e-15 in the ladder MSD by Jt = 2).
  • No partial vectors, no serial reduction.

Nothing else in the step depends on the thread count, for these reasons:

  • quspin's apply runs sequentially (parallel_hint = false).
  • The leakage merge runs in fixed chunk order.
  • The norm estimates are maxima.

Tests.

  • pc_step_is_bit_identical_across_thread_counts: a rank-capped XXZ chain with Z dephasing (n = 14, cap 1200, admit 3600, 12 steps), compared at 1/2/3/8 threads. It fails on main at 2 threads.
  • pc_step_orbit_rep_is_bit_identical_across_thread_counts: the same check on the orbit-rep path (ring with collective decay, cap binding).
  • action_term_order_is_independent_of_scratch_capacity: compute_action_terms must return the same order from a fresh and from a pre-grown scratch map. It fails without the sort. All ppvm-lindblad tests and the Python lindblad tests pass, and cargo fmt and clippy --all-targets are clean.

Measurements

Apple M4 (4P + 6E cores), 82-qubit XX ladder, T = 2, δ = 0.025.

Revisions covered. The baseline is main at 1d78a44b. The determinism, cost and breakdown numbers below were measured on the first revision, a73875f6. For the current revision, a5271fcf (canonical term order and review fixes), I re-verified bit-identity across 1/2/4/8/10 threads at B = 2^12–2^14 and peak RSS at 2^16/2^18. Its wall times are still to follow: the machine was shared with other jobs while I measured.

Determinism, B = 2^12–2^16, 1/2/4/8(6)/10 threads:

  • This branch: the MSD and the full density profile are bit-identical across all thread counts. a73875f6 was also identical to main at 1 thread; a5271fcf differs from it by ≤ 3.6e-15 in MSD.
  • main: max|ΔMSD| vs its own 1-thread run is up to 1.6e-2 (B = 2^13), 5.0e-3 (2^14) and 1.9e-3 (2^16).

Cost, 6 threads, two interleaved repetitions:

B build wall [s] peak RSS [MiB]
2^16 main 9.65 / 10.14 336 / 332
2^16 this PR 9.69 / 10.05 316 / 303
2^18 main 60.6 / 65.8 674 / 673
2^18 this PR 58.1 / 63.4 711 / 710

Expm breakdown at 2^18 / 6 threads (n = 786k, nnz = 4.4M):

main this PR
time per matvec 5.94 ms 2.35 ms (2.5× faster)
cache build, whole run 25.3 s 27.3 s (+8%)
matvecs, whole run 10.5 s 4.2 s

At 2^16 the time per matvec (ms) scales with threads as follows:

threads 1 2 4 6 10
main 1.12 0.76 0.68 0.85 1.05
this PR 1.01 0.68 0.43 0.46 0.42

The old matvec degrades above 4 threads; this one keeps scaling.

Costs, stated honestly:

  • The whole-step wall time changes within run-to-run noise (about 8%), because the old matvec was only 13–17% of the wall at 2^18.
  • The cache build is 5–8% slower.
  • Peak RSS is 6–9% lower at 2^16 but 5.5% higher at 2^18 / 6 threads. This most likely comes from bucket staging that the allocator does not fully reuse.
  • The memory trade should reverse at high thread counts, because the old partials scale with the thread count: about 1.6 GB per matvec at B = 2^21 and 32 threads, vs about 0.6 GB of staging here. This is not measured.

Deliberately not included

Canonical (|c|, word) tie-breaking in cap_basis / cap_map_to_room.

  • It is not needed for thread independence: once the matvec is deterministic, row order and hash iteration order are deterministic too.
  • Ties do decide the cut on the ladder, so the rule moves the MSD by about 2e-3 to 4e-2. That is the size of the old thread spread, and there is no accuracy argument for either rule.

Leakage merge. In local profiles, the leakage merge and room cap are the next serial bottleneck, about 60% of the wall at 10 threads. A separate branch addresses that.

🤖 Generated with Claude Code

…generator

CscOp::dot split the columns into rayon::current_num_threads() chunks,
accumulated each into a dense length-n partial vector and summed the
partials on one core. The summation order, and hence the round-off,
depended on the thread count; a rank-capped evolution amplifies that
round-off to the size of the truncation error. It also allocated
threads x n scalars per matvec and did an O(threads x n) serial reduction.

The action pass now stages each entry in an exactly-sized bucket of its
row range (a range size that depends on n only), and a second parallel
pass counting-sorts every range into CSR, dropping the buckets as it goes.
The matvec computes each output entry on one thread over its row in
ascending column order, so exp(dt M) b is bit-identical for every thread
count and equal to the previous single-thread result. No partial vectors,
no serial reduction.

Adds pc_step_is_bit_identical_across_thread_counts (fails on the old
matvec at 2 threads).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
@github-actions

github-actions Bot commented Oct 16, 2026 •

Copy link
Copy Markdown
PR Preview Action v1.8.1

QR code for preview link

🚀 View preview at
https://QuEraComputing.github.io/ppvm/pr-preview/pr-252/

Built to branch gh-pages at 2026-10-16 13:38 UTC.
Preview will be ready when the GitHub Pages deployment is complete.

@AlexSchuckert
AlexSchuckert marked this pull request as ready for review October 16, 2026 12:34
Copilot AI balanced review requested due to automatic review settings October 16, 2026 12:34

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🟢 Approval recommended

The cache conversion preserves column accumulation order, and the review identified no blocking defects.

0 open findings

What changed in this PR

Replaces the cached Lindblad generator’s column-based matvec with row-parallel CSR accumulation to eliminate thread-dependent summation order.

Changes:

  • Builds row-major caches and removes per-thread partial vectors and serial reduction.
  • Adds a rank-capped evolution regression test comparing bitwise results across 1/2/3/8 threads.
File Description
crates/​ppvm-lindblad/​src/​tests.rs Adds cross-thread determinism coverage.
crates/​ppvm-lindblad/​src/​mf_expm.rs Implements CSR caching and row-parallel matvec.

🧠 Review effort: Balanced


💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.

…cache

- compute_action_terms returns its terms sorted by (word_hash, Word).
  The drain order of the reused scratch FxHashMap follows its capacity,
  which carries over between calls on a rayon job and so depends on the
  thread count. Every caller sums or inserts in that order: the orbit-rep
  duplicate columns, raw/diag (hence the (m*, s) choice), and the leakage
  merge.
- Checked u32 range on the per-range CSR nnz, replacing the dropped
  "CSC block exceeds u32 entries" guard.
- Rows per range: a power of two with a 64-row floor (was max(1024, n/256)),
  so small and medium dims split into up to 256 tasks, and the range lookup
  in the build loop is a shift; the per-block counts buffer is reused via
  map_init.
- The matvec skips zero inputs again, as the column-wise product did.
- Docs: staging cost for complex T, row-major cache wording, `mat` naming.
- Tests: action_term_order_is_independent_of_scratch_capacity (fails
  without the sort) and pc_step_orbit_rep_is_bit_identical_across_thread_counts.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Copilot AI balanced review requested due to automatic review settings October 16, 2026 13:35

Copilot AI left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

🔵 Needs a closer look

The numerical-kernel rewrite needs human validation of regression results and historical bitwise-compatibility claims.

1 open finding

🧠 Review effort: Balanced

Comment thread crates/ppvm-lindblad/src/spec.rs

This branch has not been deployed

No deployments
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