Aller au contenu principal

Algoritmos cuánticos: Algoritmos cuánticos variacionales

remarque

Takashi Imamichi (24 de mayo de 2024)

Descarga el PDF de la conferencia original. Ten en cuenta que algunos fragmentos de código pueden quedar obsoletos, ya que son imágenes estáticas.

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

(Es posible que este notebook no se evalúe en el tiempo permitido del Plan Abierto. Por favor, utiliza los recursos de computación cuántica con responsabilidad.)

1. Introducción​

Este tutorial ofrece una visión general de un algoritmo híbrido cuántico-clásico, con especial énfasis en el solucionador cuántico variacional de valores propios (VQE) y el algoritmo cuántico de optimización aproximada (QAOA). El objetivo principal de estos algoritmos es abordar problemas de optimización mediante el uso de 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 circuitos cuánticos profundos. Para superar este desafío, VQE y QAOA adoptan un enfoque híbrido cuántico-clásico, que consiste en ejecutar de forma iterativa circuitos cuánticos relativamente cortos mediante computación cuántica y optimizar los parámetros de los circuitos cuánticos parametrizados objetivo mediante computación clásica.

QAOA tiene el potencial de proporcionar soluciones óptimas a los problemas objetivo a escala de utilidad, gracias a la aplicación de diversas técnicas de mitigación y supresión de errores. VQE tiene muchas aplicaciones (como la química cuántica) en las que resulta menos escalable. Sin embargo, han surgido varios enfoques relacionados con valores propios que complementan y amplían VQE, incluyendo la diagonalización de subespacios de Krylov y la diagonalización cuántica basada en muestreo (SQD). Comprender VQE es un primer paso fundamental 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 implementación de VQE y QAOA. Tutoriales posteriores explorarán temas avanzados y técnicas para escalar estos algoritmos. Para ejecutar este notebook necesitas la siguiente biblioteca en tu entorno. Si aún no la has instalado, puedes hacerlo descomentando y ejecutando la siguiente celda.

# 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. Cálculo del valor propio mínimo de un Hamiltoniano simple​

Comenzaremos aplicando VQE a un caso muy sencillo para ver cómo funciona. Calcularemos el valor propio mínimo de la matriz de Pauli ZZ con VQE. Para ello, importaremos 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 visualizamos en forma matricial.

op = SparsePauliOp("Z")
op.to_matrix()
array([[ 1.+0.j,  0.+0.j],
       [ 0.+0.j, -1.+0.j]])

Es fácil obtener los valores propios de forma clásica, así que podemos verificar nuestros resultados. Esto puede volverse difícil conforme nos acercamos a la escala de utilidad. Aquí usamos numpy.

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

Para obtener valores propios usando un algoritmo cuántico variacional, construimos un circuito con puertas que aceptan 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 of the previous code cell

Si queremos estimar el valor esperado de un operador (como ZZ), debemos usar Estimator. Si queremos observar los estados del sistema, usamos Sampler.

sampler = StatevectorSampler()
estimator = StatevectorEstimator()

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

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

