Skip to main content

Documentation Index

Fetch the complete documentation index at: https://docs.haiqu.ai/llms.txt

Use this file to discover all available pages before exploring further.

​
Haiqu.postprocess_skqd(results, num_shots, h1e, h2e, norb, nelec, skqd_options)

Submit SKQD postprocessing (SQD diagonalization) as an async job. Takes the raw probability distributions from all Krylov circuit executions and the Hamiltonian tensors, submits them to the Haiqu server for classical SQD diagonalization on a background worker, and returns a job handle. Call job.result() to block until the worker finishes and get an SKQDResult, or job.progress() for live updates. The caller is responsible for passing the correct tensors in the same basis used for circuit generation. For example, if circuits were built in momentum basis (SIAM), pass momentum-basis tensors.
  • Parameters:
    • results (list *[*dict *[*str , float ] ]) — List of probability distribution dicts (bitstring -> float) from all Krylov circuit executions, as returned by haiqu.run().result().
    • num_shots (int) — Number of shots used per circuit execution.
    • h1e — One-body Hamiltonian tensor, shape (norb, norb). Accepts numpy arrays or nested lists. Must be in the same basis as the circuits.
    • h2e — Two-body Hamiltonian tensor, shape (norb, norb, norb, norb). Accepts numpy arrays or nested lists. Must be in the same basis as the circuits.
    • norb (int) — Number of spatial orbitals.
    • nelec (tuple *[*int , int ]) — Tuple of (n_alpha, n_beta) electron counts.
    • skqd_options (SKQDOptions) — SKQDOptions with SQD parameters (samples_per_batch, num_batches, max_iterations, symmetrize_spin, configuration_recovery, seed).
  • Returns: Job handle. : Call job.result() to retrieve an SKQDResult containing the best ground-state energy, the CI subspace dimension, the CI amplitudes and determinant strings (ci_strs_a, ci_strs_b), per-orbital orbital_occupancies_alpha / _beta, and the per-iteration iteration_history. job.info exposes the SKQD output payload (same fields as SKQDResult). Use job.progress() for live status updates and help(job.result) for the full description of result and info contents.
  • Return type: SKQDJobModel

​
Examples

>>> from haiqu.sdk.skqd import hubbard_hamiltonian, SKQDOptions, build_hubbard_site_basis_krylov_circuits
>>> norb = 8
>>> h1e, h2e = hubbard_hamiltonian(norb, t=1.0, U=8.0)
>>> circuits = build_hubbard_site_basis_krylov_circuits(norb, krylov_dim=6, dt=0.2, h1e=h1e, h2e=h2e)
>>> run_job = haiqu.run(circuits, shots=1000, device_id="aer_simulator")
>>> skqd_job = haiqu.postprocess_skqd(
...     results=run_job.result(), num_shots=1000,
...     h1e=h1e, h2e=h2e,
...     norb=norb, nelec=(4, 4), skqd_options=SKQDOptions(),
... )
>>> result = skqd_job.result()
>>> print(result.energy)

​
How it works

This section traces what postprocess_skqd actually does to your hardware results, step by step, on the same small SIAM model used in the Configuration recovery tutorial (norb=4, nelec=(2, 2), 8 qubits). The setup: you’ve already built five Krylov circuits at evolution times 0,Δt,2Δt,3Δt,4Δt0, \Delta t, 2\Delta t, 3\Delta t, 4\Delta t with build_siam_momentum_basis_krylov_circuits, run each with 1000 shots on hardware (or a simulator), and now hold five probability distributions over 8-bit measurement outcomes. Each dist_i is a plain Python dict mapping bitstring strings to probabilities, e.g. 
dist_0 = {
    "00110011": 0.41,   # dominant determinant
    "01010011": 0.18,
    "00110101": 0.17,
    "01010101": 0.08,
    "00111100": 0.04,   # higher-energy, smaller weight
    "00010011": 0.02,   # out of sector — readout flipped one bit
    "10110011": 0.02,   # out of sector
    # ... 30 or so more bitstrings, each with small probability
}
You then call postprocess_skqd(results=[dist_0, dist_1, dist_2, dist_3, dist_4], ...). The next sections describe what happens on the worker between the call and result.energy.

​
Step 0: merge the five distributions into one pool

The worker doesn’t care which Krylov circuit produced which bitstring. It only cares about the union of bitstrings that came out of the Krylov subspace, weighted by how often they were measured. So the first thing it does is convert each distribution back into a multiset of samples and concatenate them. For each dist_i, every bitstring b with probability p is reproduced round(p × num_shots) times in a list, and the five lists are concatenated into a single BitArray of about 5000 samples. If the same bitstring appears in several circuits — which is typical, the dominant ground-state determinants show up at every Krylov index — its counts simply accumulate. For example, 00110011 might have probabilities 0.41,0.36,0.30,0.22,0.150.41, 0.36, 0.30, 0.22, 0.15 across the five circuits, contributing roughly 410+360+300+220+150=1440410 + 360 + 300 + 220 + 150 = 1440 samples to the pool.
The Krylov-circuit-of-origin information is intentionally discarded here. SQD only needs the empirical probability distribution over electron configurations; the role of the Krylov circuits is upstream — they prepare states whose measurement distribution spans the relevant subspace.
After the merge, the pool has roughly 5000 samples and somewhere around 60–80 unique bitstrings.

​
Iteration 1

