Quantenalgorithmen: Variationelle Quantenalgorithmen

hinweis

Takashi Imamichi (24. Mai 2024)

PDF herunterladen der Originalvorlesung. Beachte, dass einige Code-Snippets veraltet sein könnten, da es sich um statische Bilder handelt.

Die ungefähre QPU-Zeit für dieses Experiment beträgt 9 Minuten (getestet auf einem Eagle-Prozessor).

(Dieses Notebook lässt sich möglicherweise nicht in der im Open Plan erlaubten Zeit auswerten. Bitte nutze Quantencomputing-Ressourcen verantwortungsvoll.)

1. Einführung

Dieses Tutorial bietet einen Überblick über einen hybriden quantenklassischen Algorithmus, insbesondere den Variational Quantum Eigensolver (VQE) und den Quantum Approximate Optimization Algorithm (QAOA). Das primäre Ziel dieser Algorithmen ist die Lösung von Optimierungsproblemen durch den Einsatz von Quantenschaltkreisen mit parametrisierten Quantengattern.

Trotz der Fortschritte im Quantencomputing erschwert das Rauschen in aktuellen Quantengeräten die Extraktion aussagekräftiger Ergebnisse aus tiefen Quantenschaltkreisen. Um diese Herausforderung zu überwinden, verfolgen VQE und QAOA einen hybriden quantenklassischen Ansatz, der das iterative Ausführen relativ kurzer Quantenschaltkreise mittels Quantenberechnung und die Optimierung der Parameter der parametrisierten Quantenschaltkreise mittels klassischer Berechnung umfasst.

QAOA hat das Potenzial, optimale Lösungen für Zielprobleme im Utility-Maßstab bereitzustellen, dank der Anwendung verschiedener Techniken zur Fehlermitigation und Fehlerunterdrückung. VQE hat viele Anwendungen (wie Quantenchemie), bei denen er weniger skalierbar ist. Es sind jedoch eine Reihe von Eigenwert-bezogenen Ansätzen entstanden, die VQE ergänzen und erweitern, darunter Krylov-Unterraumdiagonalisierung und stichprobenbasierte Quantendiagonalisierung (SQD). Das Verständnis von VQE ist ein wichtiger erster Schritt zum Verständnis der breiten Palette hybrider klassisch-quantenmechanischer Algorithmen, die entstanden sind.

Dieses Modul beschreibt die grundlegenden Konzepte und die Implementierung von VQE und QAOA. Weitere Tutorials werden fortgeschrittene Themen und Techniken zur Skalierung dieser Algorithmen behandeln. Du benötigst die folgende Bibliothek in deiner Umgebung, um dieses Notebook auszuführen. Falls du sie noch nicht installiert hast, kannst du sie installieren, indem du die folgende Zelle auskommentierst und ausführst.

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "rustworkx": "rustworkx", "scipy": "scipy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

# % pip install 'qiskit[visualization]' qiskit-ibm-runtime

2. Berechnung des minimalen Eigenwerts eines einfachen Hamiltonoperators

Wir beginnen damit, VQE auf einen sehr einfachen Fall anzuwenden, um zu sehen, wie es funktioniert. Wir berechnen den minimalen Eigenwert der Pauli-$Z$-Matrix mit VQE. Wir beginnen mit dem Import einiger allgemeiner Pakete.

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

Wir definieren nun den Operator von Interesse und betrachten ihn in Matrixform.

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

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

Es ist einfach, die Eigenwerte klassisch zu berechnen, sodass wir unsere Arbeit überprüfen können. Dies könnte schwieriger werden, wenn wir Richtung Utility skalieren. Hier verwenden wir numpy.

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

Eigenvalues: [-1.  1.]

Um Eigenwerte mit einem variationellen Quantenalgorithmus zu erhalten, konstruieren wir einen Circuit mit Gates, die variationelle Parameter verwenden:

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

Wenn wir den Erwartungswert eines Operators (wie $Z$) schätzen möchten, sollten wir Estimator verwenden. Wenn wir die Zustände des Systems betrachten möchten, verwenden wir Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Wir können die Anzahl der Bitstrings 0 und 1 mit zufälligen Parameterwerten [1, 2, 3] mit dem Sampler berechnen.

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

