주 콘텐츠로 건너뛰기

Q-CTRL의 Optimization Solver를 이용한 고차 이진 최적화

참고

Qiskit Functions는 IBM Quantum® Premium Plan, Flex Plan, On-Prem (IBM Quantum Platform API 경유) Plan 사용자만 이용 가능한 실험적 기능입니다. 현재 미리 보기 출시 상태이며 변경될 수 있습니다.

사용 예상 시간: Heron r2 프로세서 기준 24분. (참고: 이는 예상치이며 실제 실행 시간은 다를 수 있습니다.)

배경

이 튜토리얼에서는 Q-CTRL Fire Opal의 Qiskit Function인 Optimization Solver를 사용하여 고차 이진 최적화(HOBO) 문제를 해결하는 방법을 설명합니다. 이 튜토리얼에서 다루는 예제는 3차 항을 가진 무작위 결합 156-Qubit Ising 모델의 바닥 상태 에너지를 찾는 최적화 문제입니다. Optimization Solver는 목적 함수로 정의할 수 있는 일반적인 최적화 문제에 사용할 수 있습니다.

Optimization Solver는 양자 하드웨어에서 최적화 문제를 해결하는 하드웨어 인식 구현 단계를 완전히 자동화하며, 양자 실행에 Performance Management를 활용함으로써 유틸리티 규모에서 정확한 해를 도출합니다. 전체 Optimization Solver 워크플로우 및 벤치마킹 결과에 대한 자세한 요약은 공개된 논문을 참조하세요.

이 튜토리얼은 다음 단계로 구성됩니다:

  1. 목적 함수로 문제 정의하기
  2. Fire Opal Optimization Solver를 사용하여 하이브리드 알고리즘 실행하기
  3. 결과 평가하기

요구 사항

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

  • Qiskit Functions (pip install qiskit-ibm-catalog)
  • SymPy (pip install sympy)

또한 Optimization Solver 함수에 대한 접근 권한을 받아야 합니다. 양식을 작성하여 접근을 요청하세요.

설정

먼저, 필요한 패키지와 도구를 가져옵니다.

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit-ibm-catalog sympy
# Qiskit Functions Catalog
from qiskit_ibm_catalog import QiskitFunctionsCatalog

# SymPy tools for constructing objective function
from sympy import Poly
from sympy import symbols, srepr

# Tools for plotting and evaluating results
import numpy as np
import matplotlib.pyplot as plt
from sympy import lambdify

튜토리얼 전반에 걸쳐 Qiskit Runtime과 Qiskit Functions 인증에 사용할 IBM Quantum Platform 자격 증명을 정의합니다.

# Credentials
token = "<YOUR-API_KEY>" # Use the 44-characters API_KEY you have created and saved from the IBM Quantum Platform Home dashboard
instance = "<YOUR_CRN>"

1단계: 목적 함수로 문제 정의하기

Optimization Solver는 목적 함수 또는 그래프를 입력으로 받습니다. 이 튜토리얼에서 Ising 스핀 글라스 최소화 문제는 목적 함수로 정의되며, IBM® 디바이스의 heavy-hex 토폴로지에 맞게 조정되었습니다.

이 목적 함수는 3차, 2차, 1차 항을 포함하기 때문에 HOBO 문제 클래스에 속하며, 이는 기존의 이차 비제약 이진 최적화(QUBO) 문제보다 훨씬 복잡한 것으로 알려져 있습니다.

문제 정의 구성과 Optimization Solver에서 얻은 이전 결과에 대한 자세한 논의는 이 기술 논문을 참조하세요. 이 문제는 원래 Los Alamos 국립연구소가 발표한 논문의 일부로 정의 및 평가되었으며, 156-Qubit IBM Quantum Heron 프로세서의 전체 디바이스 폭을 활용하도록 개선되었습니다.

qubit_count = 156

# Create symbolic variables to represent qubits
x = symbols([f"x[{i}]" for i in range(qubit_count)])

