Algoritmos cuánticos: algoritmos cuánticos variacionales

Note

Takashi Imamichi (24 de mayo de 2024)

Descargue el pdf de la conferencia original. Tenga en cuenta que algunos fragmentos de código podrían quedar obsoletos, ya que se trata de imágenes estáticas.

El tiempo aproximado de QPU para ejecutar este experimento es de 9 minutos (probado en un procesador Eagle).

(este cuaderno podría no evaluarse en el tiempo permitido en el Plan Abierto. Por favor, utilice sabiamente los recursos informáticos cuánticos)

1. Introducción

Este tutorial proporciona una visión general de un algoritmo híbrido cuántico-clásico, centrándose específicamente en el eigensolver cuántico variacional (VQE) y el algoritmo de optimización aproximada cuántica (QAOA). El objetivo principal de estos algoritmos es abordar problemas de optimización empleando circuitos cuánticos con puertas cuánticas parametrizadas.

A pesar de los avances en computación cuántica, la presencia de ruido en los dispositivos cuánticos actuales dificulta la extracción de resultados significativos de los circuitos cuánticos profundos. Para superar este reto, VQE y QAOA adoptan un enfoque híbrido cuántico-clásico, que implica la ejecución iterativa de circuitos cuánticos relativamente cortos mediante computación cuántica y la optimización de los parámetros de los circuitos cuánticos parametrizados objetivo mediante computación clásica.

Gracias a la aplicación de diversas técnicas de mitigación y supresión de errores, QAOA puede ofrecer soluciones óptimas a los problemas planteados a escala de servicios públicos. VQE tiene muchas aplicaciones (como la química cuántica) en las que es menos escalable. Pero han surgido una serie de enfoques relacionados con los valores propios para complementar y aumentar la VQE, incluida la diagonalización del subespacio de Krylov y la diagonalización cuántica basada en el muestreo (SQD). Comprender la VQE es un primer paso importante para entender la amplia gama de algoritmos híbridos clásico-cuánticos que han surgido.

Este módulo describe los conceptos fundamentales y la aplicación de VQE y QAOA. Otros tutoriales explorarán temas avanzados y técnicas para ampliar estos algoritmos.

Necesita la siguiente biblioteca en su entorno para ejecutar este cuaderno. Si aún no lo has instalado, puedes instalarlo des-comentando y ejecutando la siguiente celda.

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

2. Cálculo del valor propio mínimo de un hamiltoniano simple

Empezaremos aplicando el VQE a un caso muy sencillo, para ver cómo funciona. Calcularemos el valor propio mínimo de la matriz de Pauli $Z$ con VQE. Empezaremos importando algunos paquetes generales.

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

Ahora definimos el operador de interés y lo vemos en forma de matriz.

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

Output:

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

Es fácil obtener los valores propios clásicamente, por lo que podemos comprobar nuestro trabajo. Esto podría resultar difícil a medida que avanzamos hacia la utilidad. Aquí usamos numpy.

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

Output:

Eigenvalues: [-1.  1.]

Para obtener los valores propios mediante un algoritmo cuántico variacional, construimos un circuito con puertas que toman parámetros variacionales:

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

Output:

Si queremos estimar el valor de la expectativa de un operador (como $Z$ ), debemos utilizar Estimator. Si queremos ver los estados del sistema, utilizamos Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

Podemos calcular recuentos de cadenas de bits 0 y 1 con valores de parámetros aleatorios [1, 2, 3] utilizando Sampler.

# compute counts of bitstrings with random parameter values by Sampler
result = sampler.run([(qc, [1, 2, 3])]).result()
counts = result[0].data.c.get_counts()
counts

Output:

{'0': 783, '1': 241}

Sabemos que podemos calcular el valor de la expectativa de Z mediante $\langle Z \rangle = p_0 - p_1$ con las probabilidades $\{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())

Output:

0.529296875

Este circuito funcionaba, pero los valores de los parámetros elegidos no correspondían a un estado de muy baja energía (o de bajo valor propio). El valor propio obtenido es bastante superior al mínimo. El resultado es similar cuando se utiliza el estimador.

Nótese que el Estimador toma circuitos cuánticos sin mediciones.

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

Output:

array(0.54030231)

Tendremos que buscar entre los parámetros y encontrar los que arrojen el valor propio más bajo. Hacemos una función para recibir los valores de los parámetros de la forma variacional y devolver el valor de la expectativa $\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

Apliquemos la función SciPy's minimize para encontrar el valor propio mínimo de Z.

# minimize the cost function with scipy's minimize
min_result = minimize(cost, [0, 0, 0], method="COBYLA", tol=1e-8)
min_result

Output:

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

Output:

{'0': 1, '1': 1023}

2.1 Ejercicio

Calcule el valor propio mínimo de $Z \otimes Z$ con VQE.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

