주 콘텐츠로 건너뛰기

SQD를 활용한 fermionic 격자 모델의 에너지 추정 개선

이 튜토리얼에서는 단일 불순물 Anderson 모델로 알려진 fermionic 격자 Hamiltonian의 바닥 상태에 대한 근사치를 찾기 위해, 잡음이 있는 양자 샘플을 후처리하는 방법을 보여주는 Qiskit 패턴을 구현합니다. 점점 증가하는 시간 간격에 걸쳐 16-qubit Krylov 기저 상태 집합에서 취한 샘플을 처리하기 위해 샘플 기반 양자 대각화 접근 방식을 따릅니다. 이 상태들은 시간 진화의 Trotter화(Trotterization)를 사용하여 양자 장치에서 실현됩니다. 양자 잡음의 영향을 고려하기 위해 구성 복구(configuration recovery) 기법이 사용됩니다. 좋은 초기 상태와 바닥 상태의 희소성을 가정하면, 이 접근 방식은 효율적으로 수렴하는 것으로 증명되었습니다.

이 패턴은 네 단계로 설명할 수 있습니다:

  1. 1단계: 양자 문제로 매핑

    • 바닥 상태를 추정하기 위해 점점 증가하는 시간 간격에 걸쳐 Krylov 기저 상태 집합(즉, Trotter화된 시간 진화 회로)을 생성합니다

  2. 2단계: 문제 최적화

    • Backend용으로 회로를 트랜스파일합니다

  3. 3단계: 실험 실행

    • Sampler 프리미티브를 사용하여 회로에서 샘플을 추출합니다

  4. 4단계: 결과 후처리

    • 자기 일관적 구성 복구 루프
      • 입자 수와 가장 최근 반복에서 계산된 평균 궤도 점유율에 대한 사전 지식을 사용하여 전체 비트스트링 샘플 집합을 후처리합니다
      • 복구된 비트스트링으로부터 확률적으로 하위 샘플 배치를 생성합니다
      • 샘플링된 각 부분 공간에 대해 fermionic 격자 Hamiltonian을 투영하고 대각화합니다
      • 모든 배치에서 발견된 최소 바닥 상태 에너지를 저장하고 평균 궤도 점유율을 업데이트합니다

1단계: 문제를 양자 회로로 매핑​

먼저, 7개의 배스 사이트(8개 궤도에 8개의 전자)를 가진 1차원 단일 불순물 Anderson 모델(SIAM)을 설명하는 1체 및 2체 Hamiltonian을 생성합니다. 이 모델은 금속에 내장된 자기 불순물을 설명하는 데 사용됩니다.

그런 다음 양자 Krylov 부분 공간을 생성하는 데 사용되는 16-qubit Trotter 회로를 만듭니다.

# 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

다음으로, Trotter화된 양자 회로 집합으로 양자 Krylov 부분 공간을 생성합니다. 여기서는 초기(참조) 상태 생성과 Hamiltonian의 1체 및 2체 부분에 대한 시간 진화를 위한 헬퍼를 만듭니다. 이 모델과 회로 설계 방법에 대한 자세한 설명은 논문을 참조하세요.

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],
        )

양자 Krylov 부분 공간을 지정하는 d개의 시간 진화 상태를 생성합니다. 여기서는 d=8을 선택했습니다. Krylov 기저 상태를 샘플링할 때의 오차는 d가 증가함에 따라 수렴합니다. 이 문제 인스턴스의 특성상 OrbitalRotationJW를 사용하여 1체 진화를 효율적으로 컴파일할 수 있지만, 일반적으로는 Qiskit의 PauliEvolutionGate를 사용하여 전체 Hamiltonian의 Trotter화된 시간 진화를 구현할 수 있습니다.

# 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단계: 문제 최적화​

Trotter화된 회로를 생성한 후에는 대상 하드웨어에 맞춰 최적화합니다. 최적화하기 전에 사용할 하드웨어 장치를 선택해야 합니다. 실제 장치를 에뮬레이트하기 위해 qiskit_ibm_runtime의 fake 127-qubit Backend를 사용합니다. 실제 QPU에서 실행하려면 fake Backend를 실제 Backend로 교체하면 됩니다. 자세한 내용은 Qiskit IBM Runtime 문서를 확인하세요.

from qiskit_ibm_runtime.fake_provider.backends import FakeSherbrooke

backend = FakeSherbrooke()

다음으로, Qiskit을 사용하여 회로를 대상 Backend로 트랜스파일합니다.

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단계: 실험 실행​

하드웨어 실행에 맞춰 회로를 최적화한 후에는 대상 하드웨어에서 실행하고 바닥 상태 에너지 추정을 위한 샘플을 수집할 준비가 되었습니다. 여기서는 qiskit-ibm-runtime의 SamplerV2를 사용하여 ibm_sherbrooke Backend에서 취한 잡음이 있는 샘플을 시뮬레이션합니다. 그런 다음 각 Krylov 기저 상태의 계수(counts)를 단일 계수 딕셔너리로 결합하고 가장 흔히 샘플링된 상위 20개 비트스트링을 플로팅합니다.

참고: 트랜스파일된 회로로부터 샘플을 시뮬레이션하려면 Qiskit Aer가 필요합니다.

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단계: 결과 후처리​

이제 diagonalize_fermionic_hamiltonian 함수를 사용하여 SQD 알고리즘을 실행합니다. 이 함수의 인수에 대한 설명은 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

이제 결과를 플로팅합니다. 첫 번째 플롯은 몇 번의 반복 후에 정확한 바닥 상태 에너지를 얻는다는 것을 보여줍니다.

이 예제는 위의 출력 문장에서 볼 수 있듯이 전체 Hilbert 공간을 탐색할 수 있을 만큼 작습니다. 기억해야 할 점은, 전체 Hilbert 공간의 차원은 (num_orbitals choose nelec_a) * (num_orbitals choose nelec_b)라는 것입니다. 따라서 이 문제의 경우: (8 choose 4)**2 = 4900입니다. 향상된 구성 복구와 SQD 알고리즘이 한 반복에서 다음 반복으로 중요한 구성들을 전달한다는 사실 때문에 부분 공간 차원이 증가합니다. 마지막 반복에서는 전체 Hilbert 공간에서 대각화를 수행하게 됩니다.

두 번째 플롯은 모든 배치의 솔루션에 걸친 각 공간 궤도의 평균 점유율을 보여줍니다. 높은 확률로 각 궤도에 하나의 전자가 있음을 볼 수 있습니다.

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

이 번역은 2026년 5월 5일의 영어 버전을 기반으로 합니다
사이트 또는 번역 문제?콘텐츠 문제? IBM Quantum docs에서 편집