Funciones de coste

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

En primer lugar, conoceremos las primitivas Qiskit Runtime
Defina una función de coste $C(\vec\theta)$ . Se trata de una función específica del problema que define el objetivo que el optimizador debe minimizar (o maximizar)
Definición de una estrategia de medición con las primitivas Qiskit Runtime para optimizar la velocidad frente a la precisión

Diagrama que muestra los componentes clave de una función de coste, incluido el uso de primitivas como el estimador y el muestreador.

Primitivos

Todos los sistemas físicos, ya sean clásicos o cuánticos, pueden existir en diferentes estados. Por ejemplo, un coche en una carretera puede tener una masa, posición, velocidad o aceleración determinadas que caracterizan su estado. Del mismo modo, los sistemas cuánticos también pueden tener diferentes configuraciones o estados, pero difieren de los sistemas clásicos en la forma de tratar las mediciones y la evolución de los estados. Esto da lugar a propiedades únicas como la superposición y el entrelazamiento, exclusivas de la mecánica cuántica. Al igual que podemos describir el estado de un coche mediante propiedades físicas como la velocidad o la aceleración, también podemos describir el estado de un sistema cuántico mediante 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 estado propio. La medición de un observable para uno de sus estados propios ( $|\lambda\rangle$ ) nos dará como lectura el valor propio correspondiente ( $\lambda$ ).

Si se pregunta cómo medir un sistema cuántico y qué puede medir, Qiskit ofrece dos sitios web Operaciones básicas y fundamentales o tipo de datos. Qiskit dispone de las primitivas Sampler y Estimator que sirven como bloques de construcción para construir fácilmente cargas de trabajo algorítmicas complejas. que pueden ayudarle:

Sampler: Dado un estado cuántico $|\psi\rangle$ , esta primitiva obtiene la probabilidad de cada posible estado 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$ a partir 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 de base $2$ ).

El sitio Qiskit Runtime Sampler ejecuta el circuito varias veces en un dispositivo cuántico, realizando mediciones en cada ejecución y reconstruyendo la distribución de probabilidad a partir de las cadenas de bits recuperadas. Cuantas más ejecuciones (o disparos ) 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 disparos tendrá que crecer exponencialmente también para capturar una distribución de probabilidad densa. Por lo tanto, Sampler sólo es eficiente para distribuciones de probabilidad dispersas; donde el estado objetivo $|\psi\rangle$ debe ser expresable como una combinación lineal de los estados base computacionales, con el número de términos creciendo como máximo polinomialmente con el número de qubits:

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

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

La primitiva Estimator

La primitiva Estimator calcula el valor de expectativa 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 estados propios del observable $\hat{H}$ . El valor de expectativa se define entonces como la media de todos los resultados posibles $\lambda$ (es decir, los valores propios del observable) de una medición del estado $|\psi\rangle$ , ponderada 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 de expectativa de un observable no siempre es posible, ya que a menudo no conocemos su base propia. El sitio Qiskit Runtime Estimator utiliza un proceso algebraico complejo para estimar el valor de expectativa en un dispositivo cuántico real descomponiendo el observable en una combinación de otros observables cuyas bases propias sí conocemos.

En términos más sencillos, Estimator descompone cualquier observable que no sabe cómo medir en observables más simples y medibles denominados Conjunto de matrices utilizadas habitualmente en informática cuántica para representar y manipular estados cuánticos, formado por la matriz identidad y las tres matrices de Pauli (X, Y y Z)..

Cualquier operador puede expresarse como una 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 de base $4$ ), y $(\sigma_0, \sigma_1, \sigma_2, \sigma_3) := (I, X, Y, Z)$ .

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

