VQE를 이용한 하이젠베르크 체인의 바닥 상태 에너지 추정

사용량 추정치: Heron 프로세서에서 약 37분 (참고: 이는 추정치일 뿐이며 실제 실행 시간은 다를 수 있습니다.)

학습 목표

이 튜토리얼을 완료하면 다음 내용을 이해할 수 있습니다:

Qiskit을 사용하여 하이젠베르크 스핀 체인을 양자 해밀토니안으로 모델링하는 방법
SPSA 최적화기를 사용하여 양자 시스템의 바닥 상태 에너지를 추정하는 방법
Qiskit Runtime 프리미티브 및 세션을 사용하여 IBM® 양자 하드웨어에서 변분 워크플로우를 실행하는 방법

사전 요구 사항

다음 주제를 미리 익혀 두시길 권장합니다:

배경

하이젠베르크 스핀 체인은 응집 물질 물리학과 양자 자성에서 가장 널리 연구된 모델 중 하나입니다. 이 모델은 최근접 이웃 스핀이 교환 상호작용을 통해 결합된 상호작용하는 양자 스핀의 1차원 격자를 설명합니다. 외부 자기장이 있는 등방성 하이젠베르크 모델의 해밀토니안은 다음과 같이 주어집니다:

$H = \sum_{\langle i,j \rangle} \left( J_x X_i X_j + J_y Y_i Y_j + J_z Z_i Z_j \right) + \sum_{i} h_i Z_i,$

여기서 $X_i$ , $Y_i$ , $Z_i$ 는 위치 $i$ 에 작용하는 파울리 연산자이고, 합산 $\langle i,j \rangle$ 는 최근접 이웃 쌍에 대해 수행되며, $J_x = J_y = J_z = 0.5$ 는 교환 결합 상수(이 튜토리얼에서는 등방성), $h_i$ 는 위치에 따른 외부 자기장을 나타냅니다. 이 튜토리얼에서는 자기장 값이 $[-1, 1]$ 구간에서 무작위로 샘플링됩니다. 아래 구현에서 "최근접 이웃" 쌍의 집합은 처음 $N$ 개의 큐비트 간 하드웨어 백엔드의 네이티브 커플링에 의해 결정되므로, 디바이스 토폴로지에 따라 엄격한 선형 체인을 형성하지 않을 수 있습니다.

이 해밀토니안의 바닥 상태 에너지를 이해하는 것은 물리학에서 근본적으로 중요합니다. 바닥 상태는 양자 상전이, 얽힘 구조, 자기 정렬에 대한 정보를 담고 있습니다. 고전적으로는 스핀 수가 증가할수록 정확한 바닥 상태 에너지를 계산하기 어렵습니다. $N$ 개의 스핀에 대해 힐베르트 공간의 차원이 $2^N$ 으로 지수적으로 확장되기 때문입니다. 이 때문에 양자 시뮬레이션의 자연스러운 후보가 됩니다.

변분 양자 고유값 풀이기(VQE)는 해밀토니안의 바닥 상태 에너지를 추정하도록 설계된 하이브리드 양자-고전 알고리즘입니다. 양자 컴퓨터에서 파라미터화된 양자 상태 $|\psi(\theta)\rangle$ (앤사츠라고 함)를 준비하고 기댓값 $\langle \psi(\theta) | H | \psi(\theta) \rangle$ 을 측정합니다. 그런 다음 고전 최적화기가 파라미터 $\theta$ 를 반복적으로 조정하여 이 에너지를 최소화하며, 이는 측정된 에너지가 항상 진정한 바닥 상태 에너지의 상한임을 보장하는 변분 원리를 활용합니다.

이 튜토리얼에서는 단일 큐비트 회전 레이어와 얽힘 게이트로 구성된 Qiskit 회로 라이브러리의 efficient_su2 앤사츠를 사용합니다. 최적화는 동시 섭동 확률적 근사(SPSA) 알고리즘을 사용하여 수행되며, 파라미터 수에 관계없이 반복당 두 번의 함수 평가만으로 기울기를 추정하기 때문에 노이즈가 있는 양자 하드웨어에 적합합니다.

요구 사항

이 튜토리얼을 시작하기 전에 다음 항목이 설치되어 있는지 확인하세요:

시각화 지원이 포함된 Qiskit SDK v2.0 이상
Qiskit Runtime v0.44 이상 (pip install qiskit-ibm-runtime)

설정

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-ibm-runtime

