SOSSA hamiltonian builder and circuit mapper

[Yingrong Chen (YingrongChen)]

Implements Sum-of-Squares Spectral Amplification (SOSSA) following Low et al. 2025, PRX 15, with the qubitized walk and unary-iteration machinery it needs.

Stacked on #621

This PR is the top layer of a native GitHub stack (#686):

main  <-  #621  Alias sampling and QROM state preparation
          #539  SOSSA hamiltonian builder and circuit mapper   <- this PR

Alias sampling was merged into this branch and this PR's base retargeted onto #621, so
the diff here is now the SOSSA-only remainder: 44 files / +6,887, down from 70 files
/ +9,639 against main. Review #621 first; GitHub will cascade the merge.

What this PR adds

SOSSA block encoding and walk (SOSSAWalk.qs), the double-factorized Hamiltonian unitary builder, the SOSSA controlled circuit mapper, and resource estimation support. The alias-sampling / QROM / SelectSwap state preparation this builds on is #621, the layer below.

Data flow

%%{init: {"flowchart": {"padding": 30, "nodeSpacing": 60, "rankSpacing": 55}}}%%
flowchart TD
    subgraph S1["Stage 1 — qubit mapping"]
        HAM["Hamiltonian ⟨FactorizedHamiltonianContainer⟩<br/>U : R*B*N &nbsp; W : R*B*C &nbsp; WB : R,C<br/>signs : R, h1, E_core, E_gap"]
        HAM --> QM["&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;SOSSAQubitMapper&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;"]
        QM --> SOSC["QubitOperator ⟨SOSSAContainer⟩<br/>one_body : RotatedPaulis angles, coeffs, X/Y<br/>two_body : RotatedPaulis angles, coeffs, Z<br/>num_positive_one_body_terms, energy_shift"]
    end

    subgraph S2["Stage 2 — unitary builder"]
        SOSC --> SB["&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;SOSSABuilder&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;"]
        SB --> WALK["UnitaryRepresentation ⟨SOSSAWalkContainer⟩<br/>outer_prepare : Wavefunction ∝ c<br/>outer_prepare_probabilities ∝ c²<br/>SOSSAInnerPrepare : cond. coeffs + free-rider bits<br/>SOSSASelect : Givens angles<br/>Λ = ½Σc², power, energy_shift"]
    end

    subgraph S3["Stage 3 — circuit mapper"]
        WALK --> CM["&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;SOSSAMapper&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;"]
        CM --> CIRC["Circuit<br/>MakeControlledSOSSAWalkOp / …Circuit"]
    end

    classDef data fill:#dbeafe,stroke:#2563eb,color:#0b1b34
    classDef algo fill:#fef3c7,stroke:#d97706,color:#3b2a06
    class HAM,SOSC,WALK,CIRC data
    class QM,SB,CM algo

Loading

Testing

End-to-end energy recovery is asserted for H2 in both the iterative and unary QPE paths. Both use the quadratically suppressed tolerance Λ(1 − cos(2π/N)) rather than the linear Λ·2π/N — SOSSA's ground state sits at φ ≈ 1/2 where dE/dφ = 0, so the linear bound is roughly 15× too loose to be meaningful. Each test also asserts that the amplified bound is at least 5× tighter than the linear one, so the distinction can't silently regress.

[Copilot]

🔵 Needs a closer look

The new examples/sossa_qre.ipynb contains runtime-breaking algorithm-name and attribute-access issues that would prevent the notebook (and its new test) from executing successfully.

Review details

Suppressed comments (3)

examples/sossa_qre.ipynb:189

• SOSSAWalkContainer exposes the block-encoding normalization as container.normalization, not under container.metadata (which is FactorizedHamiltonianMetadata). Using .metadata.normalization will raise AttributeError at runtime.

    "    hf_config = Configuration.canonical_hf_configuration(n_alpha, n_beta, num_orbitals)\n",

examples/sossa_qre.ipynb:426

• metadata does not contain normalization; normalization is a property on the SOSSA walk container itself. Also, the qubit mapper variant name should be "sos" (not "sossa"). As written, this cell will fail at runtime.

    "print(f\"System qubits                       = {n2_operator.get_container().num_qubits}\")\n",

examples/sossa_qre.ipynb:185

• The SOS mapper is registered as "sos" (SOSQubitMapper). Using "sossa" here will fail algorithm lookup and prevent the notebook from executing.

    "\n",
    "    if num_queries is None:\n",
    "        sossa_walk = create(\"hamiltonian_unitary_builder\", \"sossa\").run(operator)\n",
    "        lambda_sos = sossa_walk.get_container().metadata.normalization\n",

• Files reviewed: 32/32 changed files
• Comments generated: 0 new
• Review effort level: Lite

[Copilot]

🔵 Needs a closer look

The new examples/sossa_qre.ipynb references a non-existent qubit-mapper variant ("sossa") and reads normalization from the wrong object (metadata), which will break notebook execution and the added notebook test.

Review details

Suppressed comments (4)

examples/sossa_qre.ipynb:179

• The notebook uses create("qubit_mapper", "sossa"), but the SOS factorized-Hamiltonian qubit mapper is registered under the variant name "sos" (see SOSQubitMapper.name()), so this cell will raise at runtime and the notebook test will fail.

    "    operator = create(\"qubit_mapper\", \"sossa\").run(\n",

examples/sossa_qre.ipynb:185

• SOSSAWalkContainer stores the block-encoding normalization on container.normalization, not on container.metadata. Using metadata.normalization will raise AttributeError when executing the notebook.

    "        lambda_sos = sossa_walk.get_container().metadata.normalization\n",

examples/sossa_qre.ipynb:419

• This cell also calls create("qubit_mapper", "sossa"), but the qubit-mapper variant name is "sos"; leaving this as-is will break notebook execution.

    "n2_operator = create(\"qubit_mapper\", \"sossa\").run(\n",

examples/sossa_qre.ipynb:426

• create("hamiltonian_unitary_builder", "sossa").run(...).get_container().metadata returns FactorizedHamiltonianMetadata, which does not include the block-encoding normalization. The notebook then reads n2_metadata.normalization and will fail; keep the walk container around and read normalization from it.

    "n2_metadata = create(\"hamiltonian_unitary_builder\", \"sossa\").run(n2_operator).get_container().metadata\n",
    "\n",
    "print(f\"Block-encoding normalization  Lambda = {n2_metadata.normalization:.6f} Hartree\")\n",
    "print(f\"Sum-of-squares shift         E_SOS  = {n2_metadata.energy_shift:.6f} Hartree\")\n",
    "print(f\"System qubits                       = {n2_operator.get_container().num_qubits}\")\n",

• Files reviewed: 32/32 changed files
• Comments generated: 0 new
• Review effort level: Lite

[Copilot]

🔵 Needs a closer look

The SOSSA QRE notebook currently contains malformed indentation that will cause SyntaxError during execution, and there are phase-gradient precision consistency issues in the SOSSA circuit mapper that can break valid configurations.

Review details

Suppressed comments (6)

Previously missed (4) — in code that hasn't changed since the last review.

python/src/qdk_chemistry/algorithms/circuit_mapper/sossa_mapper.py:105

• When outer_prepare is qrom, the outer PREPARE must share the same phase-gradient width as SELECT. Currently only allocate_phase_gradient is forced off; rotation_bit_precision is left at the QROM state's default (often 10), which can disagree with SOSSAMapper.rotation_bit_precision and trigger a layout/metadata mismatch.
  python/src/qdk_chemistry/algorithms/circuit_mapper/sossa_mapper.py:222
• _num_phase_gradient_qubits reads the QROM state's own rotation_bit_precision setting, which can diverge from SOSSAMapper.rotation_bit_precision (used for SELECT). Since the walk shares a single phase-gradient register, the width should be derived from a single source of truth (e.g., the mapper setting) to avoid spurious mismatches.
  python/src/qdk_chemistry/algorithms/qubit_mapper/sos.py:40
• This mapper ignores the MajoranaMapping passed to run(), but QubitMapper’s documented contract is that encoding is determined by the provided mapping. Unconditionally warning and silently forcing Jordan–Wigner can produce surprising results for callers and noisy logs; it should instead validate that the requested mapping is supported (and matches the Hamiltonian’s mode count) and then proceed without warning.
  python/tests/test_helpers.py:283
• These basis vectors are normalized but not orthogonal; calling them “orthogonal” is misleading for readers of the test helper.

examples/sossa_qre.ipynb:183

• This code cell line closes the AlgorithmRef call without the indentation needed inside the surrounding if circuit_mapper is None: block, which will raise a SyntaxError when the notebook is executed.

    ")\n",

examples/sossa_qre.ipynb:191

• This code cell line closes the builder = create(...) call without the indentation needed at the function scope, which will raise a SyntaxError when the notebook is executed.

    ")\n",

• Files reviewed: 32/32 changed files
• Comments generated: 0 new
• Review effort level: Lite

[Copilot]

🟡 Changes recommended

There are API/runtime correctness issues in newly added SOS/Pauli-LCU container code paths that should be fixed before approval.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Review details

Suppressed comments (1)

Previously missed (1) — in code that hasn't changed since the last review.

python/src/qdk_chemistry/data/qubit_operator/containers/pauli_lcu.py:143

• PauliLCUContainer stores coefficients without coercing it to a NumPy array. Several methods in this class (e.g., to_json() uses coeffs.real/coeffs.imag, hashing uses _hash_array) assume an np.ndarray, so passing a Python list (or other sequence) will break at runtime even though the constructor type hints accept np.ndarray and earlier behavior allowed array-like inputs.

• Files reviewed: 32/32 changed files
• Comments generated: 1
• Review effort level: Lite

[Copilot]

🟡 Changes recommended

SOSSA inner-PREPARE coefficient construction currently drops the sign of the SF identity-term coefficient, which breaks the phase/sign encoding needed for correct SOSSA generators.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Review details

Suppressed comments (1)

Previously missed (1) — in code that hasn't changed since the last review.

python/src/qdk_chemistry/algorithms/hamiltonian_unitary_builder/block_encoding/sossa.py:169

• The docstring says the inner conditional coefficients include the "absolute identity weight", but the Q# alias-sampling PREPARE uses coefficient sign to encode a phase. Documenting the identity term as absolute is misleading and suggests (incorrectly) that its sign is irrelevant.

• Files reviewed: 29/29 changed files
• Comments generated: 1
• Review effort level: Lite

[Rushi Gong (RushiGong)]

spinSF is in |+⟩ inside the reflection (H at :397, reflected at :402), and this Z fires exactly when isSF = 0 — so every D1/Q1 term picks up ⟨+|Z|+⟩ = 0 and the one-body Hamiltonian drops out of the circuit entirely. Driving the block encoding through the simulator on the H2 fixture, the measured block differs from the exact H_gap/Λ − I by 0.85, with only the two spin-free indices surviving and the one-body pair carrying 86% of the outer weight gone. Should the coherent Majorana sum sit on an inner register instead — and does MajoranaOp (:641-661) need fixing alongside, since it reuses one spin qubit for both the register swap and the X-vs-iY choice?

[Rushi Gong (RushiGong)]

This sandwich is exp(iθ Y_j X_{j+1}), which rotates only the X-type Majoranas and mixes the |00⟩/|11⟩ sector, so it doesn't preserve particle number — I measured a violation of 0.972 against an operator norm of 1, producing |⟨1111|B|0000⟩| = 0.13 where the exact target is identically 0. There's also a transposition: within { chain } apply gives C† · Op · C, which uses row 0 of the single-particle matrix rather than column 0, so the rotated mode is the reflection of u. Is the comment at :499-502 describing an intended design here, or documenting this as a symptom?

[Rushi Gong (RushiGong)]

The inner register is uncomputed before the reflection (SOSSAWalk.qs:394-400), so the block-encoded generator is Σ_b |β_b|² U_b — feeding raw w_b makes the spin-free generator quadratic in the weights instead of linear, and AliasSamplingStatePrep.qs:155 squares them too. Should these be sign-stripped √|w_b| normalised to the row one-norm that :160 already computes? Relatedly, where does sign(w_b) reach SELECT — both backends conjugate the phase away inside PREP† … PREP and MajoranaOp carries no per-b sign, though H2's all-positive, single-dominant weights hide both issues.

[Rushi Gong (RushiGong)]

Exact-zero eigenvalues fall through both masks, so one_body gets fewer than N rows — but outer_prep_dim (block_encoding/sossa.py:141) and _compute_free_rider_data (:254) still reserve N one-body slots, so every spin-free amplitude decodes as Q1 with the wrong rank. SOSSAWalk.TestSelectDQ reproduces it as QSharpError: index out of range: 1, and it's reachable from the public constructor — a pure two-body payload gives eigenvalues of exactly 0.0, with validation checking dimensions but not rank (and no shipped fixture is rank-deficient, so the tests can't catch it). Would you rather pad the dropped modes with √λ = 0 generators, or derive x_o from the actual one-body row count?