Wir wissen, dass wir den Erwartungswert von Z durch $\langle Z \rangle = p_0 - p_1$ mit den Wahrscheinlichkeiten $\{0: p_0, 1: p_1\}$ berechnen können.

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

0.529296875

Dieser Circuit hat funktioniert, aber die gewählten Parameterwerte entsprachen keinem Zustand mit sehr niedriger Energie (oder niedrigem Eigenwert). Der erhaltene Eigenwert ist deutlich höher als das Minimum. Das Ergebnis ist ähnlich bei Verwendung des Estimators.

Beachte, dass Estimator Quantenschaltkreise ohne Messungen entgegennimmt.

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

array(0.54030231)

Wir müssen den Parameterraum durchsuchen und die Werte finden, die den niedrigsten Eigenwert liefern. Wir erstellen eine Funktion, die Parameterwerte der variationellen Form empfängt und den Erwartungswert $\langle Z \rangle$ zurückgibt.

# 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

Wenden wir SciPys minimize-Funktion an, um den minimalen Eigenwert von Z zu finden.

# 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 Übung

Berechne den minimalen Eigenwert von $Z \otimes Z$ mit 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

Lösungen der Übung

Wir definieren den Operator von Interesse und betrachten ihn in Matrixform.

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

Um Eigenwerte mit einem variationellen Quantenalgorithmus zu erhalten, konstruieren wir einen Circuit mit Gates, die variationelle Parameter verwenden:

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

Wenn wir den Erwartungswert eines Operators (wie $Z \otimes Z$) schätzen möchten, würden wir Estimator verwenden. Wenn wir die Zustände des Systems betrachten möchten, verwenden wir 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

Dieser Circuit hat funktioniert, aber die gewählten Parameterwerte entsprachen keinem Zustand mit sehr niedriger Energie (oder niedrigem Eigenwert). Der erhaltene Eigenwert ist deutlich höher als das Minimum. Das Ergebnis ist ähnlich bei Verwendung des Estimators.

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

Wir müssen den Parameterraum durchsuchen und die Werte finden, die den niedrigsten Eigenwert liefern.

# 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

Wir haben einen Eigenwert erhalten, der dem von numpy gelieferten Minimum extrem nahe kommt.

# 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. Quantenoptimierung mit Qiskit Patterns

In diesem How-to lernen wir Qiskit Patterns und Quanten-Approximationsoptimierung kennen. Ein Qiskit Pattern ist ein intuitiver, wiederholbarer Satz von Schritten zur Implementierung eines Quantencomputing-Workflows: Wir werden die Patterns im Kontext der kombinatorischen Optimierung anwenden und zeigen, wie das max-cut-Problem mit dem Quantum Approximate Optimization Algorithm (QAOA) gelöst wird, einer hybriden (quantenklassischen) iterativen Methode.

Beachte, dass dieser QAOA-Teil auf "Part 1: Small-scale QAOA" des Tutorials Quantum approximate optimization algorithm basiert. Sieh dir das Tutorial an, um zu erfahren, wie man es hochskaliert.

3.1 (Kleinskaliges) Qiskit Pattern für Optimierung

Dieser Teil verwendet ein kleinskaliges max-cut-Problem, um die Schritte zu veranschaulichen, die zur Lösung eines Optimierungsproblems mit einem Quantencomputer erforderlich sind. Das max-cut-Problem ist ein Optimierungsproblem, das schwer zu lösen ist (genauer gesagt ein NP-hartes Problem), mit einer Reihe verschiedener Anwendungen in Clustering, Netzwerkwissenschaft und statistischer Physik. Dieses Tutorial betrachtet einen Graphen von Knoten, die durch Kanten verbunden sind, und zielt darauf ab, die Knoten in zwei Mengen aufzuteilen, indem Kanten "geschnitten" werden, sodass die Anzahl der geschnittenen Kanten maximiert wird. Um etwas Kontext zu geben, bevor dieses Problem auf einen Quantenalgorithmus abgebildet wird, kannst du besser verstehen, wie das max-cut-Problem zu einem klassischen kombinatorischen Optimierungsproblem wird, indem du zunächst die Minimierung einer Funktion $f(x)$ betrachtest

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

wobei die Eingabe $x$ ein Vektor ist, dessen Komponenten jeweils einem Knoten eines Graphen entsprechen. Anschließend wird jede dieser Komponenten so eingeschränkt, dass sie entweder $0$ oder $1$ ist (was bedeutet, im Schnitt enthalten oder nicht enthalten zu sein). Dieses Beispiel im kleinen Maßstab verwendet einen Graphen mit $n=5$ Knoten.

