Kvantové algoritmy: Variační kvantové algoritmy

poznámka

Takashi Imamichi (24. května 2024)

Stáhni si PDF původní přednášky. Některé ukázky kódu mohou být zastaralé, protože se jedná o statické obrázky.

Přibližná doba výpočtu na QPU pro tento experiment je 9 minut (testováno na procesoru Eagle).

(Tento notebook nemusí být vyhodnocen v časovém limitu povoleném v Open Planu. Prosím, využívej kvantové výpočetní zdroje ohleduplně.)

1. Úvod

Tento tutoriál poskytuje přehled hybridního kvantově-klasického algoritmu, se zaměřením na variační kvantový eigensolver (VQE) a kvantový aproximační optimalizační algoritmus (QAOA). Hlavním cílem těchto algoritmů je řešit optimalizační problémy pomocí kvantových Circuit s parametrizovanými kvantovými Gate.

Navzdory pokrokům v kvantovém výpočetnictví přítomnost šumu v současných kvantových zařízeních ztěžuje získávání smysluplných výsledků z hlubokých kvantových Circuit. Aby se překonala tato výzva, přijímají VQE a QAOA hybridní kvantově-klasický přístup, který zahrnuje iterativní spouštění relativně krátkých kvantových Circuit pomocí kvantového výpočtu a optimalizaci parametrů cílových parametrizovaných kvantových Circuit pomocí klasického výpočtu.

QAOA má potenciál poskytovat optimální řešení cílových problémů v užitkovém měřítku, díky aplikaci různých technik zmírňování a potlačování chyb. VQE má mnoho aplikací (například v kvantové chemii), kde je méně škálovatelný. Nicméně se objevila celá řada přístupů souvisejících s vlastními hodnotami, které doplňují a rozšiřují VQE, včetně diagonalizace Krylovova podprostoru a diagonalizace kvantových vzorků (SQD). Porozumění VQE je důležitým prvním krokem k pochopení širokého spektra hybridních klasicko-kvantových algoritmů, které se objevily.

Tento modul popisuje základní koncepty a implementaci VQE a QAOA. Další tutoriály prozkoumají pokročilá témata a techniky pro škálování těchto algoritmů. Pro spuštění tohoto notebooku potřebuješ ve svém prostředí následující knihovnu. Pokud ji ještě nemáš nainstalovanou, můžeš ji nainstalovat odkomentováním a spuštěním následující buňky.

# 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. Výpočet minimální vlastní hodnoty jednoduchého Hamiltoniánu

Začneme aplikací VQE na velmi jednoduchý případ, abychom viděli, jak funguje. Vypočítáme minimální vlastní hodnotu matice Pauli $Z$ pomocí VQE. Začneme importem některých obecných balíčků.

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

Nyní definujeme operátor, který nás zajímá, a zobrazíme ho ve formě matice.

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

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

Vlastní hodnoty je snadné získat klasicky, takže si můžeme zkontrolovat výsledky. S rostoucí užitkovostí se to může stát obtížným. Zde použijeme numpy.

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

Eigenvalues: [-1.  1.]

Abychom získali vlastní hodnoty pomocí variačního kvantového algoritmu, zkonstruujeme Circuit s Gate, které přijímají variační parametry:

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

Pokud chceme odhadnout střední hodnotu operátoru (jako je $Z$), použijeme Estimator. Pokud chceme zobrazit stavy systému, použijeme Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Pomocí Sampler můžeme vypočítat počty bitových řetězců 0 a 1 s náhodnými hodnotami parametrů [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}