Para cada $\hat{P}_{k}$ , el Estimator ejecuta el circuito correspondiente en un dispositivo cuántico varias veces, mide el estado de salida en la base computacional y calcula la probabilidad $p_{kj}$ de obtener cada salida posible $j$ . A continuación, busca el valor propio $\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 de la expectativa de $4^n$ Paulis es poco práctico (es decir, crece exponencialmente), Estimator sólo puede ser eficiente cuando una gran cantidad de $w_k$ son cero (es decir, descomposición Pauli dispersa en lugar de densa ). Formalmente decimos que, para que este cálculo sea eficientemente resoluble, el número de términos distintos de cero tiene que crecer como mucho polinómicamente 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 la probabilidad El proceso de tomar varias medidas de una o varias cosas. también tiene que ser eficiente como se explica para Sampler, lo que significa que

\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 single-qubit $|+\rangle := H|0\rangle = \frac{1}{\sqrt{2}}(|0\rangle + |1\rangle)$ , y observable

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

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

Como no sabemos cómo medir este observable, no podemos calcular su valor de expectativa directamente, y tenemos que volver a expresarlo como $\langle\hat{H}\rangle_+ = 2\langle X \rangle_+ - \langle Z \rangle_+$ . Se puede demostrar que se evalúa con el mismo resultado observando que $\langle+|X|+\rangle = 1$ , y $\langle+|Z|+\rangle = 0$ .

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

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

Output:

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

Output:

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

Output:

Ahora podemos realizar el cálculo manualmente utilizando Sampler y comprobar los resultados en 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}")

Output:

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 estados propios 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}

Dado que no conocemos los valores propios ni los estados propios del observable objetivo $\hat{H}$ , primero tenemos que considerar su diagonalización. Dado que $\hat{H}$ es Una hermitiana es una matriz cuadrada que es igual a su propia transpuesta conjugada, o un operador lineal que es autoadjunto., existe una transformación unitaria $V$ tal que $\hat{H}=V^\dagger \Lambda V,$ donde $\Lambda$ es la matriz diagonal de valores propios, por lo 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}

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

Es muy importante tener en cuenta que las probabilidades se toman del estado $V |\psi\rangle$ en lugar de $|\psi\rangle$ . Por este motivo, la matriz $V$ es absolutamente necesaria.

Quizá se pregunte cómo obtener la matriz $V$ y los valores propios $\Lambda$ . Si ya tuvieras los valores propios, no habría necesidad de utilizar un ordenador cuántico, ya que el objetivo de los algoritmos variacionales es encontrar estos valores propios de $\hat{H}$ .

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

Éstos son los operadores que se pueden utilizar:

\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 vamos a reescribir $\hat{H}$ con respecto a los Paulis y las 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, la 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}$ , tales que: $\hat{P_k}=V_k^\dagger \Lambda_k V_k.$

Funciones de coste

En general, las funciones de coste se utilizan para describir el objetivo de un problema y el rendimiento de un estado de prueba con respecto a ese objetivo. Esta definición puede aplicarse a varios ejemplos de química, aprendizaje automático, finanzas, optimización, etc.

Veamos un ejemplo sencillo para encontrar el estado fundamental de un sistema. Nuestro objetivo es minimizar el valor de la expectativa 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 utilizar el Estimator para evaluar el valor de la expectativa y pasar este valor a un optimizador para minimizar. Si la optimización tiene éxito, devolverá un conjunto de valores óptimos de los parámetros $\vec\theta^*$ , a partir de los cuales podremos construir el estado de solución propuesto $|\psi(\vec\theta^*)\rangle$ y calcular el valor de expectativa observado como $C(\vec\theta^*)$ .

Observe que sólo podremos minimizar la función de coste para el conjunto limitado de estados que estamos considerando. Esto nos lleva a dos posibilidades distintas:

Nuestro ansatz no define el estado de la solución en todo el espacio de búsqueda : Si este es el caso, nuestro optimizador nunca encontrará la solución, y tenemos que experimentar con otros ansatzes que podrían ser capaces de representar nuestro espacio de búsqueda con mayor precisión.
Nuestro optimizador es incapaz de encontrar esta solución válida : La optimización puede ser definida globalmente y definida localmente. Más adelante veremos qué significa esto.

