주 콘텐츠로 건너뛰기

비용 함수

이 강의에서는 비용 함수를 평가하는 방법을 배웁니다:

  • 먼저 Qiskit Runtime primitives에 대해 알아봅니다.
  • 비용 함수 C(θ)C(\vec\theta)를 정의합니다. 이것은 최적화 도구가 최소화(또는 최대화)해야 할 문제의 목표를 정의하는 문제별 함수입니다.
  • Qiskit Runtime primitives를 사용하여 속도와 정확도를 최적화하는 측정 전략 정의

 

Estimator 및 Sampler와 같은 primitives 사용을 포함한 비용 함수의 핵심 구성 요소를 보여주는 다이어그램.

Primitives

모든 물리적 시스템은 고전적이든 양자적이든 서로 다른 상태에 존재할 수 있습니다. 예를 들어, 도로 위의 자동차는 질량, 위치, 속도, 가속도 등의 물리량으로 그 상태를 특정할 수 있습니다. 마찬가지로 양자 시스템도 서로 다른 구성이나 상태를 가질 수 있지만, 측정 방식과 상태 진화 측면에서 고전 시스템과 다릅니다. 이로 인해 중첩얽힘과 같이 양자역학 고유의 특성이 나타납니다. 자동차의 상태를 속도나 가속도와 같은 물리적 특성으로 설명하듯, 양자 시스템의 상태는 수학적 객체인 *관측량(observable)*으로 설명할 수 있습니다.

양자역학에서 상태는 정규화된 복소 열벡터, 즉 켓(ket) (ψ|\psi\rangle)으로 표현되고, 관측량은 켓에 작용하는 에르미트 선형 연산자(H^=H^\hat{H}=\hat{H}^{\dagger})입니다. 관측량의 고유벡터(λ|\lambda\rangle)는 *고유상태(eigenstate)*라고 합니다. 고유상태(λ|\lambda\rangle) 중 하나에서 관측량을 측정하면 해당 고유값(λ\lambda)이 측정 결과로 나옵니다.

양자 시스템을 어떻게, 무엇을 측정할 수 있는지 궁금하다면, Qiskit이 제공하는 두 가지 primitives가 도움이 됩니다:

  • Sampler: 양자 상태 ψ|\psi\rangle가 주어지면, 가능한 각 계산 기저 상태(computational basis state)의 확률을 구합니다.
  • Estimator: 양자 관측량 H^\hat{H}와 상태 ψ|\psi\rangle가 주어지면, H^\hat{H}의 기댓값을 계산합니다.

The Sampler primitive

Sampler primitive는 상태 ψ|\psi\rangle를 준비하는 양자 Circuit이 주어졌을 때, 계산 기저의 각 가능한 상태 k|k\rangle를 얻을 확률을 계산합니다. 즉,

pk=kψ2kZ2n{0,1,,2n1},p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},

여기서 nn은 Qubit의 수이고, kk는 가능한 모든 출력 이진 문자열 {0,1}n\{0,1\}^n(즉, 2진법 정수)의 정수 표현입니다.

Qiskit Runtime Sampler는 양자 장치에서 Circuit을 여러 번 실행하고, 각 실행마다 측정을 수행한 뒤, 복원된 비트열에서 확률 분포를 재구성합니다. 실행 횟수(또는 shots)가 많을수록 결과가 더 정확해지지만, 그만큼 더 많은 시간과 양자 리소스가 필요합니다.

그러나 가능한 출력의 수는 Qubit 수 nn에 따라 지수적으로 증가(2n2^n)하므로, 밀집(dense) 확률 분포를 포착하려면 shots 수도 지수적으로 증가해야 합니다. 따라서 Sampler희소(sparse) 확률 분포에서만 효율적입니다. 즉, 대상 상태 ψ|\psi\rangle가 Qubit 수에 대해 최대 다항식으로 증가하는 항의 수로 계산 기저 상태의 선형 결합으로 표현 가능해야 합니다:

ψ=kPoly(n)wkk.|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.

Sampler는 Circuit의 일부 구간에서 확률을 가져오도록 구성할 수도 있으며, 이는 전체 가능한 상태의 부분 집합을 나타냅니다.

The Estimator primitive

