Квантові алгоритми: варіаційні квантові алгоритми

примітка

Takashi Imamichi (24 травня 2024 р.)

Завантажити PDF оригінальної лекції. Зверни увагу, що деякі фрагменти коду можуть бути застарілими, оскільки це статичні зображення.

Приблизний час роботи QPU для цього експерименту становить 9 хвилин (тестувалося на процесорі Eagle).

(цей ноутбук може не завершити виконання в межах часу, дозволеного на тарифному плані Open Plan. Будь ласка, використовуй ресурси квантових обчислень раціонально.)

1. Вступ

Цей посібник надає огляд гібридного квантово-класичного алгоритму, зокрема зосереджуючись на варіаційному квантовому власному розв'язувачі (VQE) та квантовому наближеному алгоритмі оптимізації (QAOA). Основна мета цих алгоритмів — розв'язувати задачі оптимізації за допомогою квантових Circuit із параметризованими квантовими Gate.

Незважаючи на прогрес у квантових обчисленнях, наявність шумів у сучасних квантових пристроях ускладнює отримання змістовних результатів із глибоких квантових Circuit. Щоб подолати цю проблему, VQE та QAOA використовують гібридний квантово-класичний підхід, який передбачає ітеративне виконання відносно коротких квантових Circuit за допомогою квантових обчислень та оптимізацію параметрів цільових параметризованих квантових Circuit за допомогою класичних обчислень.

QAOA має потенціал забезпечувати оптимальні рішення для цільових задач на рівні утилітарного масштабу завдяки застосуванню різних технік пом'якшення та придушення помилок. VQE має багато застосувань (наприклад, у квантовій хімії), де він є менш масштабованим. Проте з'явилася низка підходів, пов'язаних із власними значеннями, які доповнюють і розширюють VQE, зокрема діагоналізація підпростору Крилова та квантова діагоналізація на основі вибірки (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 до дуже простого випадку, щоб побачити, як він працює. Ми обчислимо мінімальне власне значення матриці Паулі $Z$ за допомогою VQE. Почнемо з імпорту деяких загальних пакетів.

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 із Gate, що приймають варіаційні параметри:

# 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")

Якщо ми хочемо оцінити математичне сподівання оператора (наприклад, $Z$), слід використовувати Estimator. Якщо ж ми хочемо розглянути стани системи, використовуємо Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Ми можемо обчислити кількість бітових рядків 0 та 1 із випадковими значеннями параметрів [1, 2, 3] за допомогою Sampler.

# 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}

Ми знаємо, що математичне сподівання Z можна обчислити як $\langle Z \rangle = p_0 - p_1$ з імовірностями $\{0: p_0, 1: p_1\}$.

# 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)

Нам потрібно перебрати параметри й знайти ті, що дають найменше власне значення. Ми створюємо функцію, яка приймає значення параметрів варіаційної форми та повертає математичне сподівання $\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

Застосуймо функцію minimize із SciPy, щоб знайти мінімальне власне значення 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 Вправа

Обчисли мінімальне власне значення $Z \otimes Z$ за допомогою VQE.

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

Щоб отримати власні значення за допомогою варіаційного квантового алгоритму, будуємо Circuit із Gate, що приймають варіаційні параметри:

# 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")