En definitiva, estaremos realizando un bucle de optimización clásico pero confiando la evaluación de la función de coste a un ordenador cuántico. Desde esta perspectiva, se podría pensar en la optimización como un esfuerzo puramente clásico en el que llamamos a algunos Dispositivo hipotético o componente de software que realiza una función específica, pero cuyo funcionamiento interno se desconoce. Un usuario sólo conoce las entradas y salidas de la caja negra, y no sabe cómo ésta procesa las entradas para producir las salidas. cada vez que el optimizador necesita evaluar la función de coste.

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

Output:

Primero lo haremos con un simulador: StatevectorEstimator. Esto suele ser aconsejable para la depuración, pero seguiremos inmediatamente la ejecución de depuración con un 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 última generación.

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

Output:

[-0.58744589]

Ahora procederemos a ejecutarlo en un ordenador cuántico real. Observe los cambios de sintaxis. Los pasos relacionados con el pass_manager se tratará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 Qiskit Runtime. Iniciar una sesión 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 espera son largos y/o se necesitan muchas iteraciones. Sólo los socios de la red IBM Quantum® pueden utilizar las sesiones Runtime. Si no tiene acceso a las sesiones, puede reducir el número de iteraciones que envía en un momento dado y guardar los parámetros más recientes para utilizarlos en futuras ejecuciones. Si envía demasiadas iteraciones o se encuentra con tiempos de espera demasiado largos, puede encontrarse con el código de error 1217, que se refiere a largos retrasos 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)

Obsérvese que los valores obtenidos de 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 del 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 tal manera que se maximice el número de aristas entre ambos conjuntos. En términos más formales, 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 del corte máximo consiste en dividir los vértices en dos subconjuntos disjuntos, $S$ y $T$ , de tal manera que se maximice el número de aristas con un extremo en $S$ y el otro en $T$ .

Podemos aplicar el método «max-cut» para resolver diversos problemas, entre los que se incluyen: la agrupación en clústeres, el diseño de redes, las transiciones de fase, etc. Empezaremos creando un gráfico 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"
)

Output:

Este problema puede expresarse como un problema de optimización binario. 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$ es uno de los grupos que etiquetaremos como $1$ y $0$ si está en el otro grupo, que etiquetaremos como $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$ . Dado que el grafo es no dirigido, $w_{ij}=w_{ji}$ . Entonces podemos formular nuestro problema como la maximización de la siguiente función de coste:

\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 un ordenador cuántico, vamos a expresar la función de coste como el valor esperado de un observable. Sin embargo, los observables que admite Qiskit consisten nativamente en operadores de Pauli, que tienen valores propios $1$ y $-1$ en lugar de $0$ y $1$ . Por eso vamos a hacer el siguiente cambio de variable:

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

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}

Así que la nueva función de coste 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 un ordenador cuántico es encontrar mínimos (normalmente la energía más baja) en lugar de máximos, así que en lugar 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 coste a minimizar cuyas variables pueden tener los valores $-1$ y $1$ , podemos hacer la siguiente analogía con la de 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$ que actúa 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 añadir después el término independiente:

\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 nodos conectados por una arista (recuérdese 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 se puede construir fácilmente utilizando el circuito QAOAAnsatz de la biblioteca de circuitos 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")

Output:

# Sum the weights, and divide by 2

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

Output:

Offset: -2.5

Con el Estimador de Tiempo de Ejecución tomando directamente un Hamiltoniano y un ansatz parametrizado, y devolviendo la energía necesaria, La función de coste para una instancia QAOA es bastante simple:

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)

Output:

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)

Output:

1.1120776913677988

Volveremos sobre este ejemplo en Aplicaciones para explorar cómo aprovechar un optimizador para iterar a través del espacio de búsqueda. En términos generales, esto incluye:

Utilizar un optimizador para encontrar los parámetros óptimos
Vinculación de parámetros óptimos al ansatz para hallar los valores propios
Trasladando los valores propios a la definición de nuestro problema

Estrategia de medición: velocidad frente a precisión