Sabemos que podemos calcular el valor esperado de Z mediante ⟨Z⟩=p0−p1\langle Z \rangle = p_0 - p_1 con las probabilidades {0:p0,1:p1}\{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

Este circuito funcionó, pero los valores de parámetros elegidos no corresponden a un estado de baja energía (o bajo valor propio). El valor propio obtenido es bastante superior al mínimo. El resultado es similar al usar Estimator.

Ten en cuenta que Estimator trabaja con circuitos cuánticos sin mediciones.

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

Necesitaremos buscar entre los parámetros y encontrar aquellos que produzcan el valor propio más bajo. Creamos una función que recibe los valores de los parámetros de la forma variacional y devuelve el valor esperado ⟨Z⟩\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 minimize de SciPy 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
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 Ejercicio​

Calcula el valor propio mínimo de Z⊗ZZ \otimes Z con 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

Soluciones del ejercicio​

Definimos el operador de interés y lo visualizamos en forma matricial.

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

Para obtener valores propios usando un algoritmo cuántico variacional, construimos un circuito con puertas que aceptan 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 of the previous code cell

Si queremos estimar el valor esperado de un operador (como Z⊗ZZ \otimes Z), usaríamos Estimator. Si queremos observar los estados del sistema, usamos 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

Este circuito funcionó, pero los valores de parámetros elegidos no corresponden a un estado de baja energía (o bajo valor propio). El valor propio obtenido es bastante superior al mínimo. El resultado es similar al usar Estimator.

# verify the expectation value of ZZ with Estimator
result = estimator.run([(qc_estimator, z2, [1, 2, 3, 4, 5, 6])]).result()
result[0].data.evs
array(-0.35316516)

Necesitaremos buscar entre los parámetros y encontrar aquellos que produzcan 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
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 proporciona 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. Optimización cuántica con patrones de Qiskit​

En este tutorial aprenderás sobre los patrones de Qiskit y la optimización aproximada cuántica. Un patrón de Qiskit es un conjunto intuitivo y repetible de pasos para implementar un flujo de trabajo de computación cuántica: "Función de Qiskit" Aplicaremos los patrones al contexto de la optimización combinatoria y mostraremos cómo resolver el problema Max-Cut usando el Algoritmo de Optimización Aproximada Cuántica (QAOA), un método iterativo híbrido (cuántico-clásico).

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

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

Esta parte usará un problema Max-Cut de pequeña escala para ilustrar los pasos necesarios para resolver un problema de optimización usando una computadora cuántica. El problema Max-Cut es un problema de optimización difícil de resolver (más específicamente, es un problema NP-hard) con diversas aplicaciones en agrupamiento, ciencias de redes y física estadística. En este tutorial se considera un grafo de nodos conectados por aristas y se busca particionar los nodos en dos conjuntos "cortando" aristas, de modo que se maximice el número de aristas cortadas. "Maxcut" Para dar algo de contexto antes de mapear este problema a un algoritmo cuántico, puedes entender mejor cómo el problema Max-Cut se convierte en un problema de optimización combinatoria clásica considerando primero la minimización de una función f(x)f(x)

min⁡x∈{0,1}nf(x),\min_{x\in \{0, 1\}^n}f(x),

donde la entrada xx es un vector cuyos componentes corresponden a cada nodo del grafo. Luego, se restringe cada componente a ser 00 o 11 (que representan si el nodo pertenece o no al corte). Este caso de ejemplo a pequeña escala usa un grafo con n=5n=5 nodos.

Puedes escribir una función para un par de nodos i,ji,j que indique si la arista correspondiente (i,j)(i,j) está en el corte. Por ejemplo, la función xi+xj−2xixjx_i + x_j - 2 x_i x_j es 1 solo si uno de los valores xix_i o xjx_j es 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∈{0,1}n∑(i,j)xi+xj−2xixj,\max_{x\in \{0, 1\}^n} \sum_{(i,j)} x_i + x_j - 2 x_i x_j,

lo cual puede reescribirse como una minimización de la forma

min⁡x∈{0,1}n∑(i,j)2xixj−xi−xj.\min_{x\in \{0, 1\}^n} \sum_{(i,j)} 2 x_i x_j - x_i - x_j.

El mínimo de f(x)f(x) en este caso ocurre cuando el número de aristas atravesadas por el corte es máximo. Como puedes ver, aún no hay nada relacionado con la computación cuántica. Es necesario reformular este problema en algo que una computadora cuántica pueda entender. Inicializa el problema creando un grafo con n=5n=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)

Salida de la celda de código anterior

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

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

  1. Utilizar una serie de reformulaciones matemáticas para representar este problema usando la notación de problemas de Optimización Binaria Cuadrática sin Restricciones (QUBO).
  2. Reescribir el problema de optimización como un Hamiltoniano cuyo estado base corresponde a la solución que minimiza la función de costo.
  3. Crear un circuito cuántico que prepare el estado base de este Hamiltoniano mediante un proceso similar al recocido cuántico.

Nota: En la metodología QAOA, el objetivo final es tener un operador (Hamiltoniano) que represente la función de costo del algoritmo híbrido, así como un circuito parametrizado (Ansatz) que represente estados cuánticos con soluciones candidatas al problema. Se puede muestrear a partir de estos estados candidatos y luego evaluarlos usando la función de costo.

Grafo → problema de optimización​

El primer paso del mapeo es un cambio de notación. Lo siguiente expresa el problema en notación QUBO:

min⁡x∈{0,1}nxTQx,\min_{x\in \{0, 1\}^n}x^T Q x,

donde QQ es una matriz n×nn\times n de números reales, nn corresponde al número de nodos del grafo, xx es el vector de variables binarias introducidas anteriormente, y xTx^T indica la transpuesta del vector xx.

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​

