Begrenzung der Unterraumdimension

In dieser Schritt-für-Schritt-Anleitung zeigen wir die Wirkung der Unterraumdimension auf die selbstkonsistente Konfigurationswiederherstellung.

A priori wissen wir nicht, welche Unterraumdimension nötig ist, um ein angestrebtes Genauigkeitsniveau zu erreichen. Wir wissen jedoch, dass eine größere Unterraumdimension die Genauigkeit der Methode erhöht. Daher können wir die Genauigkeit der Vorhersagen als Funktion der Unterraumdimension untersuchen. Lege das Molekül und seine Eigenschaften fest.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "pyscf": "pyscf", "qiskit_addon_sqd": "qiskit-addon-sqd"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

import warnings

import pyscf
import pyscf.cc
import pyscf.mcscf

warnings.filterwarnings("ignore")

# Specify molecule properties
open_shell = False
spin_sq = 0

# Build N2 molecule
mol = pyscf.gto.Mole()
mol.build(
    atom=[["N", (0, 0, 0)], ["N", (1.0, 0, 0)]],
    basis="6-31g",
    symmetry="Dooh",
)

# Define active space
n_frozen = 2
active_space = range(n_frozen, mol.nao_nr())

# Get molecular integrals
scf = pyscf.scf.RHF(mol).run()
num_orbitals = len(active_space)
n_electrons = int(sum(scf.mo_occ[active_space]))
num_elec_a = (n_electrons + mol.spin) // 2
num_elec_b = (n_electrons - mol.spin) // 2
cas = pyscf.mcscf.CASCI(scf, num_orbitals, (num_elec_a, num_elec_b))
mo = cas.sort_mo(active_space, base=0)
hcore, nuclear_repulsion_energy = cas.get_h1cas(mo)
eri = pyscf.ao2mo.restore(1, cas.get_h2cas(mo), num_orbitals)

# Compute exact energy
exact_energy = cas.run().e_tot

converged SCF energy = -108.835236570775
CASCI E = -109.046671778080  E(CI) = -32.8155692383188  S^2 = 0.0000000

Erzeuge einige zufällige Bitstrings als Stellvertreter für QPU-Samples.

import numpy as np
from qiskit_addon_sqd.counts import generate_bit_array_uniform

# Create a seed to control randomness throughout this workflow
rng = np.random.default_rng(24)

# Generate random samples
bit_array = generate_bit_array_uniform(
    10_000, num_orbitals * 2, rand_seed=rng
)

Rufe SQD mit zunehmenden Batch-Größen auf.

from qiskit_addon_sqd.fermion import diagonalize_fermionic_hamiltonian

list_samples_per_batch = [50, 200, 400, 600]

# SQD options
max_iterations = 5

# Eigenstate solver options
num_batches = 10
max_davidson_cycles = 200

energies = []
subspace_dimensions = []

for samples_per_batch in list_samples_per_batch:
    result = diagonalize_fermionic_hamiltonian(
        hcore,
        eri,
        bit_array,
        samples_per_batch=samples_per_batch,
        norb=num_orbitals,
        nelec=(num_elec_a, num_elec_b),
        num_batches=num_batches,
        max_iterations=max_iterations,
        symmetrize_spin=True,
        seed=rng,
    )
    energies.append(result.energy)
    subspace_dimensions.append(np.prod(result.sci_state.amplitudes.shape))

Dieses Diagramm zeigt, dass eine größere Unterraumdimension zu genaueren Ergebnissen führt.

import matplotlib.pyplot as plt

# Data for energies plot
x1 = subspace_dimensions
y1 = np.array(energies) + nuclear_repulsion_energy

fig, axs = plt.subplots(1, 1, figsize=(12, 6))

# Plot energies
axs.plot(x1, y1, marker=".", markersize=20, label="Estimated")
axs.set_xticks(x1)
axs.set_xticklabels(x1)
axs.axhline(y=exact_energy, color="red", linestyle="--", label="Exact")
axs.set_title("Approximated Ground State Energy vs subspace dimension")
axs.set_xlabel("Subspace dimension")
axs.set_ylabel("Energy (Ha)")
axs.legend()

plt.tight_layout()
plt.show()

Eine explizite Grenze mit max_dim setzen

Im obigen Sweep wurde die Unterraumdimension nur indirekt gesteuert. samples_per_batch legt fest, wie viele Bitstrings in jeden Batch gezogen werden, aber die Dimension des resultierenden Unterraums hängt davon ab, wie viele verschiedene Spin-Alpha- und Spin-Beta-Strings diese Bitstrings enthalten, und der Pool der für das Subsampling verfügbaren Konfigurationen wächst von einer Iteration der Konfigurationswiederherstellung zur nächsten. Deshalb sind die Unterraumdimensionen auf der horizontalen Achse oben unregelmäßige Zahlen, die wir nicht selbst gewählt haben.