Como se ha mencionado, estamos utilizando un ordenador cuántico ruidoso como oráculo de caja negra, donde el ruido puede hacer que los valores recuperados no sean deterministas, dando lugar a fluctuaciones aleatorias que, a su vez, perjudicarán -o incluso impedirán por completo- la convergencia de ciertos optimizadores a una solución propuesta. Se trata de un problema general que debemos abordar a medida que exploramos progresivamente la utilidad cuántica y avanzamos hacia la ventaja cuántica:

Gráfico que muestra cómo varía el coste de simulación con la complejidad del circuito. Con un ordenador clásico crece exponencialmente. Con la mitigación de errores cuánticos, debería haber un punto de cruce en el que eso se convierta en una ventaja. La corrección cuántica de errores permite un crecimiento lineal del coste de simulación y sin duda supondrá una ventaja.

Podemos utilizar las opciones de supresión y mitigación de errores de Qiskit Runtime Primitive para hacer frente al ruido y maximizar la utilidad de los ordenadores cuánticos actuales.

Supresión de errores

La supresión de errores se refiere 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 gestión de errores que suele suponer una pequeña fase de preprocesamiento Costes adicionales introducidos por las nuevas técnicas, en relación con una aplicación de base. en el tiempo de ejecución total. La sobrecarga implica la transpilación de circuitos para que se ejecuten en hardware cuántico mediante:

Expresar el circuito utilizando las puertas nativas disponibles en un sistema cuántico
Asignación de qubits virtuales a qubits físicos
Añadir SWAPs en función de los requisitos de conectividad
Optimización de las compuertas 1Q y 2Q
Añadir desacoplamiento dinámico a los qubits inactivos para evitar los efectos de la decoherencia.

Las primitivas permiten utilizar 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 diferentes métodos de construcción de circuitos para mejorar los resultados, pero para la mayoría de los casos, recomendamos configurar optimization_level=3.

Visualizaremos el valor de la optimización creciente en el proceso de transpilación observando un circuito de ejemplo con un comportamiento ideal simple.

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

Output:

El circuito anterior puede producir valores de expectativa 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 utilizar un simulador para demostrar la utilidad de una transpilación optimizada. Más adelante volveremos a utilizar hardware real para demostrar la utilidad de la mitigación de errores. Usaremos QiskitRuntimeService para obtener un backend real (en este caso, ibm_brisbane), y usaremos AerSimulator para simular ese backend, incluyendo su comportamiento frente al 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 utilizar un gestor de pases para transpilar el circuito a la "arquitectura de conjunto de instrucciones" 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 objetivo del backend, lo que significa que deben estar escritos en términos del ISA del backend, es decir, el conjunto de instrucciones que el dispositivo puede entender y ejecutar. Estas restricciones vienen definidas por factores como las puertas base nativas del dispositivo, su conectividad qubit y, cuando proceda, sus especificaciones de temporización de impulsos y otras instrucciones.

Tenga en cuenta que en el presente caso, vamos a hacer esto dos veces: una vez con optimization_level = 0, y una vez con él fijado en 3. Cada vez utilizaremos la primitiva Estimator para estimar los valores de la expectativa del observable en diferentes 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 representar gráficamente los resultados, y vemos que la precisión del cálculo era bastante buena incluso sin optimización, pero mejoró definitivamente al aumentar la optimización al nivel 3. Tenga en cuenta que en circuitos más profundos y complicados, es probable que la diferencia entre los niveles de optimización de 0 y 3 sea más significativa. Se trata de un circuito muy sencillo utilizado 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()

Output:

Mitigación de errores

La mitigación de errores se refiere a las técnicas que permiten a los usuarios reducir los errores de circuito mediante la modelización del ruido del dispositivo en el momento de la ejecución. Por lo general, esto conlleva 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 mediante el uso del modelo generado.

La opción resilience_level de la primitiva Qiskit Runtime especifica la cantidad de resiliencia a construir contra errores. Los niveles más altos generan resultados más precisos a costa de tiempos de procesamiento más largos debido a la sobrecarga del muestreo cuántico. Los niveles de resistencia pueden utilizarse para configurar la compensación entre coste y precisión al aplicar la mitigación de errores a su consulta primitiva.