Luego puedes reformular el problema QUBO como un Hamiltoniano (en este caso, una matriz que representa la energía de un sistema):

HC=∑ijQijZiZj+∑ibiZi.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 el problema QAOA puede reescribirse de esta manera, primero reemplaza las variables binarias xix_i por un nuevo conjunto de variables zi∈{−1,1}z_i\in\{-1, 1\} mediante

xi=1−zi2.x_i = \frac{1-z_i}{2}.

Aquí puedes ver que si xix_i es 00, entonces ziz_i debe ser 11. Cuando los xix_i se sustituyen por los ziz_i en el problema de optimización (xTQxx^TQx), se puede obtener una formulación equivalente.

xTQx=∑ijQijxixj=14∑ijQij(1−zi)(1−zj)=14∑ijQijzizj−14∑ij(Qij+Qji)zi+n24.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 bi=−∑j(Qij+Qji)b_i=-\sum_{j}(Q_{ij}+Q_{ji}), eliminamos el prefactor y el término constante n2n^2, llegamos a las dos formulaciones equivalentes del mismo problema de optimización.

minx∈{0,1}nxTQx⟺min⁡z∈{−1,1}nzTQz+bTzmin_{x\in\{0,1\}^n} x^TQx\Longleftrightarrow \min_{z\in\{-1,1\}^n}z^TQz + b^Tz

Aquí, bb depende de QQ. Ten en cuenta que para obtener zTQz+bTzz^TQz + b^Tz descartamos el factor 1/4 y un desplazamiento constante de n2n^2, los cuales no juegan ningún papel en la optimización.

Ahora, para obtener una formulación cuántica del problema, se promueven las variables ziz_i a una matriz de Pauli ZZ, como una matriz 2×22\times 2 de la forma

Zi=(100−1).Z_i = \begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}.

Cuando sustituyes estas matrices en el problema de optimización anterior, obtienes el siguiente Hamiltoniano

HC=∑ijQijZiZj+∑ibiZi.H_C=\sum_{ij}Q_{ij}Z_iZ_j + \sum_i b_iZ_i.

Recuerda también que las matrices ZZ están embebidas en el espacio computacional de la computadora cuántica, es decir, un espacio de Hilbert de tamaño 2n×2n2^n\times 2^n. Por lo tanto, debes entender términos como ZiZjZ_iZ_j como el producto tensorial Zi⊗ZjZ_i\otimes Z_j embebido en el espacio de Hilbert 2n×2n2^n\times 2^n. Por ejemplo, en un problema con cinco variables de decisión, el término Z1Z3Z_1Z_3 se entiende como I⊗Z3⊗I⊗Z1⊗II\otimes Z_3\otimes I\otimes Z_1\otimes I donde II es la matriz identidad 2×22\times 2.

Este Hamiltoniano se denomina el Hamiltoniano de la función de costo. Tiene la propiedad de que su estado base corresponde a la solución que minimiza la función de costo f(x)f(x). Por lo tanto, para resolver tu problema de optimización ahora necesitas preparar el estado base de HCH_C (o un estado con alta superposición con él) en la computadora cuántica. Luego, muestrear desde este estado producirá, con alta probabilidad, la solución de min⁡ f(x)\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

Hamiltoniano → circuito cuántico​

El Hamiltoniano HCH_C contiene la definición cuántica de tu problema. Ahora puedes crear un circuito cuántico que ayude a muestrear buenas soluciones desde la computadora cuántica. El QAOA está inspirado en el recocido cuántico y aplica capas alternadas de operadores en el circuito cuántico.

La idea general es comenzar en el estado base de un sistema conocido, H⊗n∣0⟩H^{\otimes n}|0\rangle como se indicó arriba, y luego guiar el sistema hacia el estado base del operador de costo que nos interesa. Esto se hace aplicando los operadores exp⁡{−iγkHC}\exp\{-i\gamma_k H_C\} y exp⁡{−iβkHm}\exp\{-i\beta_k H_m\} con ángulos γ1,...,γp\gamma_1,...,\gamma_p y β1,...,βp \beta_1,...,\beta_p~.

El circuito cuántico que generas está parametrizado por γi\gamma_i y βi\beta_i, de modo que puedes probar diferentes valores de γi\gamma_i y βi\beta_i y muestrear desde el estado resultante. "Diagrama del circuito QAOA" En este caso usaremos un ejemplo con 1 capa QAOA que contiene dos parámetros: γ1\gamma_1 y β1\beta_1.