[Rushi Gong (RushiGong)]

Beyond the exact-zero case above, would a screening threshold help here — eigensolver noise at ~1e-17 passes > 0.0 and becomes a real generator riding an essentially arbitrary nullspace eigenvector. On an N=54, rank-10 probe that produced 44 spurious generators and 2,332 meaningless Givens angles. The package already defaults to 1e-12 in qdk_qubit_mapper.py:39-40, though any fix still has to preserve the N outer slots, so it couples to the layout question.

[Rushi Gong (RushiGong)]

This densifies a 24-qubit dump — 16.7M amplitudes, ~256 MiB of numpy plus ~128 MiB of list pointers — and post-processes it in a Python loop I measured at 14.873 s. :165-173 does 25 qubits (~512 MiB) and :436-443 does 25 qubits twice, scaling cleanly as 2^n from 0.506 s at 20 qubits. None of the three are marked @pytest.mark.slow, so they all land in the default CI run — is that intentional?

[Rushi Gong (RushiGong)]

This docstring says the class verifies Toffoli scaling and that walk power multiplies Toffolis linearly, but the only test (:506-542) asserts a qubit-count bound. There's no Toffoli assertion, and power is never varied. Trim the docstring, or is there a missing test?

[Rushi Gong (RushiGong)]

pauli_lcu.py:465-469 writes {real, imag, dtype} and restores the dtype on load at :518-521, while this one writes {real, imag} and hardcodes complex. Harmless in practice since SOS coefficients are always complex128, but it's a near-duplicate helper with a different wire format inside the same package. Reuse the existing one, or match its shape?