import numpy as np
import matplotlib.pyplot as plt
from typing import Sequence

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp
from qiskit.primitives import BaseEstimatorV2
from qiskit.circuit.library import XGate
from qiskit.circuit.library import efficient_su2
from qiskit.transpiler import PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.transpiler.passes.scheduling import (
    ALAPScheduleAnalysis,
    PadDynamicalDecoupling,
)
from qiskit_ibm_runtime import QiskitRuntimeService, Session, EstimatorV2

def visualize_results(results):
    plt.plot(results["cost_history"], lw=2)
    plt.xlabel("Number of function evaluations")
    plt.ylabel("Energy")
    plt.show()

소규모 예제

이 섹션에서는 Qiskit 패턴의 각 단계를 소규모로 살펴보며, 워크플로우를 구성하면서 핵심 구성 요소를 설명합니다.

1단계: 고전적 입력을 양자 문제로 매핑

입력: 스핀 수
출력: 하이젠베르크 체인을 모델링하는 Ansatz와 해밀토니안

10개의 스핀으로 구성된 하이젠베르크 체인을 모델링하는 Ansatz와 해밀토니안을 구성합니다. 이 단계에서는 가장 여유 있는 백엔드의 커플링 맵을 기반으로 10스핀 하이젠베르크 해밀토니안을 구성하고 efficient_su2 앤사츠를 준비합니다.

num_spins = 10
ansatz = efficient_su2(num_qubits=num_spins, reps=2)

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=num_spins, simulator=False
)

coupling = backend.target.build_coupling_map()
reduced_coupling = coupling.reduce(list(range(num_spins)))

edge_list = reduced_coupling.graph.edge_list()
ham_list = []

for edge in edge_list:
    ham_list.append(("ZZ", edge, 0.5))
    ham_list.append(("YY", edge, 0.5))
    ham_list.append(("XX", edge, 0.5))

for qubit in reduced_coupling.physical_qubits:
    ham_list.append(("Z", [qubit], np.random.random() * 2 - 1))

hamiltonian = SparsePauliOp.from_sparse_list(ham_list, num_qubits=num_spins)

ansatz.draw("mpl", style="iqp")

이전 코드 셀의 출력

2단계: 양자 하드웨어 실행을 위한 문제 최적화

입력: 추상 Circuit, 관측 가능량
출력: 선택된 QPU에 최적화된 대상 Circuit 및 관측 가능량

Qiskit의 generate_preset_pass_manager 함수를 사용하여 선택된 QPU에 맞게 Circuit의 최적화 루틴을 자동으로 생성합니다. 사전 설정 Pass Manager 중 가장 높은 최적화 수준을 제공하는 optimization_level=3을 선택합니다. 또한 디코히어런스 오류를 억제하기 위해 ALAPScheduleAnalysis 및 PadDynamicalDecoupling 스케줄링 패스를 포함합니다.

target = backend.target
pm = generate_preset_pass_manager(optimization_level=3, target=target)
pm.scheduling = PassManager(
    [
        ALAPScheduleAnalysis(durations=target.durations()),
        PadDynamicalDecoupling(
            durations=target.durations(),
            dd_sequence=[XGate(), XGate()],
            pulse_alignment=target.pulse_alignment,
        ),
    ]
)
isa_ansatz = pm.run(ansatz)
isa_observable = hamiltonian.apply_layout(isa_ansatz.layout)
isa_ansatz.draw("mpl", scale=0.6, style="iqp", fold=-1, idle_wires=False)

이전 코드 셀의 출력

3단계: Qiskit 프리미티브를 사용하여 실행

입력: 대상 Circuit 및 관측 가능량
출력: 최적화 결과

Circuit 파라미터를 최적화하여 시스템의 추정 바닥 상태 에너지를 최소화합니다. 최적화 중 비용 함수를 평가하기 위해 Qiskit Runtime의 Estimator 프리미티브를 사용합니다.

2단계에서 백엔드에 맞게 Circuit을 최적화했으므로, 최적화된 Circuit을 전달하고 skip_transpilation=True를 설정하여 Runtime 서버에서의 트랜스파일을 생략할 수 있습니다. 이 데모에서는 qiskit-ibm-runtime 프리미티브를 사용하여 QPU에서 실행합니다. qiskit statevector 기반 프리미티브로 실행하려면 Qiskit Runtime 프리미티브를 사용하는 코드 블록을 주석 처리된 블록으로 교체하세요.

