주 콘텐츠로 건너뛰기

양자 알고리즘: 변분 양자 알고리즘

참고

Takashi Imamichi (2024년 5월 24일)

원본 강의의 PDF를 다운로드하세요. 일부 코드 스니펫은 정적 이미지이므로 더 이상 사용되지 않을 수 있습니다.

이 실험을 실행하는 데 소요되는 QPU 시간은 약 9분입니다 (Eagle 프로세서 기준).

(이 노트북은 Open Plan에서 허용된 시간 내에 실행되지 않을 수 있습니다. 양자 컴퓨팅 자원을 현명하게 사용해 주세요.)

1. 소개

이 튜토리얼은 하이브리드 양자-고전 알고리즘에 대한 개요를 제공하며, 특히 변분 양자 고유값 해석기(VQE)와 양자 근사 최적화 알고리즘(QAOA)에 초점을 맞춥니다. 이 알고리즘들의 주요 목적은 매개변수화된 양자 게이트를 사용하는 양자 Circuit을 통해 최적화 문제를 해결하는 것입니다.

양자 컴퓨팅의 발전에도 불구하고, 현재 양자 장치의 잡음 존재로 인해 깊은 양자 Circuit에서 의미 있는 결과를 추출하기 어렵습니다. 이 과제를 극복하기 위해 VQE와 QAOA는 하이브리드 양자-고전 방식을 채택하며, 이는 양자 계산을 통해 비교적 짧은 양자 Circuit을 반복적으로 실행하고 고전 계산을 통해 목표 매개변수화 양자 Circuit의 매개변수를 최적화하는 과정을 포함합니다.

QAOA는 다양한 오류 완화 및 억제 기술의 적용 덕분에 유틸리티 규모에서 목표 문제에 대한 최적 해를 제공할 잠재력을 가지고 있습니다. VQE는 확장성이 낮은 양자 화학과 같은 다양한 응용 분야를 가지고 있습니다. 그러나 VQE를 보완하고 보강하기 위해 Krylov 부분 공간 대각화와 샘플링 기반 양자 대각화(SQD) 등 여러 고유값 관련 방식들이 등장했습니다. VQE를 이해하는 것은 지금까지 등장한 광범위한 고전-양자 하이브리드 알고리즘을 이해하는 데 중요한 첫 번째 단계입니다.

이 모듈은 VQE와 QAOA의 기본 개념과 구현을 설명합니다. 이후 튜토리얼에서는 이 알고리즘들을 확장하기 위한 고급 주제와 기술을 다룰 예정입니다. 이 노트북을 실행하려면 환경에 다음 라이브러리가 필요합니다. 아직 설치하지 않은 경우, 다음 셀의 주석을 해제하고 실행하여 설치할 수 있습니다.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-ibm-runtime rustworkx scipy
# % pip install 'qiskit[visualization]' qiskit-ibm-runtime

2. 간단한 해밀토니안의 최솟값 고유값 계산

VQE가 어떻게 작동하는지 확인하기 위해 매우 간단한 경우에 적용해 보겠습니다. VQE로 Pauli ZZ 행렬의 최솟값 고유값을 계산합니다. 먼저 일반 패키지를 가져오는 것으로 시작합니다.

import numpy as np
from qiskit.circuit import ParameterVector, QuantumCircuit
from qiskit.primitives import StatevectorEstimator, StatevectorSampler
from qiskit.quantum_info import SparsePauliOp
from scipy.optimize import minimize

이제 관심 있는 연산자를 정의하고 행렬 형태로 확인합니다.

op = SparsePauliOp("Z")
op.to_matrix()
array([[ 1.+0.j,  0.+0.j],
[ 0.+0.j, -1.+0.j]])

고전적으로 고유값을 쉽게 구할 수 있으므로 결과를 검증할 수 있습니다. 유틸리티 규모로 확장될수록 이것이 어려워질 수 있습니다. 여기서는 numpy를 사용합니다.

# compute eigenvalues with numpy
result = np.linalg.eigh(op.to_matrix())
print("Eigenvalues:", result.eigenvalues)
Eigenvalues: [-1.  1.]

변분 양자 알고리즘을 사용하여 고유값을 구하려면 변분 매개변수를 받는 게이트로 Circuit을 구성합니다:

# define a variational form
param = ParameterVector("a", 3)
qc = QuantumCircuit(1, 1)
qc.u(param[0], param[1], param[2], 0)
qc_estimator = qc.copy()
qc.measure(0, 0)
qc.draw("mpl")

Output of the previous code cell

연산자(예: ZZ)의 기댓값을 추정하려면 Estimator를 사용해야 합니다. 시스템의 상태를 확인하려면 Sampler를 사용합니다.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Sampler를 사용하여 랜덤 매개변수 값 [1, 2, 3]으로 비트 문자열 0과 1의 카운트를 계산할 수 있습니다.

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3])]).result()
counts = result[0].data.c.get_counts()
counts
{'0': 783, '1': 241}

