Funciones de costo

En esta lección aprenderemos a evaluar una función de costo:

• Primero, conoceremos las primitivas de Qiskit Runtime
• Definiremos una función de costo $C(\vec\theta)$. Esta es una función específica del problema que define el objetivo que el optimizador debe minimizar (o maximizar)
• Definiremos una estrategia de medición con las primitivas de Qiskit Runtime para equilibrar velocidad y precisión

 

Primitivas

Todos los sistemas físicos, ya sean clásicos o cuánticos, pueden existir en diferentes estados. Por ejemplo, un automóvil en una carretera puede tener cierta masa, posición, velocidad o aceleración que caracterizan su estado. De manera similar, los sistemas cuánticos también pueden tener diferentes configuraciones o estados, pero se diferencian de los sistemas clásicos en la forma en que tratamos las mediciones y la evolución del estado. Esto da lugar a propiedades únicas como la superposición y el entrelazamiento, exclusivas de la mecánica cuántica. Así como podemos describir el estado de un automóvil usando propiedades físicas como la velocidad o la aceleración, también podemos describir el estado de un sistema cuántico usando observables, que son objetos matemáticos.

En mecánica cuántica, los estados se representan mediante vectores columna complejos normalizados, o kets ($|\psi\rangle$), y los observables son operadores lineales hermitianos ($\hat{H}=\hat{H}^{\dagger}$) que actúan sobre los kets. Un vector propio ($|\lambda\rangle$) de un observable se conoce como eigenestado. Medir un observable para uno de sus eigenestados ($|\lambda\rangle$) nos dará el eigenvalor correspondiente ($\lambda$) como resultado.

Si te preguntas cómo medir un sistema cuántico y qué puedes medir, Qiskit ofrece dos primitivas que pueden ayudarte:

• Sampler: Dado un estado cuántico $|\psi\rangle$, esta primitiva obtiene la probabilidad de cada posible estado de la base computacional.
• Estimator: Dado un observable cuántico $\hat{H}$ y un estado $|\psi\rangle$, esta primitiva calcula el valor esperado de $\hat{H}$.

La primitiva Sampler

La primitiva Sampler calcula la probabilidad de obtener cada posible estado $|k\rangle$ de la base computacional, dado un circuito cuántico que prepara el estado $|\psi\rangle$. Calcula

$$p_k = |\langle k | \psi \rangle|^2 \quad \forall k \in \mathbb{Z}_2^n \equiv \{0,1,\cdots,2^n-1\},$$

donde $n$ es el número de qubits, y $k$ la representación entera de cualquier posible cadena binaria de salida $\{0,1\}^n$ (es decir, enteros en base $2$).

El Sampler de Qiskit Runtime ejecuta el circuito múltiples veces en un dispositivo cuántico, realiza mediciones en cada ejecución y reconstruye la distribución de probabilidad a partir de las cadenas de bits obtenidas. Cuantas más ejecuciones (o shots) realice, más precisos serán los resultados, pero esto requiere más tiempo y recursos cuánticos.

Sin embargo, dado que el número de posibles salidas crece exponencialmente con el número de qubits $n$ (es decir, $2^n$), el número de shots también deberá crecer exponencialmente para capturar una distribución de probabilidad densa. Por lo tanto, Sampler solo es eficiente para distribuciones de probabilidad dispersas; donde el estado objetivo $|\psi\rangle$ debe poder expresarse como combinación lineal de los estados de la base computacional, con el número de términos creciendo a lo sumo polinomialmente con el número de qubits:

$$|\psi\rangle = \sum^{\text{Poly}(n)}_k w_k |k\rangle.$$

El Sampler también se puede configurar para recuperar probabilidades de una subsección del circuito, lo que representa un subconjunto del total de estados posibles.

La primitiva Estimator

La primitiva Estimator calcula el valor esperado de un observable $\hat{H}$ para un estado cuántico $|\psi\rangle$; donde las probabilidades del observable pueden expresarse como $p_\lambda = |\langle\lambda|\psi\rangle|^2$, siendo $|\lambda\rangle$ los eigenestados del observable $\hat{H}$. El valor esperado se define entonces como el promedio de todos los posibles resultados $\lambda$ (es decir, los eigenvalores del observable) de una medición del estado $|\psi\rangle$, ponderado por las probabilidades correspondientes:

$$\langle\hat{H}\rangle_\psi := \sum_\lambda p_\lambda \lambda = \langle \psi | \hat{H} | \psi \rangle$$

Sin embargo, calcular el valor esperado de un observable no siempre es posible, ya que a menudo no conocemos su eigenbasis. El Estimator de Qiskit Runtime utiliza un proceso algebraico complejo para estimar el valor esperado en un dispositivo cuántico real, descomponiendo el observable en una combinación de otros observables cuya eigenbasis sí conocemos.

En términos más simples, Estimator descompone cualquier observable que no sabe cómo medir en observables más simples y medibles llamados operadores de Pauli. Cualquier operador puede expresarse como combinación de $4^n$ operadores de Pauli.

$$\hat{P}_k := \sigma_{k_{n-1}}\otimes \cdots \otimes \sigma_{k_0} \quad \forall k \in \mathbb{Z}_4^n \equiv \{0,1,\cdots,4^n-1\}, \\$$

tal que

$$\hat{H} = \sum^{4^n-1}_{k=0} w_k \hat{P}_k$$