이 튜토리얼에서는 기울기 기반 최적화기인 동시 섭동 확률적 근사(SPSA)를 사용합니다. 다음으로 SPSA에 대한 간략한 소개와 Qiskit v2.0을 사용하여 SPSA를 구현하는 코드를 제공합니다.

SPSA 소개

동시 섭동 확률적 근사(SPSA) [1]는 각 반복에서 단 두 번의 함수 호출만으로 전체 기울기 벡터를 근사하는 최적화 알고리즘입니다. $p$ 개의 파라미터를 최적화해야 하는 비용 함수 $f:\mathbb{R}^p\rightarrow \mathbb{R}$ 와 반복의 $i$ 번째 단계에서의 파라미터 벡터 $x_i\in \mathbb{R}^p$ 가 있을 때, 기울기를 계산하기 위해 크기 $p$ 의 무작위 벡터 $\Delta_i$ 를 생성합니다. 이 벡터의 각 요소 $\Delta_{ij}$ , $\forall$ $j\in \{1,2,...,p\}$ 는 $\{-1, 1\}$ 에서 균등하게 샘플링됩니다. 이후 무작위 벡터 $\Delta_i$ 의 각 요소에 작은 값 $c_i$ 를 곱하여 무작위 섭동을 만듭니다. 기울기는 다음과 같이 추정됩니다.

$[\nabla f(x_i)]_j \approx \frac{f(x_i + c_i \Delta_i) - f(x_i - c_i \Delta_i)}{2c_i\Delta_{ij}}.$

직관적으로, 기울기 추정 시 무작위 섭동이 적용되므로 노이즈로 인한 $f$ 의 정확한 값에서의 작은 편차를 허용하고 보정할 수 있습니다. 실제로 SPSA는 노이즈에 특히 강건하기로 알려져 있으며, 각 반복에서 단 두 번의 하드웨어 호출만 필요합니다. 따라서 변분 알고리즘 구현에 매우 선호되는 최적화기 중 하나입니다.

이 튜토리얼에서 $i$ 번째 반복의 하이퍼파라미터 $a_i$ 와 $c_i$ 는 다음과 같이 계산됩니다.

$a_i = \frac{a}{(A + i + 1)^\alpha} \quad \text{and} \quad c_i = \frac{c}{(i+1)^\gamma},$

여기서 상수 값은 $A = 30$ , $\alpha = 0.9$ , $a = 0.3$ , $c = 0.1$ , $\gamma = 0.4$ 로 [2]에서 가져옵니다. SPSA에서 좋은 성능을 끌어내려면 하이퍼파라미터의 적절한 조정이 필요합니다.

def spsa(
    fun, x0, args=(), A=30, alpha=0.9, a=0.3, c=0.1, gamma=0.4, maxiter=100
):
    nparams = len(x0)
    x = np.copy(x0)

    for i in range(maxiter):
        a_i = a / (A + i + 1) ** alpha
        c_i = c / (i + 1) ** gamma
        delta_i = np.random.choice([-1, 1], nparams)

        # two hardware calls
        eval_1 = fun(x + c_i * delta_i, *args)
        eval_2 = fun(x - c_i * delta_i, *args)

        # compute the gradient and update the parameters
        grad = (eval_1 - eval_2) / (2 * c_i) * np.reciprocal(delta_i)
        x = x - a_i * grad

    return x

def cost_func(
    params: Sequence,
    ansatz: QuantumCircuit,
    hamiltonian: SparsePauliOp,
    estimator: BaseEstimatorV2,
    cost_history_dict: dict,
) -> float:
    """Ground state energy evaluation."""
    energy = (
        estimator.run([(ansatz, hamiltonian, [params])]).result()[0].data.evs
    )

    cost_history_dict["iters"] += 1
    cost_history_dict["prev_vector"] = list(params)
    cost_history_dict["cost_history"].append(float(energy[0]))

    print(
        f"Fx Iters. done: {cost_history_dict['iters']} [Current cost: {round(energy[0], 5)}]",
        end="\r",
    )

    return energy

def solve(x0, isa_ansatz, isa_observable, maxiter=150):
    cost_history_dict = {
        "prev_vector": None,
        "iters": 0,
        "cost_history": [],
        "y_min": None,
    }

    # Evaluate the problem using a QPU via Qiskit IBM Runtime
    with Session(backend=backend) as session:
        estimator = EstimatorV2(mode=session)
        estimator.skip_transpilation = True
        estimator.options.environment.job_tags = ["TUT_HSVQE"]
        x_opt = spsa(
            cost_func,
            x0=x0,
            args=(isa_ansatz, isa_observable, estimator, cost_history_dict),
            maxiter=maxiter,
        )

        y_min = cost_func(
            x_opt, isa_ansatz, isa_observable, estimator, cost_history_dict
        )

    return y_min, cost_history_dict

