Fórmulas multiproducto para reducir el error de Trotter

Estimación de tiempo de ejecución: cuatro minutos en un procesador Heron r2 (NOTA: Se trata únicamente de una estimación). (El tiempo de ejecución puede variar.)

Resultados del aprendizaje

• Cómo las fórmulas multiproducto (MPF) reducen el error de Trotter en la simulación hamiltoniana mediante la combinación de los valores esperados de múltiples circuitos poco profundos
• Cuándo los MPF resultan más ventajosos que las fórmulas estándar de los productos y cuándo no son la herramienta adecuada
• Cómo calcular los coeficientes MPF estáticos y dinámicos utilizando el qiskit_addon_mpf paquete
• Cómo ejecutar un flujo de trabajo MPF de principio a fin en un hardware d IBM Quantum®, incluyendo la transpilación, la corrección de errores y el posprocesamiento

Requisitos previos

• Métodos de compilación para circuitos de simulación hamiltonianos : introducción a los circuitos de Trotter (fórmula del producto) en Qiskit.
• Fórmulas de productos en Qiskit, en concreto las SuzukiTrotter clases de síntesis y LieTrotter .
• Qiskit primitives y la interfaz «Estimator ».

En segundo plano

¿Qué son las fórmulas multiproducto?

Al simular sistemas cuánticos en un ordenador cuántico, una tarea fundamental consiste en aproximar el operador de evolución temporal $e^{-iHt}$ para un hamiltoniano $H$. El enfoque estándar utiliza fórmulas de producto (PF), también conocidas como descomposiciones de Trotter-Suzuki. Estos descomponen $H = \sum_{a=1}^d F_a$ en términos cuyos operadores unitarios individuales $e^{-iF_a t}$ son eficientes de implementar, y a continuación aproximan la evolución completa como un producto ordenado de estos operadores unitarios más sencillos.

La fórmula del producto de primer orden (Lie-Trotter) es:

lo que da lugar a un error cuadrático: $S_1(t) = e^{-iHt} + \mathcal{O}(t^2)$. Las fórmulas simétricas de orden superior $S_{2\chi}(t)$, donde $\chi$ indica el orden de la fórmula del producto simétrico (véase la ref. [1] ), convergen más rápidamente, como se muestra en $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$, pero a costa de circuitos más complejos por paso.

Para reducir el error en un orden fijo $\chi$, normalmente se divide el tiempo total de evolución $t$ en $k$ pasos de Trotter más pequeños. Cada paso aproxima $e^{-iHt/k}$ mediante una fórmula de producto y los pasos se encadenan:

Para una fórmula simétrica de orden $2\chi$, el error residual de Trotter varía entonces proporcionalmente a $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$. Por lo tanto, al aumentar $k$ se reduce rápidamente el error de Trotter, pero también se aumenta linealmente la profundidad del circuito, lo que, en un hardware con ruido, se traduce en un mayor ruido acumulado en las puertas. Esta tensión entre el error de Trotter (que favorece valores más grandes de $k$ ) y el ruido del hardware (que favorece valores más pequeños de $k$ ) es precisamente lo que las fórmulas multiproducto están diseñadas para resolver. Ten en cuenta que las MPF consisten en combinar los resultados de diferentes opciones de $k$ en un orden fijo $\chi$; no modifican el orden de la fórmula del producto subyacente.

Las fórmulas multiproducto (MPF) [1] construyen una combinación lineal ponderada de los valores esperados obtenidos a partir de varios circuitos de Trotter menos profundos, cada uno de los cuales utiliza un número diferente de pasos de Trotter $k_1, k_2, \ldots, k_r$ (un conjunto de recuentos de pasos $r$ ):

donde $\langle A \rangle_{k_j}(t)$ es el valor esperado de un observable $A$ en el instante $t$, estimado a partir de un circuito de Trotter con $k_j$ pasos, y los coeficientes $\{x_j\}_{j=1}^r$ se eligen de tal forma que los términos principales del error de Trotter en la combinación se anulen. Volveremos sobre esta expresión en el paso 4, donde la evaluaremos explícitamente para combinar nuestros resultados de Trotter. El aspecto práctico clave es que el circuito más profundo del MPF solo necesita un $k_{\max}$ es pasos, lo cual es mucho menor que el único $k$ que se necesitaría para alcanzar directamente el mismo error efectivo de Trotter. Los circuitos menos profundos hacen que el enfoque MPF sea más adecuado para hardware ruidoso.

¿Cómo se determinan los coeficientes?

Existen dos familias de coeficientes MPF:

Los coeficientes estáticos son independientes del hamiltoniano, del estado inicial y del tiempo de evolución. Se obtienen resolviendo un sistema lineal $Ax = b$ que garantiza la cancelación de los términos principales del error de Trotter. Para un conjunto de pasos de Trotter $\{k_j\}_{j=1}^r$ utilizado con una fórmula de producto simétrico de orden $2\chi$, al desarrollar el error de Trotter en potencias inversas de $k_j$ se obtienen ecuaciones de restricción de la forma:

donde los exponentes enteros $\{\eta_n\}$ son los órdenes de los términos sucesivos del error de Trotter para la fórmula de producto elegida. Para un PF simétrico de orden $2\chi$, el error principal en $\left[S_{2\chi}(t/k)\right]^k$ varía proporcionalmente a $1/k^{2\chi}$, con correcciones posteriores en $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — por lo que los exponentes son $\eta_n = 2\chi + 2n$. En el caso de los PF no simétricos, contribuyen tanto las potencias impares como las pares y $\eta_n = 2\chi + n$. Véase la referencia [1] para la derivación completa. La primera ecuación del sistema anterior garantiza la imparcialidad (el MPF reproduce el valor exacto de la esperanza en el límite « $k_j \to \infty$ »), y las ecuaciones restantes $r-1$ cancelan sucesivamente los primeros términos de error de Trotter $r-1$. Cuando la norma $L_1$ resultante $\|x\|_1$ es demasiado grande (lo que amplifica el ruido de muestreo), se puede resolver, en su lugar, un problema de optimización aproximada que limite $\|x\|_1$ al tiempo que minimice $\|Ax - b\|$.

Los coeficientes dinámicos [2], [3] dependen además del hamiltoniano, del estado inicial y del tiempo de evolución $t$. Estos minimizan la distancia, medida en la norma de Frobenius, entre el estado real evolucionado en el tiempo y la aproximación MPF:

donde $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ es la matriz de Gram de solapamientos entre los estados evolucionados según Trotter para diferentes recuentos de pasos $k_i, k_j$, y $L_i(t) = \mathrm{Tr}[\rho(t)\,\rho_{k_i}(t)]$ mide el solapamiento con el estado exacto (aproximado). En este tutorial, estas magnitudes se calculan de forma eficiente utilizando métodos de redes tensoriales, concretamente los backends de TeNPy-based en qiskit_addon_mpf.

Cuándo utilizar los MPF

Los planes de pensiones (MPF) resultan más ventajosos cuando:

• La profundidad del circuito es el cuello de botella. Si el ruido del hardware limita la profundidad a la que se puede trabajar, utiliza los MPF para conseguir una mayor precisión efectiva de Trotter en circuitos menos profundos.
• Necesitas valores de expectativa precisos, no una preparación completa del estado. Las MPF operan a nivel de los valores esperados: combinan números clásicos, no estados cuánticos. Por lo tanto, son ideales para la estimación de observables cuando se utiliza la primitiva «Estimator».
• Combinas un número moderado de pasos de trotto. Por lo general, basta con combinar $r = 3$ – $5$, con diferentes recuentos de pasos $k_j$, para anular varios términos de error de Trotter principales, al tiempo que se mantiene $\|x\|_1$ a un nivel manejable.

Cuándo los planes de pensiones de empleo (MPF) podrían no ser de ayuda

• Tiempos de evolución muy cortos. Cuando $t$ es lo suficientemente pequeño como para que una sola fórmula de Trotter de bajo orden ya sea precisa, no es necesario el esfuerzo adicional que supone ejecutar varios circuitos.
• Tareas de preparación para el examen estatal. Los MPF producen un valor esperado corregido, no un estado cuántico corregido. Si necesitas el estado real a lo largo del tiempo (por ejemplo, como entrada para otra subrutina cuántica), los MPF no son aplicables.
• Recuentos de pasos de trote que incumplen el régimen de convergencia. La derivación del coeficiente estático expande cada « $\left[S_{2\chi}(t/k_j)\right]^{k_j}$ » individual como una serie en $t/k_j$; esta expansión solo converge bien cuando $t/k_{\min} \lesssim 1$. Si se elige un valor demasiado pequeño para $k_{\min}$ en el caso de un $t$ dado, el circuito más superficial queda muy fuera del régimen perturbativo, los términos de error de orden superior que el MPF deja sin cancelar se vuelven grandes y la cancelación puede requerir coeficientes elevados. La norma « $L_1$ » $\|x\|_1$ es el criterio de diagnóstico práctico: cuando $\|x\|_1 \gg 1$, la sobrecarga de muestreo $\propto \|x\|_1^2$ podría superar la reducción del error de Trotter. Consulta la guía sobre cómo elegir los peldaños Trotter para obtener más información.

Contenido de este tutorial

Este tutorial explica paso a paso un flujo de trabajo completo de MPF en dos fases. En primer lugar, un ejemplo de simulador a pequeña escala (cadena de Heisenberg de 10 qubits) muestra cómo plantear el problema, calcular los coeficientes MPF estáticos y dinámicos, y comparar los valores esperados resultantes con los obtenidos mediante la diagonalización exacta. A continuación, un ejemplo de hardware a gran escala (cadena XXZ de 50 qubits) muestra cómo realizar la transpilación, ejecutarlo en un hardwar IBM Quantum o con mitigación de errores y procesar posteriormente los resultados utilizando los coeficientes MPF. A lo largo de todo el proceso, utilizamos el qiskit_addon_mpf paquete junto con las herramientas estándar de Qiskit.

Requisitos

Antes de empezar este tutorial, asegúrate de que tienes instalado lo siguiente:

• Qiskit SDK v2.0 o posterior, con soporte para visualización
• Qiskit Runtime v0.22 o posterior (pip install qiskit-ibm-runtime)
• Simulador de Qiskit Aer (pip install qiskit-aer)
• Complemento MPF Qiskit con el backend « TeNPy » (pip install "qiskit-addon-mpf[tenpy]")
• Utilidades del complemento de Qiskit (pip install qiskit-addon-utils)
• SciPy (pip install scipy)

Configuración

A continuación, recopilamos en una sola celda todas las importaciones de paquetes utilizadas a lo largo de este tutorial. XXPlusYYGateTambién definimos una CollectAndCollapse pasada del transpilador que fusiona las rotaciones adyacentes rxx y ryy en una sola. Esta pasada se aplica tanto durante la construcción del circuito en el paso 1 (para mantener bajo el número de puertas) como, de forma indirecta, cuando extraemos la estructura por capas para el MPF dinámico en el paso 4 (el TeNPy e espera puertas de dos qubits, no pares de rotaciones no fusionadas).

import warnings

import numpy as np
import matplotlib.pyplot as plt
from functools import partial
from copy import deepcopy

from qiskit import QuantumCircuit
from qiskit.quantum_info import Pauli, SparsePauliOp, Statevector
from qiskit.synthesis import SuzukiTrotter
from qiskit.transpiler import CouplingMap, PassManager
from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager
from qiskit.circuit.library import XXPlusYYGate
from qiskit.transpiler.passes.optimization.collect_and_collapse import (
    CollectAndCollapse,
    collect_using_filter_function,
    collapse_to_operation,
)

from qiskit_aer import AerSimulator
from qiskit_ibm_runtime import EstimatorV2 as Estimator, QiskitRuntimeService

from qiskit_addon_utils.problem_generators import (
    generate_xyz_hamiltonian,
    generate_time_evolution_circuit,
)
from qiskit_addon_utils.slicing import slice_by_depth
from qiskit_addon_mpf.static import setup_static_lse
from qiskit_addon_mpf.dynamic import setup_dynamic_lse
from qiskit_addon_mpf.costs import (
    setup_exact_problem,
    setup_sum_of_squares_problem,
    setup_frobenius_problem,
)
from qiskit_addon_mpf.backends.tenpy_layers import (
    LayerModel,
    LayerwiseEvolver,
)
from qiskit_addon_mpf.backends.tenpy_tebd import MPOState, MPS_neel_state

from scipy.linalg import expm

# Suppress TeNPy's `unit_cell_width` future-API warning. The default
# (`unit_cell_width=len(sites)`) is correct for Chain lattices, which is what
# `CouplingMap.from_line(...)` produces here, so the warning is informational.
warnings.filterwarnings(
    "ignore",
    message=r".*unit_cell_width.*",
    category=UserWarning,
)


# --- Helper: collect XX + YY rotations into a single gate ---
def filter_function(node):
    return node.op.name in {"rxx", "ryy"}


collect_function = partial(
    collect_using_filter_function,
    filter_function=filter_function,
    split_blocks=True,
    min_block_size=1,
)