Víme, že střední hodnotu Z můžeme vypočítat jako $\langle Z \rangle = p_0 - p_1$ s pravděpodobnostmi $\{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

Tento Circuit fungoval, ale zvolené hodnoty parametrů neodpovídaly stavu s velmi nízkou energií (nebo nízkou vlastní hodnotou). Získaná vlastní hodnota je podstatně vyšší než minimum. Výsledek je podobný při použití Estimatoru.

Všimni si, že Estimator přijímá kvantové Circuit bez měření.

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

array(0.54030231)

Budeme muset prohledávat parametry a najít ty, které dávají nejnižší vlastní hodnotu. Vytvoříme funkci, která přijme hodnoty parametrů variační formy a vrátí střední hodnotu $\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

Použijme funkci minimize z SciPy k nalezení minimální vlastní hodnoty 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 Cvičení

Vypočítej minimální vlastní hodnotu $Z \otimes Z$ 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

Řešení cvičení

Definujeme operátor, který nás zajímá, a zobrazíme ho v maticové podobě.

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

Abychom mohli získat vlastní hodnoty pomocí variačního kvantového algoritmu, sestavíme Circuit s Gate, které přijímají variační parametry:

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

Pokud chceme odhadnout střední hodnotu operátoru (například $Z \otimes Z$), použijeme Estimator. Pokud chceme sledovat stavy systému, použijeme 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

Tento Circuit fungoval, ale zvolené hodnoty parametrů neodpovídaly stavu s velmi nízkou energií (nebo nízkou vlastní hodnotou). Získaná vlastní hodnota je podstatně vyšší než minimum. Výsledek je podobný i při použití Estimatoru.

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

Budeme muset prohledat parametry a najít ty, které přinesou nejnižší vlastní hodnotu.

# 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

Získali jsme vlastní hodnotu velmi blízkou minimu, které nám poskytl 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. Kvantová optimalizace s Qiskit patterns

V tomto návodu se naučíš pracovat s Qiskit patterns a kvantovou aproximační optimalizací. Qiskit pattern je intuitivní, opakovatelná sada kroků pro implementaci pracovního postupu kvantového výpočtu: Tyto patterns aplikujeme v kontextu kombinatorické optimalizace a ukážeme si, jak řešit problém Max-Cut pomocí Quantum Approximate Optimization Algorithm (QAOA) — hybridní (kvantově-klasické) iterativní metody.

Všimni si, že tato část o QAOA vychází z „Části 1: Malá QAOA" tutoriálu Quantum approximate optimization algorithm. Podívej se do tutoriálu, kde se dozvíš, jak ho škálovat.

3.1 (Malá) Qiskit pattern pro optimalizaci

Tato část využije malý problém Max-Cut k ilustraci kroků potřebných k řešení optimalizačního problému pomocí kvantového počítače. Problém Max-Cut je optimalizační problém, který je obtížné řešit (konkrétně jde o NP-těžký problém) s řadou různých aplikací v oblasti shlukování, síťové vědy a statistické fyziky. Tento tutoriál uvažuje graf uzlů propojených hranami a cílem je rozdělit uzly do dvou množin „přeříznutím" hran tak, aby byl maximalizován počet přeříznutých hran. Abychom získali kontext před mapováním tohoto problému na kvantový algoritmus, můžeš lépe pochopit, jak se problém Max-Cut stává klasickým kombinatorickým optimalizačním problémem, a to tak, že nejprve uvažuješ minimalizaci funkce $f(x)$

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

kde vstup $x$ je vektor, jehož složky odpovídají jednotlivým uzlům grafu. Každou z těchto složek pak omezíš na hodnotu $0$ nebo $1$ (která představuje zahrnutí nebo nezahrnutí do řezu). Tento malý příklad využívá graf s $n=5$ uzly.

Mohl bys napsat funkci dvojice uzlů $i,j$, která indikuje, zda odpovídající hrana $(i,j)$ leží v řezu. Například funkce $x_i + x_j - 2 x_i x_j$ je rovna 1 pouze tehdy, když je právě jedno z $x_i$ nebo $x_j$ rovno 1 (což znamená, že hrana je v řezu), a jinak je nulová. Problém maximalizace hran v řezu lze formulovat jako

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

což lze přepsat jako minimalizaci ve tvaru

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

Minimum $f(x)$ nastane v tomto případě tehdy, když je počet hran procházejících řezem maximální. Jak vidíš, zatím tu není nic, co by se týkalo kvantového výpočtu. Je potřeba tento problém přeformulovat do podoby, které kvantový počítač rozumí. Inicializuj svůj problém vytvořením grafu s $n=5$ uzly.

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. Mapování klasických vstupů na kvantový problém

Prvním krokem vzoru je mapování klasického problému (grafu) na kvantové Circuity a operátory. K tomu je třeba provést tři hlavní kroky:

1. Využít sérii matematických přeformulací, aby bylo možné tento problém reprezentovat pomocí notace problémů Quadratic Unconstrained Binary Optimization (QUBO).
2. Přepsat optimalizační problém jako Hamiltonian, jehož základní stav odpovídá řešení, které minimalizuje účelovou funkci.
3. Vytvořit kvantový Circuit, který připraví základní stav tohoto Hamiltoniánu procesem podobným kvantovému žíhání.

Poznámka: V metodologii QAOA chceš nakonec mít operátor (Hamiltonian), který reprezentuje účelovou funkci našeho hybridního algoritmu, a také parametrizovaný Circuit (Ansatz), který reprezentuje kvantové stavy s kandidátními řešeními problému. Můžeš vzorkovat z těchto kandidátních stavů a poté je vyhodnotit pomocí účelové funkce.

Graf → optimalizační problém

Prvním krokem mapování je změna notace. Následující vyjádření formuluje problém v QUBO notaci:

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

kde $Q$ je matice $n\times n$ reálných čísel, $n$ odpovídá počtu uzlů v grafu, $x$ je vektor binárních proměnných zavedených výše a $x^T$ označuje transpozici vektoru $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

Optimalizační problém → Hamiltonian

Problém QUBO pak lze přeformulovat jako Hamiltonian (zde matice reprezentující energii systému):

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

Kroky přeformulace z QAOA problému na Hamiltonian

Abychom ukázali, jak lze QAOA problém přepsat tímto způsobem, nejprve nahraď binární proměnné $x_i$ novou sadou proměnných $z_i\in\{-1, 1\}$ pomocí

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

Zde vidíš, že pokud je $x_i$ rovno $0$, pak $z_i$ musí být $1$. Když se $x_i$ nahradí $z_i$ v optimalizačním problému ($x^TQx$), lze získat ekvivalentní formulaci.

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

Pokud nyní definujeme $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$, odstraníme předfaktor a konstantní člen $n^2$, dospějeme ke dvěma ekvivalentním formulacím téhož optimalizačního problému.

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

Zde $b$ závisí na $Q$. Všimni si, že abychom získali $z^TQz + b^Tz$, vynechali jsme faktor 1/4 a konstantní posun $n^2$, které nehrají v optimalizaci žádnou roli.

Nyní, abychom získali kvantovou formulaci problému, povýšíme proměnné $z_i$ na Pauliho matici $Z$, tedy matici $2\times 2$ ve tvaru

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

Když tyto matice dosadíš do výše uvedeného optimalizačního problému, obdržíš následující Hamiltonian

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

Připomeň si také, že matice $Z$ jsou vloženy do výpočetního prostoru kvantového počítače, tedy do Hilbertova prostoru velikosti $2^n\times 2^n$. Proto je třeba chápat výrazy jako $Z_iZ_j$ jako tenzorový součin $Z_i\otimes Z_j$ vložený do Hilbertova prostoru $2^n\times 2^n$. Například v problému s pěti rozhodovacími proměnnými je člen $Z_1Z_3$ chápán jako $I\otimes Z_3\otimes I\otimes Z_1\otimes I$, kde $I$ je jednotková matice $2\times 2$.

Tento Hamiltonian se nazývá Hamiltonian účelové funkce. Má tu vlastnost, že jeho základní stav odpovídá řešení, které minimalizuje účelovou funkci $f(x)$. Abys tedy vyřešil svůj optimalizační problém, potřebuješ nyní připravit základní stav $H_C$ (nebo stav s velkou mírou překryvu s ním) na kvantovém počítači. Vzorkování z tohoto stavu pak s vysokou pravděpodobností poskytne řešení $\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 → kvantový Circuit

Hamiltonian $H_C$ obsahuje kvantovou definici tvého problému. Nyní můžeš vytvořit kvantový Circuit, který pomůže vzorkovat dobrá řešení z kvantového počítače. QAOA je inspirováno kvantovým žíháním a v kvantovém Circuitu aplikuje střídající se vrstvy operátorů.

Obecná myšlenka spočívá v tom, že začneš v základním stavu známého systému, $H^{\otimes n}|0\rangle$ výše, a poté systém nasměruješ do základního stavu operátoru nákladů, který tě zajímá. To se provede aplikací operátorů $\exp\{-i\gamma_k H_C\}$ a $\exp\{-i\beta_k H_m\}$ s úhly $\gamma_1,...,\gamma_p$ a $\beta_1,...,\beta_p~$.

Kvantový Circuit, který vygeneruješ, je parametrizován hodnotami $\gamma_i$ a $\beta_i$, takže můžeš vyzkoušet různé hodnoty $\gamma_i$ a $\beta_i$ a vzorkovat z výsledného stavu. V tomto případě vyzkoušíme příklad s 1 vrstvou QAOA, která obsahuje dva parametry: $\gamma_1$ a $\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. Optimalizace Circuitů pro spuštění na kvantovém hardwaru

Výše uvedený Circuit obsahuje řadu abstrakcí užitečných pro přemýšlení o kvantových algoritmech, které však na hardwaru nelze přímo spustit. Aby bylo možné spustit Circuit na QPU, musí projít řadou operací tvořících krok transpilace nebo optimalizace Circuitu v rámci vzoru.

Knihovna Qiskit nabízí řadu transpilačních průchodů, které pokrývají širokou škálu transformací Circuitů. Je třeba zajistit, aby byl tvůj Circuit optimalizován pro daný účel.

Transpilace může zahrnovat několik kroků, například:

• Počáteční mapování Qubitů v Circuitu (například rozhodovacích proměnných) na fyzické Qubity zařízení.
• Rozbalení instrukcí v kvantovém Circuitu na hardwarově nativní instrukce, kterým Backend rozumí.
• Směrování Qubitů v Circuitu, které spolu interagují, na fyzické Qubity, jež jsou vzájemně sousední.
• Potlačení chyb přidáním jednoQubitových Gate pro potlačení šumu pomocí dynamického oddělení.

Více informací o transpilaci najdeš v naší dokumentaci.

Následující kód transformuje a optimalizuje abstraktní Circuit do formátu připraveného ke spuštění na jednom ze zařízení dostupných přes cloud pomocí služby Qiskit IBM® Runtime.

Všimni si, že své programy můžeš lokálně otestovat v „režimu lokálního testování" před jejich odesláním na skutečné kvantové počítače. Více informací o režimu lokálního testování najdeš v dokumentaci.

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. Spuštění pomocí primitiv Qiskit

V pracovním postupu QAOA jsou optimální parametry QAOA nalezeny v iterační optimalizační smyčce, která spouští řadu vyhodnocení Circuitu a pomocí klasického optimalizátoru hledá optimální parametry $\beta_k$ a $\gamma_k$. Tato spouštěcí smyčka se provádí v následujících krocích:

1. Definování počátečních parametrů
2. Vytvoření nové Session obsahující optimalizační smyčku a primitivum použité k vzorkování Circuitu
3. Jakmile je nalezena optimální sada parametrů, spustí se Circuit naposledy, aby bylo získáno výsledné rozdělení použité v kroku následného zpracování.

Definování Circuitu s počátečními parametry

Začneme libovolně zvolenými parametry.

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

Definování Backendu a spouštěcího primitiva

K interakci s Backendy IBM® použij primitiva Qiskit Runtime. Dvě dostupná primitiva jsou Sampler a Estimator a volba primitiva závisí na tom, jaký typ měření chceš na kvantovém počítači provést. Pro minimalizaci $H_C$ použij Estimator, protože měření nákladové funkce je jednoduše střední hodnota $\langle H_C \rangle$.

Spuštění

Primitivy nabízejí různé režimy spouštění pro plánování úloh na kvantových zařízeních a pracovní postup QAOA běží iterativně v rámci Session. Funkci nákladů založenou na Sampleru lze zapojit do minimalizační rutiny SciPy a nalézt optimální 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

Optimalizátor dokázal snížit náklady a nalézt lepší parametry pro Circuit.

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

Jakmile najdeš optimální parametry pro Circuit, můžeš je přiřadit a vzorkovat výsledné rozdělení získané s optimalizovanými parametry. Právě zde by měla být použita primitiva Sampler, protože jde o rozdělení pravděpodobnosti měření bitových řetězců, které odpovídají optimálnímu řezu grafu.

Poznámka: To znamená připravit kvantový stav $\psi$ v počítači a pak ho změřit. Měření zkolabuje stav do jediného výpočetního bázového stavu – například 010101110000... – který odpovídá kandidátskému řešení $x$ našeho původního optimalizačního problému ($\max f(x)$ nebo $\min f(x)$ podle zadání).

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. Post-processing, vrácení výsledku v klasickém formátu

Krok post-processingu interpretuje výstup vzorkování a vrací řešení původního problému. V tomto případě tě zajímá bitový řetězec s nejvyšší pravděpodobností, protože ten určuje optimální řez. Symetrie problému umožňují čtyři možná řešení a proces vzorkování vrátí jedno z nich s o něco vyšší pravděpodobností, ale v níže zobrazeném rozdělení vidíš, že čtyři bitové řetězce jsou výrazně pravděpodobnější než ostatní.

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

Vizualizace nejlepšího řezu

Z optimálního bitového řetězce pak můžeš vizualizovat tento řez na původním grafu.

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)

A vypočítat hodnotu řezu. Řešení není optimální kvůli šumu (hodnota řezu optimálního řešení je 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

Tím je malý tutoriál k QAOA uzavřen. Naučíš se, jak přizpůsobit QAOA na úrovni utility, v části „Part 2: scale it up!" tutoriálu Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'