Algorytmy kwantowe: Wariacyjne algorytmy kwantowe

uwaga

Takashi Imamichi (24 maja 2024)

Pobierz plik pdf oryginalnego wykładu. Należy pamiętać, że niektóre fragmenty kodu mogą być już przestarzałe, ponieważ są to statyczne obrazy.

Przybliżony czas pracy QPU potrzebny do uruchomienia tego eksperymentu to 9 minut (testowane na procesorze Eagle).

(ten notatnik może nie zostać wykonany w czasie dozwolonym w planie Open Plan. Prosimy rozsądnie korzystać z zasobów obliczeń kwantowych.)

1. Wprowadzenie

Ten samouczek przedstawia przegląd hybrydowego algorytmu kwantowo-klasycznego, ze szczególnym uwzględnieniem wariacyjnego kwantowego rozwiązywacza własnego (VQE) oraz kwantowego przybliżonego algorytmu optymalizacyjnego (QAOA). Głównym celem tych algorytmów jest rozwiązywanie problemów optymalizacyjnych poprzez wykorzystanie obwodów kwantowych z parametryzowanymi bramkami kwantowymi.

Pomimo postępów w obliczeniach kwantowych, obecność szumu w dzisiejszych urządzeniach kwantowych sprawia, że wyodrębnianie znaczących wyników z głębokich obwodów kwantowych jest trudne. Aby pokonać to wyzwanie, VQE i QAOA przyjmują hybrydowe podejście kwantowo-klasyczne, które polega na iteracyjnym wykonywaniu stosunkowo krótkich obwodów kwantowych za pomocą obliczeń kwantowych i optymalizowaniu parametrów docelowych parametryzowanych obwodów kwantowych za pomocą obliczeń klasycznych.

QAOA ma potencjał, aby dostarczać optymalne rozwiązania problemów docelowych na skalę użyteczną, dzięki zastosowaniu różnych technik łagodzenia i tłumienia błędów. VQE ma wiele zastosowań (takich jak chemia kwantowa), w których jest mniej skalowalny. Pojawiło się jednak wiele podejść związanych z wartościami własnymi, które uzupełniają i rozszerzają VQE, w tym diagonalizacja w podprzestrzeni Krylov oraz kwantowa diagonalizacja oparta na próbkowaniu (SQD). Zrozumienie VQE jest ważnym pierwszym krokiem w zrozumieniu szerokiego zakresu hybrydowych algorytmów klasyczno-kwantowych, które się pojawiły.

Ten moduł opisuje podstawowe koncepcje i implementację VQE oraz QAOA. Kolejne samouczki będą omawiać zaawansowane tematy i techniki skalowania tych algorytmów. Do uruchomienia tego notatnika wymagana jest w środowisku następująca biblioteka. Jeśli nie została jeszcze zainstalowana, można ją zainstalować odkomentowując i uruchamiając poniższą komórkę.

# 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. Obliczanie minimalnej wartości własnej prostego Hamiltonianu

Zaczniemy od zastosowania VQE do bardzo prostego przypadku, aby zobaczyć, jak to działa. Obliczymy minimalną wartość własną macierzy Pauli $Z$ za pomocą VQE. Zaczniemy od zaimportowania kilku ogólnych pakietów.

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

Teraz zdefiniujemy interesujący nas operator i wyświetlimy go w postaci macierzowej.

op = SparsePauliOp("Z")
op.to_matrix()

array([[ 1.+0.j,  0.+0.j],
       [ 0.+0.j, -1.+0.j]])

Łatwo jest uzyskać wartości własne klasycznie, więc możemy zweryfikować naszą pracę. Może to stać się trudne w miarę zbliżania się do skali użytecznej. Tutaj używamy numpy.

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

Eigenvalues: [-1.  1.]

Aby uzyskać wartości własne za pomocą wariacyjnego algorytmu kwantowego, konstruujemy obwód z bramkami, które przyjmują parametry wariacyjne:

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

Jeśli chcemy oszacować wartość oczekiwaną operatora (takiego jak $Z$), powinniśmy użyć Estimator. Jeśli chcemy przyjrzeć się stanom systemu, używamy Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Za pomocą Samplera możemy obliczyć liczebności ciągów bitów 0 i 1 z losowymi wartościami parametrów [1, 2, 3].

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