Estimator primitive는 양자 상태 ψ|\psi\rangle에 대한 관측량 H^\hat{H}의 기댓값을 계산합니다. 관측량의 확률은 pλ=λψ2p_\lambda = |\langle\lambda|\psi\rangle|^2로 표현될 수 있으며, 여기서 λ|\lambda\rangle는 관측량 H^\hat{H}의 고유상태입니다. 기댓값은 상태 ψ|\psi\rangle의 측정에서 나타날 수 있는 모든 결과 λ\lambda(즉, 관측량의 고유값)의 평균으로, 해당 확률을 가중치로 사용하여 정의됩니다:

H^ψ:=λpλλ=ψH^ψ\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle

그러나 관측량의 기댓값을 항상 계산할 수 있는 것은 아닙니다. 왜냐하면 고유 기저를 모르는 경우가 많기 때문입니다. Qiskit Runtime Estimator는 복잡한 대수적 과정을 통해, 우리가 고유 기저를 알고 있는 다른 관측량들의 조합으로 관측량을 분해하여 실제 양자 장치에서 기댓값을 추정합니다.

간단히 말하면, Estimator는 측정 방법을 모르는 관측량을 파울리 연산자라고 하는 더 단순하고 측정 가능한 관측량으로 분해합니다. 어떤 연산자든 4n4^n개의 파울리 연산자의 조합으로 표현할 수 있습니다.

P^k:=σkn1σk0kZ4n{0,1,,4n1},\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\

이로부터

H^=k=04n1wkP^k\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k

여기서 nn은 Qubit의 수, kkn1k0k \equiv k_{n-1} \cdots k_0이며 klZ4{0,1,2,3}k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}(즉, 4진법 정수), (σ0,σ1,σ2,σ3):=(I,X,Y,Z)(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)입니다.

이 분해를 수행한 후, Estimator는 원래 Circuit으로부터 각 관측량 P^k\hat{P}_k에 대한 새로운 Circuit VkψV_k|\psi\rangle을 도출하여 파울리 관측량을 계산 기저에서 효과적으로 대각화하고 측정합니다. VkV_k는 미리 알 수 있기 때문에 파울리 관측량을 쉽게 측정할 수 있는 반면, 일반적인 다른 관측량에 대해서는 그렇지 않습니다.

P^k\hat{P}_{k}에 대해, Estimator는 해당 Circuit을 양자 장치에서 여러 번 실행하고, 출력 상태를 계산 기저에서 측정하며, 가능한 각 출력 jj를 얻을 확률 pkjp_{kj}를 계산합니다. 그런 다음 각 출력 jj에 대응하는 PkP_k의 고유값 λkj\lambda_{kj}를 찾아 wkw_k를 곱하고, 모든 결과를 합산하여 주어진 상태 ψ|\psi\rangle에 대한 관측량 H^\hat{H}의 기댓값을 구합니다.

H^ψ=k=04n1wkj=02n1pkjλkj,\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},

4n4^n개의 파울리에 대한 기댓값 계산은 비현실적(즉, 지수적 증가)이므로, Estimatorwkw_k 중 많은 수가 0일 때(즉, 밀집 분해가 아닌 희소 파울리 분해)에만 효율적으로 동작할 수 있습니다. 공식적으로, 이 계산이 효율적으로 풀 수 있으려면 0이 아닌 항의 수가 Qubit 수 nn에 대해 최대 다항식으로 증가해야 합니다: H^=kPoly(n)wkP^k.\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.

Sampler에서 설명한 바와 같이 확률 샘플링도 효율적이어야 한다는 암묵적인 가정이 있습니다. 즉,

H^ψ=kPoly(n)wkjPoly(n)pkjλkj.\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.

기댓값 계산 안내 예제

단일 Qubit 상태 +:=H0=12(0+1)|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)와 관측량

H^=(1221)=2XZ\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}

을 가정하고, 이론적인 기댓값 H^+=+H^+=2\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2를 구해봅니다.

이 관측량을 직접 측정하는 방법을 모르기 때문에 기댓값을 바로 계산할 수 없으므로, H^+=2X+Z+\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+로 재표현해야 합니다. +X+=1\langle+|X|+\rangle = 1이고 +Z+=0\langle+|Z|+\rangle = 0임을 이용하면 동일한 결과가 나옴을 확인할 수 있습니다.

X+\langle X \rangle_+Z+\langle Z \rangle_+를 직접 계산하는 방법을 살펴봅니다. XXZZ는 교환하지 않으므로(즉, 동일한 고유 기저를 공유하지 않음) 동시에 측정할 수 없습니다. 따라서 보조 Circuit이 필요합니다:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx
from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
aux_circ = original_circuit.copy()
aux_circ.barrier()
if str(pauli) == "X":
aux_circ.h(0)
elif str(pauli) == "Y":
aux_circ.sdg(0)
aux_circ.h(0)
else:
aux_circ.id(0)
aux_circ.measure_all()
aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