확률 {0:p0,1:p1}\{0: p_0, 1: p_1\}으로 Z=p0p1\langle Z \rangle = p_0 - p_1을 통해 Z의 기댓값을 계산할 수 있음을 알고 있습니다.

# compute the expectation value of Z based on the counts
(counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
0.529296875

이 Circuit은 작동했지만 선택된 매개변수 값이 매우 낮은 에너지(또는 낮은 고유값) 상태에 해당하지 않았습니다. 얻은 고유값은 최솟값보다 상당히 높습니다. Estimator를 사용할 때도 결과는 비슷합니다.

Estimator는 측정 없는 양자 Circuit을 받는다는 점에 유의하세요.

result = estimator.run([(qc_estimator, op, [1, 2, 3])]).result()
result[0].data.evs
array(0.54030231)

매개변수를 탐색하여 가장 낮은 고유값을 산출하는 매개변수를 찾아야 합니다. 변분 형식의 매개변수 값을 받아 기댓값 Z\langle Z \rangle을 반환하는 함수를 만듭니다.

# define a cost function to look for the minimum eigenvalue of Z
def cost(x):
result = sampler.run([(qc, x)]).result()
counts = result[0].data.c.get_counts()
expval = (counts.get("0", 0) - counts.get("1", 0)) / sum(counts.values())
# the following line shows the trajectory of the optimization
print(expval, counts)
return expval

SciPy의 minimize 함수를 적용하여 Z의 최솟값 고유값을 찾아보겠습니다.

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result
1.0 {'0': 1024}
0.494140625 {'0': 765, '1': 259}
0.466796875 {'0': 751, '1': 273}
0.564453125 {'0': 801, '1': 223}
-0.4296875 {'1': 732, '0': 292}
-0.984375 {'1': 1016, '0': 8}
-0.8984375 {'1': 972, '0': 52}
-0.990234375 {'1': 1019, '0': 5}
-0.892578125 {'1': 969, '0': 55}
-0.986328125 {'1': 1017, '0': 7}
-0.861328125 {'1': 953, '0': 71}
-1.0 {'1': 1024}
-0.982421875 {'1': 1015, '0': 9}
-0.99609375 {'1': 1022, '0': 2}
-0.986328125 {'1': 1017, '0': 7}
-1.0 {'1': 1024}
-0.990234375 {'1': 1019, '0': 5}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-0.99609375 {'1': 1022, '0': 2}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.998046875 {'1': 1023, '0': 1}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-0.998046875 {'1': 1023, '0': 1}
-0.994140625 {'1': 1021, '0': 3}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
-1.0 {'1': 1024}
message: Optimization terminated successfully.
success: True
status: 1
fun: -1.0
x: [ 3.182e+00 1.338e+00 1.664e-01]
nfev: 63
maxcv: 0.0
# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()
{'0': 1, '1': 1023}

2.1 연습 문제

VQE를 사용하여 ZZZ \otimes Z의 최솟값 고유값을 계산하세요.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())
SparsePauliOp(['ZZ'],
coeffs=[1.+0.j])
[[ 1.+0.j 0.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j -1.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j -1.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j 0.+0.j 1.+0.j]]
# compute eigenvalues with numpy
# define a variational form
# qc = ...
# compute counts of bitstrings with a random parameter values by Sampler
# result = sampler.run(...)
# result
# compute the expectation value of ZZ based on the counts
# verify the expectation value of ZZ with Estimator
# define a cost function to look for the minimum eigenvalue of ZZ
# def cost(x):
# expval = ...
# return expval
# minimize the cost function with scipy's minimize
# min_result = minimize(cost, [...], method="COBYLA", tol=1e-8)
# min_result
# check counts of bitstrings with the optimal parameter values
# result = sampler.run(qc, min_result.x).result()
# result

연습 문제 풀이

관심 대상 연산자를 정의하고 행렬 형태로 확인합니다.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())
SparsePauliOp(['ZZ'],
coeffs=[1.+0.j])
[[ 1.+0.j 0.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j -1.+0.j 0.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j -1.+0.j 0.+0.j]
[ 0.+0.j 0.+0.j 0.+0.j 1.+0.j]]

변분 양자 알고리즘으로 고유값을 구하기 위해, 변분 매개변수를 받는 Gate들로 구성된 Circuit을 구성합니다.

# define a variational form
param = ParameterVector("a", 6)
qc = QuantumCircuit(2, 2)
qc.u(param[0], param[1], param[2], 0)
qc.u(param[3], param[4], param[5], 1)
qc_estimator = qc.copy()
qc.measure([0, 1], [0, 1])
qc.draw("mpl")