np.random.seed(42)
num_params = ansatz.num_parameters
params = 2 * np.pi * np.random.random(num_params)

여기서 maxiter = 50으로 설정합니다. 각 반복에서 기울기를 계산하기 위해 함수를 두 번 호출하므로, 총 함수 호출 횟수는 $2 \times \text{maxiter}$ 가 됩니다. 더 나은 에너지 추정을 위해 maxiter를 더 높은 값으로 늘릴 수 있습니다.

maxiter = 50
spsa_min, spsa_history = solve(
    params, isa_ansatz, isa_observable, maxiter=maxiter
)

Fx Iters. done: 101 [Current cost: -3.03843]

4단계: 후처리 및 원하는 고전적 형식으로 결과 반환

입력: 최적화 중 바닥 상태 에너지 추정값
출력: 추정된 바닥 상태 에너지

print(f"Estimated ground state energy: {spsa_min}")

Estimated ground state energy: [-3.03842968]

results = {
    "spsa": spsa_history,
}

visualize_results(spsa_history)

이전 코드 셀의 출력

대규모 하드웨어 예제

이 튜토리얼에는 대규모 하드웨어 예제가 포함되어 있지 않습니다. 큐비트 수가 증가하면 VQE는 척박한 고원(barren plateau) 현상으로 인해 심각한 어려움을 겪습니다. 시스템 크기가 커질수록 비용 함수의 기울기가 지수적으로 소멸하여, 대형 회로에서의 최적화가 사실상 불가능해집니다. 하드웨어 노이즈와 결합되면, 더 큰 스핀 체인으로 VQE를 확장해도 재현 가능한 결과를 얻기 어렵습니다. 이러한 한계를 극복하는 접근 방식은 아래 다음 단계 섹션을 참조하세요.

도전 과제

하이젠베르크 체인에 대한 VQE 구현이 완성되었으니 다음을 시도해 보세요:

앤사츠 깊이 실험: efficient_su2의 reps 파라미터를 수정해 보세요(예: reps=1과 reps=3 시도). 앤사츠 깊이가 추정 바닥 상태 에너지와 수렴 속도에 어떤 영향을 미치나요? 어느 시점에서 수익 감소 또는 불안정성이 관찰되나요?
SPSA 하이퍼파라미터 튜닝: 학습률 스케줄 파라미터(a, c, alpha, gamma, A)를 조정하고 수렴에 어떤 영향을 미치는지 관찰하세요. 여기서 사용된 기본값보다 더 빠르게 수렴하는 설정을 찾을 수 있나요?
커플링 토폴로지 비교: 백엔드의 네이티브 커플링 맵을 사용하는 대신, 간단한 최근접 이웃 선형 체인을 구성하고 결과를 비교해 보세요. 물리적 하드웨어의 연결성이 트랜스파일된 회로 깊이와 최종 에너지 추정에 어떤 영향을 미치나요?

참고 문헌

[1] Spall, J. C. (2002). Implementation of the simultaneous perturbation algorithm for stochastic optimization. IEEE Transactions on Aerospace and Electronic Systems, 34(3), 817-823.

[2] Sahin, M. Emre, et al. (2025). Qiskit Machine Learning: an open-source library for quantum machine learning tasks at scale on quantum hardware and classical simulators. arXiv:2505.17756.

학습 목표​

사전 요구 사항​

배경​

요구 사항​

설정​

소규모 예제​

1단계: 고전적 입력을 양자 문제로 매핑​

2단계: 양자 하드웨어 실행을 위한 문제 최적화​

3단계: Qiskit 프리미티브를 사용하여 실행​

SPSA 소개​

4단계: 후처리 및 원하는 고전적 형식으로 결과 반환​

대규모 하드웨어 예제​

도전 과제​

참고 문헌​

다음 단계​

학습 목표

사전 요구 사항

배경

요구 사항

설정

소규모 예제

1단계: 고전적 입력을 양자 문제로 매핑

2단계: 양자 하드웨어 실행을 위한 문제 최적화

3단계: Qiskit 프리미티브를 사용하여 실행

SPSA 소개

4단계: 후처리 및 원하는 고전적 형식으로 결과 반환

대규모 하드웨어 예제

도전 과제

참고 문헌

다음 단계