Of the ~5000 merged samples, the ones with the wrong electron count get discarded. With 8 qubits at ~1% per-qubit readout-error rate, roughly 84% of samples land in-sector (have exactly two ones in each half). So we go from ~5000 to ~4200 samples and from ~80 unique bitstrings to ~36 (or fewer, if some valid configurations weren’t sampled). 
| Bitstring   | Total count (pool) | In-sector? | Surviving count |
|-------------|-------------------:|:----------:|----------------:|
| `00110011`  | 1440               | yes        | 1440            |
| `01010011`  | 620                | yes        | 620             |
| `00110101`  | 580                | yes        | 580             |
| `01010101`  | 280                | yes        | 280             |
| `00111100`  | 110                | yes        | 110             |
| …           | …                  | …          | …               |
| `00010011`  | 95                 | no         | 0 (dropped)     |
| `10110011`  | 80                 | no         | 0 (dropped)     |

With `configuration_recovery=True` this step would instead *repair* the out-of-sector bitstrings by flipping bits toward the current orbital occupancies. On iteration 1 there are no occupancies yet, so postselection is the only option in either case.
The worker draws num_batches × samples_per_batch = 5 × 100 = 500 bitstrings from the in-sector pool, with replacement, weighted by count. The five batches are independent random draws and almost always differ from each other. 
Each batch becomes the input to one independent diagonalization. Running five independent batches and keeping the best is a variance-reduction trick — it gives both a tighter upper bound on the energy and a sense of how stable the answer is.
For each batch, the worker splits each bitstring into its alpha half (right four bits) and beta half (left four bits), takes the unique alpha and beta substrings, and forms the CI subspace as their Cartesian product. 
For our small SIAM, the (2, 2) sector contains $\binom{4}{2}^2 = 36$ determinants total. With 100 samples per batch and only six possible alpha substrings, each batch typically covers the full sector after a few dozen draws. So in this example each batch's subspace is the full 36-determinant sector and the diagonalization returns the *exact* ground-state energy of the sector, $E = -8.3089$.

On bigger systems the subspace is only a tiny corner of the full sector; the five batches each cover *different* corners, each gives a different (looser) energy estimate, and the lowest one wins.
From the winning batch’s ground state, the worker reads off the per-orbital occupancies — for our SIAM these come out to [0.99,0.50,0.50,0.01][0.99, 0.50, 0.50, 0.01] for both spin sectors — and stores them in current_occupancies. 
It also picks out every CI string whose amplitude in the ground state exceeds `carryover_threshold` (default $10^{-4}$) and saves them in `carryover_strings_a` and `carryover_strings_b` ([fermion.py:399-410](https://github.com/Qiskit/qiskit-addon-sqd/blob/main/qiskit_addon_sqd/fermion.py#L399)). On the next iteration, these strings will be prepended to whatever the new random batches happen to draw, so they can never be lost to subsampling randomness.

​
Iteration 2: what changes

The same four-step sequence runs again, but with two pieces of memory now in play.
  1. The bitstring pool may shift. With configuration_recovery=False (the default) the pool is the same in-sector samples as iteration 1. With configuration_recovery=True, the worker re-runs recover_configurations against the new occupancies [0.99,0.50,0.50,0.01][0.99, 0.50, 0.50, 0.01], repairing the previously out-of-sector samples back into the (2, 2) sector and adding them to the pool. See Configuration recovery for the exact bit-flip rule.
  2. Subsamples are fresh. The RNG state has advanced, so even with the same pool the five batches draw different bitstrings.
  3. Carryover is prepended. Each batch’s CI string list is now (include + carryover_from_iter_1 + new_samples) deduplicated (fermion.py:360-361). The high-weight determinants from iteration 1 are guaranteed to be in iteration 2’s basis regardless of what the random subsamples bring in.
For our 4-orbital SIAM this doesn’t show — the subspace already saturated at 36 determinants in iteration 1, so iteration 2 produces the same energy. But on bigger systems, where the full sector is enormous and each batch covers only a small corner, iteration 2’s basis is strictly larger than iteration 1’s. That is exactly what produces the subspace-dimension growth across iterations you see in result.iteration_history — for example, a typical 360→400→485→485360 \to 400 \to 485 \to 485 trajectory: each round adds the strings carryover preserved plus whatever new high-weight strings the fresh random subsamples surface, until new draws stop adding anything new and the subspace plateaus.

​
Convergence

The loop checks two conditions at the end of each iteration (fermion.py:384-394):
  • the energy difference from the previous iteration is below energy_tol, and
  • the maximum change in any orbital occupancy is below occupancies_tol.
If both hold, it stops. Otherwise it iterates again, up to max_iterations. The function returns the best result seen across all iterations, not necessarily the last one. For our small SIAM, convergence happens at iteration 2 — energy and occupancies haven’t moved since iteration 1.

​
What the result looks like

print(result.energy)              # -8.3089158728  (matches exact diagonalization)
print(result.subspace_dimension)  # 36  (the full (2, 2) sector)
print(len(result.iteration_history))  # 2

print(result.orbital_occupancies_alpha)
# array([0.99, 0.50, 0.50, 0.01])

# Per-iteration record
for i, it in enumerate(result.iteration_history):
    print(i, it["energy"], it["subspace_dimension"], it["runtime"])
# 0 -8.308915 36 0.017
# 1 -8.308915 36 0.014
iteration_history is the place to look when diagnosing convergence — if energy keeps falling at the last iteration, you ran out of max_iterations before convergence and should raise it; if subspace dimension keeps growing without the energy improving, your batches are too small or carryover_threshold is too low.

​
See also

  • Configuration recovery — what changes when configuration_recovery=True, with worked bit-flip examples.
  • Yu, J. et al. Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization. arXiv:2501.09702 (2025).
  • Implementation: “qiskit_addon_sqd.fermion.diagonalize_fermionic_hamiltonian <[https://github.com/Qiskit/qiskit-addon-sqd/blob/main/qiskit_addon_sqd/fermion.py](https://github.com/Qiskit/qiskit-addon-sqd/blob/main/qiskit_addon_sqd/fermion.py)>\__.