Якщо ми хочемо оцінити математичне сподівання оператора (наприклад, $Z \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 — це інтуїтивний, повторюваний набір кроків для реалізації робочого процесу квантових обчислень: Ми застосуємо ці патерни до задач комбінаторної оптимізації та покажемо, як розв'язати задачу Max-Cut за допомогою Квантового наближеного алгоритму оптимізації (QAOA) — гібридного (квантово-класичного) ітераційного методу.

Зверни увагу: ця частина про QAOA базується на «Частині 1: малоформатний QAOA» туторіалу Квантовий наближений алгоритм оптимізації. Переглянь туторіал, щоб дізнатися, як масштабувати його.

3.1 (Малоформатний) патерн Qiskit для оптимізації

У цій частині ми використаємо малоформатну задачу Max-Cut, щоб проілюструвати кроки, необхідні для розв'язання оптимізаційної задачі на квантовому комп'ютері. Задача Max-Cut — це важка для вирішення оптимізаційна задача (точніше, NP-важка), що має низку практичних застосувань у кластеризації, мережевій науці та статистичній фізиці. У цьому туторіалі розглядається граф із вузлами, з'єднаними ребрами, і ставиться мета розбити вузли на дві множини, «розрізаючи» ребра таким чином, щоб кількість розрізаних ребер була максимальною. Щоб зрозуміти, як відобразити цю задачу на квантовий алгоритм, спершу розглянемо, як задача Max-Cut стає класичною комбінаторною оптимізаційною задачею, мінімізуючи функцію $f(x)$

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

де вхідне $x$ — вектор, компоненти якого відповідають кожному вузлу графу. Кожна з цих компонент обмежена значеннями $0$ або $1$ (що позначає, чи включено вузол у розріз). У цьому малоформатному прикладі граф містить $n=5$ вузлів.

Можна записати функцію для пари вузлів $i,j$, яка показує, чи входить відповідне ребро $(i,j)$ у розріз. Наприклад, функція $x_i + x_j - 2 x_i x_j$ дорівнює 1 лише тоді, коли рівно одне з $x_i$ або $x_j$ дорівнює 1 (тобто ребро входить у розріз), і 0 в іншому разі. Задачу максимізації кількості розрізаних ребер можна сформулювати як

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

що можна переписати у вигляді мінімізації

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

Мінімум $f(x)$ у цьому випадку досягається, коли кількість ребер, що перетинаються розрізом, є максимальною. Як ти бачиш, поки що нічого квантового тут немає. Нам потрібно переформулювати задачу у форму, зрозумілу квантовому комп'ютеру. Ініціалізуємо задачу, створивши граф із $n=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)

3.2 Крок 1. Відображення класичних вхідних даних на квантову задачу

Перший крок патерну — відобразити класичну задачу (граф) на квантові Circuit та оператори. Для цього потрібно виконати три основні дії:

1. Застосувати низку математичних перетворень, щоб записати задачу у форматі QUBO (квадратична безобмежена двійкова оптимізація).
2. Переформулювати оптимізаційну задачу у вигляді Гамільтоніана, основний стан якого відповідає розв'язку, що мінімізує функцію витрат.
3. Побудувати квантовий Circuit, який готуватиме основний стан цього Гамільтоніана через процес, схожий на квантовий відпал.

Примітка: У методології QAOA ти зрештою хочеш мати оператор (Гамільтоніан), який представляє функцію витрат гібридного алгоритму, а також параметризований Circuit (Ansatz), що представляє квантові стани-кандидати на розв'язок задачі. Ти можеш семплювати з цих станів-кандидатів, а потім оцінювати їх за допомогою функції витрат.

Граф → оптимізаційна задача

Перший крок відображення — це зміна нотації. Задача записується у форматі QUBO:

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

де $Q$ — матриця розміру $n\times n$ з дійсних чисел, $n$ — кількість вузлів графу, $x$ — вектор двійкових змінних, введених вище, а $x^T$ позначає транспонований вектор $x$.

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 як Гамільтоніан (тут — матриця, що представляє енергію системи):

$$H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.$$

Кроки переформулювання від задачі QAOA до Гамільтоніана

Щоб продемонструвати, як задача QAOA переписується в такій формі, спершу замінимо двійкові змінні $x_i$ на нові змінні $z_i\in\{-1, 1\}$ за допомогою

$$x_i = \frac{1-z_i}{2}.$$

Тут видно: якщо $x_i = 0$, то $z_i = 1$. Після підстановки $z_i$ замість $x_i$ в оптимізаційну задачу ($x^TQx$) отримуємо еквівалентну формулювання.

$$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}.$$

Якщо визначимо $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, відкинемо передмножник і константний член $n^2$, отримаємо два еквівалентних формулювання однієї й тієї ж оптимізаційної задачі.

$$min_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz$$

Тут $b$ залежить від $Q$. Зауваж, що для отримання $z^TQz + b^Tz$ ми відкинули множник 1/4 і константне зміщення $n^2$, оскільки вони не впливають на оптимізацію.

Тепер, щоб отримати квантове формулювання задачі, підвищимо $z_i$ до матриці Паулі $Z$ — матриці $2\times 2$ вигляду

$$Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.$$

Підставивши ці матриці в оптимізаційну задачу вище, отримаємо такий Гамільтоніан

$$H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.$$

Також нагадаємо, що матриці $Z$ вбудовані в обчислювальний простір квантового комп'ютера — простір Гільберта розміру $2^n\times 2^n$. Тому такі вирази, як $Z_iZ_j$, слід розуміти як тензорний добуток $Z_i\otimes Z_j$, вбудований у простір Гільберта розміру $2^n\times 2^n$. Наприклад, для задачі з п'ятьма змінними рішення термін $Z_1Z_3$ означає $I\otimes Z_3\otimes I\otimes Z_1\otimes I$, де $I$ — одинична матриця $2\times 2$.

Цей Гамільтоніан називається Гамільтоніаном функції витрат. Його особливість: основний стан відповідає розв'язку, що мінімізує функцію витрат $f(x)$. Отже, щоб розв'язати оптимізаційну задачу, тепер потрібно підготувати основний стан $H_C$ (або стан із великим перекриттям з ним) на квантовому комп'ютері. Тоді семплювання з цього стану з високою імовірністю дасть розв'язок задачі $\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

Гамільтоніан $H_C$ містить квантове визначення твоєї задачі. Тепер можна побудувати квантовий Circuit, який допоможе семплювати гарні розв'язки з квантового комп'ютера. QAOA натхненний квантовим відпалом і застосовує поперемінні шари операторів у квантовому Circuit.