Output of the previous code cell

연산자(ZZZ \otimes Z 등)의 기댓값을 추정하려면 Estimator를 사용하고, 시스템의 상태를 관찰하려면 Sampler를 사용합니다.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()
# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3, 4, 5, 6])]).result()
counts = result[0].data.c.get_counts()
counts
{'10': 661, '11': 203, '01': 47, '00': 113}
# compute the expectation value of ZZ based on the counts
(
counts.get("00", 0)
- counts.get("01", 0)
- counts.get("10", 0)
+ counts.get("11", 0)
) / sum(counts.values())
-0.3828125

이 Circuit은 잘 작동하였지만, 선택한 매개변수 값이 에너지(또는 고유값)가 매우 낮은 상태에 해당하지는 않습니다. 얻어진 고유값은 최솟값보다 상당히 높습니다. Estimator를 사용한 결과도 유사합니다.

# verify the expectation value of ZZ with Estimator
result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()
result[0].data.evs
array(-0.35316516)

매개변수를 탐색하여 가장 낮은 고유값을 산출하는 값을 찾아야 합니다.

# define a cost function to look for the minimum eigenvalue of ZZ
def cost(x):
result = sampler.run([(qc, x)]).result()
counts = result[0].data.c.get_counts()
expval = (
counts.get("00", 0)
- counts.get("01", 0)
- counts.get("10", 0)
+ counts.get("11", 0)
) / sum(counts.values())
print(expval, counts)
return expval
# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0, 0, 0, 0], method="COBYLA", tol=1e-8)
min_result
1.0 {'00': 1024}
0.578125 {'00': 808, '01': 216}
0.5234375 {'00': 780, '01': 244}
0.548828125 {'00': 793, '01': 231}
0.3515625 {'00': 637, '10': 164, '11': 55, '01': 168}
0.3359375 {'00': 638, '11': 46, '10': 174, '01': 166}
0.283203125 {'00': 602, '10': 181, '01': 186, '11': 55}
-0.087890625 {'01': 414, '00': 184, '10': 143, '11': 283}
0.236328125 {'10': 27, '11': 623, '01': 364, '00': 10}
-0.0625 {'11': 261, '01': 403, '00': 219, '10': 141}
0.248046875 {'01': 366, '11': 628, '00': 11, '10': 19}
-0.0625 {'10': 145, '11': 254, '01': 399, '00': 226}
0.228515625 {'01': 373, '11': 609, '00': 20, '10': 22}
0.0546875 {'11': 376, '10': 273, '01': 211, '00': 164}
-0.447265625 {'01': 731, '10': 10, '11': 267, '00': 16}
-0.71484375 {'01': 871, '11': 99, '00': 47, '10': 7}
-0.46484375 {'01': 741, '00': 253, '10': 9, '11': 21}
-0.87890625 {'01': 962, '00': 39, '11': 23}
-0.640625 {'00': 176, '01': 837, '11': 8, '10': 3}
-0.88671875 {'01': 966, '00': 41, '11': 17}
-0.994140625 {'01': 1021, '11': 3}
-0.91796875 {'01': 982, '11': 35, '00': 7}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.939453125 {'01': 993, '00': 31}
-0.990234375 {'01': 1019, '11': 5}
-0.90234375 {'01': 974, '00': 21, '11': 29}
-0.98046875 {'01': 1014, '11': 10}
-0.994140625 {'01': 1021, '00': 3}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.98828125 {'01': 1018, '11': 6}
-0.990234375 {'01': 1019, '11': 4, '00': 1}
-0.994140625 {'01': 1021, '11': 2, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '00': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '00': 1, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '00': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.99609375 {'01': 1022, '11': 1, '00': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-0.994140625 {'01': 1021, '00': 3}
-0.998046875 {'01': 1023, '00': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '00': 1}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.998046875 {'01': 1023, '11': 1}
-0.99609375 {'01': 1022, '11': 2}
-1.0 {'01': 1024}
-0.998046875 {'01': 1023, '11': 1}
message: Optimization terminated successfully.
success: True
status: 1
fun: -0.998046875
x: [ 3.167e+00 6.940e-01 1.033e+00 -2.894e-02 8.933e-01
1.885e+00]
nfev: 128
maxcv: 0.0
message: Optimization terminated successfully.
success: True
status: 1
fun: -0.99609375
x: [ 3.098e+00 -5.402e-01 1.091e+00 -1.004e-02 3.615e-01
6.913e-01]
nfev: 115
maxcv: 0.0

numpy에서 알려준 최솟값과 매우 가까운 고유값을 얻었습니다.

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()
{'01': 1024}

3. Qiskit 패턴을 활용한 양자 최적화

이 실습에서는 Qiskit 패턴과 양자 근사 최적화에 대해 배웁니다. Qiskit 패턴이란 양자 컴퓨팅 워크플로를 구현하기 위한 직관적이고 반복 가능한 일련의 단계입니다: "Qiskit function" 조합 최적화 문제에 이 패턴을 적용하여, 하이브리드(양자-고전) 반복 방법인 **양자 근사 최적화 알고리즘(QAOA)**을 사용해 Max-Cut 문제를 푸는 방법을 보여주겠습니다.

이 QAOA 부분은 양자 근사 최적화 알고리즘 튜토리얼의 "Part 1: Small-scale QAOA"를 기반으로 합니다. 이를 확장하는 방법은 해당 튜토리얼을 참조하세요.

3.1 (소규모) 최적화를 위한 Qiskit 패턴

이 부분에서는 소규모 Max-Cut 문제를 사용하여 양자 컴퓨터로 최적화 문제를 푸는 데 필요한 단계를 설명합니다. Max-Cut 문제는 풀기 어려운 최적화 문제로(구체적으로는 NP-hard 문제), 클러스터링, 네트워크 과학, 통계 물리학 등 다양한 분야에 응용됩니다. 이 튜토리얼에서는 노드가 엣지로 연결된 그래프를 다루며, 자르는 엣지의 수를 최대화하도록 노드를 두 집합으로 분할하는 것을 목표로 합니다. "Maxcut" 이 문제를 양자 알고리즘에 매핑하기 전에, Max-Cut 문제가 어떻게 고전적인 조합 최적화 문제가 되는지 먼저 함수 f(x)f(x)의 최소화를 고려해 이해할 수 있습니다.

minx{0,1}nf(x),\min_{x\in \{0, 1\}^n}f(x),

여기서 입력 xx는 각 성분이 그래프의 각 노드에 대응하는 벡터입니다. 그런 다음 각 성분을 00 또는 11로 제한합니다(이는 컷에 포함되거나 포함되지 않음을 나타냅니다). 이 소규모 예제는 n=5n=5개의 노드를 가진 그래프를 사용합니다.

노드 쌍 i,ji, j에 대해, 대응하는 엣지 (i,j)(i,j)가 컷에 포함되는지를 나타내는 함수를 작성할 수 있습니다. 예를 들어, 함수 xi+xj2xixjx_i + x_j - 2 x_i x_jxix_i 또는 xjx_j 중 하나만 1일 때(즉, 해당 엣지가 컷에 있을 때)만 1이 되고, 그 외에는 0입니다. 컷의 엣지 수를 최대화하는 문제는 다음과 같이 정식화할 수 있습니다.

maxx{0,1}n(i,j)xi+xj2xixj,\max_{x\in \{0, 1\}^n} \sum_{(i,j)} x_i + x_j - 2 x_i x_j,

이는 다음과 같은 최소화 형태로 다시 쓸 수 있습니다.

minx{0,1}n(i,j)2xixjxixj.\min_{x\in \{0, 1\}^n} \sum_{(i,j)} 2 x_i x_j - x_i - x_j.

이 경우 f(x)f(x)의 최솟값은 컷이 통과하는 엣지의 수가 최대일 때입니다. 보면 알 수 있듯이, 아직 양자 컴퓨팅과 관련된 내용은 없습니다. 이 문제를 양자 컴퓨터가 이해할 수 있는 형태로 재정식화해야 합니다. n=5n=5개의 노드를 가진 그래프를 생성하여 문제를 초기화합니다.

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import rustworkx as rx
from rustworkx.visualization import mpl_draw
n = 5

graph = rx.PyGraph()
graph.add_nodes_from(range(1, n + 1))
edge_list = [
(0, 1, 1.0),
(0, 2, 1.0),
(1, 2, 1.0),
(1, 3, 1.0),
(2, 4, 1.0),
(3, 4, 1.0),
]
graph.add_edges_from(edge_list)
pos = rx.spring_layout(graph, seed=2)
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str)

Output of the previous code cell

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

패턴의 첫 번째 단계는 고전적인 문제(그래프)를 양자 Circuit연산자로 매핑하는 것입니다. 이를 위해 다음 세 가지 주요 단계를 수행합니다:

  1. 일련의 수학적 재정식화를 활용하여 이 문제를 이차 제약 없는 이진 최적화(QUBO) 문제 표기법으로 표현합니다.
  2. 최적화 문제를 기저 상태가 비용 함수를 최소화하는 해에 대응하는 해밀토니안으로 재작성합니다.
  3. 양자 어닐링과 유사한 과정을 통해 이 해밀토니안의 기저 상태를 준비할 양자 Circuit을 생성합니다.

참고: QAOA 방법론에서 최종적으로는 하이브리드 알고리즘의 비용 함수를 나타내는 연산자(해밀토니안)와, 문제에 대한 후보 해를 가진 양자 상태를 나타내는 매개변수화된 Circuit(Ansatz)이 필요합니다. 이 후보 상태에서 샘플링하고 비용 함수를 사용해 평가할 수 있습니다.

그래프 → 최적화 문제

매핑의 첫 번째 단계는 표기법을 변경하는 것입니다. 다음은 QUBO 표기법으로 문제를 표현한 것입니다:

minx{0,1}nxTQx,\min_{x\in \{0, 1\}^n}x^T Q x,

여기서 QQ는 실수로 이루어진 n×nn\times n 행렬, nn은 그래프의 노드 수, xx는 위에서 소개한 이진 변수 벡터, xTx^T는 벡터 xx의 전치를 나타냅니다.

Problem name: maxcut

Minimize
2*x_1*x_2 + 2*x_1*x_3 + 2*x_2*x_3 + 2*x_2*x_4 + 2*x_3*x_5 + 2*x_4*x_5 - 2*x_1
- 3*x_2 - 3*x_3 - 2*x_4 - 2*x_5

Subject to
No constraints

Binary variables (5)
x_1 x_2 x_3 x_4 x_5

최적화 문제 → 해밀토니안

그런 다음 QUBO 문제를 해밀토니안(시스템의 에너지를 나타내는 행렬)으로 재정식화할 수 있습니다:

HC=ijQijZiZj+ibiZi.H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

QAOA 문제에서 해밀토니안으로의 재정식화 단계

QAOA 문제를 이런 방식으로 재작성하는 방법을 보여 드리기 위해, 먼저 이진 변수 xix_i를 다음 식을 통해 새로운 변수 zi{1,1}z_i\in\{-1, 1\}로 대체합니다:

xi=1zi2.x_i = \frac{1-z_i}{2}.

여기서 xix_i00이면 ziz_i11이 되어야 합니다. 최적화 문제(xTQxx^TQx)에서 xix_iziz_i로 치환하면 동등한 정식화를 얻을 수 있습니다.

xTQx=ijQijxixj=14ijQij(1zi)(1zj)=14ijQijzizj14ij(Qij+Qji)zi+n24.x^TQx=\sum_{ij}Q_{ij}x_ix_j \\ =\frac{1}{4}\sum_{ij}Q_{ij}(1-z_i)(1-z_j) \\=\frac{1}{4}\sum_{ij}Q_{ij}z_iz_j-\frac{1}{4}\sum_{ij}(Q_{ij}+Q_{ji})z_i + \frac{n^2}{4}.

이제 bi=j(Qij+Qji)b_i=-\sum_{j}(Q_{ij}+Q_{ji})로 정의하고, 전치인수와 상수항 n2n^2을 제거하면, 동일한 최적화 문제의 두 가지 동등한 정식화에 도달합니다.

minx{0,1}nxTQxminz{1,1}nzTQz+bTzmin_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz

여기서 bbQQ에 의존합니다. zTQz+bTzz^TQz + b^Tz를 얻기 위해 최적화에 영향을 미치지 않는 1/4 인수와 상수 오프셋 n2n^2을 생략했음에 유의하세요.

이제 문제의 양자적 정식화를 얻기 위해, ziz_i 변수를 다음과 같은 형태의 2×22\times 2 행렬인 파울리 ZZ 행렬로 승격합니다:

Zi=(1001).Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.

이 행렬을 위의 최적화 문제에 대입하면 다음 해밀토니안을 얻습니다:

HC=ijQijZiZj+ibiZi.H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

ZZ 행렬은 양자 컴퓨터의 계산 공간, 즉 크기 2n×2n2^n\times 2^n의 힐베르트 공간에 내장되어 있음을 기억하세요. 따라서 ZiZjZ_iZ_j와 같은 항은 2n×2n2^n\times 2^n 힐베르트 공간에 내장된 텐서 곱 ZiZjZ_i\otimes Z_j로 이해해야 합니다. 예를 들어, 결정 변수가 5개인 문제에서 Z1Z3Z_1Z_3 항은 IZ3IZ1II\otimes Z_3\otimes I\otimes Z_1\otimes I를 의미하며, 여기서 II2×22\times 2 항등 행렬입니다.

이 해밀토니안을 비용 함수 해밀토니안이라고 합니다. 이 해밀토니안은 기저 상태가 비용 함수 f(x)f(x)를 최소화하는 해에 대응한다는 특성을 갖습니다. 따라서 최적화 문제를 풀기 위해서는 이제 양자 컴퓨터에서 HCH_C의 기저 상태(또는 기저 상태와 높은 중첩을 가진 상태)를 준비해야 합니다. 그러면 이 상태에서 샘플링하면 높은 확률로 min f(x)\min~f(x)의 해를 얻을 수 있습니다.

def build_max_cut_operator(graph: rx.PyGraph) -> tuple[SparsePauliOp, float]:
sp_list = []
constant = 0
for s, t in graph.edge_list():
w = graph.get_edge_data(s, t)
sp_list.append(("ZZ", [s, t], w / 2))
constant -= 1 / 2
return SparsePauliOp.from_sparse_list(
sp_list, num_qubits=graph.num_nodes()
), constant
cost_hamiltonian, constant = build_max_cut_operator(graph)
print("Cost Function Hamiltonian:", cost_hamiltonian)
print("Constant:", constant)
Cost Function Hamiltonian: SparsePauliOp(['IIIZZ', 'IIZIZ', 'IIZZI', 'IZIZI', 'ZIZII', 'ZZIII'],
coeffs=[0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j, 0.5+0.j])
Constant: -3.0

해밀토니안 → 양자 Circuit

해밀토니안 HCH_C에는 문제의 양자적 정의가 담겨 있습니다. 이제 양자 컴퓨터에서 좋은 해를 샘플링하는 데 도움이 될 양자 Circuit을 만들 수 있습니다. QAOA는 양자 어닐링에서 영감을 받아 양자 Circuit에서 연산자를 번갈아 가며 레이어로 적용합니다.

일반적인 아이디어는 알려진 시스템의 기저 상태인 Hn0H^{\otimes n}|0\rangle에서 시작하여, 관심 있는 비용 연산자의 기저 상태로 시스템을 유도하는 것입니다. 이는 각도 γ1,...,γp\gamma_1,...,\gamma_pβ1,...,βp \beta_1,...,\beta_p~를 가진 연산자 exp{iγkHC}\exp\{-i\gamma_k H_C\}exp{iβkHm}\exp\{-i\beta_k H_m\}를 적용함으로써 이루어집니다.

생성된 양자 Circuit은 γi\gamma_iβi\beta_i매개변수화되어 있어, 다양한 γi\gamma_iβi\beta_i 값을 시도하고 결과 상태에서 샘플링할 수 있습니다. "QAOA circuit diagram" 이 경우, 두 개의 매개변수 γ1\gamma_1β1\beta_1을 포함하는 1개의 QAOA 레이어를 사용한 예제를 시도해 보겠습니다.

from qiskit.circuit.library import QAOAAnsatz
circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

Output of the previous code cell

circuit.decompose(reps=3).draw("mpl", fold=-1)

Output of the previous code cell

circuit.parameters
ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 Step 2. 양자 하드웨어 실행을 위한 Circuit 최적화

위의 Circuit에는 양자 알고리즘을 사고하는 데 유용한 여러 추상화 요소가 포함되어 있지만, 실제 하드웨어에서는 바로 실행할 수 없습니다. QPU에서 실행하려면 Circuit이 패턴의 트랜스파일(transpilation) 또는 Circuit 최적화 단계를 구성하는 일련의 연산을 거쳐야 합니다.

Qiskit 라이브러리는 다양한 Circuit 변환을 지원하는 일련의 **트랜스파일 패스(transpilation passes)**를 제공합니다. 목적에 맞게 Circuit이 최적화되어 있는지 확인해야 합니다.

트랜스파일에는 다음과 같은 여러 단계가 포함될 수 있습니다.

  • Circuit 내 Qubit(의사결정 변수 등)을 디바이스의 물리적 Qubit에 초기 매핑합니다.
  • 양자 Circuit의 명령어를 Backend가 이해하는 하드웨어 네이티브 명령어로 **언롤링(unrolling)**합니다.
  • 서로 상호작용하는 Qubit를 인접한 물리적 Qubit에 배치하도록 **라우팅(routing)**합니다.
  • 동적 디커플링(dynamical decoupling)을 위한 단일 Qubit Gate를 추가하여 노이즈를 억제하는 **오류 억제(error suppression)**를 수행합니다.

트랜스파일에 대한 자세한 내용은 문서를 참고하세요.

아래 코드는 Qiskit IBM® Runtime 서비스를 통해 클라우드에서 접근 가능한 디바이스 중 하나에서 실행할 수 있는 형식으로 추상 Circuit을 변환하고 최적화합니다.

실제 양자 컴퓨터로 전송하기 전에 "로컬 테스트 모드(local testing mode)"를 사용하여 프로그램을 로컬에서 테스트할 수 있습니다. 로컬 테스트 모드에 대한 자세한 내용은 문서를 참고하세요.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Use a quantum device
service = QiskitRuntimeService()
backend = service.least_busy(min_num_qubits=127)
# backend = service.backend("ibm_kingston")

# You can test your programs locally with a fake backend (local testing mode)
# backend = FakeBrisbane()

print(backend)

# Create pass manager for transpilation
pm = generate_preset_pass_manager(optimization_level=3, backend=backend)

candidate_circuit = pm.run(circuit)
candidate_circuit.draw("mpl", fold=False, idle_wires=False)
service = QiskitRuntimeService(channel="ibm_quantum_platform")
<IBMBackend('ibm_strasbourg')>

Output of the previous code cell

3.4 Step 3. Qiskit 프리미티브를 사용한 실행

QAOA 워크플로에서는 일련의 Circuit 평가를 실행하고 고전적 최적화기를 사용하여 최적의 βk\beta_kγk\gamma_k 매개변수를 찾는 반복 최적화 루프를 통해 최적의 QAOA 매개변수를 구합니다. 이 실행 루프는 다음 단계를 통해 수행됩니다.

  1. 초기 매개변수를 정의합니다.
  2. 최적화 루프와 Circuit 샘플링에 사용할 프리미티브를 포함하는 새 Session을 인스턴스화합니다.
  3. 최적의 매개변수 집합을 찾으면, 사후 처리 단계에서 사용할 최종 분포를 얻기 위해 Circuit을 마지막으로 한 번 더 실행합니다.

초기 매개변수로 Circuit 정의하기

임의로 선택한 매개변수부터 시작합니다.

initial_gamma = np.pi
initial_beta = np.pi / 2
init_params = [initial_gamma, initial_beta]

Backend 및 실행 프리미티브 정의하기

IBM® Backend와 상호작용하려면 Qiskit Runtime 프리미티브를 사용하세요. 두 가지 프리미티브는 Sampler와 Estimator이며, 양자 컴퓨터에서 실행할 측정 유형에 따라 프리미티브를 선택합니다. HCH_C의 최솟값 구하기에는 비용 함수 측정이 단순히 HC\langle H_C \rangle의 기댓값이므로 Estimator를 사용합니다.

실행하기

프리미티브는 양자 디바이스에서 워크로드를 스케줄링하기 위한 다양한 실행 모드를 제공하며, QAOA 워크플로는 Session 안에서 반복적으로 실행됩니다.

&quot;execution mode&quot;

Sampler 기반 비용 함수를 SciPy 최소화 루틴에 연결하여 최적의 매개변수를 찾을 수 있습니다.

def cost_func_estimator(params, ansatz, hamiltonian, estimator):
# transform the observable defined on virtual qubits to
# an observable defined on all physical qubits
isa_hamiltonian = hamiltonian.apply_layout(ansatz.layout)

pub = (ansatz, isa_hamiltonian, params)
job = estimator.run([pub])

results = job.result()[0]
cost = results.data.evs

objective_func_vals.append(cost)

return cost
from qiskit_ibm_runtime import Session, EstimatorV2
from scipy.optimize import minimize

objective_func_vals = [] # Global variable
with Session(backend=backend) as session:
# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `session=`
estimator = EstimatorV2(mode=session)
estimator.options.default_shots = 1000

# Set simple error suppression/mitigation options
estimator.options.dynamical_decoupling.enable = True
estimator.options.dynamical_decoupling.sequence_type = "XY4"
estimator.options.twirling.enable_gates = True
estimator.options.twirling.num_randomizations = "auto"

result = minimize(
cost_func_estimator,
init_params,
args=(candidate_circuit, cost_hamiltonian, estimator),
method="COBYLA",
tol=1e-2,
)
print(result)
message: Optimization terminated successfully.
success: True
status: 1
fun: -0.6557925874481715
x: [ 2.873e+00 9.414e-01]
nfev: 21
maxcv: 0.0

최적화기가 비용을 줄이고 Circuit에 더 나은 매개변수를 찾을 수 있었습니다.

plt.figure(figsize=(12, 6))
plt.plot(objective_func_vals)
plt.xlabel("Iteration")
plt.ylabel("Cost")
plt.show()

Output of the previous code cell

최적의 매개변수를 찾은 후에는 이 매개변수를 Circuit에 할당하고, 최적화된 매개변수로 얻은 최종 분포를 샘플링할 수 있습니다. 여기서는 그래프의 최적 컷에 해당하는 비트스트링 측정값의 확률 분포가 필요하므로 Sampler 프리미티브를 사용해야 합니다.

참고: 이는 양자 컴퓨터에서 양자 상태 ψ\psi를 준비한 후 측정하는 것을 의미합니다. 측정을 수행하면 상태가 단일 계산 기저 상태(예: 010101110000...)로 붕괴되며, 이는 초기 최적화 문제(maxf(x)\max f(x) 또는 minf(x)\min f(x))에 대한 후보 해 xx에 해당합니다.

optimized_circuit = candidate_circuit.assign_parameters(result.x)
optimized_circuit.draw("mpl", fold=False, idle_wires=False)

Output of the previous code cell

from qiskit_ibm_runtime import SamplerV2

# If using qiskit-ibm-runtime<0.24.0, change `mode=` to `backend=`
sampler = SamplerV2(mode=backend)

# Set simple error suppression/mitigation options
sampler.options.dynamical_decoupling.enable = True
sampler.options.dynamical_decoupling.sequence_type = "XY4"
sampler.options.twirling.enable_gates = True
sampler.options.twirling.num_randomizations = "auto"

pub = (optimized_circuit,)
job = sampler.run([pub], shots=int(1e4))
counts_int = job.result()[0].data.meas.get_int_counts()
counts_bin = job.result()[0].data.meas.get_counts()
shots = sum(counts_int.values())
final_distribution_int = {key: val / shots for key, val in counts_int.items()}
final_distribution_bin = {key: val / shots for key, val in counts_bin.items()}
print(final_distribution_int)
{12: 0.0652, 31: 0.0089, 4: 0.0085, 13: 0.0731, 26: 0.0256, 28: 0.0246, 17: 0.0405, 25: 0.0591, 20: 0.031, 15: 0.0221, 8: 0.017, 21: 0.0371, 14: 0.0461, 16: 0.0229, 19: 0.0723, 23: 0.0199, 22: 0.0478, 18: 0.0708, 24: 0.0165, 6: 0.0525, 7: 0.0155, 5: 0.0245, 3: 0.0231, 29: 0.0121, 30: 0.0062, 10: 0.0363, 1: 0.0097, 9: 0.042, 27: 0.0094, 11: 0.0349, 0: 0.0129, 2: 0.0119}

3.5 4단계. 후처리 및 결과를 고전 형식으로 반환

후처리 단계에서는 샘플링 출력을 해석하여 원래 문제에 대한 해를 반환합니다. 이 경우, 최적 컷을 결정하는 가장 높은 확률을 가진 비트스트링에 관심이 있습니다. 문제의 대칭성으로 인해 네 가지 가능한 해가 존재하며, 샘플링 과정은 그 중 하나를 약간 더 높은 확률로 반환하게 됩니다. 아래에 그려진 분포에서 네 개의 비트스트링이 나머지보다 뚜렷하게 높은 확률을 가지는 것을 확인할 수 있습니다.

# auxiliary functions to sample most likely bitstring
def to_bitstring(integer, num_bits):
result = np.binary_repr(integer, width=num_bits)
return [int(digit) for digit in result]

keys = list(final_distribution_int.keys())
values = list(final_distribution_int.values())
most_likely = keys[np.argmax(np.abs(values))]
most_likely_bitstring = to_bitstring(most_likely, len(graph))
most_likely_bitstring.reverse()

print("Result bitstring:", most_likely_bitstring)
Result bitstring: [1, 0, 1, 1, 0]
import matplotlib.pyplot as plt

matplotlib.rcParams.update({"font.size": 10})
final_bits = final_distribution_bin
values = np.abs(list(final_bits.values()))
top_4_values = sorted(values, reverse=True)[:4]
positions = []
for value in top_4_values:
positions.append(np.where(values == value)[0])
fig = plt.figure(figsize=(11, 6))
ax = fig.add_subplot(1, 1, 1)
plt.xticks(rotation=45)
plt.title("Result Distribution")
plt.xlabel("Bitstrings (reversed)")
plt.ylabel("Probability")
ax.bar(list(final_bits.keys()), list(final_bits.values()), color="tab:grey")
for p in positions:
ax.get_children()[p[0].item()].set_color("tab:purple")
plt.show()

Output of the previous code cell

최적 컷 시각화

최적 비트스트링으로부터 원래 그래프에 해당 컷을 시각화할 수 있습니다.

colors = ["tab:grey" if i == 0 else "tab:purple" for i in most_likely_bitstring]
mpl_draw(graph, node_size=600, pos=pos, with_labels=True, labels=str, node_color=colors)

Output of the previous code cell

그리고 컷의 값을 계산합니다. 노이즈로 인해 해가 최적이 아닐 수 있습니다 (최적 해의 컷 값은 5입니다).

from typing import Sequence

def evaluate_sample(x: Sequence[int], graph: rx.PyGraph) -> float:
assert len(x) == len(
list(graph.nodes())
), "The length of x must coincide with the number of nodes in the graph."
return sum(
x[u] * (1 - x[v]) + x[v] * (1 - x[u]) for u, v in list(graph.edge_list())
)

cut_value = evaluate_sample(most_likely_bitstring, graph)
print("The value of the cut is:", cut_value)
The value of the cut is: 5

이것으로 소규모 QAOA 튜토리얼을 마칩니다. 양자 근사 최적화 알고리즘 튜토리얼의 "2부: 규모 확장하기!"에서 유틸리티 규모의 QAOA를 적용하는 방법을 배울 수 있습니다.

# Check Qiskit version
import qiskit

qiskit.__version__
'2.0.2'