def collapse_to_xx_plus_yy(block):
    param = 0.0
    for node in block.data:
        param += node.operation.params[0]
    return XXPlusYYGate(param)


collapse_function = partial(
    collapse_to_operation,
    collapse_function=collapse_to_xx_plus_yy,
)

pm = PassManager()
pm.append(CollectAndCollapse(collect_function, collapse_function))

Ejemplo de simulador a pequeña escala

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

Comenzamos con un modelo de Heisenberg de 10 qubits en una línea, utilizando el estado de Néel $\vert 0101\ldots01 \rangle$ como estado inicial. El hamiltoniano es:

donde « $J$ » es la intensidad del acoplamiento entre vecinos más cercanos. Medimos el correlador ZZ $Z_{L/2-1} Z_{L/2}$ en un par de qubits situados en el centro de la cadena, y utilizamos los pasos de Trotter $k_j = [1, 2, 4]$ con una fórmula de producto de segundo orden.

L = 10

# Generate coupling map and Hamiltonian
coupling_map = CouplingMap.from_line(L, bidirectional=False)

hamiltonian = generate_xyz_hamiltonian(
    coupling_map,
    coupling_constants=(1.0, 1.0, 1.0),
    ext_magnetic_field=(0.0, 0.0, 0.0),
)
print(hamiltonian)

Output:

SparsePauliOp(['IIIIIIIXXI', 'IIIIIIIYYI', 'IIIIIIIZZI', 'IIIIIXXIII', 'IIIIIYYIII', 'IIIIIZZIII', 'IIIXXIIIII', 'IIIYYIIIII', 'IIIZZIIIII', 'IXXIIIIIII', 'IYYIIIIIII', 'IZZIIIIIII', 'IIIIIIIIXX', 'IIIIIIIIYY', 'IIIIIIIIZZ', 'IIIIIIXXII', 'IIIIIIYYII', 'IIIIIIZZII', 'IIIIXXIIII', 'IIIIYYIIII', 'IIIIZZIIII', 'IIXXIIIIII', 'IIYYIIIIII', 'IIZZIIIIII', 'XXIIIIIIII', 'YYIIIIIIII', 'ZZIIIIIIII'],
              coeffs=[1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j,
 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j, 1.+0.j])

# Observable: ZZ on the middle pair of qubits
observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)
print(observable)

Output:

SparsePauliOp(['IIIIZZIIII'],
              coeffs=[1.+0.j])

# MPF parameters
mpf_trotter_steps = [1, 2, 4]
order = 2
symmetric = False

trotter_times = np.arange(0.5, 1.55, 0.1)
exact_evolution_times = np.arange(trotter_times[0], 1.55, 0.05)

Construir circuitos de Trotter

Creamos los circuitos aplicando las evoluciones temporales aproximadas de Trotter para cada instante y cada número de pasos de Trotter. El CollectAndCollapse paso definido en la sección «Configuración» agrupa las rotaciones XX y YY en puertas únicas de tipo XX+YY, con el fin de facilitar una simulación más eficiente de la red tensorial posteriormente.

# Initial Neel state preparation
initial_state_circ = QuantumCircuit(L)
initial_state_circ.x([i for i in range(L) if i % 2 != 0])


all_circs = []
for total_time in trotter_times:
    mpf_trotter_circs = [
        generate_time_evolution_circuit(
            hamiltonian,
            time=total_time,
            synthesis=SuzukiTrotter(reps=num_steps, order=order),
        )
        for num_steps in mpf_trotter_steps
    ]

    mpf_trotter_circs = pm.run(
        mpf_trotter_circs
    )  # Collect XX and YY into XX + YY

    mpf_circuits = [
        initial_state_circ.compose(circuit) for circuit in mpf_trotter_circs
    ]
    all_circs.append(mpf_circuits)

mpf_circuits[-1].draw("mpl", fold=-1)

Output:

Paso 2: Optimizar el problema para la ejecución en hardware cuántico

Para el ejemplo a pequeña escala, nos centramos en el simulador Aer. Antes de que los circuitos estén listos para ejecutarse, se producen dos transformaciones:

1. Recopilación de puertas a nivel de simulación hamiltoniana. XXPlusYYGateEn la celda «Configuración» hemos creado un CollectAndCollapse paso que fusiona las rotaciones adyacentes rxx y ryy en una sola. Ya aplicamos esta pasada cuando creamos los circuitos de Trotter en el paso 1 (la pm.run(...) llamada). Esto reduce el número de puertas de dos qubits y da lugar a una estructura más adecuada para la simulación mediante redes tensoriales, con vistas al cálculo posterior de los coeficientes dinámicos.

2. Reducción a la ISA del simulador. A continuación, ejecutamos el gestor de pasadas predefinido de Qiskit para optimization_level=3 adaptar cada circuito de Trotter a la arquitectura del conjunto de instrucciones (ISA) del simulador.

aer_sim = AerSimulator()
pm_sim = generate_preset_pass_manager(backend=aer_sim, optimization_level=3)

isa_circs_all_times = [
    pm_sim.run([deepcopy(c) for c in mpf_circuits])
    for mpf_circuits in all_circs
]

Paso 3: Ejecutar utilizando Qiskit primitives

Para el ejemplo a pequeña escala, aplicamos los circuitos de Trotter reducidos a ISA a la EstimatorV2 primitiva respaldada por Aer. De este modo, obtenemos un valor de referencia sin ruido para cada par « $(k_j, t)$ »; estos son los valores « $\langle A \rangle_{k_j}(t)$ » que el MPF combinará en el paso 4. Analizamos los periodos de evolución para poder trazar posteriormente la curva completa de la serie temporal de cada fórmula de producto individual y de la MPF.

estimator = Estimator(mode=aer_sim)

mpf_expvals_all_times, mpf_stds_all_times = [], []
for isa_circuits in isa_circs_all_times:
    result = estimator.run(
        [(circuit, observable) for circuit in isa_circuits], precision=0.005
    ).result()
    mpf_expvals_all_times.append([res.data.evs for res in result])
    mpf_stds_all_times.append([res.data.stds for res in result])

Paso 4: Procesamiento posterior y devolución del resultado en el formato clásico deseado

El paso 4 es donde se construye realmente el MPF. Aunque aquí se calculan los coeficientes $x_j$ (y, en el caso de la variante dinámica, este cálculo puede ser muy exigente), conceptualmente constituyen una fórmula clásica para combinar las mediciones cuánticas del paso 3 en un único valor esperado corregido; por lo tanto, consideramos que todo el proceso de cálculo de coeficientes y combinación forma parte del posprocesamiento.