이전 코드 셀의 출력

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

이전 코드 셀의 출력

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

이전 코드 셀의 출력

이제 Sampler를 사용하여 수동으로 계산을 수행하고 Estimator에서 결과를 확인할 수 있습니다:

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np

## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
data_pub = job.result()[index].data
bitstrings = data_pub.meas.get_bitstrings()
counts = data_pub.meas.get_counts()
quasi_dist = QuasiDistribution(
{outcome: freq / shots for outcome, freq in counts.items()}
)

# Use the probabilities and known eigenvalues of Pauli operators to estimate the expectation value.
val = 0

if str(pauli) == "X":
val += -1 * quasi_dist.get(1, 0)
val += 1 * quasi_dist.get(0, 0)

if str(pauli) == "Y":
val += -1 * quasi_dist.get(1, 0)
val += 1 * quasi_dist.get(0, 0)

if str(pauli) == "Z":
val += 1 * quasi_dist.get(0, 0)
val += -1 * quasi_dist.get(1, 0)

expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
print(f" >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f" >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
*H.paulis,
H,
] # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
if obs is not H:
print(f" >> Expected value of {str(obs)}: {expval:.5f}")
else:
print(f" >> Total expected value: {expval:.5f}")
Sampler results:
>> Expected value of X: 1.00000
>> Expected value of Z: 0.00420
>> Total expected value: 1.99580
Estimator results:
>> Expected value of X: 1.00000
>> Expected value of Z: 0.00000
>> Total expected value: 2.00000

수학적 엄밀성 (선택 사항)

H^\hat{H}의 고유상태 기저에 대해 ψ|\psi\rangleψ=λaλλ|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle로 표현하면 다음이 성립합니다:

ψH^ψ=(λaλλ)H^(λaλλ)=λλaλaλλH^λ=λλaλaλλλλ=λλaλaλλδλ,λ=λaλ2λ=λpλλ\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}

대상 관측량 H^\hat{H}의 고유값이나 고유상태를 모르기 때문에 먼저 대각화를 고려해야 합니다. H^\hat{H}에르미트(Hermitian)이므로, H^=VΛV\hat{H}=V^\dagger \Lambda V를 만족하는 유니터리 변환 VV가 존재합니다. 여기서 Λ\Lambda는 대각 고유값 행렬로, jΛk=0\langle j | \Lambda | k \rangle = 0 (jkj\neq k), jΛj=λj\langle j | \Lambda | j \rangle = \lambda_j입니다.

따라서 기댓값을 다음과 같이 다시 쓸 수 있습니다:

ψH^ψ=ψVΛVψ=ψV(j=02n1jj)Λ(k=02n1kk)Vψ=j=02n1k=02n1ψVjjΛkkVψ=j=02n1ψVjjΛjjVψ=j=02n1jVψ2λj\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}

시스템이 상태 ϕ=Vψ|\phi\rangle = V |\psi\rangle에 있을 때 j| j\rangle를 측정할 확률이 pj=jϕ2p_j = |\langle j|\phi \rangle|^2임을 감안하면, 위의 기댓값은 다음과 같이 표현됩니다:

ψH^ψ=j=02n1pjλj.\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.

확률이 ψ|\psi\rangle가 아닌 VψV |\psi\rangle 상태에서 구해진다는 점에 유의하는 것이 매우 중요합니다. 이것이 바로 행렬 VV가 반드시 필요한 이유입니다. 행렬 VV와 고유값 Λ\Lambda를 어떻게 구하는지 궁금할 수 있습니다. 고유값을 이미 알고 있다면, 변분 알고리즘의 목표가 바로 H^\hat{H}의 고유값을 찾는 것이므로 양자 컴퓨터를 사용할 필요가 없을 것입니다.

다행히 우회 방법이 있습니다: 2n×2n2^n \times 2^n 행렬은 nn개의 파울리 행렬과 항등 행렬의 텐서 곱 4n4^n개의 선형 결합으로 나타낼 수 있으며, 이들은 모두 VVΛ\Lambda가 알려진 에르미트이자 유니터리 행렬입니다. 이것이 바로 Runtime의 Estimator가 내부적으로 Operator 객체를 SparsePauliOp으로 분해하여 수행하는 작업입니다.

사용 가능한 연산자는 다음과 같습니다:

OperatorσVΛIσ0=(1001)V0=IΛ0=I=(1001)Xσ1=(0110)V1=H=12(1111)Λ1=σ3=(1001)Yσ2=(0ii0)V2=HS=12(1111)(100i)=12(1i1i)Λ2=σ3=(1001)Zσ3=(1001)V3=IΛ3=σ3=(1001)\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}

이제 H^\hat{H}를 파울리와 항등 행렬로 다시 표현해봅니다:

H^=kn1=03...k0=03wkn1...k0σkn1...σk0=k=04n1wkP^k,\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,

여기서 k=l=0n14lklkn1...k0k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0 (kn1,...,k0{0,1,2,3}k_{n-1},...,k_0\in \{0,1,2,3\}, 즉 4진법)이고, P^k:=σkn1...σk0\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}입니다:

ψH^ψ=k=04n1wkj=02n1jVkψ2jΛkj=k=04n1wkj=02n1pkjλkj,\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}

여기서 Vk:=Vkn1...Vk0V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}이고 Λk:=Λkn1...Λk0\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}로, Pk^=VkΛkVk\hat{P_k}=V_k^\dagger \Lambda_k V_k를 만족합니다.

비용 함수

일반적으로 비용 함수는 문제의 목표와 시험 상태(trial state)가 그 목표에 얼마나 잘 부합하는지를 나타내는 데 사용됩니다. 이 정의는 화학, 머신러닝, 금융, 최적화 등 다양한 분야에 적용할 수 있습니다.

간단한 예로, 시스템의 기저 상태(ground state)를 찾는 경우를 생각해 봅시다. 목표는 에너지를 나타내는 관측량(Hamiltonian H^\hat{\mathcal{H}})의 기댓값을 최소화하는 것입니다:

minθψ(θ)H^ψ(θ)\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle

Estimator를 사용하여 기댓값을 계산하고, 이 값을 최적화기(optimizer)에 전달하여 최소화할 수 있습니다. 최적화가 성공하면 최적 파라미터 값 θ\vec\theta^*이 반환되고, 이를 통해 제안된 해 상태 ψ(θ)|\psi(\vec\theta^*)\rangle를 구성하고 관측된 기댓값을 C(θ)C(\vec\theta^*)로 계산할 수 있습니다.

우리가 고려하는 제한된 상태 집합에 대해서만 비용 함수를 최소화할 수 있다는 점에 주의하세요. 이는 두 가지 가능성으로 이어집니다:

  • 앤사츠(ansatz)가 탐색 공간에서 해 상태를 정의하지 못하는 경우: 이 경우 최적화기는 해를 찾지 못하며, 탐색 공간을 더 정확하게 표현할 수 있는 다른 앤사츠를 실험해 봐야 합니다.
  • 최적화기가 유효한 해를 찾지 못하는 경우: 최적화는 전역(global) 최적화와 지역(local) 최적화로 나눌 수 있습니다. 이에 대해서는 이후 섹션에서 자세히 살펴볼 것입니다.

결론적으로, 우리는 고전적인 최적화 루프를 수행하되, 비용 함수의 평가는 양자 컴퓨터에 의존하게 됩니다. 이 관점에서 보면, 최적화기가 비용 함수를 평가해야 할 때마다 어떤 블랙박스 양자 오라클을 호출하는 순수한 고전적 작업으로 최적화를 생각할 수도 있습니다.

def cost_func_vqe(params, circuit, hamiltonian, estimator):
"""Return estimate of energy from estimator

Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance

Returns:
float: Energy estimate
"""
pub = (circuit, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
return cost
from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
2,
rotation_blocks=["rz", "ry"],
entanglement_blocks="cx",
entanglement="linear",
reps=1,
)
ansatz = reference_circuit.compose(variational_form)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

Output of the previous code cell

먼저 시뮬레이터인 StatevectorEstimator를 사용하여 실행합니다. 이는 디버깅 시 일반적으로 권장되는 방법이지만, 디버깅 실행 직후 실제 양자 하드웨어에서도 계산을 수행할 것입니다. 점점 더 많은 관심 문제들이 최첨단 슈퍼컴퓨팅 시설 없이는 고전적으로 시뮬레이션이 불가능해지고 있습니다.

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)
[-0.58744589]