from qiskit.circuit.library import QAOAAnsatz
circuit = QAOAAnsatz(cost_operator=cost_hamiltonian, reps=1)
circuit.measure_all()
circuit.draw("mpl")

Salida de la celda de código anterior

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

Salida de la celda de código anterior

circuit.parameters
ParameterView([ParameterVectorElement(β[0]), ParameterVectorElement(γ[0])])

3.3 Paso 2. Optimizar circuitos para ejecución en hardware cuántico​

El circuito anterior contiene una serie de abstracciones útiles para razonar sobre algoritmos cuánticos, pero que no pueden ejecutarse directamente en el hardware. Para poder correr en una QPU, el circuito debe pasar por una serie de operaciones que conforman el paso de transpilación u optimización de circuitos del patrón.

La biblioteca Qiskit ofrece una serie de pasadas de transpilación que cubren una amplia gama de transformaciones de circuitos. Debes asegurarte de que tu circuito esté optimizado para tu propósito.

La transpilación puede involucrar varios pasos, como:

  • Mapeo inicial de los qubits del circuito (como las variables de decisión) a qubits físicos del dispositivo.
  • Descomposición de las instrucciones del circuito cuántico en instrucciones nativas del hardware que el backend entiende.
  • Enrutamiento de los qubits del circuito que interactúan hacia qubits físicos adyacentes entre sí.
  • Supresión de errores mediante la adición de compuertas de un qubit para suprimir el ruido con desacoplamiento dinámico.

Más información sobre transpilación está disponible en nuestra documentación.

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

Ten en cuenta que puedes probar tus programas localmente usando el "modo de prueba local" antes de enviarlos a computadoras cuánticas reales. Más información sobre el modo de prueba local está disponible 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)
service = QiskitRuntimeService(channel="ibm_quantum_platform")
<IBMBackend('ibm_strasbourg')>

Salida de la celda de código anterior

3.4 Paso 3. Ejecutar usando las primitivas de Qiskit​

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

  1. Definir los parámetros iniciales
  2. Instanciar una nueva Session que contenga el bucle de optimización y la primitiva usada para muestrear el circuito
  3. Una vez encontrado el conjunto óptimo de parámetros, ejecutar el circuito una última vez para obtener una distribución final que se usará en el paso de post-procesamiento.

Definir el circuito con parámetros iniciales​

Comenzamos con parámetros elegidos de manera arbitraria.

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

Definir el backend y la primitiva de ejecución​

Usa 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 quieras realizar en la computadora cuántica. Para la minimización de HCH_C, usa el Estimator, ya que la medición de la función de costo es simplemente el valor esperado de ⟨HC⟩\langle H_C \rangle.

Ejecutar​

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 de forma iterativa en una sesión. &quot;Modo de ejecución&quot; Puedes conectar la función de costo basada en el Sampler a la rutina de minimización de 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)
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 fue capaz de reducir el costo 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()

Salida de la celda de código anterior

Una vez que has encontrado los parámetros óptimos para el circuito, puedes asignarlos y muestrear la distribución final obtenida con los parámetros optimizados. Aquí es donde se debe usar la primitiva Sampler, ya que es la distribución de probabilidad de las mediciones de cadenas de bits la que corresponde al corte óptimo del grafo.

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

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

Salida de la celda de código anterior

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 Paso 4. Post-procesar y devolver el resultado en formato clásico​

El paso de post-procesamiento interpreta la salida del muestreo para devolver una solución al problema original. En este caso, nos 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 más alta; sin embargo, en la distribución graficada a continuación puedes ver que cuatro de las cadenas de bits son notablemente 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)
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()

Salida de la celda de código anterior

Visualizar el mejor corte​

A partir de la cadena de bits óptima, puedes visualizar este corte sobre el grafo 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)

Salida de la celda de código anterior

Y calcular el valor del corte. La solución no es óptima debido al ruido (el valor del 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)
The value of the cut is: 5

Con esto concluye el tutorial de QAOA a pequeña escala. Aprenderás cómo adaptar el QAOA a escala de utilidad en "Parte 2: ¡escálalo!" del tutorial Algoritmo de optimización aproximada cuántica.

# Check Qiskit version
import qiskit

qiskit.__version__
'2.0.2'