Al aplicar cualquier técnica de mitigación de errores, esperamos que la Desviación sistemática de las magnitudes medidas, generalmente causada por errores. de nuestros resultados se reduzca con respecto al sesgo anterior no mitigado. En algunos casos, el sesgo puede incluso desaparecer. Sin embargo, esto tiene un coste. A medida que reduzcamos el sesgo de nuestras cantidades estimadas, aumentará la variabilidad estadística (es decir, la varianza), que podemos tener en cuenta aumentando aún más el número de disparos por circuito en nuestro proceso de muestreo. Esto introducirá una sobrecarga mayor que la necesaria para reducir el sesgo, por lo que no se hace por defecto. Podemos optar fácilmente por este comportamiento ajustando el número de disparos por circuito en options.executions.shots, como se muestra en el ejemplo siguiente.

Un diagrama que muestra distribuciones más amplias o más estrechas, como en la compensación de sesgo /var iance.

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

Extinción del error de lectura giratoria (T-REx)

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

Flujo de trabajo general:

Adquisición de datos para el estado cero con cambio aleatorio de bits (Pauli X antes de la medición)
Adquirir datos para el estado deseado (ruidoso) con cambios aleatorios de bits (Pauli X antes de la medición)
Calcula la función especial para cada conjunto de datos y divídela.

Esquema de los circuitos de medición y calibración del T-REX.

Podemos configurarlo con options.resilience_level = 1, como se muestra en el siguiente ejemplo.

Extrapolación sin ruido

La extrapolación de ruido cero (ZNE) funciona amplificando primero el ruido en el circuito que está preparando el estado cuántico deseado, obteniendo mediciones para varios niveles diferentes de ruido y utilizando esas mediciones para inferir el resultado sin ruido.

Flujo de trabajo general:

Amplificar el ruido del circuito para varios factores de ruido
Ejecuta todos los circuitos amplificados por ruido
Extrapolar hasta el límite de ruido cero

Un diagrama que muestra los pasos de la ZNE. El ruido se amplifica artificialmente por diferentes factores. A continuación, los valores se extrapolan a lo que deberían ser con ruido cero.

Podemos configurarlo con options.resilience_level = 2. Podemos optimizar esto aún más explorando diversas opciones de noise_factors, noise_amplifiers, y extrapolators, pero esto queda fuera del alcance de este curso. Te animamos a que pruebes estas opciones tal y como se describen aquí.

Cada método lleva asociada su propia sobrecarga: un compromiso entre el número de cálculos cuánticos necesarios (tiempo) y la precisión de nuestros 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 explica cómo calcular un valor de expectativa utilizando la mitigación y supresión de errores en Qiskit Runtime. Podemos utilizar exactamente el mismo circuito y observable que antes, pero esta vez manteniendo el nivel de optimización fijo en el nivel 2, y ahora ajustando la resiliencia o la(s) técnica(s) de mitigación de errores utilizada(s). Este proceso de mitigación de errores se produce varias veces a lo largo de un 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 trazar los valores de expectativa resultantes en función del ángulo de fase para los tres niveles de mitigación de errores utilizados. Con gran dificultad, se puede ver 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 complicados.

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

Output:

Resumen

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

Crear una función de costes
Cómo aprovechar las primitivas Qiskit Runtime para mitigar y suprimir el ruido
Cómo definir una estrategia de medición para optimizar la velocidad frente a la precisión

Aquí está nuestra carga de trabajo variacional de alto nivel:

Un diagrama que muestra el circuito cuántico con unitarios que preparan el estado de referencia y el estado variacional, seguido de mediciones. Se utilizan para evaluar la función de coste.

Nuestra función de costes se ejecuta durante cada iteración del bucle de optimización. La siguiente lección explorará cómo el optimizador clásico utiliza nuestra evaluación de la 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())

Output:

1.1.0
0.23.0

¿Le ha resultado útil esta página?

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