Загальна ідея: почати в основному стані відомої системи, $H^{\otimes n}|0\rangle$, а потім скерувати систему до основного стану оператора витрат, що нас цікавить. Це досягається застосуванням операторів $\exp\{-i\gamma_k H_C\}$ та $\exp\{-i\beta_k H_m\}$ із кутами $\gamma_1,...,\gamma_p$ та $\beta_1,...,\beta_p~$.

Згенерований квантовий Circuit параметризований значеннями $\gamma_i$ та $\beta_i$, тому можна перебирати різні значення $\gamma_i$ і $\beta_i$ та семплювати з отриманого стану. У цьому прикладі розглянемо 1 шар QAOA, що містить два параметри: $\gamma_1$ та $\beta_1$.

from qiskit.circuit.library import QAOAAnsatz

circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

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

circuit.parameters

ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 Крок 2. Оптимізація Circuit для виконання на квантовому залізі

Circuit вище містить низку абстракцій, корисних для обговорення квантових алгоритмів, але неможливих для запуску на залізі. Щоб запустити схему на QPU, вона повинна пройти через ряд операцій, що складають крок транспіляції або оптимізації Circuit у патерні.

Бібліотека Qiskit пропонує серію проходів транспіляції, що охоплюють широкий спектр перетворень Circuit. Потрібно переконатися, що твій Circuit оптимізований під конкретну мету.

Транспіляція може включати кілька кроків, зокрема:

• Початкове відображення кубітів у Circuit (наприклад, змінних рішення) на фізичні Qubit пристрою.
• Розгортання інструкцій квантового Circuit до апаратно-нативних інструкцій, які розуміє Backend.
• Маршрутизація кубітів у Circuit, що взаємодіють, до фізично сусідніх Qubit.
• Придушення помилок шляхом додавання однокубітних Gate для придушення шумів методом динамічного розчеплення.

Більше інформації про транспіляцію — у нашій документації.

Наведений нижче код перетворює та оптимізує абстрактний Circuit у формат, готовий до виконання на одному з пристроїв, доступних через хмару, використовуючи сервіс Qiskit IBM® Runtime.

Зверни увагу: програми можна тестувати локально у «режимі локального тестування» перед відправкою на справжні квантові комп'ютери. Більше інформації про режим локального тестування — у документації.

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')>

3.4 Крок 3. Виконання за допомогою примітивів Qiskit

У робочому процесі QAOA оптимальні параметри QAOA знаходяться в ітераційному циклі оптимізації, який виконує серію обчислень Circuit та використовує класичний оптимізатор для пошуку оптимальних параметрів $\beta_k$ та $\gamma_k$. Цей цикл виконання складається з таких кроків:

1. Визначити початкові параметри
2. Ініціалізувати нову Session, що містить цикл оптимізації та примітив для семплювання Circuit
3. Після знаходження оптимального набору параметрів — виконати Circuit ще раз для отримання фінального розподілу, який використовуватиметься на кроці постобробки.

Визначення Circuit із початковими параметрами

Почнемо з довільно обраних параметрів.

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

Визначення Backend та примітиву виконання

Для взаємодії з Backend IBM® використовуй примітиви Qiskit Runtime. Два примітиви — це Sampler та Estimator, і вибір примітиву залежить від типу вимірювання, яке ти хочеш виконати на квантовому комп'ютері. Для мінімізації $H_C$ використовуй Estimator, оскільки вимірювання функції витрат — це просто очікуване значення $\langle H_C \rangle$.

Запуск

Примітиви пропонують різноманітні режими виконання для планування навантажень на квантових пристроях, а робочий процес QAOA виконується ітераційно в сесії. Функцію витрат на основі 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()

Після знаходження оптимальних параметрів Circuit можна призначити ці параметри та семплювати фінальний розподіл, отриманий з оптимізованими параметрами. Саме тут слід використовувати примітив Sampler, оскільки він дає розподіл імовірностей вимірювань бітових рядків, що відповідають оптимальному розрізу графу.

Примітка: Це означає підготовку квантового стану $\psi$ у комп'ютері та подальше його вимірювання. Вимірювання згорне стан до єдиного обчислювального базисного стану — наприклад, 010101110000... — що відповідає стану-кандидату $x$ для нашої вихідної оптимізаційної задачі ($\max f(x)$ або $\min f(x)$ залежно від задачі).

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

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()

Візуалізація найкращого розрізу

З оптимального бітового рядка можна візуалізувати цей розріз на вихідному графі.

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)

І обчислимо значення розрізу. Результат не є оптимальним через шуми (значення розрізу оптимального розв'язку дорівнює 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 завершується. Як адаптувати QAOA до задач утилітарного масштабу, дізнаєшся у «Частині 2: масштабуємо!» туторіалу Квантовий наближений алгоритм оптимізації.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'