Перейти до основного вмісту

Покращення оцінки енергії фермійонної ґраткової моделі за допомогою SQD

У цьому посібнику ми реалізуємо шаблон Qiskit, що показує, як постобробляти зашумлені квантові зразки для знаходження наближення до основного стану фермійонного ґраткового гамільтоніана, відомого як модель одиничної домішки Андерсона. Ми дотримуватимемося підходу квантової діагоналізації на основі вибірки для обробки зразків, взятих із набору базисних станів Крилова з 16 кубітами за зростаючими часовими інтервалами. Ці стани реалізуються на квантовому пристрої за допомогою тротеризації часової еволюції. Щоб врахувати вплив квантового шуму, використовується метод відновлення конфігурацій. За умови гарного початкового стану та розрідженості основного стану, доведено, що цей підхід ефективно сходиться.

Шаблон можна описати у чотири кроки:

  1. Крок 1: Відображення на квантову задачу

    • Генерація набору базисних станів Крилова (тобто тротеризованих схем часової еволюції) за зростаючими часовими інтервалами для оцінки основного стану

  2. Крок 2: Оптимізація задачі

    • Транспіляція Circuit для Backend

  3. Крок 3: Виконання експериментів

    • Отримання зразків із Circuit за допомогою примітиву Sampler

  4. Крок 4: Постобробка результатів

    • Самоузгоджений цикл відновлення конфігурацій
      • Постобробка повного набору бітрядкових зразків з використанням апріорних знань про число частинок і середню орбітальну заповненість, обчислену на останній ітерації
      • Імовірнісне створення пакетів підвибірок із відновлених бітрядків
      • Проєкція та діагоналізація фермійонного ґраткового гамільтоніана в кожному відібраному підпросторі
      • Збереження мінімальної знайденої енергії основного стану по всіх пакетах і оновлення середньої орбітальної заповненості

Крок 1: Відображення задачі на квантову схему​

Спочатку ми створимо одно- та двочастинкові гамільтоніани, що описують одновимірну модель одиничної домішки Андерсона (SIAM) з 7 вузлами ванни (8 електронів у 8 орбіталях). Ця модель використовується для опису магнітних домішок, вбудованих у метали.

Потім ми створимо тротеризовані Circuit з 16 кубітами, які використовуються для генерації квантового підпростору Крилова.

# Added by doQumentation — required packages for this notebook
!pip install -q ffsim matplotlib numpy qiskit qiskit-addon-sqd qiskit-ibm-runtime scipy
import numpy as np

n_bath = 7  # number of bath sites

V = 1  # hybridization amplitude
t = 1  # bath hopping amplitude
U = 10  # Impurity onsite repulsion
eps = -U / 2  # Chemical potential for the impurity

# Place the impurity on the first qubit
impurity_index = 0

# One body matrix elements in the "position" basis
h1e = -t * np.diag(np.ones(n_bath), k=1) - t * np.diag(np.ones(n_bath), k=-1)
h1e[impurity_index, impurity_index + 1] = -V
h1e[impurity_index + 1, impurity_index] = -V
h1e[impurity_index, impurity_index] = eps

# Two body matrix elements in the "position" basis
h2e = np.zeros((n_bath + 1, n_bath + 1, n_bath + 1, n_bath + 1))
h2e[impurity_index, impurity_index, impurity_index, impurity_index] = U

Далі ми генеруємо квантовий підпростір Крилова за допомогою набору тротеризованих квантових схем. Тут ми створюємо допоміжні функції для генерації початкового (еталонного) стану, а також часової еволюції для одно- та двочастинкових частин гамільтоніана. Для детальнішого опису цієї моделі та принципу побудови схем звернись до статті.

import ffsim
import scipy
from qiskit import QuantumCircuit, QuantumRegister
from qiskit.circuit.library import CPhaseGate, XGate, XXPlusYYGate