Para evaluar en qué medida el MPF reproduce la dinámica real, calculamos en primer lugar los valores esperados exactos a lo largo del tiempo elevando directamente el hamiltoniano a una potencia exponencial. Esto solo es factible porque $L = 10$; en el ejemplo de hardware a gran escala que se muestra a continuación, tendremos que recurrir, en su lugar, a estimaciones basadas en redes tensoriales.

exact_expvals = []
for t in exact_evolution_times:
    exp_H = expm(-1j * t * hamiltonian.to_matrix())
    initial_state = Statevector(initial_state_circ).data
    time_evolved_state = exp_H @ initial_state

    exact_obs = (
        time_evolved_state.conj()
        @ observable.to_matrix()
        @ time_evolved_state
    ).real
    exact_expvals.append(exact_obs)

Coeficientes MPF estáticos

Los MPF estáticos utilizan coeficientes $x_j$ que son independientes del tiempo de evolución, del hamiltoniano y del estado inicial. Establecemos el sistema lineal $Ax = b$ descrito en la sección «Antecedentes» y resolvemos los coeficientes. La matriz $A$ viene determinada por el número de pasos de Trotter $k_j$, el orden $\chi$ de la fórmula del producto y si la fórmula es simétrica (lo que determina los exponentes $\eta_n$ ).

Para nuestro ejemplo a pequeña escala utilizamos un modelo de tipo « $k_j = [1, 2, 4]$ » con una fórmula de Suzuki-Trotter de orden no simétrico — $2\chi=2$ — (por lo que $\chi=1$ y $\eta_n = 2 + n$, lo que da $\eta_0 = 2,\, \eta_1 = 3$ ). El sistema queda así:

La primera fila garantiza la imparcialidad ( $\sum_j x_j = 1$ ); la segunda y la tercera fila eliminan, respectivamente, los términos de error de Trotter de primer orden $1/k^2$ y de orden siguiente $1/k^3$.

Configurar el LSE

Utilizamos setup_static_lse desde qiskit_addon_mpf.static para construir la matriz $A$ y el vector del lado derecho $b$ descritos anteriormente. La matriz $A$ depende no solo de $k_j$, sino también de la fórmula del producto que elijamos —en concreto, de su orden $\chi$ y de si es simétrica—. El symmetric indicador controla el patrón del exponente $\eta_n$ (las fórmulas simétricas solo producen términos de error de Trotter de potencia par; véase la ref. [1] ). Cabe señalar que, tal y como se muestra en la ref. [2], establecer symmetric=True no es estrictamente necesario, incluso cuando la PF subyacente es simétrica: la LSE no simétrica sigue siendo válida (aunque impone restricciones adicionales innecesarias).

Para nuestro ejemplo, ya hemos establecido order = 2 y symmetric = False en el paso 1.

lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)

Comprueba la matriz $A$ y el vector $b$ que se han construido para confirmar que coinciden con el sistema descrito anteriormente.

lse.A

Output:

array([[1.      , 1.      , 1.      ],
       [1.      , 0.25    , 0.0625  ],
       [1.      , 0.125   , 0.015625]])

lse.b

Output:

array([1., 0., 0.])

Una vez obtenida la ecuación de LSE, resolvemos los coeficientes estáticos $x_j$ mediante lse.solve() (esta es la solución directa $x = A^{-1}b$ ).

mpf_coeffs = lse.solve()
print(
    f"The static coefficients associated with the ansatze are: {mpf_coeffs}"
)

Output:

The static coefficients associated with the ansatze are: [ 0.04761905 -0.57142857  1.52380952]

Optimizar para $x$ utilizando un modelo exacto

Como alternativa al cálculo de $x = A^{-1}b$, puedes utilizar setup\_exact\_model para construir una instancia de cvxpy.Problem que utilice el LSE como restricciones y cuya solución óptima dé como resultado $x$.

model_exact, coeffs_exact = setup_exact_problem(lse)
model_exact.solve()
print(coeffs_exact.value)

Output:

[ 0.04761905 -0.57142857  1.52380952]

print(
    "L1 norm of the exact coefficients:",
    np.linalg.norm(coeffs_exact.value, ord=1),
)

Output:

L1 norm of the exact coefficients: 2.1428571428556378

Optimizar para $x$ utilizando un modelo aproximado

Podría darse el caso de que la norma « $L_1$ » para el conjunto elegido de valores de « $k_j$ » se considere demasiado elevada. Si ese es el caso y no puedes elegir otro conjunto de valores de $k_j$, puedes utilizar una solución aproximada que limite la norma $L_1$ a un umbral elegido, al tiempo que minimiza $\|Ax - b\|$. Consulta la guía sobre «Cómo utilizar el modelo aproximado ».

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=1.5
)
model_approx.solve()
print(coeffs_approx.value)
print(
    "L1 norm of the approximate coefficients:",
    np.linalg.norm(coeffs_approx.value, ord=1),
)

Output:

[-1.10294118e-03 -2.48897059e-01  1.25000000e+00]
L1 norm of the approximate coefficients: 1.5

Coeficientes MPF dinámicos

El MPF estático anula los términos de error de Trotter de una forma independiente del hamiltoniano y del estado, por lo que no produce necesariamente el menor error de aproximación posible para un hamiltoniano y un estado inicial dados. Por el contrario, el MPF dinámico (Ref. [2], [3] ) determina coeficientes dependientes del tiempo $x_i(t)$ que minimizan la distancia de la norma de Frobenius $\|\rho(t) - \mu^D(t)\|_F^2$ en cada instante $t$. Tal y como se muestra en la sección «Antecedentes», esto requiere la matriz de solapamiento $M_{ij}(t)$ entre los estados evolucionados según Trotter y el solapamiento $L_i(t)$ con el estado exacto; ambos se estiman utilizando backends de redes tensoriales ( TeNPy ) en qiskit_addon_mpf.

Para configurar el LSE dinámico necesitamos tres elementos:

1. Una fábrica de evolutores aproximada que el complemento ejecutará para cada $k_j$ con el fin de generar $\rho_{k_j}(t)$ como MPS/MPO. Lo construimos a partir de la estructura por capas del circuito de Trotter de orden $2$ (una capa por slice_by_depth), envuelto como un LayerwiseEvolver con parámetros de truncamiento de tipo « TeNPy ».
2. Una fábrica de evolutores exactos que produce un $\rho(t)$ de referencia de alta precisión. Utilizamos un circuito de Suzuki-Trotter de cuarto orden con paso de tiempo pequeño (dt=0.1, order=4) como aproximación a la evolución exacta.
3. Una fábrica de identidades y un MPS de estado inicial que sirven de base para la simulación « TeNPy ».