이제 실제 양자 컴퓨터에서 실행하겠습니다. 문법 변경에 주의하세요. pass_manager와 관련된 단계는 다음 예제에서 더 자세히 설명할 것입니다. 변분 알고리즘에서 특히 중요한 단계 중 하나는 Qiskit Runtime Session을 사용하는 것입니다. Session을 시작하면 파라미터가 업데이트될 때마다 새 대기열에서 기다리지 않고 변분 알고리즘의 여러 반복을 실행할 수 있습니다. 이는 대기열 시간이 길거나 많은 반복이 필요한 경우 중요합니다. IBM Quantum® Network 파트너만 Runtime Session을 사용할 수 있습니다. Session에 접근할 수 없는 경우, 한 번에 제출하는 반복 횟수를 줄이고 가장 최근 파라미터를 저장하여 이후 실행에 사용할 수 있습니다. 너무 많은 반복을 제출하거나 대기열 시간이 너무 길어지면 오류 코드 1217이 발생할 수 있으며, 이는 작업 제출 사이의 긴 지연을 나타냅니다.

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
QiskitRuntimeService,
Session,
EstimatorOptions,
EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)
cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

위의 두 계산에서 얻은 값이 매우 유사하다는 점에 유의하세요. 결과를 개선하는 기법들은 아래에서 더 자세히 설명합니다.

비물리적 시스템에 대한 매핑 예시

최대 절단(Max-Cut) 문제는 그래프의 꼭짓점을 두 개의 서로소인 집합으로 나누어 두 집합 사이의 간선 수를 최대화하는 조합 최적화 문제입니다. 더 공식적으로 말하면, 무방향 그래프 G=(V,E)G=(V,E)가 주어졌을 때, 여기서 VV는 꼭짓점의 집합이고 EE는 간선의 집합이며, Max-Cut 문제는 꼭짓점을 두 개의 서로소인 부분 집합 SSTT로 분할하여 한 끝점은 SS에, 다른 끝점은 TT에 있는 간선의 수를 최대화하는 것을 요구합니다.

Max-Cut은 클러스터링, 네트워크 설계, 상전이 등 다양한 문제를 해결하는 데 적용할 수 있습니다. 먼저 문제 그래프를 생성해 봅시다:

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Output of the previous code cell

이 문제는 이진 최적화 문제로 표현할 수 있습니다. 각 노드 0i<n0 \leq i < n에 대해 (여기서 nn은 그래프의 노드 수이며 이 경우 n=4n=4), 이진 변수 xix_i를 고려합니다. 이 변수는 노드 ii11로 레이블된 그룹 중 하나에 있으면 값 11을, 00으로 레이블된 다른 그룹에 있으면 00을 가집니다. 또한 wijw_{ij}를 노드 ii에서 노드 jj로 가는 간선의 가중치(인접 행렬 ww(i,j)(i,j) 원소)로 표기합니다. 그래프가 무방향이므로 wij=wjiw_{ij}=w_{ji}입니다. 그러면 다음 비용 함수를 최대화하는 문제로 공식화할 수 있습니다:

C(x)=i,j=0nwijxi(1xj)=i,j=0nwijxii,j=0nwijxixj=i,j=0nwijxii=0nj=0i2wijxixj\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}

이 문제를 양자 컴퓨터로 풀기 위해, 비용 함수를 관측량의 기댓값으로 표현해야 합니다. 그러나 Qiskit이 기본적으로 지원하는 관측량은 0011 대신 고유값 111-1을 가지는 Pauli 연산자로 구성됩니다. 따라서 다음과 같이 변수를 변환합니다:

여기서 x=(x0,x1,,xn1)\vec{x}=(x_0,x_1,\cdots ,x_{n-1})입니다. 인접 행렬 ww를 사용하면 모든 간선의 가중치에 편리하게 접근할 수 있습니다. 이를 통해 비용 함수를 다음과 같이 구합니다:

zi=12xixi=1zi2z_i = 1-2x_i \rightarrow x_i = \frac{1-z_i}{2}

이는 다음을 의미합니다:

xi=0zi=1xi=1zi=1.\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}

따라서 최대화하고자 하는 새로운 비용 함수는 다음과 같습니다:

C(z)=i,j=0nwij(1zi2)(11zj2)=i,j=0nwij4i,j=0nwij4zizj=i=0nj=0iwij2i=0nj=0iwij2zizj\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}

또한, 양자 컴퓨터는 최댓값 대신 최솟값(보통 가장 낮은 에너지)을 찾는 경향이 있으므로, C(z)C(\vec{z})를 최대화하는 대신 다음을 최소화합니다:

C(z)=i=0nj=0iwij2zizji=0nj=0iwij2-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}