[Rushi Gong (RushiGong)]

Berry2024 is the only bibliography change in this PR, but seems like nothing cites it.

[Rushi Gong (RushiGong)]

This comment writes the identity coefficient as W^{rc} I, but the code puts identity_weights there — the paper's w_B^{(rc)} — while the local w0 at :101 is the actual W^{(rc)}. The two symbols are swapped in the one place where the distinction decides a sign. Renaming w0 to something like w_rc and correcting the comment would help, especially given the np.abs question on that same weight.

[Copilot]

🟡 Changes recommended

There are at least two correctness issues that can break SOSSA invariants (generator count vs N, and a Q# inner-PREPARE bit-width edge case) that should be fixed before approval.

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Review details

Suppressed comments (1)

Previously missed (1) — in code that hasn't changed since the last review.

python/src/qdk_chemistry/utils/qsharp/src/SOSSAWalk.qs:718

• MakeInnerPrepareAliasSampling computes nIndexBits as BitSizeI(nCoeffs - 1) without guarding nCoeffs<=1. If B=0 (so B+1=1) or a degenerate coefficient row is passed, this can yield 0 index bits and make slices like innerReg[0..nIndexBits - 1] invalid. The direct inner-prepare path already uses the guarded BitSizeI((if nCoeffs > 1 { nCoeffs } else { 2 }) - 1) pattern; the alias-sampling path should match it.

• Files reviewed: 29/29 changed files
• Comments generated: 1
• Review effort level: Lite

[Copilot]

🟡 Changes recommended

There are verified correctness issues in the SOS/SOSSA coefficient construction that can alter the encoded operator (identity-term handling and sign preservation).

Once you've addressed the issues Copilot identified, you can request another Copilot review.

Review details

Suppressed comments (1)

python/src/qdk_chemistry/algorithms/hamiltonian_unitary_builder/block_encoding/sossa.py:179

• The inner-PREPARE coefficient table keeps signed coefficients so ConditionalAliasSamplingPrepareWithFreeRider can apply the sign as a phase (see Q# doc comment in SOSSAWalk.qs). Taking np.abs of the identity-column (b == B) silently drops the sign of that coefficient, which changes the SF generator (the identity term participates in cross-terms inside the square).

Keep the signed identity coefficient instead of its magnitude.

        sf = np.asarray(sossa.two_body.coeffs)
        sf_rows = np.zeros((sf.shape[0], b_plus_1))
        if sf.size:
            sf_rows[:, :-1] = sf[:, :-1].real
            sf_rows[:, -1] = np.abs(sf[:, -1])

• Files reviewed: 29/29 changed files
• Comments generated: 1
• Review effort level: Lite