La celda siguiente crea la fábrica de evolutores aproximados.

# Create approximate time-evolution circuits
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)  # collect XX and YY

# Find layers in the circuit
layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)

# Create tensor network models
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

# Create the time-evolution object
approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

dt=0.1Aproximamos el estado exacto evolucionado en el tiempo mediante una fórmula de Suzuki-Trotter de cuarto orden, utilizando un paso de tiempo pequeño. Los parámetros de truncamiento de « TeNPy » pueden afectar a la precisión, por lo que es importante probar diferentes valores.

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {
            "chi_max": 64,
            "svd_min": 1e-8,
            "trunc_cut": None,
        },
        "max_delta_t": 2,
    },
)

Por último, definimos un identity_factory que da lugar al estado MPO inicial y preparamos el estado inicial de Néel como un MPS que se ajusta a la red utilizada por el modelo de Trotter en capas.

def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)


mps_initial_state = MPS_neel_state(models[0].lat)

Una vez establecidas las fábricas, calculamos ahora los coeficientes dinámicos en cada instante de evolución. Para cada $t$, setup_dynamic_lse se calculan las matrices de solapamiento pertinentes mediante TeNPy, y setup_frobenius_problem se devuelve un valor cvxpy.Problem que minimiza el coste de la norma de Frobenius. El solucionador devuelve los coeficientes $x_j(t)$ adaptados a ese momento; los recopilamos en mpf_dynamic_coeffs_list. Si el solucionador falla para un « $t$ » determinado, recurrimos a coeficientes nulos para que el bucle continúe.

mpf_dynamic_coeffs_list = []
for t in trotter_times:
    print(f"Computing dynamic coefficients for time={t}")
    lse = setup_dynamic_lse(
        mpf_trotter_steps,
        t,
        identity_factory,
        exact_factory,
        approx_factory,
        mps_initial_state,
    )
    problem, coeffs = setup_frobenius_problem(lse)
    try:
        problem.solve()
        mpf_dynamic_coeffs_list.append(coeffs.value)
    except Exception as error:
        mpf_dynamic_coeffs_list.append(np.zeros(len(mpf_trotter_steps)))
        print(error, "Calculation Failed for time", t)
    print("")

Output:

Computing dynamic coefficients for time=0.5

Computing dynamic coefficients for time=0.6

Computing dynamic coefficients for time=0.7

Computing dynamic coefficients for time=0.7999999999999999

Computing dynamic coefficients for time=0.8999999999999999

Computing dynamic coefficients for time=0.9999999999999999

Computing dynamic coefficients for time=1.0999999999999999

Computing dynamic coefficients for time=1.1999999999999997

Computing dynamic coefficients for time=1.2999999999999998

Computing dynamic coefficients for time=1.4

Computing dynamic coefficients for time=1.4999999999999998

Combinar los valores esperados de Trotter con los coeficientes del MPF

Ahora calculamos $\langle A \rangle_{\text{MPF}}(t) = \sum_j x_j \, \langle A \rangle_{k_j}(t)$ para cada conjunto de coeficientes (estático-exacto, estático-aproximado y dinámico), propagamos los errores estándar por circuito y representamos gráficamente las series temporales resultantes en función de la curva de diagonalización exacta.

sym = {1: "^", 2: "s", 4: "p"}
# Get expectation values at all times for each Trotter step
for k, step in enumerate(mpf_trotter_steps):
    trotter_curve, trotter_curve_error = [], []
    for trotter_expvals, trotter_stds in zip(
        mpf_expvals_all_times, mpf_stds_all_times
    ):
        trotter_curve.append(trotter_expvals[k])
        trotter_curve_error.append(trotter_stds[k])

    plt.errorbar(
        trotter_times,
        trotter_curve,
        yerr=trotter_curve_error,
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

# Get expectation values at all times for the static MPF with exact coeffs
exact_mpf_curve, exact_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_exact.value, trotter_stds)
            ]
        )
    )
    exact_mpf_curve_error.append(mpf_std)
    exact_mpf_curve.append(trotter_expvals @ coeffs_exact.value)

plt.errorbar(
    trotter_times,
    exact_mpf_curve,
    yerr=exact_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Exact",
    color="purple",
)


# Get expectation values at all times for the static MPF with approximate coeffs
approx_mpf_curve, approx_mpf_curve_error = [], []
for trotter_expvals, trotter_stds in zip(
    mpf_expvals_all_times, mpf_stds_all_times
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(coeffs_approx.value, trotter_stds)
            ]
        )
    )
    approx_mpf_curve_error.append(mpf_std)
    approx_mpf_curve.append(trotter_expvals @ coeffs_approx.value)

plt.errorbar(
    trotter_times,
    approx_mpf_curve,
    yerr=approx_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Static MPF - Approx",
    color="orange",
)


# Get expectation values at all times for the dynamic MPF
dynamic_mpf_curve, dynamic_mpf_curve_error = [], []
for trotter_expvals, trotter_stds, dynamic_coeffs in zip(
    mpf_expvals_all_times, mpf_stds_all_times, mpf_dynamic_coeffs_list
):
    mpf_std = np.sqrt(
        sum(
            [
                (coeff**2) * (std**2)
                for coeff, std in zip(dynamic_coeffs, trotter_stds)
            ]
        )
    )
    dynamic_mpf_curve_error.append(mpf_std)
    dynamic_mpf_curve.append(trotter_expvals @ dynamic_coeffs)

plt.errorbar(
    trotter_times,
    dynamic_mpf_curve,
    yerr=dynamic_mpf_curve_error,
    markersize=4,
    marker="o",
    label="Dynamic MPF",
    color="pink",
)


# Exact expectation values
plt.plot(
    exact_evolution_times,
    exact_expvals,
    color="red",
    linestyle="--",
    label="Exact time-evolution",
)