이제 변수가 1-111의 값을 가질 수 있는 최소화할 비용 함수를 얻었으므로, Pauli ZZ와 다음과 같이 유사성을 도출할 수 있습니다:

ziZi=In1...Zi...I0z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}

다시 말해, 변수 ziz_i는 Qubit ii에 작용하는 ZZ Gate에 해당합니다. 또한:

Zixn1x0=zixn1x0xn1x0Zixn1x0=ziZ_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i

따라서 고려할 관측량은 다음과 같습니다:

H^=i=0nj=0iwij2ZiZj\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j

여기에 나중에 독립 항을 더해야 합니다:

offset=i=0nj=0iwij2\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}

연산자는 간선으로 연결된 노드에 Z 연산자가 있는 항들의 선형 결합입니다 (0번째 Qubit이 가장 오른쪽에 위치): IIZZ+IZIZ+IZZI+ZIIZ+ZZIIIIZZ + IZIZ + IZZI + ZIIZ + ZZII. 연산자가 구성되면 QAOA 알고리즘의 앤사츠는 Qiskit Circuit 라이브러리의 QAOAAnsatz Circuit을 사용하여 쉽게 구축할 수 있습니다.

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
[("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

Output of the previous code cell

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")
Offset: -2.5

Runtime Estimator가 Hamiltonian과 매개변수화된 앤사츠를 직접 받아 필요한 에너지를 반환하므로, QAOA 인스턴스의 비용 함수는 매우 간단합니다:

def cost_func(params, ansatz, hamiltonian, estimator):
"""Return estimate of energy from estimator

Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (Estimator): Estimator primitive instance

Returns:
float: Energy estimate
"""
pub = (ansatz, hamiltonian, params)
cost = estimator.run([pub]).result()[0].data.evs
# cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
return cost
import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)
1.473098768180865
# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)
cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)
1.1120776913677988

이 예시는 이후 응용(Applications) 섹션에서 탐색 공간을 반복하기 위해 최적화기를 활용하는 방법을 살펴볼 때 다시 다룰 것입니다. 일반적으로 다음과 같은 과정을 포함합니다:

  • 최적화기를 활용하여 최적 파라미터 찾기
  • 최적 파라미터를 앤사츠에 바인딩하여 고유값 찾기
  • 고유값을 문제 정의로 변환하기

측정 전략: 속도 대 정확도

앞서 언급했듯이, 우리는 노이즈가 있는 양자 컴퓨터를 블랙박스 오라클로 사용하고 있습니다. 노이즈는 가져온 값을 비결정적으로 만들어 무작위 변동을 일으키고, 이는 결국 특정 옵티마이저가 제안된 해로 수렴하는 것을 방해하거나 완전히 막을 수 있습니다. 이는 양자 유용성을 점진적으로 탐색하고 양자 우위를 향해 나아가면서 반드시 해결해야 할 일반적인 문제입니다.

회로 복잡도에 따른 시뮬레이션 비용의 변화를 나타낸 그래프. 고전 컴퓨터를 사용하면 지수적으로 증가합니다. 양자 오류 경감을 적용하면 이점이 생기는 교차점이 존재해야 합니다. 양자 오류 수정은 시뮬레이션 비용의 선형 성장을 가능하게 하여 확실한 이점을 제공합니다.

Qiskit Runtime Primitive의 오류 억제(error suppression) 및 오류 경감(error mitigation) 옵션을 활용하면 노이즈를 처리하고 오늘날의 양자 컴퓨터 활용도를 극대화할 수 있습니다.

오류 억제 (Error Suppression)

오류 억제는 오류를 최소화하기 위해 컴파일 과정에서 Circuit을 최적화하고 변환하는 기법을 의미합니다. 이는 기본적인 오류 처리 기법으로, 전체 런타임에 일부 고전적 전처리 오버헤드가 발생합니다. 이 오버헤드에는 양자 하드웨어에서 Circuit을 실행하기 위한 트랜스파일 과정이 포함됩니다.

  • 양자 시스템에서 사용 가능한 네이티브 Gate로 Circuit 표현
  • 가상 Qubit을 물리적 Qubit에 매핑
  • 연결성 요구 사항에 따른 SWAP 추가
  • 1Q 및 2Q Gate 최적화
  • 결어긋남(decoherence) 영향을 방지하기 위해 유휴 Qubit에 동적 디커플링(dynamical decoupling) 추가