Das Argument max_dim von diagonalize_fermionic_hamiltonian begrenzt die Dimension direkt. Es beschränkt die Anzahl der CI-Strings, die in jedem Spinsektor behalten werden, und da der Unterraum vom kartesischen Produkt der beiden Sektoren aufgespannt wird, ist seine Dimension das Produkt der beiden Sektordimensionen:

• max_dim=d beschränkt beide Sektoren auf d Strings, sodass die Unterraumdimension höchstens d ** 2 beträgt.

• max_dim=(d_a, d_b) beschränkt den Spin-Alpha- und den Spin-Beta-Sektor getrennt, sodass die Unterraumdimension höchstens d_a * d_b beträgt. Bei symmetrize_spin=True müssen die beiden Grenzen gleich sein.

• max_dim=None, der Standardwert, setzt keine Grenze.

Wenn ein Batch mehr Strings enthält, als die Grenze zulässt, werden die unwichtigsten verworfen: Zuerst werden die über include_configurations übergebenen Strings behalten, danach die aus der vorherigen Iteration übernommenen Strings (siehe carryover_threshold), und die übrigen Plätze werden mit den gesampelten Strings aufgefüllt, und zwar in absteigender Reihenfolge ihrer Samplehäufigkeit.

Da max_dim die Größe der größten Diagonalisierung begrenzt, die vom Eigenzustands-Solver verlangt wird, ist es ein bequemer Weg, den Solver innerhalb eines Laufzeit- und Speicherbudgets zu halten: Wähle samples_per_batch großzügig und lass max_dim entscheiden, wie groß der Unterraum werden darf.

Im Folgenden wiederholen wir die Untersuchung, diesmal mit festem samples_per_batch und zunehmendem max_dim.

list_max_dim = [100, 200, 400, 800]

# Use a generous number of samples per batch, so that max_dim is what limits
# the dimension of the subspace
samples_per_batch = 600

energies_max_dim = []
subspace_dimensions_max_dim = []

for max_dim in list_max_dim:
    result = diagonalize_fermionic_hamiltonian(
        hcore,
        eri,
        bit_array,
        samples_per_batch=samples_per_batch,
        norb=num_orbitals,
        nelec=(num_elec_a, num_elec_b),
        num_batches=num_batches,
        max_iterations=max_iterations,
        symmetrize_spin=True,
        max_dim=max_dim,
        seed=rng,
    )
    energies_max_dim.append(result.energy)
    subspace_dimensions_max_dim.append(
        int(np.prod(result.sci_state.amplitudes.shape))
    )

print(f"max_dim:            {list_max_dim}")
print(f"subspace dimension: {subspace_dimensions_max_dim}")

max_dim:            [100, 200, 400, 800]
subspace dimension: [10000, 40000, 160000, 640000]

Die Unterraumdimension beträgt nun in jedem Lauf genau max_dim ** 2: Sie wurde gewählt statt entdeckt. Die gemeinsame Darstellung beider Sweeps zeigt, dass die Genauigkeit von der Unterraumdimension selbst bestimmt wird, unabhängig davon, über welches Argument man sie erreicht hat. Die beiden Kurven müssen jedoch nicht exakt übereinstimmen. Die Läufe mit max_dim ziehen 600 Bitstrings pro Batch und behalten dann die am häufigsten gesampelten Strings, sodass sie bei einer gegebenen Dimension eine bessere Menge an Konfigurationen behalten können als ein Lauf, der diese Dimension von vornherein durch das Sampling weniger Bitstrings erreicht hat.

fig, axs = plt.subplots(1, 1, figsize=(12, 6))

axs.plot(
    subspace_dimensions,
    np.array(energies) + nuclear_repulsion_energy,
    marker=".",
    markersize=20,
    label="Varying samples_per_batch",
)
axs.plot(
    subspace_dimensions_max_dim,
    np.array(energies_max_dim) + nuclear_repulsion_energy,
    marker="^",
    markersize=12,
    label="Varying max_dim",
)
axs.axhline(y=exact_energy, color="red", linestyle="--", label="Exact")
axs.set_xscale("log")
axs.ticklabel_format(axis="y", useOffset=False)
axs.set_title("Approximated Ground State Energy vs subspace dimension")
axs.set_xlabel("Subspace dimension")
axs.set_ylabel("Energy (Ha)")
axs.legend()

plt.tight_layout()
plt.show()