Output:

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

Soluciones del ejercicio

Definimos el operador de interés y lo vemos en forma de matriz.

z2 = SparsePauliOp("ZZ")
print(z2)
print(z2.to_matrix())

Output:

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

Para obtener los valores propios mediante un algoritmo cuántico variacional, construimos un circuito con puertas que toman parámetros variacionales:

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

Output:

Si queremos estimar el valor de la expectativa de un operador (como $Z \otimes Z$ ), utilizaríamos Estimator. Si queremos ver los estados del sistema, utilizamos 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

Output:

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

Output:

-0.3828125

# 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

Output:

array(-0.35316516)

Tendremos que buscar entre los parámetros y encontrar los que arrojen el valor propio más bajo.

# 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

Output:

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

Obtuvimos un valor propio extremadamente cercano al mínimo que nos daba numpy.

# check counts of bitstrings with the optimal parameters
result = sampler.run([(qc, min_result.x)]).result()
result[0].data.c.get_counts()

Output:

{'01': 1024}

3. Optimización cuántica con patrones Qiskit

En este tutorial aprenderemos sobre los patrones Qiskit y la optimización aproximada cuántica. Un patrón Qiskit es un conjunto intuitivo y repetible de pasos para implementar un flujo de trabajo de computación cuántica:

Aplicaremos estos patrones al ámbito de la optimización combinatoria y mostraremos cómo resolver el problema del corte máximo utilizando el algoritmo de optimización cuántica aproximada (QAOA), un método iterativo híbrido (cuántico-clásico).

Tenga en cuenta que esta parte de QAOA se basa en la "Parte 1: QAOA a pequeña escala" del tutorial Algoritmo de optimización aproximada cuántica. Consulta el tutorial para aprender a escalarlo.

3.1 Patrón Qiskit (a pequeña escala) para la optimización

En esta sección se utilizará un problema de «corte máximo» a pequeña escala para ilustrar los pasos necesarios para resolver un problema de optimización mediante un ordenador cuántico.

El problema del corte máximo es un problema de optimización difícil de resolver (más concretamente, es un problema NP-difícil) que tiene diversas aplicaciones en la agrupación de datos, la ciencia de redes y la física estadística. Este tutorial analiza un grafo formado por nodos conectados por aristas y tiene como objetivo dividir los nodos en dos conjuntos «cortando» aristas, de tal manera que se maximice el número de aristas cortadas.

Para situar el problema en su contexto antes de trasladarlo a un algoritmo cuántico, se puede comprender mejor cómo el problema del corte máximo se convierte en un problema de optimización combinatoria clásica si se parte de la minimización de una función $f(x)$

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

donde la entrada $x$ es un vector cuyos componentes corresponden a cada nodo de un grafo. A continuación, restrinja cada uno de estos componentes para que sean $0$ o $1$ (que representan estar incluido o no en el corte). Este ejemplo a pequeña escala utiliza un grafo con $n=5$ nodos.

Se podría escribir una función de un par de nodos $i,j$ que indique si la arista correspondiente $(i,j)$ está en el corte. Por ejemplo, la función $x_i + x_j - 2 x_i x_j$ es 1 sólo si una de las dos $x_i$ o $x_j$ son 1 (lo que significa que la arista está en el corte) y cero en caso contrario. El problema de maximizar las aristas en el corte puede formularse como

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

que puede reescribirse como una minimización de la forma

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

El mínimo de $f(x)$ en este caso es cuando el número de aristas atravesadas por el corte es máximo. Como puede ver, aún no hay nada relacionado con la computación cuántica. Necesitas reformular este problema en algo que un ordenador cuántico pueda entender.

Inicialice su problema creando un grafo con $n=5$ nodos.

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)

Output:

3.2 Paso 1. Asignar entradas clásicas a un problema cuántico

El primer paso del patrón consiste en mapear el problema clásico (grafo) en circuitos y operadores cuánticos. Para ello, hay que seguir tres pasos principales:

Utilizar una serie de reformulaciones matemáticas, para representar este problema utilizando la notación de problemas de Optimización Cuadrática Binaria No Restringida (QUBO).
Reescribir el problema de optimización como un Hamiltoniano para el que el estado base corresponde a la solución que minimiza la función de coste.
Cree un circuito cuántico que prepare el estado fundamental de este Hamiltoniano mediante un proceso similar al recocido cuántico.

Nota: En la metodología QAOA, en última instancia se desea tener un operador (Hamiltoniano ) que represente la función de coste de nuestro algoritmo híbrido, así como un circuito parametrizado (Ansatz ) que represente estados cuánticos con soluciones candidatas al problema. Puede tomar muestras de estos estados candidatos y luego evaluarlas utilizando la función de coste.

Gráfico → problema de optimización

El primer paso del mapeo es un cambio de notación, A continuación se expresa el problema en notación QUBO:

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