Primitive는 optimization_level 옵션을 설정하고 고급 트랜스파일 옵션을 선택함으로써 오류 억제 기법을 사용할 수 있습니다. 이후 과정에서는 결과를 개선하기 위한 다양한 Circuit 구성 방법을 살펴볼 예정이지만, 대부분의 경우 optimization_level=3으로 설정하는 것을 권장합니다.

간단한 이상적 동작을 가진 예시 Circuit을 통해 트랜스파일 과정에서 최적화 수준을 높이는 것이 얼마나 가치 있는지 시각적으로 확인해 보겠습니다.

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

이전 코드 셀의 출력

위의 Circuit은 [0,2π][0,2\pi]와 같은 적절한 구간에 걸쳐 위상(phase)을 입력할 경우, 주어진 관측 가능량의 정현파 기댓값을 생성합니다.

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

시뮬레이터를 활용해 최적화된 트랜스파일의 유용성을 살펴봐요. 오류 경감의 유용성을 실제 하드웨어로 시연하는 것은 아래에서 다시 다룹니다. QiskitRuntimeService를 사용해 실제 Backend(여기서는 ibm_brisbane)를 가져오고, AerSimulator로 해당 Backend의 노이즈 동작을 포함한 시뮬레이션을 수행합니다.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

이제 pass manager를 사용해 Circuit을 Backend의 "명령어 집합 아키텍처(instruction set architecture, ISA)"로 트랜스파일할 수 있습니다. 이는 Qiskit Runtime의 새로운 요구 사항입니다. Backend에 제출되는 모든 Circuit은 해당 Backend의 타겟 제약 조건을 준수해야 합니다. 즉, Backend의 ISA — 디바이스가 이해하고 실행할 수 있는 명령어 집합 — 로 작성되어야 합니다. 이 타겟 제약 조건은 디바이스의 네이티브 기저 Gate, Qubit 연결성, 그리고 관련이 있는 경우 펄스 및 기타 명령어 타이밍 사양 등의 요소에 의해 정의됩니다.

현재 경우에는 이 작업을 두 번 수행합니다. 한 번은 optimization_level = 0으로, 한 번은 3으로 설정합니다. 매번 Estimator Primitive를 사용해 다양한 위상 값에서 관측 가능량의 기댓값을 추정합니다.

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc
# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

마지막으로 결과를 플롯할 수 있으며, 최적화 없이도 계산 정밀도가 상당히 양호했지만 최적화 수준을 3으로 높이면 확실히 개선되는 것을 확인할 수 있습니다. 더 깊고 복잡한 Circuit에서는 최적화 수준 0과 3의 차이가 더 크게 나타날 수 있습니다. 이는 장난감 모델로 사용된 매우 단순한 Circuit입니다.

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

이전 코드 셀의 출력

오류 완화

오류 완화란 실행 시점에 장치 노이즈를 모델링하여 회로 오류를 줄이는 기법을 말합니다. 일반적으로 이 과정에서는 모델 학습과 관련된 양자 전처리 오버헤드와, 생성된 모델을 이용해 원시 결과의 오류를 완화하는 고전 후처리 오버헤드가 발생합니다.

Qiskit Runtime 프리미티브의 resilience_level 옵션은 오류에 대한 복원력 수준을 지정합니다. 레벨이 높을수록 양자 샘플링 오버헤드로 인해 처리 시간이 길어지는 대신 더 정확한 결과를 얻을 수 있습니다. 복원력 레벨을 사용하면 프리미티브 쿼리에 오류 완화를 적용할 때 비용과 정확도 간의 균형을 조정할 수 있습니다.

오류 완화 기법을 적용하면 이전의 완화되지 않은 결과에 비해 편향(bias)이 줄어들 것으로 기대합니다. 경우에 따라서는 편향이 완전히 사라지기도 합니다. 그러나 이는 비용을 수반합니다. 추정값의 편향을 줄이면 통계적 변동성(즉, 분산)이 증가하는데, 이는 샘플링 과정에서 회로당 샷 수를 늘려 보완할 수 있습니다. 이로 인해 편향 감소에 필요한 것 이상의 추가 오버헤드가 발생하므로 기본값으로는 적용되지 않습니다. 아래 예시처럼 options.executions.shots에서 회로당 샷 수를 조정하면 이 동작을 쉽게 활성화할 수 있습니다.

A diagram showing broader or narrowing distributions as in the bias/variance tradeoff.

이 과정에서는 Qiskit Runtime 프리미티브가 전체 구현 세부 사항 없이도 수행할 수 있는 오류 완화를 설명하기 위해 이러한 오류 완화 모델을 개략적으로 살펴봅니다.