n_modes = n_bath + 1
nelec = (n_modes // 2, n_modes // 2)

dt = 0.2
Utar = scipy.linalg.expm(-1j * dt * h1e)

# The reference state
def initial_state(q_circuit, norb, nocc):
    """Prepare an initial state."""
    for i in range(nocc):
        q_circuit.append(XGate(), [i])
        q_circuit.append(XGate(), [norb + i])
    rot = XXPlusYYGate(np.pi / 2, -np.pi / 2)

    for i in range(3):
        for j in range(nocc - i - 1, nocc + i, 2):
            q_circuit.append(rot, [j, j + 1])
            q_circuit.append(rot, [norb + j, norb + j + 1])
    q_circuit.append(rot, [j + 1, j + 2])
    q_circuit.append(rot, [norb + j + 1, norb + j + 2])

# The one-body time evolution
free_fermion_evolution = ffsim.qiskit.OrbitalRotationJW(n_modes, Utar)

# The two-body time evolution
def append_diagonal_evolution(dt, U, impurity_qubit, num_orb, q_circuit):
    """Append two-body time evolution to a quantum circuit."""
    if U != 0:
        q_circuit.append(
            CPhaseGate(-dt / 2 * U),
            [impurity_qubit, impurity_qubit + num_orb],
        )

Генеруємо d станів із часовою еволюцією, що визначають квантовий підпростір Крилова. Тут ми обрали d=8. Похибка від вибірки базисних станів Крилова зменшується зі збільшенням d. Зауваж, що особливості цього примірника задачі дозволяють ефективно компілювати однотілову еволюцію за допомогою OrbitalRotationJW; однак загалом можна використовувати PauliEvolutionGate від Qiskit для реалізації тротеризованої часової еволюції повного гамільтоніана.

# Generate the initial state
qubits = QuantumRegister(2 * n_modes, name="q")
init_state = QuantumCircuit(qubits)
initial_state(init_state, n_modes, n_modes // 2)
init_state.draw("mpl", scale=0.4, fold=-1)

d = 8  # Number of Krylov basis states
circuits = []
for i in range(d):
    circ = init_state.copy()
    circuits.append(circ)
    for _ in range(i):
        append_diagonal_evolution(dt, U, impurity_index, n_modes, circ)
        circ.append(free_fermion_evolution, qubits)
        append_diagonal_evolution(dt, U, impurity_index, n_modes, circ)
    circ.measure_all()
circuits[0].draw("mpl", scale=0.4, fold=-1)

Quantum circuit diagram

circuits[-1].draw("mpl", scale=0.4, fold=-1)

Quantum circuit diagram

Крок 2: Оптимізація задачі​

Після створення тротеризованих схем ми оптимізуємо їх для цільового апаратного забезпечення. Перед оптимізацією потрібно вибрати апаратний пристрій. Ми використаємо фіктивний Backend на 127 кубітів із qiskit_ibm_runtime для емуляції реального пристрою. Щоб запустити на реальному QPU, просто замін фіктивний Backend на реальний. Докладніше дивись у документації Qiskit IBM Runtime.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

Далі ми транспілюємо Circuit до цільового Backend за допомогою Qiskit.

from qiskit.transpiler import generate_preset_pass_manager

# The circuit needs to be transpiled to the AerSimulator target
pass_manager = generate_preset_pass_manager(optimization_level=3, backend=backend)
isa_circuits = pass_manager.run(circuits)

Крок 3: Виконання експериментів​

Після оптимізації Circuit для виконання на апаратному забезпеченні ми готові запустити їх на цільовому пристрої та зібрати зразки для оцінки енергії основного стану. Тут ми використовуємо SamplerV2 із qiskit-ibm-runtime для симуляції зашумлених зразків від Backend ibm_sherbrooke. Потім ми об'єднуємо підрахунки від кожного з базисних станів Крилова в єдиний словник підрахунків і будуємо графік 20 найчастіше вибираних бітрядків.

Примітка: Qiskit Aer необхідний для симуляції зразків із транспільованих Circuit.

from qiskit.primitives import BitArray
from qiskit.visualization import plot_histogram
from qiskit_ibm_runtime import SamplerV2 as Sampler

# Sample from the circuits
noisy_sampler = Sampler(backend, options={"simulator": {"seed_simulator": 24}})
job = noisy_sampler.run(isa_circuits, shots=500)

# Combine the counts from the individual Trotter circuits
bit_array = BitArray.concatenate_shots([result.data.meas for result in job.result()])

plot_histogram(bit_array.get_counts(), number_to_keep=20)

Quantum circuit diagram

Крок 4: Постобробка результатів​

Тепер ми запускаємо алгоритм SQD за допомогою функції diagonalize_fermionic_hamiltonian. Пояснення аргументів цієї функції дивись у документації API.

from qiskit_addon_sqd.fermion import SCIResult, diagonalize_fermionic_hamiltonian

# List to capture intermediate results
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):
        print(f"\tSubsample {i}")
        print(f"\t\tEnergy: {result.energy}")
        print(f"\t\tSubspace dimension: {np.prod(result.sci_state.amplitudes.shape)}")

rng = np.random.default_rng(24)
result = diagonalize_fermionic_hamiltonian(
    h1e,
    h2e,
    bit_array,
    samples_per_batch=300,
    norb=n_modes,
    nelec=nelec,
    num_batches=3,
    max_iterations=10,
    symmetrize_spin=True,
    callback=callback,
    seed=rng,
)
Iteration 1
	Subsample 0
		Energy: -13.257128325607997
		Subspace dimension: 3969
	Subsample 1
		Energy: -13.257128325607997
		Subspace dimension: 3969
	Subsample 2
		Energy: -13.257128325607997
		Subspace dimension: 3969
Iteration 2
	Subsample 0
		Energy: -13.319666127542039
		Subspace dimension: 4096
	Subsample 1
		Energy: -13.420534292304595
		Subspace dimension: 4624
	Subsample 2
		Energy: -9.136171594591085
		Subspace dimension: 4624
Iteration 3
	Subsample 0
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 1
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 2
		Energy: -13.422491814612833
		Subspace dimension: 4900
Iteration 4
	Subsample 0
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 1
		Energy: -13.422491814612833
		Subspace dimension: 4900
	Subsample 2
		Energy: -13.422491814612833
		Subspace dimension: 4900

Тепер будуємо графіки результатів. Перший графік показує, що після кількох ітерацій ми отримуємо точну енергію основного стану.

Цей приклад достатньо малий, щоб ми могли дослідити весь простір Гільберта, що видно з вище наведених виведень. Пам'ятай, що повний простір Гільберта має розмірність (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b). Отже, для цієї задачі: (8 choose 4)**2 = 4900. Розмірності підпросторів збільшуються завдяки покращеному відновленню конфігурацій, а також тому, що алгоритм SQD переносить важливі конфігурації з однієї ітерації до наступної. До останньої ітерації ми діагоналізуємо в усьому просторі Гільберта.

Другий графік показує середню заповненість кожної просторової орбіталі по всіх рішеннях пакетів. Ми бачимо, що з високою імовірністю кожна орбіталь містить один електрон.

import matplotlib.pyplot as plt

exact_energy = -13.422491814605827
min_es = [min(result, key=lambda res: res.energy).energy for result in result_history]
min_id, min_e = min(enumerate(min_es), key=lambda x: x[1])

# Data for energies plot
x1 = range(len(result_history))
yt1 = list(np.arange(-13.5, -13.1, 0.1))
ytl = [f"{i:.1f}" for i in yt1]

# Data for avg spatial orbital occupancy
y2 = np.sum(result.orbital_occupancies, axis=0)
x2 = range(len(y2))

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

# Plot energies
axs[0].plot(x1, min_es, label="SQD energy", marker="o")
axs[0].set_xticks(x1)
axs[0].set_xticklabels(x1)
axs[0].set_yticks(yt1)
axs[0].set_yticklabels(ytl)
axs[0].axhline(y=exact_energy, color="#BF5700", linestyle="--", label="Exact energy")
axs[0].set_title("Approximated Ground State Energy vs SQD Iterations")
axs[0].set_xlabel("Iteration Index", fontdict={"fontsize": 12})
axs[0].set_ylabel("Energy", fontdict={"fontsize": 12})
axs[0].legend()

# Plot orbital occupancy
axs[1].bar(x2, y2, width=0.8)
axs[1].set_xticks(x2)
axs[1].set_xticklabels(x2)
axs[1].set_title("Avg Occupancy per Spatial Orbital")
axs[1].set_xlabel("Orbital Index", fontdict={"fontsize": 12})
axs[1].set_ylabel("Avg Occupancy", fontdict={"fontsize": 12})

print(f"Exact energy: {exact_energy:.5f}")
print(f"SQD energy: {min_e:.5f}")
print(f"Absolute error: {abs(min_e - exact_energy):.5f}")
plt.tight_layout()
plt.show()
Exact energy: -13.42249
SQD energy: -13.42249
Absolute error: 0.00000

Plot output

Цей переклад базується на англійській версії від 5 трав. 2026 р.
Проблема із сайтом чи перекладом?Проблема з вмістом? Редагувати в IBM Quantum docs