# # Define a polynomial representing a spin glass model
spin_glass_poly = Poly(
-4 * x[0] * x[1]
- 8 * x[1] * x[2] * x[3]
+ 8 * x[1] * x[2]
+ 4 * x[1] * x[3]
- 4 * x[2]
+ 8 * x[3] * x[4] * x[5]
- 4 * x[3] * x[5]
- 8 * x[3] * x[16] * x[23]
+ 4 * x[3] * x[23]
- 2 * x[3]
- 4 * x[4]
- 8 * x[5] * x[6] * x[7]
+ 8 * x[5] * x[6]
+ 4 * x[5] * x[7]
- 2 * x[5]
+ 8 * x[6] * x[7]
- 4 * x[6]
- 8 * x[7] * x[8] * x[9]
+ 4 * x[7] * x[9]
- 8 * x[7] * x[17] * x[27]
+ 4 * x[7] * x[27]
- 6 * x[7]
+ 8 * x[8] * x[9]
+ 8 * x[9] * x[10] * x[11]
- 4 * x[9] * x[11]
- 2 * x[9]
- 8 * x[10] * x[11]
+ 4 * x[10]
- 8 * x[11] * x[12] * x[13]
+ 4 * x[11] * x[13]
- 8 * x[11] * x[18] * x[31]
+ 8 * x[11] * x[18]
+ 4 * x[11] * x[31]
- 2 * x[11]
+ 8 * x[12] * x[13]
+ 8 * x[13] * x[14] * x[15]
- 4 * x[13] * x[15]
- 2 * x[13]
- 8 * x[14] * x[15]
+ 4 * x[14]
- 8 * x[15] * x[19] * x[35]
+ 8 * x[15] * x[19]
+ 4 * x[15] * x[35]
- 2 * x[15]
+ 8 * x[16] * x[23]
+ 8 * x[17] * x[27]
- 4 * x[17]
+ 8 * x[18] * x[31]
- 8 * x[18]
+ 8 * x[19] * x[35]
- 8 * x[19]
+ 4 * x[20] * x[21]
- 4 * x[20]
- 8 * x[21] * x[22] * x[23]
+ 8 * x[21] * x[22]
+ 4 * x[21] * x[23]
- 8 * x[21] * x[36] * x[41]
+ 4 * x[21] * x[41]
- 4 * x[21]
+ 8 * x[22] * x[23]
- 8 * x[22]
+ 8 * x[23] * x[24] * x[25]
- 4 * x[23] * x[25]
- 10 * x[23]
- 8 * x[24] * x[25]
+ 8 * x[25] * x[26] * x[27]
- 8 * x[25] * x[26]
- 4 * x[25] * x[27]
+ 8 * x[25] * x[37] * x[45]
- 8 * x[25] * x[37]
- 4 * x[25] * x[45]
+ 14 * x[25]
- 8 * x[26] * x[27]
+ 4 * x[26]
+ 8 * x[27] * x[28] * x[29]
- 4 * x[27] * x[29]
- 2 * x[27]
- 8 * x[28] * x[29]
- 8 * x[29] * x[30] * x[31]
+ 4 * x[29] * x[31]
+ 8 * x[29] * x[38] * x[49]
- 8 * x[29] * x[38]
- 4 * x[29] * x[49]
+ 6 * x[29]
+ 8 * x[30] * x[31]
- 4 * x[30]
- 8 * x[31] * x[32] * x[33]
+ 4 * x[31] * x[33]
- 6 * x[31]
+ 8 * x[33] * x[34] * x[35]
- 4 * x[33] * x[35]
- 8 * x[33] * x[39] * x[53]
+ 8 * x[33] * x[39]
+ 4 * x[33] * x[53]
- 6 * x[33]
- 8 * x[34] * x[35]
+ 2 * x[35]
+ 8 * x[36] * x[41]
- 8 * x[37] * x[45]
+ 4 * x[37]
- 8 * x[38] * x[49]
+ 4 * x[38]
+ 4 * x[40] * x[41]
- 8 * x[41] * x[42] * x[43]
+ 4 * x[41] * x[43]
- 8 * x[41]
+ 8 * x[42] * x[43]
- 4 * x[42]
- 8 * x[43] * x[44] * x[45]
+ 8 * x[43] * x[44]
+ 4 * x[43] * x[45]
- 8 * x[43] * x[56] * x[63]
+ 4 * x[43] * x[63]
- 6 * x[43]
- 4 * x[44]
- 8 * x[45] * x[46] * x[47]
+ 4 * x[45] * x[47]
+ 2 * x[45]
+ 4 * x[46]
- 8 * x[47] * x[48] * x[49]
+ 8 * x[47] * x[48]
+ 4 * x[47] * x[49]
- 8 * x[47] * x[57] * x[67]
+ 4 * x[47] * x[67]
- 2 * x[47]
- 4 * x[48]
- 8 * x[49] * x[50] * x[51]
+ 8 * x[49] * x[50]
+ 4 * x[49] * x[51]
- 2 * x[49]
+ 8 * x[50] * x[51]
- 8 * x[50]
- 8 * x[51] * x[52] * x[53]
+ 8 * x[51] * x[52]
+ 4 * x[51] * x[53]
- 8 * x[51] * x[58] * x[71]
+ 4 * x[51] * x[71]
- 6 * x[51]
+ 8 * x[52] * x[53]
- 8 * x[52]
+ 8 * x[53] * x[54] * x[55]
- 8 * x[53] * x[54]
- 4 * x[53] * x[55]
- 2 * x[53]
+ 4 * x[54]
- 8 * x[55] * x[59] * x[75]
+ 4 * x[55] * x[75]
- 2 * x[55]
+ 8 * x[56] * x[63]
+ 8 * x[57] * x[67]
- 4 * x[57]
+ 8 * x[58] * x[71]
+ 8 * x[59] * x[75]
- 4 * x[59]
+ 4 * x[60] * x[61]
+ 8 * x[61] * x[62] * x[63]
- 4 * x[61] * x[63]
+ 8 * x[61] * x[76] * x[81]
- 8 * x[61] * x[76]
- 4 * x[61] * x[81]
- 8 * x[63] * x[64] * x[65]
+ 8 * x[63] * x[64]
+ 4 * x[63] * x[65]
- 6 * x[63]
+ 8 * x[65] * x[66] * x[67]
- 8 * x[65] * x[66]
- 4 * x[65] * x[67]
- 8 * x[65] * x[77] * x[85]
+ 4 * x[65] * x[85]
+ 2 * x[65]
+ 4 * x[66]
- 8 * x[67] * x[68] * x[69]
+ 8 * x[67] * x[68]
+ 4 * x[67] * x[69]
- 10 * x[67]
+ 8 * x[68] * x[69]
- 4 * x[68]
+ 8 * x[69] * x[70] * x[71]
- 4 * x[69] * x[71]
- 8 * x[69] * x[78] * x[89]
+ 4 * x[69] * x[89]
- 6 * x[69]
+ 8 * x[71] * x[72] * x[73]
- 8 * x[71] * x[72]
- 4 * x[71] * x[73]
+ 2 * x[71]
- 8 * x[72] * x[73]
+ 8 * x[72]
- 8 * x[73] * x[74] * x[75]
+ 8 * x[73] * x[74]
+ 4 * x[73] * x[75]
- 8 * x[73] * x[79] * x[93]
+ 8 * x[73] * x[79]
+ 4 * x[73] * x[93]
- 6 * x[73]
+ 8 * x[74] * x[75]
- 4 * x[74]
- 10 * x[75]
+ 4 * x[76]
+ 8 * x[78] * x[89]
- 4 * x[78]
- 4 * x[79]
- 4 * x[80] * x[81]
+ 4 * x[80]
- 8 * x[81] * x[82] * x[83]
+ 8 * x[81] * x[82]
+ 4 * x[81] * x[83]
+ 8 * x[82] * x[83]
- 8 * x[82]
- 8 * x[83] * x[84] * x[85]
+ 4 * x[83] * x[85]
- 8 * x[83] * x[96] * x[103]
+ 4 * x[83] * x[103]
- 2 * x[83]
- 8 * x[85] * x[86] * x[87]
+ 8 * x[85] * x[86]
+ 4 * x[85] * x[87]
- 6 * x[85]
+ 8 * x[86] * x[87]
- 4 * x[86]
- 8 * x[87] * x[88] * x[89]
+ 4 * x[87] * x[89]
+ 8 * x[87] * x[97] * x[107]
- 8 * x[87] * x[97]
- 4 * x[87] * x[107]
+ 2 * x[87]
+ 4 * x[88]
- 8 * x[89] * x[90] * x[91]
+ 8 * x[89] * x[90]
+ 4 * x[89] * x[91]
- 10 * x[89]
+ 8 * x[90] * x[91]
- 8 * x[90]
- 8 * x[91] * x[92] * x[93]
+ 4 * x[91] * x[93]
- 8 * x[91] * x[98] * x[111]
+ 8 * x[91] * x[98]
+ 4 * x[91] * x[111]
- 10 * x[91]
+ 8 * x[92] * x[93]
- 4 * x[92]
- 8 * x[93] * x[94] * x[95]
+ 4 * x[93] * x[95]
- 6 * x[93]
+ 8 * x[95] * x[99] * x[115]
- 8 * x[95] * x[99]
- 4 * x[95] * x[115]
+ 2 * x[95]
+ 4 * x[96]
- 8 * x[97] * x[107]
+ 4 * x[97]
- 4 * x[98]
- 8 * x[99] * x[115]
+ 4 * x[99]
- 4 * x[100] * x[101]
+ 8 * x[101] * x[102] * x[103]
- 8 * x[101] * x[102]
- 4 * x[101] * x[103]
- 8 * x[101] * x[116] * x[121]
+ 8 * x[101] * x[116]
+ 4 * x[101] * x[121]
+ 4 * x[101]
- 8 * x[103] * x[104] * x[105]
+ 4 * x[103] * x[105]
+ 2 * x[103]
+ 8 * x[105] * x[106] * x[107]
- 4 * x[105] * x[107]
- 8 * x[105] * x[117] * x[125]
+ 4 * x[105] * x[125]
+ 2 * x[105]
- 8 * x[106] * x[107]
+ 4 * x[106]
+ 8 * x[107] * x[108] * x[109]
- 4 * x[107] * x[109]
+ 6 * x[107]
- 4 * x[108]
+ 8 * x[109] * x[110] * x[111]
- 4 * x[109] * x[111]
- 8 * x[109] * x[118] * x[129]
+ 4 * x[109] * x[129]
+ 2 * x[109]
- 8 * x[110] * x[111]
+ 4 * x[110]
- 8 * x[111] * x[112] * x[113]
+ 8 * x[111] * x[112]
+ 4 * x[111] * x[113]
+ 2 * x[111]
+ 8 * x[112] * x[113]
- 8 * x[112]
- 8 * x[113] * x[114] * x[115]
+ 4 * x[113] * x[115]
- 8 * x[113] * x[119] * x[133]
+ 4 * x[113] * x[133]
- 2 * x[113]
+ 6 * x[115]
- 4 * x[116]
+ 4 * x[118]
+ 4 * x[119]
+ 4 * x[120] * x[121]
- 8 * x[121] * x[122] * x[123]
+ 4 * x[121] * x[123]
- 4 * x[121]
+ 4 * x[122]
- 8 * x[123] * x[124] * x[125]
+ 4 * x[123] * x[125]
- 8 * x[123] * x[136] * x[143]
+ 4 * x[123] * x[143]
- 2 * x[123]
+ 8 * x[124] * x[125]
- 4 * x[124]
+ 8 * x[125] * x[126] * x[127]
- 8 * x[125] * x[126]
- 4 * x[125] * x[127]
+ 2 * x[125]
- 8 * x[127] * x[128] * x[129]
+ 8 * x[127] * x[128]
+ 4 * x[127] * x[129]
+ 8 * x[127] * x[137] * x[147]
- 8 * x[127] * x[137]
- 4 * x[127] * x[147]
- 2 * x[127]
+ 8 * x[129] * x[130] * x[131]
- 4 * x[129] * x[131]
+ 2 * x[129]
- 4 * x[130]
- 8 * x[131] * x[132] * x[133]
+ 4 * x[131] * x[133]
- 8 * x[131] * x[138] * x[151]
+ 4 * x[131] * x[151]
- 2 * x[131]
+ 8 * x[133] * x[134] * x[135]
- 4 * x[133] * x[135]
+ 2 * x[133]
- 8 * x[134] * x[135]
+ 4 * x[134]
- 8 * x[135] * x[139] * x[155]
+ 8 * x[135] * x[139]
+ 4 * x[135] * x[155]
+ 2 * x[135]
+ 8 * x[136] * x[143]
- 4 * x[136]
+ 4 * x[138]
+ 8 * x[139] * x[155]
- 4 * x[139]
- 4 * x[140] * x[141]
- 8 * x[141] * x[142] * x[143]
+ 8 * x[141] * x[142]
+ 4 * x[141] * x[143]
+ 8 * x[142] * x[143]
- 8 * x[142]
- 8 * x[143] * x[144] * x[145]
+ 8 * x[143] * x[144]
+ 4 * x[143] * x[145]
- 14 * x[143]
+ 8 * x[144] * x[145]
- 8 * x[144]
- 8 * x[145] * x[146] * x[147]
+ 8 * x[145] * x[146]
+ 4 * x[145] * x[147]
- 6 * x[145]
+ 8 * x[146] * x[147]
- 4 * x[146]
- 8 * x[147] * x[148] * x[149]
+ 8 * x[147] * x[148]
+ 4 * x[147] * x[149]
- 6 * x[147]
- 4 * x[148]
- 8 * x[149] * x[150] * x[151]
+ 8 * x[149] * x[150]
+ 4 * x[149] * x[151]
- 6 * x[149]
+ 8 * x[151] * x[152] * x[153]
- 4 * x[151] * x[153]
+ 2 * x[151]
+ 8 * x[153] * x[154] * x[155]
- 8 * x[153] * x[154]
- 4 * x[153] * x[155]
+ 2 * x[153]
- 8 * x[154] * x[155]
+ 4 * x[154]
- 2 * x[155]
+ 46,
x,
domain="ZZ",
)