Du könntest eine Funktion eines Knotenpaars $i,j$ schreiben, die angibt, ob die entsprechende Kante $(i,j)$ im Schnitt liegt. Zum Beispiel ist die Funktion $x_i + x_j - 2 x_i x_j$ nur dann 1, wenn entweder $x_i$ oder $x_j$ gleich 1 ist (was bedeutet, dass die Kante im Schnitt liegt), und sonst null. Das Problem der Maximierung der Kanten im Schnitt kann formuliert werden als

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

was als Minimierung der Form umgeschrieben werden kann

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

Das Minimum von $f(x)$ liegt in diesem Fall vor, wenn die Anzahl der vom Schnitt durchquerten Kanten maximal ist. Wie du sehen kannst, gibt es bisher nichts, was mit Quantencomputing zu tun hat. Du musst dieses Problem in etwas umformulieren, das ein Quantencomputer verstehen kann. Initialisiere dein Problem, indem du einen Graphen mit $n=5$ Knoten erstellst.

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 Schritt 1. Klassische Eingaben auf ein Quantenproblem abbilden

Der erste Schritt des Patterns besteht darin, das klassische Problem (Graph) auf Quanten-Circuits und Operatoren abzubilden. Dazu sind drei Hauptschritte erforderlich:

1. Nutze eine Reihe mathematischer Umformulierungen, um dieses Problem in der Notation der Quadratischen Unbeschränkten Binären Optimierung (QUBO) darzustellen.

2. Schreibe das Optimierungsproblem als Hamiltonoperator um, dessen Grundzustand der Lösung entspricht, die die Kostenfunktion minimiert.

3. Erstelle einen Quantenschaltkreis, der den Grundzustand dieses Hamiltonoperators durch einen Prozess ähnlich dem Quanten-Annealing präpariert.

Hinweis: Bei der QAOA-Methodik möchtest du letztendlich einen Operator (Hamiltonoperator) haben, der die Kostenfunktion unseres hybriden Algorithmus darstellt, sowie einen parametrisierten Circuit (Ansatz), der Quantenzustände mit Kandidatenlösungen für das Problem darstellt. Du kannst aus diesen Kandidatenzuständen sampeln und sie dann mit der Kostenfunktion auswerten.

Graph → Optimierungsproblem

Der erste Schritt der Abbildung ist eine Notationsänderung. Das Folgende drückt das Problem in QUBO-Notation aus:

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

wobei $Q$ eine $n\times n$-Matrix reeller Zahlen ist, $n$ der Anzahl der Knoten in deinem Graphen entspricht, $x$ der Vektor der binären Variablen ist, die oben eingeführt wurden, und $x^T$ die Transponierte des Vektors $x$ angibt.

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

Optimierungsproblem → Hamiltonoperator

Du kannst dann das QUBO-Problem als einen Hamiltonoperator (hier eine Matrix, die die Energie eines Systems darstellt) umformulieren:

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

Umformulierungsschritte vom QAOA-Problem zum Hamiltonoperator

Um zu demonstrieren, wie das QAOA-Problem auf diese Weise umgeschrieben werden kann, ersetze zunächst die binären Variablen $x_i$ durch einen neuen Satz von Variablen $z_i\in\{-1, 1\}$ mittels

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

Hier siehst du: Wenn $x_i$ gleich $0$ ist, muss $z_i$ gleich $1$ sein. Wenn die $x_i$ durch die $z_i$ im Optimierungsproblem ($x^TQx$) substituiert werden, erhält man eine äquivalente Formulierung.

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

Wenn wir nun $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$ definieren, den Vorfaktor und den konstanten Term $n^2$ entfernen, gelangen wir zu den beiden äquivalenten Formulierungen desselben Optimierungsproblems.

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

Hier hängt $b$ von $Q$ ab. Beachte, dass wir zur Gewinnung von $z^TQz + b^Tz$ den Faktor 1/4 und einen konstanten Offset von $n^2$ weggelassen haben, die bei der Optimierung keine Rolle spielen.

Um nun eine Quantenformulierung des Problems zu erhalten, hebe die $z_i$-Variablen zu einer Pauli-$Z$-Matrix an, also einer $2\times 2$-Matrix der Form

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