Wiemy, że możemy obliczyć wartość oczekiwaną Z jako $\langle Z \rangle = p_0 - p_1$ przy prawdopodobieństwach $\{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

Ten obwód zadziałał, ale wybrane wartości parametrów nie odpowiadały stanowi o bardzo niskiej energii (lub o niskiej wartości własnej). Uzyskana wartość własna jest znacznie wyższa niż minimum. Wynik jest podobny podczas używania estimatora.

Należy zauważyć, że Estimator przyjmuje obwody kwantowe bez pomiarów.

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

array(0.54030231)

Będziemy musieli przeszukać parametry i znaleźć te, które dają najniższą wartość własną. Tworzymy funkcję, która przyjmuje wartości parametrów formy wariacyjnej i zwraca wartość oczekiwaną $\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

Zastosujmy funkcję minimize z biblioteki SciPy, aby znaleźć minimalną wartość własną 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 Ćwiczenie

Oblicz minimalną wartość własną $Z \otimes Z$ za pomocą 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

Rozwiązania ćwiczenia

Definiujemy interesujący nas operator i przedstawiamy go w postaci macierzowej.

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

Aby uzyskać wartości własne przy użyciu wariacyjnego algorytmu kwantowego, konstruujemy obwód z bramkami, które przyjmują parametry wariacyjne:

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

Jeżeli chcemy oszacować wartość oczekiwaną operatora (takiego jak $Z \otimes Z$), użyjemy Estimator. Jeżeli chcemy przyjrzeć się stanom systemu, używamy 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

Ten obwód zadziałał, ale wybrane wartości parametrów nie odpowiadały stanowi o bardzo niskiej energii (lub niskiej wartości własnej). Uzyskana wartość własna jest znacznie wyższa od minimum. Wynik jest podobny przy użyciu estymatora.

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

Musimy przeszukać parametry i znaleźć te, które dają najniższą wartość własną.

# 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

Uzyskaliśmy wartość własną niezwykle bliską minimum podanemu przez 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. Optymalizacja kwantowa z wykorzystaniem wzorców Qiskit

W tym poradniku poznamy wzorce Qiskit oraz kwantową przybliżoną optymalizację. Wzorzec Qiskit to intuicyjny, powtarzalny zestaw kroków do implementacji przepływu pracy w obliczeniach kwantowych: Zastosujemy te wzorce w kontekście optymalizacji kombinatorycznej i pokażemy, jak rozwiązać problem Max-Cut przy użyciu Quantum Approximate Optimization Algorithm (QAOA) — hybrydowej (kwantowo-klasycznej) metody iteracyjnej.

Zwróć uwagę, że ta część dotycząca QAOA jest oparta na "Part 1: Small-scale QAOA" z tutoriala Quantum approximate optimization algorithm. Zajrzyj do tutoriala, aby dowiedzieć się, jak skalować to rozwiązanie.

3.1 Wzorzec Qiskit dla optymalizacji (w małej skali)

W tej części wykorzystamy problem Max-Cut w małej skali, aby zilustrować kroki wymagane do rozwiązania problemu optymalizacyjnego z użyciem komputera kwantowego. Problem Max-Cut to problem optymalizacyjny, który jest trudny do rozwiązania (konkretniej, jest to problem NP-trudny) i ma wiele różnych zastosowań w klasteryzacji, naukach o sieciach i fizyce statystycznej. Ten tutorial rozważa graf węzłów połączonych krawędziami i ma na celu podział węzłów na dwa zbiory poprzez "przecinanie" krawędzi tak, aby liczba przeciętych krawędzi była maksymalna. Aby zapewnić pewien kontekst przed odwzorowaniem tego problemu na algorytm kwantowy, możesz lepiej zrozumieć, w jaki sposób problem Max-Cut staje się klasycznym problemem optymalizacji kombinatorycznej, rozważając najpierw minimalizację funkcji $f(x)$

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

gdzie wejście $x$ jest wektorem, którego komponenty odpowiadają poszczególnym węzłom grafu. Następnie nałóż ograniczenie, aby każdy z tych komponentów był równy $0$ albo $1$ (co reprezentuje to, czy element jest włączony do przecięcia, czy nie). Ten przykład w małej skali wykorzystuje graf z $n=5$ węzłami.

Można zapisać funkcję pary węzłów $i,j$, która wskazuje, czy odpowiadająca im krawędź $(i,j)$ znajduje się w przecięciu. Na przykład funkcja $x_i + x_j - 2 x_i x_j$ przyjmuje wartość 1 tylko wtedy, gdy dokładnie jedno z $x_i$ lub $x_j$ jest równe 1 (co oznacza, że krawędź jest w przecięciu), a w przeciwnym razie wartość 0. Problem maksymalizacji liczby krawędzi w przecięciu można sformułować jako

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

co można przepisać jako problem minimalizacji w postaci

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

Minimum $f(x)$ w tym przypadku występuje wtedy, gdy liczba krawędzi przeciętych przez cięcie jest maksymalna. Jak widać, nie ma tu jeszcze nic związanego z obliczeniami kwantowymi. Musisz przeformułować ten problem w coś, co komputer kwantowy potrafi zrozumieć. Zainicjuj swój problem, tworząc graf z $n=5$ węzłami.

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 Krok 1. Odwzorowanie klasycznych danych wejściowych na problem kwantowy

Pierwszym krokiem wzorca jest odwzorowanie klasycznego problemu (grafu) na kwantowe obwody i operatory. Aby to zrobić, należy wykonać trzy główne kroki:

1. Wykorzystać serię matematycznych przekształceń, aby przedstawić ten problem w notacji problemów Quadratic Unconstrained Binary Optimization (QUBO).
2. Przepisać problem optymalizacji jako Hamiltonian, którego stan podstawowy odpowiada rozwiązaniu minimalizującemu funkcję kosztu.
3. Utworzyć obwód kwantowy, który przygotuje stan podstawowy tego Hamiltonianu za pomocą procesu podobnego do wyżarzania kwantowego.

Uwaga: W metodologii QAOA ostatecznie chcesz mieć operator (Hamiltonian), który reprezentuje funkcję kosztu naszego algorytmu hybrydowego, a także sparametryzowany obwód (Ansatz), który reprezentuje stany kwantowe z kandydatami na rozwiązania problemu. Możesz próbkować te stany kandydujące, a następnie oceniać je za pomocą funkcji kosztu.

Graf → problem optymalizacyjny

Pierwszym krokiem odwzorowania jest zmiana notacji. Poniżej wyrażono problem w notacji QUBO:

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

gdzie $Q$ jest macierzą $n\times n$ liczb rzeczywistych, $n$ odpowiada liczbie węzłów w grafie, $x$ jest wektorem zmiennych binarnych wprowadzonym powyżej, a $x^T$ oznacza transpozycję wektora $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

Problem optymalizacyjny → Hamiltonian

Możesz następnie przeformułować problem QUBO jako Hamiltonian (tutaj macierz, która reprezentuje energię układu):

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

Kroki przeformułowania z problemu QAOA do Hamiltonianu

Aby pokazać, w jaki sposób problem QAOA można zapisać w ten sposób, najpierw zastąp zmienne binarne $x_i$ nowym zbiorem zmiennych $z_i\in\{-1, 1\}$ poprzez

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

Tutaj widać, że jeśli $x_i$ wynosi $0$, to $z_i$ musi być równe $1$. Gdy $x_i$ zostaną zastąpione przez $z_i$ w problemie optymalizacyjnym ($x^TQx$), otrzymamy równoważne sformułowanie.

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

Teraz, jeśli zdefiniujemy $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, usuniemy prefaktor oraz stałą $n^2$, otrzymamy dwa równoważne sformułowania tego samego problemu optymalizacyjnego.

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

Tutaj $b$ zależy od $Q$. Zauważ, że aby uzyskać $z^TQz + b^Tz$, pominęliśmy czynnik 1/4 oraz stałe przesunięcie $n^2$, które nie odgrywają roli w optymalizacji.

Teraz, aby uzyskać kwantowe sformułowanie problemu, podnieś zmienne $z_i$ do macierzy Pauliego $Z$, takiej jak macierz $2\times 2$ postaci

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

Gdy podstawisz te macierze do powyższego problemu optymalizacji, otrzymasz następujący Hamiltonian

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

Przypomnij sobie również, że macierze $Z$ są osadzone w przestrzeni obliczeniowej komputera kwantowego, tj. w przestrzeni Hilberta o rozmiarze $2^n\times 2^n$. Dlatego wyrażenia takie jak $Z_iZ_j$ należy rozumieć jako iloczyn tensorowy $Z_i\otimes Z_j$ osadzony w przestrzeni Hilberta $2^n\times 2^n$. Na przykład w problemie z pięcioma zmiennymi decyzyjnymi wyrażenie $Z_1Z_3$ rozumiane jest jako $I\otimes Z_3\otimes I\otimes Z_1\otimes I$, gdzie $I$ jest macierzą jednostkową $2\times 2$.

Ten Hamiltonian nazywany jest Hamiltonianem funkcji kosztu. Ma on tę właściwość, że jego stan podstawowy odpowiada rozwiązaniu, które minimalizuje funkcję kosztu $f(x)$. Dlatego, aby rozwiązać problem optymalizacji, musisz teraz przygotować stan podstawowy $H_C$ (lub stan o dużym nakładaniu się z nim) na komputerze kwantowym. Następnie próbkowanie z tego stanu z dużym prawdopodobieństwem da rozwiązanie $\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

Hamiltonian → obwód kwantowy

Hamiltonian $H_C$ zawiera kwantową definicję twojego problemu. Teraz możesz stworzyć obwód kwantowy, który pomoże próbkować dobre rozwiązania z komputera kwantowego. QAOA jest inspirowany wyżarzaniem kwantowym i stosuje naprzemienne warstwy operatorów w obwodzie kwantowym.

Ogólna idea polega na tym, aby rozpocząć w stanie podstawowym znanego układu, $H^{\otimes n}|0\rangle$ powyżej, a następnie skierować układ do stanu podstawowego interesującego cię operatora kosztu. Odbywa się to poprzez zastosowanie operatorów $\exp\{-i\gamma_k H_C\}$ i $\exp\{-i\beta_k H_m\}$ z kątami $\gamma_1,...,\gamma_p$ i $\beta_1,...,\beta_p~$.

Generowany obwód kwantowy jest parametryzowany przez $\gamma_i$ i $\beta_i$, więc możesz wypróbować różne wartości $\gamma_i$ i $\beta_i$ i próbkować z uzyskanego stanu. W tym przypadku wypróbujemy przykład z 1 warstwą QAOA, która zawiera dwa parametry: $\gamma_1$ i $\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 Krok 2. Optymalizacja obwodów do wykonania na sprzęcie kwantowym

Powyższy obwód zawiera szereg abstrakcji przydatnych do rozważań o algorytmach kwantowych, ale niemożliwych do uruchomienia na sprzęcie. Aby móc uruchomić go na QPU, obwód musi przejść szereg operacji, które składają się na krok transpilacji lub optymalizacji obwodu wzorca.

Biblioteka Qiskit oferuje szereg przebiegów transpilacji, które obsługują szeroki zakres transformacji obwodów. Musisz upewnić się, że twój obwód jest zoptymalizowany do twojego celu.

Transpilacja może obejmować kilka kroków, takich jak:

• Początkowe mapowanie kubitów w obwodzie (takich jak zmienne decyzyjne) na fizyczne kubity w urządzeniu.
• Rozwijanie instrukcji w obwodzie kwantowym do instrukcji natywnych dla sprzętu, które backend rozumie.
• Routing wszelkich kubitów w obwodzie, które oddziałują na siebie, do fizycznych kubitów sąsiadujących ze sobą.
• Tłumienie błędów poprzez dodawanie bramek jednokubitowych w celu tłumienia szumu za pomocą dynamicznego rozprzęgania.

Więcej informacji o transpilacji znajduje się w naszej dokumentacji.

Poniższy kod przekształca i optymalizuje abstrakcyjny obwód do formatu gotowego do wykonania na jednym z urządzeń dostępnych przez chmurę za pomocą usługi Qiskit IBM® Runtime.

Zauważ, że możesz testować swoje programy lokalnie przy użyciu "trybu testowania lokalnego" przed wysłaniem ich do prawdziwych komputerów kwantowych. Więcej informacji o trybie testowania lokalnego znajduje się w dokumentacji.

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 Krok 3. Wykonanie przy użyciu prymitywów Qiskit

W przepływie pracy QAOA optymalne parametry QAOA są znajdowane w iteracyjnej pętli optymalizacji, która uruchamia serię ewaluacji obwodu i używa klasycznego optymalizatora do znalezienia optymalnych parametrów $\beta_k$ i $\gamma_k$. Ta pętla wykonawcza jest realizowana za pomocą następujących kroków:

1. Zdefiniuj parametry początkowe
2. Utwórz nową Session zawierającą pętlę optymalizacji i prymityw używany do próbkowania obwodu
3. Gdy zostanie znaleziony optymalny zestaw parametrów, wykonaj obwód po raz ostatni, aby uzyskać ostateczny rozkład, który zostanie użyty w kroku przetwarzania końcowego.

Zdefiniuj obwód z parametrami początkowymi

Zaczynamy od arbitralnie wybranych parametrów.

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

Zdefiniuj backend i prymityw wykonania

Użyj prymitywów Qiskit Runtime, aby komunikować się z backendami IBM®. Dwa prymitywy to Sampler i Estimator, a wybór prymitywu zależy od tego, jaki rodzaj pomiaru chcesz uruchomić na komputerze kwantowym. Do minimalizacji $H_C$ użyj Estimatora, ponieważ pomiar funkcji kosztu jest po prostu wartością oczekiwaną $\langle H_C \rangle$.

Uruchom

Prymitywy oferują różnorodne tryby wykonania do planowania obciążeń na urządzeniach kwantowych, a przepływ pracy QAOA działa iteracyjnie w sesji. Możesz podłączyć funkcję kosztu opartą na Samplerze do procedury minimalizującej SciPy, aby znaleźć optymalne parametry.

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

Optymalizator był w stanie zmniejszyć koszt i znaleźć lepsze parametry dla obwodu.

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

Gdy już znajdziesz optymalne parametry dla obwodu, możesz je przypisać i próbkować końcowy rozkład uzyskany z zoptymalizowanymi parametrami. W tym miejscu należy użyć prymitywu Sampler, ponieważ to właśnie rozkład prawdopodobieństwa pomiarów ciągów bitów odpowiada optymalnemu przecięciu grafu.

Uwaga: Oznacza to przygotowanie stanu kwantowego $\psi$ w komputerze, a następnie jego pomiar. Pomiar spowoduje kolaps stanu do pojedynczego stanu bazy obliczeniowej — na przykład 010101110000... — który odpowiada kandydującemu rozwiązaniu $x$ naszego początkowego problemu optymalizacyjnego ($\max f(x)$ lub $\min f(x)$ w zależności od zadania).

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 Krok 4. Przetwórz końcowo, zwróć wynik w formacie klasycznym

Krok przetwarzania końcowego interpretuje wynik próbkowania, aby zwrócić rozwiązanie oryginalnego problemu. W tym przypadku interesuje Cię ciąg bitów o najwyższym prawdopodobieństwie, ponieważ to on wyznacza optymalne przecięcie. Symetrie występujące w problemie umożliwiają istnienie czterech możliwych rozwiązań, a proces próbkowania zwróci jedno z nich z nieznacznie wyższym prawdopodobieństwem, jednak na poniższym wykresie rozkładu widać, że cztery ciągi bitów są wyraźnie bardziej prawdopodobne niż pozostałe.

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

Wizualizacja najlepszego przecięcia

Na podstawie optymalnego ciągu bitów możesz następnie zwizualizować to przecięcie na oryginalnym grafie.

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)

Oraz obliczyć wartość przecięcia. Rozwiązanie nie jest optymalne z powodu szumu (wartość przecięcia dla rozwiązania optymalnego wynosi 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

Na tym kończy się samouczek QAOA w małej skali. Dowiesz się, jak dostosować QAOA do skali użytkowej w części „Part 2: scale it up!" samouczka Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'