2단계: Fire Opal Optimization Solver를 사용한 하이브리드 알고리즘 실행

이제 Optimization Solver Qiskit Function을 사용하여 알고리즘을 실행합니다. 내부적으로 Optimization Solver는 문제를 하이브리드 양자 알고리즘에 매핑하고, 오류 억제를 적용하여 양자 Circuit을 실행하며, 고전적 최적화를 수행합니다.

# Authenticate to the Qiskit Functions Catalog
catalog = QiskitFunctionsCatalog(
token=token,
instance=instance,
)

# Load the function
solver = catalog.load("q-ctrl/optimization_solver")

선택한 디바이스에 최소 156개의 Qubit이 있는지 확인합니다.

# Specify the target backend name
backend_name = "<CHOOSE_A_BACKEND>"

Solver는 목적 함수의 문자열 표현을 입력으로 받습니다.

# Convert the objective function to string format
spin_glass_poly_as_str = srepr(spin_glass_poly)
# Run the problem
spin_glass_job = solver.run(
problem=spin_glass_poly_as_str,
run_options={"backend_name": backend_name},
)

익숙한 Qiskit Serverless API를 사용하여 Qiskit Function 워크로드의 상태를 확인할 수 있습니다:

# Get job status
spin_glass_job.status()

Solver는 솔루션과 관련 메타데이터(솔루션 비트스트링, 반복 횟수, 변수와 비트스트링 간의 매핑 등)가 담긴 딕셔너리를 반환합니다. Solver의 입출력에 대한 전체 정의는 문서를 참고하세요.