Wenn du diese Matrizen in das obige Optimierungsproblem einsetzt, erhältst du den folgenden Hamiltonoperator

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

Erinnere dich auch daran, dass die $Z$-Matrizen in den Berechnungsraum des Quantencomputers eingebettet sind, also in einen Hilbertraum der Größe $2^n\times 2^n$. Daher solltest du Terme wie $Z_iZ_j$ als das Tensorprodukt $Z_i\otimes Z_j$ verstehen, eingebettet in den $2^n\times 2^n$-Hilbertraum. Zum Beispiel wird bei einem Problem mit fünf Entscheidungsvariablen der Term $Z_1Z_3$ als $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ verstanden, wobei $I$ die $2\times 2$-Einheitsmatrix ist.

Dieser Hamiltonoperator wird als Kostenfunktions-Hamiltonoperator bezeichnet. Er hat die Eigenschaft, dass sein Grundzustand der Lösung entspricht, die die Kostenfunktion $f(x)$ minimiert. Um dein Optimierungsproblem zu lösen, musst du daher den Grundzustand von $H_C$ (oder einen Zustand mit hoher Überlappung damit) auf dem Quantencomputer vorbereiten. Das Sampeln aus diesem Zustand liefert dann mit hoher Wahrscheinlichkeit die Lösung für $\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

Hamiltonoperator → Quantenschaltkreis

Der Hamiltonoperator $H_C$ enthält die Quantendefinition deines Problems. Jetzt kannst du einen Quantenschaltkreis erstellen, der beim Sampeln guter Lösungen vom Quantencomputer hilft. QAOA ist vom Quanten-Annealing inspiriert und wendet abwechselnde Schichten von Operatoren im Quantenschaltkreis an.

Die allgemeine Idee ist, im Grundzustand eines bekannten Systems zu starten, $H^{\otimes n}|0\rangle$ oben, und dann das System in den Grundzustand des Kostenoperators zu steuern, an dem du interessiert bist. Dies geschieht durch Anwenden der Operatoren $\exp\{-i\gamma_k H_C\}$ und $\exp\{-i\beta_k H_m\}$ mit den Winkeln $\gamma_1,...,\gamma_p$ und $\beta_1,...,\beta_p~$.

Der erzeugte Quantenschaltkreis ist parametrisiert durch $\gamma_i$ und $\beta_i$, sodass du verschiedene Werte von $\gamma_i$ und $\beta_i$ ausprobieren und aus dem resultierenden Zustand sampeln kannst. In diesem Fall probieren wir ein Beispiel mit 1 QAOA-Schicht, die zwei Parameter enthält: $\gamma_1$ und $\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 Schritt 2. Circuits für die Ausführung auf Quantenhardware optimieren

Der obige Circuit enthält eine Reihe von Abstraktionen, die nützlich sind, um über Quantenalgorithmen nachzudenken, aber nicht direkt auf der Hardware ausgeführt werden können. Um auf einer QPU laufen zu können, muss der Circuit eine Reihe von Operationen durchlaufen, die den Transpilations- oder **Circuit-Optimierungs-**Schritt des Patterns ausmachen.

Die Qiskit-Bibliothek bietet eine Reihe von Transpilations-Passes, die eine breite Palette von Circuit-Transformationen abdecken. Du musst sicherstellen, dass dein Circuit für deinen Zweck optimiert ist.

Transpilation kann mehrere Schritte umfassen, wie zum Beispiel:

• Initiales Mapping der Qubits im Circuit (wie Entscheidungsvariablen) auf physische Qubits auf dem Gerät.

• Unrolling der Anweisungen im Quantenschaltkreis in die hardware-nativen Anweisungen, die das Backend versteht.

• Routing aller Qubits im Circuit, die interagieren, auf physische Qubits, die zueinander benachbart sind.

• Fehlerunterdrückung durch Hinzufügen von Ein-Qubit-Gates zur Rauschunterdrückung mit dynamischer Entkopplung.

Weitere Informationen zur Transpilation findest du in unserer Dokumentation.

Der folgende Code transformiert und optimiert den abstrakten Circuit in ein Format, das für die Ausführung auf einem der über die Cloud zugänglichen Geräte bereit ist, unter Verwendung des Qiskit IBM® Runtime Service.