plt.title(f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ vs time")
plt.xlabel("Time")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

Output:

El gráfico anterior ilustra la relación entre el error de Trotter y el error de muestreo.

• Error de Trotter. Las fórmulas de cada producto (marcadores grises) se desvían cada vez más de la curva exacta a medida que pasa el tiempo. El circuito « $k=1$ » presenta la mayor desviación y es el menos profundo, pero también se encuentra ya en el régimen en el que « $t/k \gtrsim 1$ », por lo que el término de error principal « $1/k^{2}$ » es grande. Las combinaciones MPF (marcadores de color) anulan varios de estos términos de error de Trotter principales, por lo que siguen la curva exacta mucho más fielmente que cualquier circuito « $k_j$ » por sí solo. La desviación restante refleja los términos de Trotter de orden superior que el MPF no cancela: un MPF estático de orden $2$, $r=3$ solo elimina los dos primeros órdenes de error, y a gran $t/k_{\min}$ la cola no cancelada acaba dominando; por lo tanto, el MPF no garantiza que los circuitos muy poco profundos mantengan su precisión en momentos arbitrarios.

• Error de muestreo. Las barras de error más anchas en las curvas del MPF son una consecuencia directa de la combinación lineal: al propagar los errores estándar independientes por circuito $\sigma_{k_j}$ se obtiene una varianza total $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$. Por lo tanto, cuanto mayor sea el $\|x\|_2$ (y, en la práctica, el $\|x\|_1$, que es lo que controlamos), más mediciones se necesitarán para alcanzar una incertidumbre objetivo determinada. Esta es la compensación que subyace a la opción «solucionador aproximado» en «Antecedentes»: limitamos $\|x\|_1$ para que esta sobrecarga sea manejable. Es fundamental señalar que, a diferencia del error de Trotter, el error de muestreo se reduce con el « $1/\sqrt{N_{\text{shots}}}$ », por lo que siempre puede minimizarse realizando más mediciones.

En el ejemplo de hardware a gran escala que se muestra a continuación, el ruido del hardware se introduce como una fuente de error adicional en cada « $\langle A \rangle_{k_j}$ », que, de forma similar, se amplifica mediante los coeficientes del MPF. En esa sección veremos cómo interactúa la mitigación de errores con los MPF.

Ejemplo de hardware a gran escala

En esta sección ampliamos la escala del problema más allá de lo que es posible simular con exactitud. Reproducimos algunos de los resultados presentados en la ref. [3], utilizando una cadena XXZ de 50 qubits en el intervalo de tiempo $t = 3$. Seguimos el mismo procedimiento de cuatro pasos que en el ejemplo a pequeña escala, pero en esta ocasión nos centramos en hardware cuántico real con mitigación de errores. Al igual que en la plantilla, cada paso está marcado directamente en el código, y un mismo paso puede abarcar varias celdas cuando conviene examinar sus resultados intermedios.

La representación sigue el ejemplo a pequeña escala: definir un hamiltoniano, elegir los parámetros de Trotter, calcular los coeficientes MPF (estáticos y dinámicos) y construir circuitos. Las diferencias principales son:

• Un hamiltoniano XXZ en 50 sitios con acoplamientos aleatorios extraídos de $\mathcal{U}(0.5, 1.5)$ (Ref. [3] ).
• Una fórmula de Trotter simétrica de segundo orden con « $k_j = [3, 4, 6]$ » (por lo que $\chi=1$, symmetric=True).
• Un único tiempo de evolución fijo $t = 3$. Con $k_{\min}=3$, esto da como resultado $t/k_{\min}=1$, lo que mantiene los componentes superficiales dentro del régimen de convergencia de Trotter, donde es válido el modelo de error dominante en el que se basa el MPF.
• Una comparación adicional de un solo circuito ejecutada con pasos de Trotter $k = 10$, utilizada como línea base. Elegimos « $k = 10$ » porque su profundidad de dos qubits en el hardware es mayor que la del componente MPF más profundo ( $k_{\max}=6$ ) más la sobrecarga que supone ejecutar múltiples circuitos MPF —lo suficientemente profunda como para estar limitada por el ruido, que es el régimen en el que se espera que la combinación de MPF supere al circuito único de referencia—. Se trata de una comparación de «un solo circuito profundo» frente a la combinación MPF, no de un circuito que tenga como objetivo el error de Trotter efectivo del MPF (lo cual requeriría muchos más pasos).

Ten en cuenta que, aunque aquí todavía nos encontramos en el paso 1 (mapeo y construcción del circuito), en esta celda también calculamos previamente los coeficientes dinámicos junto con los estáticos. Los coeficientes dinámicos dependen de $H$ y $t$, pero no de las mediciones cuánticas, por lo que pueden calcularse en cualquier momento antes del paso 4. Lo hacemos ahora para tener toda la configuración específica del MPF en un solo lugar.

# -------------------------Step 1-------------------------
L = 50
coupling_map = CouplingMap.from_line(L, bidirectional=False)

# XXZ Hamiltonian with random couplings (Ref. [3])
np.random.seed(0)
even_edges = list(coupling_map.get_edges())[::2]
odd_edges = list(coupling_map.get_edges())[1::2]

Js = np.random.uniform(0.5, 1.5, size=L)
hamiltonian = SparsePauliOp(Pauli("I" * L))
for i, edge in enumerate(even_edges + odd_edges):
    hamiltonian += SparsePauliOp.from_sparse_list(
        [
            ("XX", (edge), 2 * Js[i]),
            ("YY", (edge), 2 * Js[i]),
            ("ZZ", (edge), 4 * Js[i]),
        ],
        num_qubits=L,
    )

observable = SparsePauliOp.from_sparse_list(
    [("ZZ", (L // 2 - 1, L // 2), 1.0)], num_qubits=L
)

total_time = 3
mpf_trotter_steps = [3, 4, 6]
order = 2
symmetric = True

# Static coefficients
lse = setup_static_lse(mpf_trotter_steps, order=order, symmetric=symmetric)
mpf_coeffs = lse.solve()
print(f"Static coefficients: {mpf_coeffs}")
print(f"L1 norm: {np.linalg.norm(mpf_coeffs, ord=1)}")

model_approx, coeffs_approx = setup_sum_of_squares_problem(
    lse, max_l1_norm=2.0
)
model_approx.solve()
print(f"Approximate coefficients: {coeffs_approx.value}")
print(f"L1 norm (approx): {np.linalg.norm(coeffs_approx.value, ord=1)}")

# -------------------------Dynamic coefficients-------------------------
single_2nd_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=order)
)
single_2nd_order_circ = pm.run(single_2nd_order_circ)

layers = slice_by_depth(single_2nd_order_circ, max_slice_depth=1)
models = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz") for layer in layers
]

approx_factory = partial(
    LayerwiseEvolver,
    layers=models,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 4,
    },
)

single_4th_order_circ = generate_time_evolution_circuit(
    hamiltonian, time=1.0, synthesis=SuzukiTrotter(reps=1, order=4)
)
single_4th_order_circ = pm.run(single_4th_order_circ)
exact_model_layers = [
    LayerModel.from_quantum_circuit(layer, conserve="Sz")
    for layer in slice_by_depth(single_4th_order_circ, max_slice_depth=1)
]