# Poll for results
result = spin_glass_job.result()
# Get the final bitstring distribution and set the number of shots
distribution = result["final_bitstring_distribution"]

3단계: 결과 평가

# Get the solution ground state energy
print(f"Minimum ground state energy: {result["solution_bitstring_cost"]}")
Minimum ground state energy: -242.0

Solver는 고전적 최적화 소프트웨어를 통해 검증된 정확한 솔루션을 찾았습니다. 이 유틸리티 규모 문제의 복잡성을 고전적으로 풀기 위해서는 IBM ILOG CPLEX Optimization Studio (CPLEX)Gurobi Optimization과 같은 고급 최적화 소프트웨어가 필요합니다. 결과의 품질을 시각적으로 분석하기 위해, 비트스트링과 확률로부터 비용 값을 계산하여 결과를 플롯할 수 있습니다. 비교를 위해, 무작위로 샘플링된 비트스트링의 분포(고전적 "무차별 대입" 솔루션에 해당)와 함께 결과를 나란히 표시합니다. 알고리즘이 일관되게 더 낮은 비용을 찾는다면, 양자 알고리즘이 최적화 문제를 효과적으로 풀고 있음을 나타냅니다.

def plot_cost_histogram(
costs, probabilities, distribution, qubit_count, bitstring_cost
):
"""Plots a histogram comparing the cost distributions of Q-CTRL Solver and random sampling."""

