> ## 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.

# Postprocessing

#### 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*](helpers.md#haiqu.sdk.skqd.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

```python theme={null}
>>> 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](configuration_recovery.md) tutorial
(`norb=4`, `nelec=(2, 2)`, 8 qubits).

The setup: you’ve already built five Krylov circuits at evolution times
$0, \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.

```python theme={null}
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.15$ across the five
circuits, contributing roughly
$410 + 360 + 300 + 220 + 150 = 1440$ samples to the pool.

<Note>
  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.
</Note>

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).

```default theme={null}
| 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.

```default theme={null}
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.

```default theme={null}
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]$ for both spin sectors — and stores them
in `current_occupancies`.

```default theme={null}
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]$, repairing the previously
   out-of-sector samples back into the (2, 2) sector and adding them to
   the pool. See [Configuration recovery](configuration_recovery.md)
   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](https://github.com/Qiskit/qiskit-addon-sqd/blob/main/qiskit_addon_sqd/fermion.py#L360)).
   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 \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](https://github.com/Qiskit/qiskit-addon-sqd/blob/main/qiskit_addon_sqd/fermion.py#L384)):

* 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

```python theme={null}
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](configuration_recovery.md) — 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](https://arxiv.org/abs/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)>\`\_\_.