Quickstart

This guide demonstrates a minimal working example of the qiskit-addon-sqd package. We use sample-based quantum diagonalization (SQD) to refine noisy samples taken from a quantum circuit that approximates the ground state of a chemistry Hamiltonian. Ground state energies are approximated by diagonalizing the Hamiltonian in subspaces formed by the refined samples.

For more explanatory end-to-end workflows that run on quantum hardware, check out the tutorials:

• Sample-based quantum diagonalization: chemistry Hamiltonian workflow
• Sample-based Krylov quantum diagonalization: fermionic lattice model workflow

1. Prepare the inputs for SQD

First, we define the chemistry Hamiltonian and quantum circuit that approximates its ground state. The system studied in this example is N₂ at equilibrium geometry in STO-3G basis with two core orbitals frozen. This system has 8 spatial orbitals and 10 electrons, which means the quantum circuit is defined on 16 qubits, where qubits 0-7 represent spin-up orbitals and qubits 8-15 represent spin-down orbitals. Since this is a closed-shell system with 10 electrons, we know each physically-valid sample should have hamming weight 10, where the hamming weight of both the left and right halves of the bitstring is 5.

from pyscf import ao2mo, gto, mcscf, scf

# Construct the Hamiltonian, pull out integrals and other information about the system
mol = gto.M(
    atom="N 0 0 0; N 0 0 1.09768", basis="sto-3g", spin=0, charge=0, verbose=0
)
mf = scf.RHF(mol).run()
cas = mcscf.CASCI(mf, ncas=8, nelecas=10)
hcore, nuclear_repulsion_energy = cas.get_h1eff()
eri = cas.get_h2eff()
num_orbitals = cas.ncas
eri = ao2mo.restore(1, eri, num_orbitals)
n_alpha, n_beta = cas.nelecas
nelec = (n_alpha, n_beta)

print(
    f"Spatial orbitals: {num_orbitals}\nQubits: {num_orbitals * 2}\nElectrons (alpha, beta): {nelec}"
)

Output:

Spatial orbitals: 8
Qubits: 16
Electrons (alpha, beta): (5, 5)

2. Generate noisy samples

In a real workflow you would design an ansatz circuit, optimize it for hardware, and sample it on the QPU. Here we skip that and draw 10,000 uniformly random bitstrings instead. While our samples contain no information about the actual ground state, we will see that the SQD protocol can still recover the true ground state by iteratively refining these samples, diagonalizing, and using the energies at each iteration to inform the next round of bitstring refinement. We can eventually recover the exact ground state in this manner because the subspace containing the ground state is small enough to calculate exactly. Specifically, since we know each half of each bitstring (8 bits) should have hamming weight 5, we know the full Hilbert subspace of the system is size $\binom{8}{5}^2=3136$.

import numpy as np
from qiskit_addon_sqd.counts import generate_bit_array_uniform

# Generate some uniformly-random bitstrings to simulate QPU samples
rng = np.random.default_rng(24)
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)
print(
    f"Generated {bit_array.num_shots} uniformly-random, {bit_array.num_bits}-qubit samples."
)

Output:

Generated 10000 uniformly-random, 16-qubit samples.

3. Post-process with SQD

The objective of the SQD algorithm is to refine noisy samples taken from a QPU and ultimately diagonalize the Hamiltonian in a subspace spanned by those refined samples. This is done iteratively as follows:

• Refine samples by postselecting only those which have correct particle number for both the spin-up and spin-down subsystems
• Diagonalize in small subspaces spanned by the refined samples and collect information about the subspaces that yield the the best (lowest) energies.
• For each iteration:
  • Refine the noisy samples by probabilistically flipping bits in samples with incorrect particle number. The probabilities are calculated using the subspace information from the most recent diagonalization.
  • Diagonalize in small subspaces spanned by the refined samples
  • Collect information about the subspaces yielding the best (lowest) energies

We can see from the output that the algorithm converges to the exact ground state energy after 6 iterations and that the diagonalization is performed over the full Hilbert subspace of dimension $\binom{8}{5}^2=3136$. Of course, diagonalizing over the full subspace is usually not feasible, and therefore, systems with ground states having sparse support are the most likely to benefit from subspace diagonalization techniques such as SQD.

For more information on the diagonalize_fermionic_hamiltonian function, visit the API docs.

from functools import partial

from pyscf import fci
from qiskit_addon_sqd.fermion import (
    SCIResult,
    diagonalize_fermionic_hamiltonian,
    solve_sci_batch,
)

exact_energy, _ = fci.direct_spin1.kernel(
    hcore, eri, int(num_orbitals), nelec
)
exact_energy += nuclear_repulsion_energy
print(f"Exact (FCI) reference energy: {exact_energy:.6f} Ha\n")

sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=200)

result_history = []


def callback(results: list[SCIResult]):
    result_history.append(results)
    iteration = len(result_history)
    print(f"Iteration {iteration}")
    for i, result in enumerate(results):
        error = result.energy + nuclear_repulsion_energy - exact_energy
        print(f"  Subsample {i}")
        print(f"    Energy: {result.energy + nuclear_repulsion_energy:.6f}")
        print(
            f"    Subspace dimension: {np.prod(result.sci_state.amplitudes.shape)}"
        )
        print(f"    Error vs exact: {error:.6f} Ha")


samples_per_batch = 20
result = diagonalize_fermionic_hamiltonian(
    hcore,
    eri,
    bit_array,
    samples_per_batch=samples_per_batch,
    norb=num_orbitals,
    nelec=nelec,
    occupancies_tol=1e-7,
    max_iterations=30,
    sci_solver=sci_solver,
    symmetrize_spin=True,
    callback=callback,
    seed=np.random.default_rng(32),
)

Output:

Exact (FCI) reference energy: -107.652521 Ha

Iteration 1
  Subsample 0
    Energy: -106.908481
    Subspace dimension: 784
    Error vs exact: 0.744040 Ha
Iteration 2
  Subsample 0
    Energy: -106.934921
    Subspace dimension: 1764
    Error vs exact: 0.717600 Ha
Iteration 3
  Subsample 0
    Energy: -106.953430
    Subspace dimension: 2209
    Error vs exact: 0.699091 Ha
Iteration 4
  Subsample 0
    Energy: -107.651946
    Subspace dimension: 2916
    Error vs exact: 0.000575 Ha
Iteration 5
  Subsample 0
    Energy: -107.652003
    Subspace dimension: 3025
    Error vs exact: 0.000518 Ha
Iteration 6
  Subsample 0
    Energy: -107.652521
    Subspace dimension: 3136
    Error vs exact: 0.000000 Ha
Iteration 7
  Subsample 0
    Energy: -107.652521
    Subspace dimension: 3136
    Error vs exact: 0.000000 Ha