Beachte, dass du deine Programme lokal im "lokalen Testmodus" testen kannst, bevor du sie an echte Quantencomputer sendest. Weitere Informationen zum lokalen Testmodus findest du in der Dokumentation.

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 Schritt 3. Ausführung mit IBM Quantum Primitives

Im QAOA-Workflow werden die optimalen QAOA-Parameter in einer iterativen Optimierungsschleife gefunden, die eine Reihe von Circuit-Auswertungen durchführt und einen klassischen Optimierer verwendet, um die optimalen $\beta_k$- und $\gamma_k$-Parameter zu finden. Diese Ausführungsschleife wird über die folgenden Schritte ausgeführt:

1. Definiere die Anfangsparameter

2. Instanziiere eine neue Session, die die Optimierungsschleife und die zum Sampeln des Circuits verwendete Primitive enthält

3. Sobald ein optimaler Satz von Parametern gefunden ist, führe den Circuit ein letztes Mal aus, um eine endgültige Verteilung zu erhalten, die im Nachbearbeitungsschritt verwendet wird.

Definiere Circuit mit Anfangsparametern

Wir beginnen mit willkürlich gewählten Parametern.

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

Definiere Backend und Ausführungs-Primitive

Verwende die IBM Quantum Primitives, um mit IBM®-Backends zu interagieren. Die beiden Primitives sind Sampler und Estimator, und die Wahl der Primitive hängt davon ab, welche Art von Messung du auf dem Quantencomputer durchführen möchtest. Für die Minimierung von $H_C$ verwende Estimator, da die Messung der Kostenfunktion einfach der Erwartungswert von $\langle H_C \rangle$ ist.

Ausführung

Die Primitives bieten verschiedene Ausführungsmodi zur Planung von Workloads auf Quantengeräten, und ein QAOA-Workflow läuft iterativ in einer Session. Du kannst die sampler-basierte Kostenfunktion in die SciPy-Minimierungsroutine einfügen, um die optimalen Parameter zu finden.

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

Der Optimierer konnte die Kosten reduzieren und bessere Parameter für den Circuit finden.

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

Sobald du die optimalen Parameter für den Circuit gefunden hast, kannst du diese Parameter zuweisen und die endgültige Verteilung sampeln, die mit den optimierten Parametern erhalten wurde. Hier sollte die Sampler-Primitive verwendet werden, da es die Wahrscheinlichkeitsverteilung der Bitstring-Messungen ist, die dem optimalen Schnitt des Graphen entspricht.

Hinweis: Das bedeutet, einen Quantenzustand $\psi$ im Computer zu präparieren und ihn dann zu messen. Eine Messung lässt den Zustand in einen einzelnen Zustand der Rechenbasis kollabieren - zum Beispiel 010101110000... - der einer Kandidatenlösung $x$ für unser ursprüngliches Optimierungsproblem ($\max f(x)$ oder $\min f(x)$ je nach Aufgabe) entspricht.

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 Schritt 4. Nachbearbeitung, Ergebnis im klassischen Format zurückgeben

Der Nachbearbeitungsschritt interpretiert die Sampling-Ausgabe, um eine Lösung für dein ursprüngliches Problem zurückzugeben. In diesem Fall interessiert dich der Bitstring mit der höchsten Wahrscheinlichkeit, da dieser den optimalen Schnitt bestimmt. Die Symmetrien im Problem erlauben vier mögliche Lösungen, und der Sampling-Prozess gibt eine davon mit einer etwas höheren Wahrscheinlichkeit zurück, aber du kannst in der unten dargestellten Verteilung sehen, dass vier der Bitstrings deutlich wahrscheinlicher sind als der Rest.

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

Besten Schnitt visualisieren

Aus dem optimalen Bitstring kannst du dann diesen Schnitt auf dem ursprünglichen Graphen visualisieren.

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)

Und berechne den Wert des Schnitts. Die Lösung ist aufgrund von Rauschen nicht optimal (der Schnittwert der optimalen Lösung beträgt 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

Damit ist das kleinskalige QAOA-Tutorial abgeschlossen. Du wirst lernen, wie du QAOA im Utility-Scale anpasst, in "Part 2: scale it up!" des Tutorials Quantum approximate optimization algorithm.

# Check Qiskit version
import qiskit

qiskit.__version__

'2.0.2'