donde $n$ es el número de qubits, $k \equiv k_{n-1} \cdots k_0$ para $k_l \in \mathbb{Z}_4 \equiv \{0, 1, 2, 3\}$ (es decir, enteros en base $4$), y $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$.

Tras realizar esta descomposición, Estimator deriva un nuevo circuito $V_k|\psi\rangle$ para cada observable $\hat{P}_k$ (a partir del circuito original), con el fin de diagonalizar efectivamente el observable de Pauli en la base computacional y medirlo. Podemos medir fácilmente los observables de Pauli porque conocemos $V_k$ de antemano, lo cual no ocurre en general con otros observables.

Para cada $\hat{P}_{k}$, el Estimator ejecuta el circuito correspondiente en un dispositivo cuántico múltiples veces, mide el estado de salida en la base computacional y calcula la probabilidad $p_{kj}$ de obtener cada posible salida $j$. Luego busca el eigenvalor $\lambda_{kj}$ de $P_k$ correspondiente a cada salida $j$, multiplica por $w_k$ y suma todos los resultados para obtener el valor esperado del observable $\hat{H}$ para el estado dado $|\psi\rangle$.

$$\langle\hat{H}\rangle_\psi = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj},$$

Dado que calcular el valor esperado de $4^n$ Paulis es poco práctico (es decir, crece exponencialmente), Estimator solo puede ser eficiente cuando una gran cantidad de $w_k$ son cero (es decir, descomposición de Pauli dispersa en lugar de densa). Formalmente decimos que, para que este cálculo sea eficientemente resoluble, el número de términos no nulos debe crecer a lo sumo polinomialmente con el número de qubits $n$: $\hat{H} = \sum^{\text{Poly}(n)}_k w_k \hat{P}_k.$

El lector puede notar la suposición implícita de que el muestreo de probabilidades también debe ser eficiente, como se explicó para Sampler, lo que significa

$$\langle\hat{H}\rangle_\psi = \sum_{k}^{\text{Poly}(n)} w_k \sum_{j}^{\text{Poly}(n)}p_{kj} \lambda_{kj}.$$

Ejemplo guiado para calcular valores esperados

Supongamos el estado de un qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$, y el observable

$$\begin{aligned} \hat{H} & = \begin{pmatrix} -1 & 2 \\ 2 & 1 \\ \end{pmatrix}\\[1mm] & = 2X - Z \end{aligned}$$

con el siguiente valor esperado teórico $\langle\hat{H}\rangle_+ = \langle+|\hat{H}|+\rangle = 2.$

Como no sabemos cómo medir este observable, no podemos calcular su valor esperado directamente y necesitamos reexpresarlo como $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$. Se puede demostrar que esto da el mismo resultado notando que $\langle+|X|+\rangle = 1$ y $\langle+|Z|+\rangle = 0$.

Veamos cómo calcular $\langle X \rangle_+$ y $\langle Z \rangle_+$ directamente. Como $X$ y $Z$ no conmutan (es decir, no comparten la misma eigenbasis), no pueden medirse simultáneamente, por lo que necesitamos los circuitos auxiliares:

# Added by doQumentation — required packages for this notebook
!pip install -q matplotlib numpy qiskit qiskit-aer qiskit-ibm-runtime rustworkx

from qiskit import QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

# The following code will work for any other initial single-qubit state and observable
original_circuit = QuantumCircuit(1)
original_circuit.h(0)

H = SparsePauliOp(["X", "Z"], [2, -1])

aux_circuits = []
for pauli in H.paulis:
    aux_circ = original_circuit.copy()
    aux_circ.barrier()
    if str(pauli) == "X":
        aux_circ.h(0)
    elif str(pauli) == "Y":
        aux_circ.sdg(0)
        aux_circ.h(0)
    else:
        aux_circ.id(0)
    aux_circ.measure_all()
    aux_circuits.append(aux_circ)

original_circuit.draw("mpl")

# Auxiliary circuit for X
aux_circuits[0].draw("mpl")

# Auxiliary circuit for Z
aux_circuits[1].draw("mpl")

Ahora podemos llevar a cabo el cálculo manualmente usando Sampler y verificar los resultados con Estimator:

from qiskit.primitives import StatevectorSampler, StatevectorEstimator
from qiskit.result import QuasiDistribution
import numpy as np

## SAMPLER
shots = 10000
sampler = StatevectorSampler()
job = sampler.run(aux_circuits, shots=shots)

# Run the sampler job and step through results
expvals = []
for index, pauli in enumerate(H.paulis):
    data_pub = job.result()[index].data
    bitstrings = data_pub.meas.get_bitstrings()
    counts = data_pub.meas.get_counts()
    quasi_dist = QuasiDistribution(
        {outcome: freq / shots for outcome, freq in counts.items()}
    )

    # Use the probabilities and known eigenvalues of Pauli operators to estimate the expectation value.
    val = 0

    if str(pauli) == "X":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Y":
        val += -1 * quasi_dist.get(1, 0)
        val += 1 * quasi_dist.get(0, 0)

    if str(pauli) == "Z":
        val += 1 * quasi_dist.get(0, 0)
        val += -1 * quasi_dist.get(1, 0)

    expvals.append(val)

# Print expectation values

print("Sampler results:")
for pauli, expval in zip(H.paulis, expvals):
    print(f"  >> Expected value of {str(pauli)}: {expval:.5f}")