# Set figure DPI for higher resolution and font size for labels
plt.rcParams["figure.dpi"] = 300
plt.rcParams.update({"font.size": 6}) # Set default font size to 6

# Define labels and colors for the plot
labels = ["Q-CTRL Solver", "Random Sampling"]
colors = ["#680CE9", "#E04542"]

# Calculate total shots (total number of bitstrings in the distribution)
shots = sum(distribution.values())

# Generate random bitstrings for comparison (random sampling)
rng = np.random.default_rng(seed=0)
random_array = rng.integers(
0, 2, size=(shots, qubit_count)
) # Generate random bitstrings (0 or 1 for each qubit)
random_bitstrings = ["".join(row.astype(str)) for row in random_array]

# Compute the cost for each random bitstring
random_costs = [bitstring_cost(k) for k in random_bitstrings]

# Set uniform probabilities for the random sampling
random_probabilities = (
np.ones(shape=(shots,)) / shots
) # Equal probability for each random bitstring

# Find the minimum and maximum costs for binning the histogram
min_cost = np.min(costs)
max_cost = np.max(random_costs)

# Create a histogram plot with a smaller figure size (4x2 inches)
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4, 2))

# Plot histograms for the Q-CTRL solver and random sampling costs
_, _, _ = ax.hist(
[costs, random_costs], # Data for the two histograms
np.arange(min_cost, max_cost, 2), # Bins for the histogram
weights=[
probabilities,
random_probabilities,
], # Probabilities for each data set
label=labels, # Labels for the legend
color=colors, # Colors for each histogram
histtype="stepfilled", # Filled step histogram
align="mid", # Align bars to the bin center
alpha=0.8, # Transparency
)