exact_factory = partial(
    LayerwiseEvolver,
    layers=exact_model_layers,
    dt=0.1,
    options={
        "preserve_norm": False,
        "trunc_params": {"chi_max": 64, "svd_min": 1e-8, "trunc_cut": None},
        "max_delta_t": 3,
    },
)


def identity_factory():
    return MPOState.initialize_from_lattice(models[0].lat, conserve=True)


mps_initial_state = MPS_neel_state(models[0].lat)

print(f"Computing dynamic coefficients for time={total_time}")
lse_dyn = setup_dynamic_lse(
    mpf_trotter_steps,
    total_time,
    identity_factory,
    exact_factory,
    approx_factory,
    mps_initial_state,
)
problem, coeffs_dyn = setup_frobenius_problem(lse_dyn)
try:
    problem.solve()
    mpf_dynamic_coeffs = coeffs_dyn.value
except Exception as error:
    mpf_dynamic_coeffs = np.zeros(len(mpf_trotter_steps))
    print(error, "Calculation Failed")

# -------------------------Step 1 (cont): Build circuits-------------------------
mpf_circuits = []
for k in mpf_trotter_steps:
    circuit = QuantumCircuit(L)
    circuit.x([i for i in range(L) if i % 2])
    trotter_circ = generate_time_evolution_circuit(
        hamiltonian,
        synthesis=SuzukiTrotter(reps=k, order=order),
        time=total_time,
    )
    circuit.compose(trotter_circ, qubits=range(L), inplace=True)
    mpf_circuits.append(circuit)

# Baseline "single deep circuit" comparison run with k=10 Trotter steps.
# Its two-qubit depth is deeper than the deepest MPF constituent (k_max=6) plus
# the overhead of running multiple circuits, pushing it into the noise-limited
# regime where MPF is expected to outperform. It does NOT target the MPF's effective
# Trotter error (which would require many more steps).
comp_circuit = QuantumCircuit(L)
comp_circuit.x([i for i in range(L) if i % 2])
trotter_circ = generate_time_evolution_circuit(
    hamiltonian,
    synthesis=SuzukiTrotter(reps=10, order=order),
    time=total_time,
)
comp_circuit.compose(trotter_circ, qubits=range(L), inplace=True)
mpf_circuits.append(comp_circuit)

Output:

Static coefficients: [ 0.42857143 -1.82857143  2.4       ]
L1 norm: 4.65714285714286
Approximate coefficients: [-0.4942491   0.40206845  1.09218065]
L1 norm (approx): 1.9884981979026675
Computing dynamic coefficients for time=3

Ahora optimizamos los circuitos para el backend elegido. optimization_level=3Utilizamos el gestor de pasadas preconfigurado de Qiskit, que selecciona automáticamente un buen conjunto de qubits físicos y asigna cada circuito a la topología del dispositivo.

# -------------------------Step 2-------------------------
service = QiskitRuntimeService()
# backend = service.least_busy(operational=True, simulator=False, min_num_qubits=L)
backend = service.backend("ibm_fez")
print(backend)

transpiler = generate_preset_pass_manager(
    optimization_level=3, backend=backend
)
transpiled_circuits = [transpiler.run(circ) for circ in mpf_circuits]

isa_observables = [
    observable.apply_layout(circ.layout) for circ in transpiled_circuits
]

Output:

<IBMBackend('ibm_fez')>

La ejecución de circuitos más complejos en hardware real requiere medidas enérgicas de mitigación de errores. Permitimos el desacoplamiento dinámico, la rotación de puertas y mediciones, la mitigación de errores de medición y la extrapolación sin ruido (ZNE). Cabe señalar que los factores de ruido ZNE que utilizamos aquí (1, 1.2, 1.4) son menores que en un escenario de circuito superficial, ya que los componentes MPF más profundos ya se encuentran cerca del umbral de ruido y unas amplificaciones de ruido elevadas los llevarían más allá del punto en el que la extrapolación ZNE resulta fiable.

Enviamos los cuatro circuitos (los tres componentes del MPF en $k_j = [3, 4, 6]$ más la línea de referencia de $k = 10$ ) en un único trabajo de Estimator.

# -------------------------Step 3-------------------------
estimator = Estimator(mode=backend)
estimator.options.default_shots = 30000

# Error suppression/mitigation
estimator.options.dynamical_decoupling.enable = True
estimator.options.twirling.enable_gates = True
estimator.options.twirling.enable_measure = True
estimator.options.twirling.num_randomizations = "auto"
estimator.options.twirling.strategy = "active-accum"
estimator.options.resilience.measure_mitigation = True
estimator.options.experimental.execution_path = "gen3-turbo"

estimator.options.resilience.zne_mitigation = True
estimator.options.resilience.zne.noise_factors = (1, 1.2, 1.4)
estimator.options.resilience.zne.extrapolator = "linear"

estimator.options.environment.job_tags = ["TUT_MPF"]

job_50 = estimator.run(
    [
        (circ, observable)
        for circ, observable in zip(transpiled_circuits, isa_observables)
    ]
)

Extraemos los valores esperados y las desviaciones estándar por circuito de los resultados del trabajo y, a continuación, los combinamos con cada conjunto de coeficientes MPF exactamente igual que en el ejemplo a pequeña escala: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$, con la varianza propagada $\sigma^2 = \sum_j x_j^2 \sigma_{k_j}^2$.

# -------------------------Step 4-------------------------
result = job_50.result()
evs = [res.data.evs for res in result]
std = [res.data.stds for res in result]

print(evs)
print(std)

Output:

[array(-0.07916195), array(-0.04479681), array(-0.2560756), array(-0.06045848)]
[array(0.04605538), array(0.10056336), array(0.14426151), array(0.04059092)]

exact_mpf_std = np.sqrt(
    sum([(coeff**2) * (std**2) for coeff, std in zip(mpf_coeffs, std[:3])])
)
print(
    "Exact static MPF expectation value: ",
    evs[:3] @ mpf_coeffs,
    "+-",
    exact_mpf_std,
)
approx_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(coeffs_approx.value, std[:3])
        ]
    )
)
print(
    "Approximate static MPF expectation value: ",
    evs[:3] @ coeffs_approx.value,
    "+-",
    approx_mpf_std,
)
dynamic_mpf_std = np.sqrt(
    sum(
        [
            (coeff**2) * (std**2)
            for coeff, std in zip(mpf_dynamic_coeffs, std[:3])
        ]
    )
)
print(
    "Dynamic MPF expectation value: ",
    evs[:3] @ mpf_dynamic_coeffs,
    "+-",
    dynamic_mpf_std,
)

Output:

Exact static MPF expectation value:  -0.5665938395816946 +- 0.3925273058119915
Approximate static MPF expectation value:  -0.25856647611537903 +- 0.164249927266166
Dynamic MPF expectation value:  -0.12667812062949296 +- 0.06059471006973169

sym = {3: "^", 4: "s", 6: "p"}
for k, step in enumerate(mpf_trotter_steps):
    plt.errorbar(
        k,
        evs[k],
        yerr=std[k],
        alpha=0.5,
        markersize=4,
        marker=sym[step],
        color="grey",
        label=f"{mpf_trotter_steps[k]} Trotter steps",
    )

plt.errorbar(
    3,
    evs[-1],
    yerr=std[-1],
    alpha=0.5,
    markersize=8,
    marker="x",
    color="blue",
    label="10 Trotter steps",
)

plt.errorbar(
    4,
    evs[:3] @ mpf_coeffs,
    yerr=exact_mpf_std,
    markersize=4,
    marker="o",
    color="purple",
    label="Static MPF",
)

plt.errorbar(
    5,
    evs[:3] @ coeffs_approx.value,
    yerr=approx_mpf_std,
    markersize=4,
    marker="o",
    color="orange",
    label="Approximate static MPF",
)

plt.errorbar(
    6,
    evs[:3] @ mpf_dynamic_coeffs,
    yerr=dynamic_mpf_std,
    markersize=4,
    marker="o",
    color="pink",
    label="Dynamic MPF",
)

exact_obs = -0.24384471447172074  # Calculated via Tensor Network calculation
plt.axhline(
    y=exact_obs, linestyle="--", color="red", label="Exact time-evolution"
)

plt.title(
    f"$\\langle Z_{{{L//2-1}}} Z_{{{L//2}}} \\rangle$ at time {total_time} for the different methods"
)
plt.xlabel("Method")
plt.ylabel("Expectation Value")
plt.legend(loc="upper center", bbox_to_anchor=(0.5, -0.2), ncol=2)
plt.grid(alpha=0.1)
plt.tight_layout()
plt.show()

Output:

Algunas observaciones sobre los resultados de hardware anteriores:

• Profundizar más no es gratuito en lo que respecta al hardware. Las curvas de referencia de un solo circuito lo dejan claro: el circuito $k = 6$ es prácticamente exacto ( $-0.256$ frente al de referencia $-0.244$ ), mientras que la curva de referencia más profunda $k = 10$ es peor ( $-0.061$, con un error de $\sim 0.18$ ), y no mejor. Una vez que el error de Trotter ya es pequeño, añadir pasos solo sirve, en su mayor parte, para hacer más profundo el circuito y acumular más ruido de puerta y decoherencia. Este es precisamente el régimen para el que se han diseñado los MPF: alcanzar la precisión de un circuito profundo utilizando únicamente componentes superficiales.

• Un MPF de norma pequeña supera al circuito único profundo. El MPF «aproximadamente estático» (con un límite máximo de $\|x\|_1 \approx 2$ ) se sitúa en $-0.259$, a $\sim 0.015$ del valor de referencia y mucho más cerca que la referencia de $k = 10$. El MPF dinámico ( $-0.127$ ) también supera con holgura ese nivel de referencia. Ambos combinan únicamente los circuitos superficiales $k_j = [3, 4, 6]$, pero obtienen una respuesta que el circuito único profundo no podía proporcionar.

• La norma del coeficiente es más importante que la optimalidad matemática. El MPF estático exacto tiene una norma de coeficientes de $\|x\|_1 = 4.66$ y es el peor estimador de todos ( $-0.567$, con un error de más de $0.3$ ): la elevada norma de coeficientes amplifica el ruido residual de la puerta, la decoherencia y el error ZNE en cada $\langle A \rangle_{k_j}$ aproximadamente en el mismo factor, lo que anula la cancelación del error de Trotter que proporciona. La limitación de la norma (el solucionador estático aproximado, $\|x\|_1 \approx 2$ ) elimina esta sobrecarga y ofrece la mejor estimación, aunque sus coeficientes ya no anulen exactamente el error principal de Trotter.

• Los circuitos individuales de poca profundidad pueden seguir siendo competitivos. El único componente « $k = 6$ » ( $-0.256$ ) es, en sí mismo, prácticamente exacto en este caso; en esta ejecución, incluso se acerca ligeramente más que el MPF «aproximate-static». El problema es que no se sabe de antemano qué $k$ concreto se encuentra en el punto óptimo de «convergente pero aún no limitado por el ruido», y la opción que parece más segura —simplemente profundizar más ( $k = 10$ ) para garantizar la convergencia de Trotter— es precisamente la que falla. El MPF ofrece una combinación basada en principios de circuitos poco profundos que no requiere adivinar la profundidad adecuada.

En la práctica, esto significa que, en el ámbito del hardware, los MPF deben combinarse con una sólida mitigación de errores en cada « $\langle A \rangle_{k_j}$ » individual; la norma del coeficiente $L_1$ debe mantenerse en un nivel moderado (utilizando el solucionador aproximado o el MPF dinámico); y los pasos de Trotter $k_j$ deben elegirse de tal forma que $t/k_{\min} \lesssim 1$ —aquí $k_{\min} = 3$ en $t = 3$ da como resultado $t/k_{\min} = 1$ —, manteniendo los componentes dentro del régimen de convergencia en el que es válido el modelo de error principal en el que se basa el MPF estático. Con esas opciones, los MPF de norma pequeña aquí se equiparan a un circuito único convergente, mientras que la línea de base «simplemente ir más profundo» no lo hace, recuperando así la ventaja de profundidad frente a precisión mostrada en la ref. [3]. Cabe señalar también que las ejecuciones individuales presentan ruido: en un envío diferente del mismo trabajo (o en un backend diferente), el orden exacto puede variar; las tendencias generales indican que los MPF de pequeño $\|x\|_1$ dan buenos resultados, que el MPF de gran $\|x\|_1$ estático exacto se ve amplificado por el ruido del hardware y que el circuito único excesivamente profundo está limitado por el ruido.

Próximos pasos

Referencias

[1] Vázquez, A. C., Egger, D. J., Ochsner, D., y Woerner, S. Fórmulas multiproducto bien condicionadas para la simulación hamiltoniana adaptada al hardware. Quantum, 7, 1067 (2023)

[2] Zhuk, S., Robertson, N. F., y Bravyi, S. Límites de error de Trotter y fórmulas dinámicas multiproducto para la simulación hamiltoniana. Physical Review Research, 6(3), 033309 (2024)

[3] Robertson, N. F., et al. Fórmulas dinámicas multiproducto mejoradas mediante redes tensoriales. arXiv:2407.17405 (2024)