total_expval = np.sum(H.coeffs * expvals).real
print(f"  >> Total expected value: {total_expval:.5f}")

# Use estimator for comparison
observables = [
    *H.paulis,
    H,
]  # Note: run for individual Paulis as well as full observable H

estimator = StatevectorEstimator()
job = estimator.run([(original_circuit, observables)])
estimator_expvals = job.result()[0].data.evs

# Print results
print("Estimator results:")
for obs, expval in zip(observables, estimator_expvals):
    if obs is not H:
        print(f"  >> Expected value of {str(obs)}: {expval:.5f}")
    else:
        print(f"  >> Total expected value: {expval:.5f}")

Sampler results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00420
  >> Total expected value: 1.99580
Estimator results:
  >> Expected value of X: 1.00000
  >> Expected value of Z: 0.00000
  >> Total expected value: 2.00000

Rigor matemático (opcional)

Expresando $|\psi\rangle$ con respecto a la base de eigenestados de $\hat{H}$, $|\psi\rangle = \sum_\lambda a_\lambda |\lambda\rangle$, se deduce:

$$\begin{aligned} \langle \psi | \hat{H} | \psi \rangle & = \bigg(\sum_{\lambda'}a^*_{\lambda'} \langle \lambda'|\bigg) \hat{H} \bigg(\sum_{\lambda} a_\lambda | \lambda\rangle\bigg)\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \langle \lambda'|\hat{H}| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \langle \lambda'| \lambda\rangle\\[1mm] & = \sum_{\lambda}\sum_{\lambda'} a^*_{\lambda'}a_{\lambda} \lambda \cdot \delta_{\lambda, \lambda'}\\[1mm] & = \sum_\lambda |a_\lambda|^2 \lambda\\[1mm] & = \sum_\lambda p_\lambda \lambda\\[1mm] \end{aligned}$$

Como no conocemos los eigenvalores ni los eigenestados del observable objetivo $\hat{H}$, primero necesitamos considerar su diagonalización. Dado que $\hat{H}$ es Hermitiano, existe una transformación unitaria $V$ tal que $\hat{H}=V^\dagger \Lambda V,$ donde $\Lambda$ es la matriz diagonal de eigenvalores, de modo que $\langle j | \Lambda | k \rangle = 0$ si $j\neq k$, y $\langle j | \Lambda | j \rangle = \lambda_j$.

Esto implica que el valor esperado puede reescribirse como:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \langle\psi|V^\dagger \Lambda V|\psi\rangle\\[1mm] & = \langle\psi|V^\dagger \bigg(\sum_{j=0}^{2^n-1} |j\rangle \langle j|\bigg) \Lambda \bigg(\sum_{k=0}^{2^n-1} |k\rangle \langle k|\bigg) V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1} \sum_{k=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |k\rangle \langle k| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}\langle\psi|V^\dagger |j\rangle \langle j| \Lambda |j\rangle \langle j| V|\psi\rangle\\[1mm] & = \sum_{j=0}^{2^n-1}|\langle j| V|\psi\rangle|^2 \lambda_j\\[1mm] \end{aligned}$$

Dado que si un sistema está en el estado $|\phi\rangle = V |\psi\rangle$ la probabilidad de medir $| j\rangle$ es $p_j = |\langle j|\phi \rangle|^2$, el valor esperado anterior puede expresarse como:

$$\langle\psi|\hat{H}|\psi\rangle = \sum_{j=0}^{2^n-1} p_j \lambda_j.$$

Es muy importante notar que las probabilidades se toman del estado $V |\psi\rangle$ en lugar de $|\psi\rangle$. Por eso la matriz $V$ es absolutamente necesaria. Quizás te preguntes cómo obtener la matriz $V$ y los eigenvalores $\Lambda$. Si ya tuvieras los eigenvalores, no habría necesidad de usar un computador cuántico, ya que el objetivo de los algoritmos variacionales es precisamente encontrar esos eigenvalores de $\hat{H}$.

Afortunadamente, existe una solución: cualquier matriz $2^n \times 2^n$ puede escribirse como combinación lineal de $4^n$ productos tensoriales de $n$ matrices de Pauli e identidades, todas las cuales son hermitianas y unitarias con $V$ y $\Lambda$ conocidos. Esto es lo que el Estimator de Runtime hace internamente al descomponer cualquier objeto Operator en un SparsePauliOp.

Estos son los operadores que se pueden usar:

$$\begin{array}{c|c|c|c} \text{Operator} & \sigma & V & \Lambda \\[1mm] \hline I & \sigma_0 = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} & V_0 = I & \Lambda_0 = I = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \\[4mm] X & \sigma_1 = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} & V_1 = H =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix} & \Lambda_1 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Y & \sigma_2 = \begin{pmatrix} 0 & -i \\ i & 0 \end{pmatrix} & V_2 = HS^\dagger =\frac{1}{\sqrt{2}} \begin{pmatrix} 1 & 1 \\ 1 & -1 \end{pmatrix}\cdot \begin{pmatrix} 1 & 0 \\ 0 & -i \end{pmatrix} = \frac{1}{\sqrt{2}} \begin{pmatrix} 1 & -i \\ 1 & i \end{pmatrix}\quad & \Lambda_2 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \\[4mm] Z & \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} & V_3 = I & \Lambda_3 = \sigma_3 = \begin{pmatrix} 1 & 0 \\ 0 & -1 \end{pmatrix} \end{array}$$