donde $Q$ es una matriz $n\times n$ de números reales, $n$ corresponde al número de nodos de su grafo, $x$ es el vector de variables binarias introducido anteriormente, y $x^T$ indica la transposición del vector $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

Problema de optimización → Hamiltoniano

A continuación, puede reformular el problema QUBO como un Hamiltoniano (en este caso, una matriz que representa la energía de un sistema):

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

Pasos de reformulación del problema QAOA al Hamiltoniano

Para demostrar cómo se puede reescribir el problema QAOA de esta forma, sustituya primero las variables binarias $x_i$ por un nuevo conjunto de variables $z_i\in\{-1, 1\}$ mediante

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

Aquí se puede ver que si $x_i$ es $0$ , entonces $z_i$ debe ser $1$ . Cuando los $x_i$ 's se sustituyen por los $z_i$ 's en el problema de optimización ( $x^TQx$ ), se puede obtener una formulación equivalente.

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

Ahora, si definimos $b_i=-\sum_{j}(Q_{ij}+Q_{ji})$ , eliminamos el prefactor y el término constante $n^2$ , llegamos a las dos formulaciones equivalentes del mismo problema de optimización.

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

Aquí, $b$ depende de $Q$ . Nótese que para obtener $z^TQz + b^Tz$ eliminamos el factor de 1/4 y un desplazamiento constante de $n^2$ que no juegan ningún papel en la optimización.

Ahora, para obtener una formulación cuántica del problema, promueva las variables $z_i$ a una matriz de Pauli $Z$ , tal como una matriz $2\times 2$ de la forma

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

Al sustituir estas matrices en el problema de optimización anterior, se obtiene el siguiente Hamiltoniano

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

Recuerde también que las matrices $Z$ están incrustadas en el espacio computacional del ordenador cuántico, es decir, un espacio de Hilbert de tamaño $2^n\times 2^n$ . Por lo tanto, debe entender términos como $Z_iZ_j$ como el producto tensorial $Z_i\otimes Z_j$ incrustado en el espacio de Hilbert $2^n\times 2^n$ . Por ejemplo, en un problema con cinco variables de decisión, el término $Z_1Z_3$ se entiende como $I\otimes Z_3\otimes I\otimes Z_1\otimes I$ donde $I$ es la matriz identidad $2\times 2$ .

Este Hamiltoniano se llama función de coste Hamiltoniano. Tiene la propiedad de que su estado terreno corresponde a la solución que minimiza la función de coste $f(x)$ . Por lo tanto, para resolver tu problema de optimización ahora necesitas preparar el estado fundamental de $H_C$ (o un estado con un alto solapamiento con él) en el ordenador cuántico. Entonces, el muestreo a partir de este estado arrojará, con una alta probabilidad, la solución a $\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)

Output:

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

Hamiltoniano → circuito cuántico

El Hamiltoniano $H_C$ contiene la definición cuántica de su problema. Ahora puedes crear un circuito cuántico que ayude a muestrear buenas soluciones del ordenador cuántico. El QAOA se inspira en el recocido cuántico y aplica capas alternas de operadores en el circuito cuántico.

La idea general es comenzar en el estado fundamental de un sistema conocido, $H^{\otimes n}|0\rangle$ arriba, y luego dirigir el sistema hacia el estado fundamental del operador de costes que le interesa. Para ello se aplican los operadores $\exp\{-i\gamma_k H_C\}$ y $\exp\{-i\beta_k H_m\}$ con ángulos $\gamma_1,...,\gamma_p$ y $\beta_1,...,\beta_p~$ .

El circuito cuántico que generas está parametrizado por $\gamma_i$ y $\beta_i$ , por lo que puedes probar diferentes valores de $\gamma_i$ y $\beta_i$ y tomar muestras del estado resultante.

En este caso probaremos un ejemplo con 1 capa QAOA que contiene dos parámetros: $\gamma_1$ y $\beta_1$ .

from qiskit.circuit.library import QAOAAnsatz

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

Output:

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

Output:

circuit.parameters

Output:

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

3.3 Paso 2. Optimizar circuitos para la ejecución de hardware cuántico

El circuito anterior contiene una serie de abstracciones útiles para pensar en algoritmos cuánticos, pero no posibles de ejecutar en el hardware. Para poder ejecutarse en una QPU, el circuito debe someterse a una serie de operaciones que conforman el paso de transpilación u optimización del circuito del patrón.

La biblioteca Qiskit ofrece una serie de pases de transpilación que se adaptan a una amplia gama de transformaciones de circuitos. Debes asegurarte de que tu circuito está optimizado para tu propósito.

La transpilación puede implicar varios pasos, como:

Asignación inicial de los qubits en el circuito (como variables de decisión) a qubits físicos en el dispositivo.
Desenrollo de las instrucciones del circuito cuántico a las instrucciones nativas del hardware que entiende el backend.
Enrutamiento de los qubits del circuito que interactúan a qubits físicos adyacentes entre sí.
Supresión de errores añadiendo puertas de un solo qubit para suprimir el ruido con desacoplamiento dinámico.

Encontrará más información sobre la transpilación en nuestra documentación.

El siguiente código transforma y optimiza el circuito abstracto en un formato listo para su ejecución en uno de los dispositivos accesibles a través de la nube mediante el servicio Qiskit IBM® Runtime.

Tenga en cuenta que puede probar sus programas localmente mediante el "modo de prueba local" antes de enviarlos a ordenadores cuánticos reales. Encontrará más información sobre el modo de prueba local en la documentación.

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)

Output:

  service = QiskitRuntimeService(channel="ibm_quantum_platform")

<IBMBackend('ibm_strasbourg')>

3.4 Paso 3. Ejecutar utilizando Qiskit primitives

En el flujo de trabajo QAOA, los parámetros QAOA óptimos se encuentran en un bucle de optimización iterativo, que ejecuta una serie de evaluaciones de circuitos y utiliza un optimizador clásico para encontrar los parámetros $\beta_k$ y $\gamma_k$ óptimos. Este bucle de ejecución se ejecuta mediante los siguientes pasos:

Definir los parámetros iniciales
Instanciar un nuevo Session que contenga el bucle de optimización y la primitiva utilizada para muestrear el circuito
Una vez encontrado un conjunto óptimo de parámetros, ejecute el circuito una última vez para obtener una distribución final que se utilizará en la etapa de postprocesamiento.

Definir circuito con parámetros iniciales

Partimos de parámetros elegidos arbitrariamente.

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

Definir backend y primitiva de ejecución

Utiliza las primitivas de « Qiskit Runtime » para interactuar con los backends de « IBM® ». Las dos primitivas son «Sampler» y «Estimator», y la elección de la primitiva depende del tipo de medición que se desee realizar en el ordenador cuántico. Para minimizar $H_C$ , utiliza el estimador, ya que la medida de la función de coste es simplemente el valor esperado de $\langle H_C \rangle$ .

Ejecute

Las primitivas ofrecen una variedad de modos de ejecución para programar cargas de trabajo en dispositivos cuánticos, y un flujo de trabajo QAOA se ejecuta iterativamente en una sesión.

Puede introducir la función de coste basada en el muestreador en la rutina de minimización SciPy para encontrar los parámetros óptimos.

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)

Output:

 message: Optimization terminated successfully.
 success: True
  status: 1
     fun: -0.6557925874481715
       x: [ 2.873e+00  9.414e-01]
    nfev: 21
   maxcv: 0.0

El optimizador consiguió reducir el coste y encontrar mejores parámetros para el circuito.

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

Output:

Una vez que haya encontrado los parámetros óptimos para el circuito, puede asignar estos parámetros y muestrear la distribución final obtenida con los parámetros optimizados. Aquí es donde debe utilizarse la primitiva Sampler, ya que es la distribución de probabilidad de las medidas de la cadena de bits la que corresponde al corte óptimo del grafo.

Nota: Esto significa preparar un estado cuántico $\psi$ en el ordenador y luego medirlo. Una medición colapsará el estado en un único estado base computacional - por ejemplo, 010101110000... - que corresponde a una solución candidata $x$ a nuestro problema de optimización inicial ( $\max f(x)$ o $\min f(x)$ dependiendo de la tarea).

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

Output:

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)

Output:

{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 Paso 4. Procesamiento posterior, devolución del resultado en formato clásico

El paso de postprocesamiento interpreta el resultado del muestreo para obtener una solución al problema original. En este caso, te interesa la cadena de bits con mayor probabilidad, ya que determina el corte óptimo. Las simetrías del problema permiten cuatro soluciones posibles, y el proceso de muestreo devolverá una de ellas con una probabilidad ligeramente mayor, pero en la distribución que se muestra a continuación se puede ver que cuatro de las cadenas de bits son claramente más probables que el resto.

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

Output:

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

Output:

Visualizar el mejor corte

A partir de la cadena de bits óptima, se puede visualizar este corte en el gráfico original.

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)

Output:

Y calcula el valor del corte. La solución no es óptima debido al ruido (el valor de corte de la solución óptima es 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)

Output:

The value of the cut is: 5

Con esto concluye el tutorial de QAOA a pequeña escala. Aprenderás a adaptar QAOA a escala de servicios públicos en la "Parte 2: ¡aumentar la escala!" del tutorial del algoritmo de optimización aproximada Quantum.

# Check Qiskit version
import qiskit

qiskit.__version__

Output:

'2.0.2'

¿Le ha resultado útil esta página?

Informe de un error, de una errata o solicite contenido en GitHub.