트월드 판독 오류 소거(T-REx)

트월드 판독 오류 소거(T-REx)는 Pauli 트월링(twirling)이라는 기법을 사용하여 양자 측정 과정에서 발생하는 노이즈를 줄입니다. 이 기법은 특정 노이즈 형태를 가정하지 않으므로 매우 범용적이고 효과적입니다.

전체 워크플로:

  1. 무작위 비트 플립(측정 전 Pauli X)을 적용한 영(zero) 상태 데이터 획득
  2. 무작위 비트 플립(측정 전 Pauli X)을 적용한 원하는 (노이즈 포함) 상태 데이터 획득
  3. 각 데이터 세트에 대한 특수 함수를 계산하고 나누기

 

A diagram showing measurement and calibration circuits for T-REX.

options.resilience_level = 1로 설정하면 되며, 아래 예시에서 확인할 수 있습니다.

제로 노이즈 외삽(ZNE)

제로 노이즈 외삽(ZNE)은 원하는 양자 상태를 준비하는 회로의 노이즈를 먼저 증폭하고, 여러 노이즈 수준에서 측정값을 획득한 뒤, 이를 이용해 노이즈가 없는 결과를 추론하는 방식으로 동작합니다.

전체 워크플로:

  1. 여러 노이즈 배율에 대해 회로 노이즈 증폭
  2. 노이즈 증폭된 모든 회로 실행
  3. 제로 노이즈 한계로 외삽

 

A diagram showing steps in ZNE. Noise is artificially amplified by different factors. Then the values are extrapolated to what they should be at zero noise.

options.resilience_level = 2로 설정하면 됩니다. noise_factors, noise_amplifiers, extrapolators를 다양하게 탐색하여 더욱 최적화할 수 있지만, 이는 이 과정의 범위를 벗어납니다. 여기에 설명된 옵션을 직접 실험해 보세요.

각 방법에는 고유한 오버헤드가 수반됩니다. 필요한 양자 계산 횟수(시간)와 결과 정확도 간의 균형 관계입니다.

MethodsR=1, T-RExR=2, ZNEAssumptionsNoneAbility to scale noiseQubit overhead11Sampling overhead2Nnoise-factorsBias0O(λNnoise-factors)\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}

Qiskit Runtime의 완화 및 억제 옵션 사용

다음은 Qiskit Runtime에서 오류 완화와 억제를 적용하면서 기댓값을 계산하는 방법입니다. 이전과 동일한 회로와 옵저버블을 사용하되, 이번에는 최적화 레벨을 2로 고정하고 복원력(resilience) 또는 사용되는 오류 완화 기법을 조정합니다. 이 오류 완화 과정은 최적화 루프의 매 반복마다 수행됩니다.

오류 완화는 시뮬레이터에서는 사용할 수 없으므로 이 부분은 실제 하드웨어에서 수행합니다.

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
Session,
EstimatorOptions,
EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
# operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

estimator_options = EstimatorOptions(resilience_level=0, default_shots=10_000)

with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

estimator_options = EstimatorOptions(resilience_level=2, default_shots=10_000)

with Session(backend=backend) as session:
estimator = Estimator(mode=session, options=estimator_options)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

이전과 마찬가지로 사용된 세 가지 오류 완화 레벨에 대한 기댓값을 위상각의 함수로 그래프로 나타낼 수 있습니다. 오류 완화가 결과를 약간 개선한다는 것을 확인할 수 있지만, 그 차이는 미미합니다. 다시 말하지만, 이 효과는 더 깊고 복잡한 회로에서 훨씬 더 두드러집니다.

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Output of the previous code cell

요약

이 강의에서는 비용 함수를 만드는 방법을 배웠습니다.

  • 비용 함수 생성
  • Qiskit Runtime 프리미티브를 활용하여 노이즈를 완화하고 억제하는 방법
  • 속도 대 정확도를 최적화하는 측정 전략 정의 방법

다음은 우리의 고수준 변분 워크로드입니다.

A diagram showing the quantum circuit with unitaries preparing the reference state and variational state, followed by measurements. These are used to evaluate the cost function.

비용 함수는 최적화 루프의 매 반복마다 실행됩니다. 다음 강의에서는 고전 옵티마이저가 비용 함수 평가값을 사용하여 새로운 파라미터를 선택하는 방법을 살펴봅니다.

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())
1.1.0
0.23.0