Así que reescribamos $\hat{H}$ con respecto a los Paulis e identidades:

$$\hat{H} = \sum_{k_{n-1}=0}^3... \sum_{k_0=0}^3 w_{k_{n-1}...k_0} \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0} = \sum_{k=0}^{4^n-1} w_k \hat{P}_k,$$

donde $k = \sum_{l=0}^{n-1} 4^l k_l \equiv k_{n-1}...k_0$ para $k_{n-1},...,k_0\in \{0,1,2,3\}$ (es decir, base $4$), y $\hat{P}_{k} := \sigma_{k_{n-1}}\otimes ... \otimes \sigma_{k_0}$:

$$\begin{aligned} \langle\psi|\hat{H}|\psi\rangle & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}|\langle j| V_k|\psi\rangle|^2 \langle j| \Lambda_k |j\rangle \\[1mm] & = \sum_{k=0}^{4^n-1} w_k \sum_{j=0}^{2^n-1}p_{kj} \lambda_{kj}, \\[1mm] \end{aligned}$$

donde $V_k := V_{k_{n-1}}\otimes ... \otimes V_{k_0}$ y $\Lambda_k := \Lambda_{k_{n-1}}\otimes ... \otimes \Lambda_{k_0}$, tal que: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Funciones de costo

En general, las funciones de costo se usan para describir el objetivo de un problema y qué tan bien está funcionando un estado de prueba con respecto a ese objetivo. Esta definición puede aplicarse a distintos contextos, como química, aprendizaje automático, finanzas, optimización, entre otros.

Consideremos un ejemplo sencillo: encontrar el estado base de un sistema. Nuestro objetivo es minimizar el valor esperado del observable que representa la energía (Hamiltoniano $\hat{\mathcal{H}}$):

$$\min_{\vec\theta} \langle\psi(\vec\theta)|\hat{\mathcal{H}}|\psi(\vec\theta)\rangle$$

Podemos usar el Estimator para evaluar el valor esperado y pasárselo a un optimizador para que lo minimice. Si la optimización tiene éxito, devolverá un conjunto de valores de parámetros óptimos $\vec\theta^*$, a partir de los cuales podremos construir el estado solución propuesto $|\psi(\vec\theta^*)\rangle$ y calcular el valor esperado observado como $C(\vec\theta^*)$.

Nota que solo podremos minimizar la función de costo para el conjunto limitado de estados que estamos considerando. Esto nos lleva a dos posibilidades distintas:

• Nuestro ansatz no define el estado solución en el espacio de búsqueda: En este caso, el optimizador nunca encontrará la solución y tendremos que experimentar con otros ansatzes que puedan representar nuestro espacio de búsqueda con mayor precisión.
• El optimizador es incapaz de encontrar una solución válida: La optimización puede definirse de forma global o local. Exploraremos qué significa esto en la sección posterior.

En definitiva, estaremos ejecutando un ciclo de optimización clásico pero apoyándonos en la evaluación de la función de costo por parte de una computadora cuántica. Desde esta perspectiva, se podría concebir la optimización como una tarea puramente clásica en la que llamamos a algún oráculo cuántico de caja negra cada vez que el optimizador necesita evaluar la función de costo.