# Set the x and y labels for the plot
ax.set_xlabel("Cost")
ax.set_ylabel("Probability")

# Add the legend to the plot
ax.legend()

# Show the plot
plt.show()
# Convert spin_glass_poly into a NumPy-compatible function
poly_as_numpy_function = lambdify(x, spin_glass_poly.as_expr(), "numpy")

# Function to compute the cost of a given bitstring using spin_glass_poly
def bitstring_cost(bitstring: str) -> float:
# Convert bitstring to a reversed list of integers (0s and 1s)
return float(
poly_as_numpy_function(*[int(b) for b in str(bitstring[::-1])])
)

# Calculate the cost of each bitstring in the distribution
costs = [bitstring_cost(k) for k, _ in distribution.items()]

# Extract probabilities from the bitstring distribution
probabilities = np.array([v for _, v in distribution.items()])
probabilities = probabilities / sum(
probabilities
) # Normalize to get probabilities

plot_cost_histogram(
costs, probabilities, distribution, qubit_count, bitstring_cost
)

이전 코드 셀의 출력

이 최적화 알고리즘의 목표는 Ising 모델의 최소 기저 상태를 찾는 것이므로, 값이 낮을수록 더 좋은 솔루션을 의미합니다. 따라서 Fire Opal Optimization Solver가 생성한 솔루션이 무작위 선택보다 훨씬 우수하다는 것을 시각적으로 명확히 확인할 수 있습니다.