Skip to content

SQD tutorial: add section on passing initial occupancies #4521

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kevinsung wants to merge 10 commits into
base: main
Choose a base branch
from
Open

SQD tutorial: add section on passing initial occupancies #4521

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
8ffa5bb
SQD tutorial: add section on passing initial occupancies
kevinsung Jan 10, 2026
cf0ba3c
add to existing workflow instead of separate section
kevinsung Jan 14, 2026
7ba019b
add cell to count valid bitstrings
kevinsung Jan 14, 2026
7449680
Merge branch 'main' into kjs/sqd-initial-occupancies
kevinsung Jan 14, 2026
7db26e0
edit variable names
kevinsung Jan 14, 2026
36494a6
re-run with boston
kevinsung Jan 14, 2026
a334951
re-run with https://github.com/Qiskit/qiskit-addon-sqd/pull/276
kevinsung Jan 15, 2026
0290050
Merge branch 'main' into kjs/sqd-initial-occupancies
kevinsung Jan 15, 2026
69c594d
run ./fix
kevinsung Jan 15, 2026
e6c44fd
print final energy and update reference energy
kevinsung Jan 15, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
205 changes: 147 additions & 58 deletions docs/tutorials/sample-based-quantum-diagonalization.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,16 +167,16 @@
"metadata": {},
"outputs": [],
"source": [
"import math\n",
"\n",
"import ffsim\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import pyscf\n",
"import pyscf.cc\n",
"import pyscf.mcscf\n",
"import ffsim\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"from qiskit import QuantumCircuit, QuantumRegister\n",
"from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
"\n",
"from qiskit_ibm_runtime import QiskitRuntimeService\n",
"from qiskit_ibm_runtime import SamplerV2 as Sampler"
]
Expand Down Expand Up @@ -224,17 +224,17 @@
"\n",
"# Get molecular integrals\n",
"scf = pyscf.scf.RHF(mol).run()\n",
"num_orbitals = len(active_space)\n",
"norb = len(active_space)\n",
"n_electrons = int(sum(scf.mo_occ[active_space]))\n",
"num_elec_a = (n_electrons + mol.spin) // 2\n",
"num_elec_b = (n_electrons - mol.spin) // 2\n",
"cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))\n",
"n_alpha = (n_electrons + mol.spin) // 2\n",
"n_beta = (n_electrons - mol.spin) // 2\n",
"cas = pyscf.mcscf.CASCI(scf, norb, (n_alpha, n_beta))\n",
"mo = cas.sort_mo(active_space, base=0)\n",
"hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)\n",
"eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), norb)\n",
"\n",
"# Store reference energy from SCI calculation performed separately\n",
"exact_energy = -109.22690201485733"
"reference_energy = -109.22802921665716"
]
},
{
Expand All @@ -255,7 +255,7 @@
"name": "stdout",
"output_type": "stream",
"text": [
"E(CCSD) = -109.2177884185543 E_corr = -0.2879500329450045\n"
"E(CCSD) = -109.2177884185545 E_corr = -0.2879500329450042\n"
]
}
],
Expand Down Expand Up @@ -285,8 +285,8 @@
"outputs": [],
"source": [
"n_reps = 1\n",
"alpha_alpha_indices = [(p, p + 1) for p in range(num_orbitals - 1)]\n",
"alpha_beta_indices = [(p, p) for p in range(0, num_orbitals, 4)]\n",
"alpha_alpha_indices = [(p, p + 1) for p in range(norb - 1)]\n",
"alpha_beta_indices = [(p, p) for p in range(0, norb, 4)]\n",
"\n",
"\n",
"ucj_op = ffsim.UCJOpSpinBalanced.from_t_amplitudes(\n",
Expand All @@ -301,14 +301,14 @@
" options=dict(maxiter=1000),\n",
")\n",
"\n",
"nelec = (num_elec_a, num_elec_b)\n",
"nelec = (n_alpha, n_beta)\n",
"\n",
"# create an empty quantum circuit\n",
"qubits = QuantumRegister(2 * num_orbitals, name=\"q\")\n",
"qubits = QuantumRegister(2 * norb, name=\"q\")\n",
"circuit = QuantumCircuit(qubits)\n",
"\n",
"# prepare Hartree-Fock state as the reference state and append it to the quantum circuit\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(num_orbitals, nelec), qubits)\n",
"circuit.append(ffsim.qiskit.PrepareHartreeFockJW(norb, nelec), qubits)\n",
"\n",
"# apply the UCJ operator to the reference state\n",
"circuit.append(ffsim.qiskit.UCJOpSpinBalancedJW(ucj_op), qubits)\n",
Expand Down Expand Up @@ -612,13 +612,13 @@
"name": "stdout",
"output_type": "stream",
"text": [
"Gate counts (w/o pre-init passes): OrderedDict({'rz': 12438, 'sx': 12169, 'cz': 3986, 'x': 573, 'measure': 52, 'barrier': 1})\n",
"Gate counts (w/ pre-init passes): OrderedDict({'sx': 7059, 'rz': 6962, 'cz': 1876, 'measure': 52, 'x': 35, 'barrier': 1})\n"
"Gate counts (w/o pre-init passes): OrderedDict({'rz': 12407, 'sx': 12155, 'cz': 3986, 'x': 581, 'measure': 52, 'barrier': 1})\n",
"Gate counts (w/ pre-init passes): OrderedDict({'sx': 7059, 'rz': 6980, 'cz': 1876, 'measure': 52, 'x': 35, 'barrier': 1})\n"
]
}
],
"source": [
"initial_layout, _ = get_zigzag_physical_layout(num_orbitals, backend=backend)\n",
"initial_layout, _ = get_zigzag_physical_layout(norb, backend=backend)\n",
"\n",
"pass_manager = generate_preset_pass_manager(\n",
" optimization_level=3, backend=backend, initial_layout=initial_layout\n",
Expand Down Expand Up @@ -683,15 +683,90 @@
},
{
"cell_type": "markdown",
"id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
"id": "a8736983",
"metadata": {},
"source": [
"Now, we estimate the ground state energy of the Hamiltonian using the `diagonalize_fermionic_hamiltonian` function. This function performs the self-consistent configuration recovery procedure to iteratively refine the noisy quantum samples to improve the energy estimate. We pass a callback function so that we can save the intermediate results for later analysis. See the [API documentation](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) for explanations of the arguments to `diagonalize_fermionic_hamiltonian`."
"A useful metric to judge the quality of the QPU output is the number of valid configurations returned. A valid configuration has the correct particle number and spin Z, which means that the right half of the bitstring has Hamming weight equal to the number of spin-up electrons, and the left half has Hamming weight equal to the number of spin-down electrons. The following cell computes the fraction of sampled configurations that are valid."
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "718f8517",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fraction of sampled configurations that are valid: 0.00065\n"
]
}
],
"source": [
"def is_valid_bitstring(\n",
" bitstring: str, norb: int, nelec: tuple[int, int]\n",
") -> bool:\n",
" n_alpha, n_beta = nelec\n",
" return (\n",
" len(bitstring) == 2 * norb\n",
" and bitstring[norb:].count(\"1\") == n_alpha\n",
" and bitstring[:norb].count(\"1\") == n_beta\n",
" )\n",
"\n",
"\n",
"bit_array = pub_result.data.meas\n",
"num_valid = sum(\n",
" is_valid_bitstring(b, norb, nelec) for b in bit_array.get_bitstrings()\n",
")\n",
"valid_fraction = num_valid / bit_array.num_shots\n",
"print(f\"Fraction of sampled configurations that are valid: {valid_fraction}\")"
]
},
{
"cell_type": "markdown",
"id": "7b126f3a",
"metadata": {},
"source": [
"While this fraction might seem fairly small, it is significantly higher than the fraction of bitstrings you would expect to be valid if the bitstrings were sampled uniformly at random, as shown by the following cell."
]
},
{
"cell_type": "code",
"execution_count": 11,
"id": "6b3e4bca",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Expected fraction of valid configurations from uniformly random bitstrings: 9.607888706852918e-07\n"
]
}
],
"source": [
"expected_fraction_random = (\n",
" math.comb(norb, n_alpha) * math.comb(norb, n_beta) / 2 ** (2 * norb)\n",
")\n",
"print(\n",
" f\"Expected fraction of valid configurations from uniformly random bitstrings: {expected_fraction_random}\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "eb704101-0fe8-4d12-b572-b1d844e35a90",
"metadata": {},
"source": [
"Now, we estimate the ground state energy of the Hamiltonian using the `diagonalize_fermionic_hamiltonian` function. This function performs the self-consistent configuration recovery procedure to iteratively refine the noisy quantum samples to improve the energy estimate. We pass a callback function so that we can save the intermediate results for later analysis. See the [API documentation](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian) for explanations of the arguments to `diagonalize_fermionic_hamiltonian`.\n",
"\n",
"Here, we use the `initial_occupancies` argument to `diagonalize_fermionic_hamiltonian` to specify the Hartree-Fock configuration as the initial guess for the orbital occupancies in the ground state. This approach is sensible for systems where the ground state has significant support on the Hartree-Fock configuration, but it might not be appropriate in other situations, though more advanced computational methods might yield better initial guesses in those cases. Specifying `initial_occupancies` also allows configuration recovery to run even if no valid configurations were sampled, as may be the case when sampling a large circuit on a noisy QPU. Without this argument, the configuration recovery would fail and raise an error if no valid configurations were provided."
]
},
{
"cell_type": "code",
"execution_count": 12,
"id": "2f32a352",
"metadata": {},
"outputs": [
Expand All @@ -701,54 +776,56 @@
"text": [
"Iteration 1\n",
"\tSubsample 0\n",
"\t\tEnergy: -108.97781410104506\n",
"\t\tSubspace dimension: 28561\n",
"\t\tEnergy: -109.04961576529212\n",
"\t\tSubspace dimension: 291600\n",
"\tSubsample 1\n",
"\t\tEnergy: -108.97781410104506\n",
"\t\tSubspace dimension: 28561\n",
"\t\tEnergy: -109.05419845387246\n",
"\t\tSubspace dimension: 298116\n",
"\tSubsample 2\n",
"\t\tEnergy: -108.97781410104506\n",
"\t\tSubspace dimension: 28561\n",
"\t\tEnergy: -109.03983430633693\n",
"\t\tSubspace dimension: 291600\n",
"Iteration 2\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.05150860754739\n",
"\t\tSubspace dimension: 287296\n",
"\t\tEnergy: -109.13341224824521\n",
"\t\tSubspace dimension: 432964\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.08152283263908\n",
"\t\tSubspace dimension: 264196\n",
"\t\tEnergy: -109.12829311613213\n",
"\t\tSubspace dimension: 434281\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.11636893067873\n",
"\t\tSubspace dimension: 284089\n",
"\t\tEnergy: -109.13563153146234\n",
"\t\tSubspace dimension: 444889\n",
"Iteration 3\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.15878555367885\n",
"\t\tSubspace dimension: 429025\n",
"\t\tEnergy: -109.15805816816365\n",
"\t\tSubspace dimension: 644809\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.16462831275786\n",
"\t\tSubspace dimension: 473344\n",
"\t\tEnergy: -109.17104671394021\n",
"\t\tSubspace dimension: 588289\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.15895026995382\n",
"\t\tSubspace dimension: 435600\n",
"\t\tEnergy: -109.17084043977475\n",
"\t\tSubspace dimension: 585225\n",
"Iteration 4\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.17784051223317\n",
"\t\tSubspace dimension: 622521\n",
"\t\tEnergy: -109.18154982106486\n",
"\t\tSubspace dimension: 776161\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.1774651326829\n",
"\t\tSubspace dimension: 657721\n",
"\t\tEnergy: -109.17757205183068\n",
"\t\tSubspace dimension: 762129\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.18085212360191\n",
"\t\tSubspace dimension: 617796\n",
"\t\tEnergy: -109.17890903245448\n",
"\t\tSubspace dimension: 760384\n",
"Iteration 5\n",
"\tSubsample 0\n",
"\t\tEnergy: -109.18636242518915\n",
"\t\tSubspace dimension: 815409\n",
"\t\tEnergy: -109.18924157389883\n",
"\t\tSubspace dimension: 966289\n",
"\tSubsample 1\n",
"\t\tEnergy: -109.18451014767518\n",
"\t\tSubspace dimension: 837225\n",
"\t\tEnergy: -109.18449630323602\n",
"\t\tSubspace dimension: 952576\n",
"\tSubsample 2\n",
"\t\tEnergy: -109.18333728638888\n",
"\t\tSubspace dimension: 857476\n"
"\t\tEnergy: -109.18632692423985\n",
"\t\tSubspace dimension: 942841\n",
"Final energy: -109.18924157389883\n",
"Final energy error: 0.03878764275833646\n"
]
}
],
Expand All @@ -773,6 +850,12 @@
"carryover_threshold = 1e-4\n",
"max_cycle = 200\n",
"\n",
"# Use the Hartree-Fock configuration as an initial guess for the orbital occupancies\n",
"initial_occupancies = (\n",
" np.array([1] * n_alpha + [0] * (norb - n_alpha)),\n",
" np.array([1] * n_beta + [0] * (norb - n_beta)),\n",
")\n",
"\n",
"# Pass options to the built-in eigensolver. If you just want to use the defaults,\n",
"# you can omit this step, in which case you would not specify the sci_solver argument\n",
"# in the call to diagonalize_fermionic_hamiltonian below.\n",
Expand All @@ -797,20 +880,26 @@
"result = diagonalize_fermionic_hamiltonian(\n",
" hcore,\n",
" eri,\n",
" pub_result.data.meas,\n",
" bit_array,\n",
" samples_per_batch=samples_per_batch,\n",
" norb=num_orbitals,\n",
" norb=norb,\n",
" nelec=nelec,\n",
" num_batches=num_batches,\n",
" energy_tol=energy_tol,\n",
" occupancies_tol=occupancies_tol,\n",
" max_iterations=max_iterations,\n",
" sci_solver=sci_solver,\n",
" symmetrize_spin=symmetrize_spin,\n",
" initial_occupancies=initial_occupancies,\n",
" carryover_threshold=carryover_threshold,\n",
" callback=callback,\n",
" seed=12345,\n",
")"
")\n",
"\n",
"final_energy = result.energy + nuclear_repulsion_energy\n",
"energy_error = final_energy - reference_energy\n",
"print(f\"Final energy: {final_energy}\")\n",
"print(f\"Final energy error: {energy_error}\")"
]
},
{
Expand All @@ -820,14 +909,14 @@
"source": [
"### Visualize the results\n",
"\n",
"The first plot shows that after a couple of iterations we estimate the ground state energy within ``~41 mH`` (chemical accuracy is typically accepted to be ``1 kcal/mol`` $\\approx$ ``1.6 mH``). The energy can be improved by allowing more iterations of configuration recovery or increasing the number of samples per batch.\n",
"The first plot shows that after a couple of iterations we estimate the ground state energy within ``~40 mH`` (chemical accuracy is typically accepted to be ``1 kcal/mol`` $\\approx$ ``1.6 mH``). The energy can be improved by allowing more iterations of configuration recovery or increasing the number of samples per batch.\n",
"\n",
"The second plot shows the average occupancy of each spatial orbital after the final iteration. We can see that both the spin-up and spin-down electrons occupy the first five orbitals with high probability in our solutions."
]
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 13,
"id": "caffd888-e89c-4aa9-8bae-4d1bb723b35e",
"metadata": {},
"outputs": [
Expand All @@ -848,7 +937,7 @@
" min(result, key=lambda res: res.energy).energy + nuclear_repulsion_energy\n",
" for result in result_history\n",
"]\n",
"e_diff = [abs(e - exact_energy) for e in min_e]\n",
"e_diff = [abs(e - reference_energy) for e in min_e]\n",
"yt1 = [1.0, 1e-1, 1e-2, 1e-3, 1e-4]\n",
"\n",
"# Chemical accuracy (+/- 1 milli-Hartree)\n",
Expand Down
Binary file modified ...sed-quantum-diagonalization/extracted-outputs/caffd888-e89c-4aa9-8bae-4d1bb723b35e-0.avif
View file
Open in desktop
Binary file not shown.
Loading