def cost_func_vqe(params, circuit, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (circuit, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    return cost

from qiskit.circuit.library import TwoLocal

observable = SparsePauliOp.from_list([("XX", 1), ("YY", -3)])

reference_circuit = QuantumCircuit(2)
reference_circuit.x(0)

variational_form = TwoLocal(
    2,
    rotation_blocks=["rz", "ry"],
    entanglement_blocks="cx",
    entanglement="linear",
    reps=1,
)
ansatz = reference_circuit.compose(variational_form)

theta_list = (2 * np.pi * np.random.rand(1, 8)).tolist()
ansatz.decompose().draw("mpl")

Primero ejecutaremos esto con un simulador: el StatevectorEstimator. Esto es lo recomendable para depurar, pero enseguida repetiremos el cálculo en hardware cuántico real. Cada vez más, los problemas de interés ya no son simulables clásicamente sin instalaciones de supercomputación de vanguardia.

estimator = StatevectorEstimator()
cost = cost_func_vqe(theta_list, ansatz, observable, estimator)
print(cost)

[-0.58744589]

Ahora procederemos a ejecutar en una computadora cuántica real. Ten en cuenta los cambios de sintaxis. Los pasos relacionados con el pass_manager se analizarán con más detalle en el siguiente ejemplo. Un paso de especial importancia en los algoritmos variacionales es el uso de una sesión de Qiskit Runtime. Abrir una sesión te permite ejecutar múltiples iteraciones de un algoritmo variacional sin tener que esperar en una nueva cola cada vez que se actualizan los parámetros. Esto es importante si los tiempos de cola son largos y/o se necesitan muchas iteraciones. Solo los socios de la red IBM Quantum® pueden usar las sesiones de Runtime. Si no tienes acceso a sesiones, puedes reducir el número de iteraciones que envías a la vez y guardar los parámetros más recientes para usarlos en ejecuciones futuras. Si envías demasiadas iteraciones o encuentras tiempos de cola demasiado largos, es posible que obtengas el código de error 1217, que hace referencia a demoras largas entre envíos de trabajos.

# Estimated usage: < 1 min. Benchmarked at 7 seconds on an Eagle processor
# Load necessary packages:

from qiskit_ibm_runtime import (
    QiskitRuntimeService,
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

# Select the least busy backend:

service = QiskitRuntimeService()
backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)
# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_observable = observable.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(theta_list, isa_ansatz, isa_observable, estimator)

session.close()
print(cost)

Ten en cuenta que los valores obtenidos en los dos cálculos anteriores son muy similares. Las técnicas para mejorar los resultados se analizarán más adelante.

Ejemplo de mapeo a sistemas no físicos

El problema de corte máximo (Max-Cut) es un problema de optimización combinatoria que consiste en dividir los vértices de un grafo en dos conjuntos disjuntos de manera que se maximice el número de aristas entre ambos conjuntos. Más formalmente, dado un grafo no dirigido $G=(V,E)$, donde $V$ es el conjunto de vértices y $E$ es el conjunto de aristas, el problema Max-Cut pide dividir los vértices en dos subconjuntos disjuntos, $S$ y $T$, de modo que se maximice el número de aristas con un extremo en $S$ y el otro en $T$.

Podemos aplicar Max-Cut para resolver distintos problemas, como agrupamiento (clustering), diseño de redes, transiciones de fase, entre otros. Empezaremos creando un grafo del problema:

import rustworkx as rx
from rustworkx.visualization import mpl_draw

n = 4
G = rx.PyGraph()
G.add_nodes_from(range(n))
# The edge syntax is (start, end, weight)
edges = [(0, 1, 1.0), (0, 2, 1.0), (0, 3, 1.0), (1, 2, 1.0), (2, 3, 1.0)]
G.add_edges_from(edges)

mpl_draw(
    G, pos=rx.shell_layout(G), with_labels=True, edge_labels=str, node_color="#1192E8"
)

Este problema puede expresarse como un problema de optimización binaria. Para cada nodo $0 \leq i < n$, donde $n$ es el número de nodos del grafo (en este caso $n=4$), consideraremos la variable binaria $x_i$. Esta variable tendrá el valor $1$ si el nodo $i$ pertenece a uno de los grupos que llamaremos $1$, y $0$ si pertenece al otro grupo, que llamaremos $0$. También denotaremos como $w_{ij}$ (elemento $(i,j)$ de la matriz de adyacencia $w$) el peso de la arista que va del nodo $i$ al nodo $j$. Como el grafo es no dirigido, $w_{ij}=w_{ji}$. Entonces podemos formular nuestro problema como la maximización de la siguiente función de costo:

$$\begin{aligned} C(\vec{x}) & =\sum_{i,j=0}^n w_{ij} x_i(1-x_j)\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i,j=0}^n w_{ij} x_ix_j\\[1mm] & = \sum_{i,j=0}^n w_{ij} x_i - \sum_{i=0}^n \sum_{j=0}^i 2w_{ij} x_ix_j \end{aligned}$$

Para resolver este problema con una computadora cuántica, expresaremos la función de costo como el valor esperado de un observable. Sin embargo, los observables que Qiskit admite de forma nativa consisten en operadores de Pauli, que tienen valores propios $1$ y $-1$ en lugar de $0$ y $1$. Por eso realizaremos el siguiente cambio de variable:

Donde $\vec{x}=(x_0,x_1,\cdots ,x_{n-1})$. Podemos usar la matriz de adyacencia $w$ para acceder cómodamente a los pesos de todas las aristas. Esto nos servirá para obtener nuestra función de costo:

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

Esto implica que:

$$\begin{array}{lcl} x_i=0 & \rightarrow & z_i=1 \\ x_i=1 & \rightarrow & z_i=-1.\end{array}$$

Entonces la nueva función de costo que queremos maximizar es:

$$\begin{aligned} C(\vec{z}) & = \sum_{i,j=0}^n w_{ij} \bigg(\frac{1-z_i}{2}\bigg)\bigg(1-\frac{1-z_j}{2}\bigg)\\[1mm] & = \sum_{i,j=0}^n \frac{w_{ij}}{4} - \sum_{i,j=0}^n \frac{w_{ij}}{4} z_iz_j\\[1mm] & = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j \end{aligned}$$

Además, la tendencia natural de una computadora cuántica es encontrar mínimos (generalmente la energía más baja) en lugar de máximos, así que en vez de maximizar $C(\vec{z})$ vamos a minimizar:

$$-C(\vec{z}) = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} z_iz_j - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

Ahora que tenemos una función de costo que minimizar cuyas variables pueden tomar los valores $-1$ y $1$, podemos establecer la siguiente analogía con la Pauli $Z$:

$$z_i \equiv Z_i = \overbrace{I}^{n-1}\otimes ... \otimes \overbrace{Z}^{i} \otimes ... \otimes \overbrace{I}^{0}$$

En otras palabras, la variable $z_i$ será equivalente a una puerta $Z$ actuando sobre el qubit $i$. Además:

$$Z_i|x_{n-1}\cdots x_0\rangle = z_i|x_{n-1}\cdots x_0\rangle \rightarrow \langle x_{n-1}\cdots x_0 |Z_i|x_{n-1}\cdots x_0\rangle = z_i$$

Entonces el observable que vamos a considerar es:

$$\hat{H} = \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2} Z_iZ_j$$

al que tendremos que sumar el término independiente después:

$$\texttt{offset} = - \sum_{i=0}^n \sum_{j=0}^i \frac{w_{ij}}{2}$$

El operador es una combinación lineal de términos con operadores Z en los nodos conectados por una arista (recuerda que el qubit 0 está más a la derecha): $IIZZ + IZIZ + IZZI + ZIIZ + ZZII$. Una vez construido el operador, el ansatz para el algoritmo QAOA puede construirse fácilmente usando el circuito QAOAAnsatz de la biblioteca de circuitos de Qiskit.

from qiskit.circuit.library import QAOAAnsatz
from qiskit.quantum_info import SparsePauliOp

hamiltonian = SparsePauliOp.from_list(
    [("IIZZ", 1), ("IZIZ", 1), ("IZZI", 1), ("ZIIZ", 1), ("ZZII", 1)]
)

ansatz = QAOAAnsatz(hamiltonian, reps=2)
# Draw
ansatz.decompose(reps=3).draw("mpl")

# Sum the weights, and divide by 2

offset = -sum(edge[2] for edge in edges) / 2
print(f"""Offset: {offset}""")

Offset: -2.5

Con el Estimator de Runtime tomando directamente un Hamiltoniano y un ansatz parametrizado, y devolviendo la energía necesaria, la función de costo para una instancia de QAOA es bastante sencilla:

def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

    Parameters:
        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (Estimator): Estimator primitive instance

    Returns:
        float: Energy estimate
    """
    pub = (ansatz, hamiltonian, params)
    cost = estimator.run([pub]).result()[0].data.evs
    #    cost = estimator.run(ansatz, hamiltonian, parameter_values=params).result().values[0]
    return cost

import numpy as np

x0 = 2 * np.pi * np.random.rand(ansatz.num_parameters)

estimator = StatevectorEstimator()
cost = cost_func_vqe(x0, ansatz, hamiltonian, estimator)
print(cost)

1.473098768180865

# Estimated usage: < 1 min, benchmarked at 6 seconds on ibm_osaka, 5-23-24
# Load some necessary packages:

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import Session, EstimatorV2 as Estimator

# Select the least busy backend:

backend = service.least_busy(
    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
)

# Or get a specific backend:
# backend = service.backend("ibm_brisbane")

# Use a pass manager to transpile the circuit and observable for the specific backend being used:

pm = generate_preset_pass_manager(backend=backend, optimization_level=1)
isa_ansatz = pm.run(ansatz)
isa_hamiltonian = hamiltonian.apply_layout(layout=isa_ansatz.layout)

# Set estimator options
estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

# Open a Runtime session:

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)
    cost = cost_func_vqe(x0, isa_ansatz, isa_hamiltonian, estimator)

# Close session after done
session.close()
print(cost)

1.1120776913677988

Volveremos a este ejemplo en la sección de Aplicaciones para explorar cómo usar un optimizador para recorrer el espacio de búsqueda. En términos generales, esto incluye:

• Usar un optimizador para encontrar los parámetros óptimos
• Vincular los parámetros óptimos al ansatz para encontrar los valores propios
• Traducir los valores propios a nuestra definición del problema

Estrategia de medición: velocidad versus precisión

Como mencionamos, estamos usando una computadora cuántica ruidosa como un oráculo de caja negra, donde el ruido puede hacer que los valores recuperados sean no deterministas, lo que genera fluctuaciones aleatorias que, a su vez, dificultarán — o incluso impedirán por completo — la convergencia de ciertos optimizadores hacia una solución propuesta. Este es un problema general que debemos abordar a medida que exploramos progresivamente la utilidad cuántica y avanzamos hacia la ventaja cuántica:

Podemos usar las opciones de supresión y mitigación de errores de los Primitivos de Qiskit Runtime para hacer frente al ruido y maximizar la utilidad de las computadoras cuánticas actuales.

Supresión de errores

La supresión de errores hace referencia a las técnicas utilizadas para optimizar y transformar un circuito durante la compilación con el fin de minimizar los errores. Se trata de una técnica básica de manejo de errores que normalmente introduce cierto sobrecoste de preprocesamiento clásico en el tiempo de ejecución total. Ese sobrecoste incluye transpilar los circuitos para ejecutarlos en hardware cuántico mediante:

• Expresar el circuito usando las puertas nativas disponibles en el sistema cuántico
• Mapear los qubits virtuales a qubits físicos
• Añadir operaciones SWAP según los requisitos de conectividad
• Optimizar las puertas de 1 y 2 qubits
• Añadir desacoplamiento dinámico a los qubits inactivos para prevenir los efectos de la decoherencia.

Los Primitivos permiten usar técnicas de supresión de errores configurando la opción optimization_level y seleccionando opciones avanzadas de transpilación. En un curso posterior, profundizaremos en distintos métodos de construcción de circuitos para mejorar los resultados, pero en la mayoría de los casos recomendamos usar optimization_level=3.

Visualizaremos el valor de aumentar la optimización en el proceso de transpilación examinando un circuito de ejemplo con un comportamiento ideal sencillo.

from qiskit.circuit import Parameter, QuantumCircuit
from qiskit.quantum_info import SparsePauliOp

theta = Parameter("theta")

qc = QuantumCircuit(2)
qc.x(1)
qc.h(0)
qc.cp(theta, 0, 1)
qc.h(0)
observables = SparsePauliOp.from_list([("ZZ", 1)])

qc.draw("mpl")

El circuito anterior puede producir valores esperados sinusoidales del observable dado, siempre que insertemos fases que abarquen un intervalo apropiado, como $[0,2\pi]$.

## Setup phases
import numpy as np

phases = np.linspace(0, 2 * np.pi, 50)

# phases need to be expressed as a list of lists in order to work
individual_phases = [[phase] for phase in phases]

Podemos usar un simulador para demostrar la utilidad de una transpilación optimizada. Más adelante volveremos a usar hardware real para mostrar la utilidad de la mitigación de errores. Usaremos QiskitRuntimeService para obtener un backend real (en este caso, ibm_brisbane) y AerSimulator para simular ese backend, incluyendo su comportamiento de ruido.

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_aer import AerSimulator

# get a real backend from the runtime service
service = QiskitRuntimeService()
backend = service.backend("ibm_brisbane")

# generate a simulator that mimics the real quantum system with the latest calibration results
backend_sim = AerSimulator.from_backend(backend)

Ahora podemos usar un pass manager para transpilar el circuito al "conjunto de instrucciones de arquitectura" o ISA del backend. Este es un nuevo requisito en Qiskit Runtime: todos los circuitos enviados a un backend deben ajustarse a las restricciones del target del backend, lo que significa que deben estar escritos en términos de la ISA del backend, es decir, el conjunto de instrucciones que el dispositivo puede entender y ejecutar. Estas restricciones del target están definidas por factores como las puertas base nativas del dispositivo, la conectividad entre qubits y, cuando corresponde, las especificaciones de pulsos y otros tiempos de instrucción.

Ten en cuenta que en este caso lo haremos dos veces: una con optimization_level = 0 y otra con el nivel fijado en 3. En cada caso usaremos la primitiva Estimator para estimar los valores esperados del observable para distintos valores de fase.

# Import estimator and specify that we are using the simulated backend:

from qiskit_ibm_runtime import EstimatorV2 as Estimator

estimator = Estimator(mode=backend_sim)

circuit = qc

# Use a pass manager to transpile the circuit and observable for the backend being simulated.
# Start with no optimization:

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=0)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

noisy_exp_values = []
pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
noisy_exp_values = cost[0]

# Repeat above steps, but now with optimization = 3:

exp_values_with_opt_es = []
pm = generate_preset_pass_manager(backend=backend_sim, optimization_level=3)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

pub = (isa_circuit, isa_observables, [individual_phases])
cost = estimator.run([pub]).result()[0].data.evs
exp_values_with_opt_es = cost[0]

Por último, podemos graficar los resultados. Vemos que la precisión del cálculo fue bastante buena incluso sin optimización, pero mejoró claramente al subir el nivel de optimización a 3. Ten en cuenta que en circuitos más profundos y complejos, la diferencia entre los niveles de optimización 0 y 3 probablemente será más significativa. Este es un circuito muy simple que sirve como modelo de juguete.

import matplotlib.pyplot as plt

plt.plot(phases, noisy_exp_values, "o", label="opt=0")
plt.plot(phases, exp_values_with_opt_es, "o", label="opt=3")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Mitigación de errores

La mitigación de errores hace referencia a técnicas que permiten reducir los errores de un circuito modelando el ruido del dispositivo en el momento de la ejecución. Por lo general, esto implica una sobrecarga de preprocesamiento cuántico relacionada con el entrenamiento del modelo y una sobrecarga de posprocesamiento clásico para mitigar los errores en los resultados brutos a partir del modelo generado.

La opción resilience_level de las primitivas de Qiskit Runtime especifica el nivel de resiliencia frente a errores que se desea construir. Los niveles más altos generan resultados más precisos a costa de tiempos de procesamiento más largos, debidos a la mayor sobrecarga de muestreo cuántico. Los niveles de resiliencia se pueden usar para configurar el equilibrio entre coste y precisión al aplicar mitigación de errores a tu consulta de primitivas.

Al implementar cualquier técnica de mitigación de errores, esperamos que el sesgo de nuestros resultados se reduzca con respecto al sesgo previo sin mitigar. En algunos casos, el sesgo puede incluso desaparecer. Sin embargo, esto tiene un coste. A medida que reducimos el sesgo en nuestras cantidades estimadas, la variabilidad estadística aumenta (es decir, la varianza), lo cual podemos compensar incrementando aún más el número de shots por circuito en nuestro proceso de muestreo. Esto introduce una sobrecarga adicional más allá de la necesaria para reducir el sesgo, por lo que no se aplica por defecto. Podemos activar fácilmente este comportamiento ajustando el número de shots por circuito en options.executions.shots, como se muestra en el ejemplo a continuación.

En este curso exploraremos estos modelos de mitigación de errores a alto nivel para ilustrar la mitigación que pueden realizar las primitivas de Qiskit Runtime, sin necesidad de conocer todos los detalles de implementación.

Extinción de errores de lectura por twirling (T-REx)

La extinción de errores de lectura por twirling (T-REx) utiliza una técnica conocida como Pauli twirling para reducir el ruido introducido durante el proceso de medición cuántica. Esta técnica no asume ninguna forma específica de ruido, lo que la hace muy general y eficaz.

Flujo de trabajo general:

1. Adquirir datos para el estado cero con bits aleatoriamente invertidos (Pauli X antes de la medición)
2. Adquirir datos para el estado deseado (ruidoso) con bits aleatoriamente invertidos (Pauli X antes de la medición)
3. Calcular la función especial para cada conjunto de datos y dividir.

 

Podemos configurar esto con options.resilience_level = 1, como se muestra en el ejemplo a continuación.

Extrapolación a ruido cero

La extrapolación a ruido cero (ZNE, por sus siglas en inglés) funciona amplificando primero el ruido en el circuito que prepara el estado cuántico deseado, obteniendo mediciones para varios niveles de ruido distintos y usando esas mediciones para inferir el resultado sin ruido.

Flujo de trabajo general:

1. Amplificar el ruido del circuito para varios factores de ruido
2. Ejecutar cada circuito con el ruido amplificado
3. Extrapolar hasta el límite de ruido cero

 

Podemos configurar esto con options.resilience_level = 2. Podemos optimizarlo aún más explorando distintos valores de noise_factors, noise_amplifiers y extrapolators, pero eso está fuera del alcance de este curso. Te animamos a experimentar con estas opciones descritas aquí.

Cada método conlleva su propia sobrecarga asociada: un equilibrio entre el número de cómputos cuánticos necesarios (tiempo) y la precisión de los resultados:

$$\begin{array}{c|c|c|c} \text{Methods} & R=1 \text{, T-REx} & R=2 \text{, ZNE} \\[1mm] \hline \text{Assumptions} & \text{None} & \text{Ability to scale noise} \\[1mm] \text{Qubit overhead} & 1 & 1 \\[1mm] \text{Sampling overhead} & 2 & N_{\text{noise-factors}} \\[1mm] \text{Bias} & 0 & \mathcal{O}(\lambda^{N_{\text{noise-factors}}}) \\[1mm] \end{array}$$

Uso de las opciones de mitigación y supresión de Qiskit Runtime

A continuación se muestra cómo calcular un valor esperado usando mitigación y supresión de errores en Qiskit Runtime. Podemos aprovechar exactamente el mismo circuito y observable que antes, pero esta vez manteniendo el nivel de optimización fijo en 2 y ajustando la resiliencia o las técnicas de mitigación de errores que se utilizan. Este proceso de mitigación de errores ocurre múltiples veces a lo largo del bucle de optimización.

Realizamos esta parte en hardware real, ya que la mitigación de errores no está disponible en simuladores.

# Estimated usage: 8 minutes, benchmarked on an Eagle processor, 5-23-24

from qiskit_ibm_runtime import QiskitRuntimeService
from qiskit_ibm_runtime import (
    Session,
    EstimatorOptions,
    EstimatorV2 as Estimator,
)

# We select the least busy backend

# Select the least busy backend
# backend = service.least_busy(
#    operational=True, min_num_qubits=ansatz.num_qubits, simulator=False
# )

# Or use a specific backend
backend = service.backend("ibm_brisbane")

# Initialize some variables to save the results from different runs:

exp_values_with_em0_es = []
exp_values_with_em1_es = []
exp_values_with_em2_es = []

# Use a pass manager to optimize the circuit and observables for the backend chosen:

pm = generate_preset_pass_manager(backend=backend, optimization_level=2)
isa_circuit = pm.run(circuit)
isa_observables = observables.apply_layout(layout=isa_circuit.layout)

# Open a session and run with no error mitigation:

estimator_options = EstimatorOptions(resilience_level=0, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em0_es = cost[0]

# Open a session and run with resilience = 1:

estimator_options = EstimatorOptions(resilience_level=1, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em1_es = cost[0]

# Open a session and run with resilience = 2:

estimator_options = EstimatorOptions(resilience_level=2, default_shots=10_000)

with Session(backend=backend) as session:
    estimator = Estimator(mode=session, options=estimator_options)

    pub = (isa_circuit, isa_observables, [individual_phases])
    cost = estimator.run([pub]).result()[0].data.evs

session.close()

exp_values_with_em2_es = cost[0]

Como antes, podemos graficar los valores esperados resultantes en función del ángulo de fase para los tres niveles de mitigación de errores utilizados. Con cierta dificultad, se puede apreciar que la mitigación de errores mejora ligeramente los resultados. De nuevo, este efecto es mucho más pronunciado en circuitos más profundos y complejos.

import matplotlib.pyplot as plt

plt.plot(phases, exp_values_with_em0_es, "o", label="unmitigated")
plt.plot(phases, exp_values_with_em1_es, "o", label="resil = 1")
plt.plot(phases, exp_values_with_em2_es, "o", label="resil = 2")
plt.plot(phases, 2 * np.sin(phases / 2) ** 2 - 1, label="ideal")
plt.ylabel("Expectation")
plt.legend()
plt.show()

Resumen

Con esta lección aprendiste a crear una función de coste:

• Crear una función de coste
• Cómo aprovechar las primitivas de Qiskit Runtime para mitigar y suprimir el ruido
• Cómo definir una estrategia de medición para optimizar velocidad versus precisión

Este es nuestro flujo de trabajo variacional de alto nivel:

Nuestra función de coste se ejecuta en cada iteración del bucle de optimización. La próxima lección explorará cómo el optimizador clásico usa la evaluación de nuestra función de coste para seleccionar nuevos parámetros.

import qiskit
import qiskit_ibm_runtime

print(qiskit.version.get_version_info())
print(qiskit_ibm_runtime.version.get_version_info())

1.1.0
0.23.0