Fórmulas multiproducto para reducir el error de Trotter

Uso estimado: Cuatro minutos en un procesador Heron r2 (NOTA: Esto es solo una estimación. Tu tiempo de ejecución puede variar.)

Resultados de aprendizaje

Cómo las fórmulas multiproducto (MPF) reducen el error de Trotter en la simulación hamiltoniana combinando valores esperados de múltiples circuitos poco profundos
Cuándo las MPF son beneficiosas frente a las fórmulas producto estándar y cuándo no son la herramienta adecuada
Cómo calcular coeficientes MPF estáticos y dinámicos usando el paquete qiskit_addon_mpf
Cómo ejecutar un flujo de trabajo MPF de principio a fin en hardware de IBM Quantum®, incluyendo transpilación, mitigación de errores y postprocesamiento

Requisitos previos

Métodos de compilación para circuitos de simulación hamiltoniana — presenta los circuitos de Trotter (fórmula producto) en Qiskit.
Fórmulas producto en Qiskit, en particular las clases de síntesis SuzukiTrotter y LieTrotter.
Primitivas de Qiskit y la interfaz Estimator.

Antecedentes

¿Qué son las fórmulas multiproducto?

Al simular sistemas cuánticos en una computadora cuántica, una tarea central es aproximar el operador de evolución temporal $e^{-iHt}$ para un hamiltoniano $H$ . El enfoque estándar utiliza fórmulas producto (PF), también conocidas como descomposiciones de Trotter-Suzuki. Estas descomponen $H = \sum_{a=1}^d F_a$ en términos cuyas unitarias individuales $e^{-iF_a t}$ son eficientes de implementar, y luego aproximan la evolución completa como un producto ordenado de estas unitarias más simples.

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

S_1(t) := \prod_{a=1}^d e^{-i F_a t},

que incurre en 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 producto simétrica (ver Ref. [1]), convergen más rápido como $e^{-iHt} + \mathcal{O}(t^{2\chi+1})$ , pero a costa de circuitos más profundos por paso.

Para reducir el error a un orden fijo $\chi$ , generalmente 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 concatenan:

e^{-iHt} \approx \left[S_{2\chi}(t/k)\right]^k.

Para una fórmula simétrica de orden $2\chi$ , el error residual de Trotter escala entonces como $\mathcal{O}\!\left(t^{2\chi+1} / k^{2\chi}\right)$ . Así, aumentar $k$ suprime rápidamente el error de Trotter — pero también profundiza el circuito de forma lineal, y en hardware con ruido eso significa más ruido de compuertas acumulado. Esta tensión entre el error de Trotter (favorece un $k$ mayor) y el ruido del hardware (favorece un $k$ menor) es precisamente lo que las fórmulas multiproducto están diseñadas para resolver. Ten en cuenta que las MPF se tratan de combinar resultados de diferentes elecciones de $k$ a un orden fijo $\chi$ — no cambian el orden de la fórmula de producto subyacente.

Las fórmulas multiproducto (MPF) [1] construyen una combinación lineal ponderada de valores esperados obtenidos de varios circuitos de Trotter más superficiales, cada uno usando un número diferente de pasos de Trotter $k_1, k_2, \ldots, k_r$ (un conjunto de $r$ recuentos de pasos):

\langle A \rangle_{\text{MPF}}(t) = \sum_{j=1}^r x_j \, \langle A \rangle_{k_j}(t),

donde $\langle A \rangle_{k_j}(t)$ es el valor esperado de un observable $A$ en el tiempo $t$ estimado a partir de un circuito de Trotter con $k_j$ pasos, y los coeficientes $\{x_j\}_{j=1}^r$ se eligen de modo que los términos principales del error de Trotter en la combinación se cancelen. Retomaremos esta expresión en el Paso 4, donde la evaluaremos explícitamente para combinar nuestros resultados de Trotter. El punto práctico clave es que el circuito más profundo en la MPF solo necesita $k_{\max}$ pasos, que es mucho menor que el único $k$ que se requeriría para alcanzar directamente el mismo error de Trotter efectivo. Los circuitos más superficiales hacen que el enfoque MPF sea más adecuado para hardware con ruido.

¿Cómo se determinan los coeficientes?

Existen dos familias de coeficientes MPF:

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

\sum_{j=1}^r x_j = 1, \quad \sum_{j=1}^r \frac{x_j}{k_j^{\eta_n}} = 0 \quad (n = 0, \ldots, r-2),

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 una PF simétrica de orden $2\chi$ , el error principal en $\left[S_{2\chi}(t/k)\right]^k$ escala como $1/k^{2\chi}$ , con correcciones subsiguientes en $1/k^{2\chi+2}, 1/k^{2\chi+4}, \ldots$ — por lo que los exponentes son $\eta_n = 2\chi + 2n$ . Para PF no simétricas, contribuyen tanto potencias impares como pares y $\eta_n = 2\chi + n$ . Consulta la Ref. [1] para la derivación completa. La primera ecuación del sistema anterior garantiza la insesgadez (la MPF reproduce el valor esperado exacto en el límite $k_j \to \infty$ ), y las restantes $r-1$ ecuaciones cancelan sucesivamente los primeros $r-1$ términos del error de Trotter. Cuando la norma $L_1$ resultante $\|x\|_1$ es demasiado grande (lo que amplifica el ruido de muestreo), puedes resolver en su lugar una optimización aproximada que limita $\|x\|_1$ mientras minimiza $\|Ax - b\|$ .

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

\|\rho(t) - \mu^D(t)\|_F^2 = 1 + \sum_{i,j} M_{ij}(t)\, x_i(t)\, x_j(t) - 2\sum_i L_i(t)\, x_i(t),

donde $M_{ij}(t) = \mathrm{Tr}[\rho_{k_i}(t)\,\rho_{k_j}(t)]$ es la matriz de Gram de solapamientos entre estados evolucionados por 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 cantidades se calculan de manera eficiente usando métodos de redes tensoriales, específicamente los backends basados en TeNPy en qiskit_addon_mpf.

Cuándo usar MPF

Las MPF son más beneficiosas cuando:

La profundidad del circuito es el cuello de botella. Si el ruido del hardware limita la profundidad que puedes ejecutar, usa MPF para lograr mayor precisión efectiva de Trotter a partir de circuitos más superficiales.
Necesitas valores esperados precisos, no una preparación completa del estado. Las MPF operan a nivel de valores esperados — combinan números clásicos, no estados cuánticos. Por eso son ideales para la estimación de observables al usar la primitiva Estimator.
Combinas un número modesto de recuentos de pasos de Trotter. Normalmente, combinar $r = 3$ – $5$ recuentos de pasos diferentes $k_j$ es suficiente para cancelar varios de los términos principales del error de Trotter mientras se mantiene $\|x\|_1$ manejable.

Cuándo las MPF podrían no ayudar

Tiempos de evolución muy cortos. Cuando $t$ es lo suficientemente pequeño como para que una fórmula de Trotter de bajo orden ya sea precisa, la sobrecarga de ejecutar múltiples circuitos es innecesaria.
Tareas de preparación de estado. Las MPF producen un valor esperado corregido, no un estado cuántico corregido. Si necesitas el estado real evolucionado en el tiempo (por ejemplo, como entrada para otra subrutina cuántica), las MPF no son aplicables.
Recuentos de pasos de Trotter que violan el régimen de convergencia. La derivación de coeficientes estáticos 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 $k_{\min}$ se elige demasiado pequeño para el $t$ dado, el circuito más superficial queda muy fuera del régimen perturbativo, los términos de error de orden superior que la MPF deja sin cancelar se vuelven grandes, y la cancelación puede requerir coeficientes grandes. La norma $L_1$ $\|x\|_1$ es el 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 pasos de Trotter para más detalles.

Qué cubre este tutorial

Este tutorial recorre un flujo de trabajo MPF de extremo a extremo en dos etapas. Primero, un ejemplo de simulador a pequeña escala (cadena de Heisenberg de 10 qubits) demuestra cómo plantear el problema, calcular coeficientes MPF estáticos y dinámicos, y comparar los valores esperados resultantes con la diagonalización exacta. Luego, un ejemplo de hardware a gran escala (cadena XXZ de 50 qubits) muestra cómo transpilar, ejecutar en hardware de IBM Quantum con mitigación de errores, y procesar los resultados usando los coeficientes MPF. A lo largo del proceso, usamos el paquete qiskit_addon_mpf junto con herramientas estándar de Qiskit.

Requisitos

Antes de comenzar este tutorial, asegúrate de tener instalado lo siguiente:

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

Configuración

A continuación, recopilamos todas las importaciones de paquetes usadas a lo largo de este tutorial en una sola celda. También definimos un pase de transpilador CollectAndCollapse que fusiona rotaciones adyacentes rxx y ryy en una única XXPlusYYGate. Este pase se aplica tanto durante la construcción del circuito en el Paso 1 (para mantener bajo el recuento de compuertas) como indirectamente cuando extraemos la estructura por capas para la MPF dinámica en el Paso 4 (TeNPy espera compuertas de dos qubits, no pares de rotaciones sin fusionar).

# Added by doQumentation — installs the packages this notebook needs if they are missing
import importlib.util

_needed = {"matplotlib": "matplotlib", "numpy": "numpy", "qiskit": "qiskit", "qiskit_addon_mpf": "qiskit-addon-mpf", "qiskit_addon_utils": "qiskit-addon-utils", "qiskit_aer": "qiskit-aer", "qiskit_ibm_runtime": "qiskit-ibm-runtime", "scipy": "scipy", "tenpy": "physics-tenpy"}
_missing = [pip for module, pip in _needed.items()
            if importlib.util.find_spec(module) is None]
# One at a time, so a package that fails to install does not block the others
for _pip in _missing:
    %pip install -q {_pip}
if not _missing:
    print("\u2713 All required packages are installed")

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: Mapear las entradas clásicas a un problema cuántico

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

\hat{\mathcal{H}}_{\text{Heis}} = J \sum_{i=1}^{L-1} \left(X_i X_{i+1} + Y_i Y_{i+1} + Z_i Z_{i+1}\right),

donde $J$ es la fuerza de acoplamiento entre vecinos más cercanos. Medimos el correlador ZZ $Z_{L/2-1} Z_{L/2}$ en un par de qubits en el medio de la cadena, y usamos 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)

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)

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 que implementan las evoluciones aproximadas de Trotter en el tiempo para cada punto temporal y cada recuento de pasos de Trotter. El pase CollectAndCollapse definido en la sección de Configuración agrupa las rotaciones XX y YY en compuertas XX+YY únicas, para preparar una simulación de redes tensoriales más eficiente más adelante.

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

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

Para el ejemplo a pequeña escala, apuntamos al simulador Aer. Ocurren dos transformaciones antes de que los circuitos estén listos para ejecutarse:

Recolección de compuertas a nivel de simulación del hamiltoniano. En la celda de Configuración construimos un pase CollectAndCollapse que fusiona rotaciones adyacentes rxx y ryy en una única XXPlusYYGate. Ya aplicamos este pase cuando construimos los circuitos de Trotter en el Paso 1 (la llamada pm.run(...)). Esto reduce el recuento de compuertas de dos qubits y produce una estructura más apta para la simulación de redes tensoriales para el cálculo de coeficientes dinámicos más adelante.
Reducción a la ISA del simulador. A continuación, ejecutamos el gestor de pases preestablecido de Qiskit con optimization_level=3 para reducir cada circuito de Trotter a la arquitectura de 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 primitivas de Qiskit

Para el ejemplo a pequeña escala, ejecutamos los circuitos de Trotter reducidos a la ISA a través de la primitiva EstimatorV2 respaldada por Aer. Al hacerlo, obtenemos un valor de referencia sin ruido para cada par $(k_j, t)$ — estos son los valores $\langle A \rangle_{k_j}(t)$ que la MPF combinará en el Paso 4. Recorremos varios tiempos de evolución para poder graficar más adelante la serie temporal completa 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: Post-procesar y devolver el resultado en el formato clásico deseado

El Paso 4 es donde se construye realmente la MPF. Aunque los coeficientes $x_j$ se calculan aquí (y, para la variante dinámica, este cálculo puede ser intensivo), conceptualmente son una receta clásica para combinar las mediciones cuánticas del Paso 3 en un único valor esperado corregido — por eso tratamos todo el flujo de trabajo de coeficientes y combinación como posprocesamiento.

Para evaluar qué tan bien la MPF sigue la dinámica verdadera, primero calculamos los valores esperados exactos evolucionados en el tiempo exponenciando directamente el hamiltoniano. Esto solo es factible porque $L = 10$ ; en el ejemplo de hardware a gran escala más adelante tendremos que depender de estimaciones de redes tensoriales en su lugar.

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

Las MPF estáticas usan coeficientes $x_j$ independientes del tiempo de evolución, el hamiltoniano y el estado inicial. Planteamos el sistema lineal $Ax = b$ descrito en la sección Antecedentes y resolvemos los coeficientes. La matriz $A$ está determinada por los recuentos de pasos de Trotter $k_j$ , el orden $\chi$ de la fórmula de producto, y si la fórmula es simétrica (lo que controla los exponentes $\eta_n$ ).

Para nuestro ejemplo a pequeña escala usamos $k_j = [1, 2, 4]$ con una fórmula de Suzuki-Trotter no simétrica de orden $2\chi=2$ (por lo que $\chi=1$ y $\eta_n = 2 + n$ , dando $\eta_0 = 2,\, \eta_1 = 3$ ). El sistema queda:

A = \begin{bmatrix} 1 & 1 & 1\\ 1 & \frac{1}{2^2} & \frac{1}{4^2} \\ 1 & \frac{1}{2^3} & \frac{1}{4^3} \\ \end{bmatrix}, \quad b = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}.

La primera fila impone la insesgadez ( $\sum_j x_j = 1$ ); la segunda y tercera filas cancelan los términos principales de error de Trotter $1/k^2$ y $1/k^3$ de siguiente orden, respectivamente.

Configurar el LSE

Usamos setup_static_lse de qiskit_addon_mpf.static para ensamblar 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 nuestra elección de fórmula de producto — en particular su orden $\chi$ y si es simétrica. El indicador symmetric controla el patrón de exponentes $\eta_n$ (las fórmulas simétricas solo producen términos de error de Trotter de potencia par; consulta la Ref. [1]). Ten en cuenta que, como se muestra en la Ref. [2], establecer symmetric=True no es estrictamente necesario incluso cuando la PF subyacente es simétrica — el LSE no simétrico sigue siendo válido (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)

Inspecciona la matriz $A$ y el vector $b$ construidos para confirmar que coinciden con el sistema escrito anteriormente.

lse.A

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

lse.b

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

Con el LSE listo, 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}"
)

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

Optimizar para $x$ utilizando un modelo exacto

Como alternativa a calcular $x = A^{-1}b$ , puedes usar setup_exact_model para construir una instancia de cvxpy.Problem que utiliza el LSE como restricciones y cuya solución óptima produzca $x$ .

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

[ 0.04761905 -0.57142857  1.52380952]

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

L1 norm of the exact coefficients: 2.1428571428556378

Optimizar para $x$ utilizando un modelo aproximado

Puede ocurrir que la norma $L_1$ para el conjunto elegido de valores $k_j$ se considere demasiado alta. Si ese es el caso y no puedes elegir un conjunto diferente de valores $k_j$ , puedes usar una solución aproximada que limite la norma $L_1$ a un umbral elegido mientras minimiza $\|Ax - b\|$ . Consulta la guía sobre Cómo usar 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),
)

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

Coeficientes MPF dinámicos

La MPF estática cancela los términos de error de Trotter de una manera agnóstica al hamiltoniano y al estado, por lo que no necesariamente produce el error de aproximación más pequeño posible para un hamiltoniano y estado inicial dados. La MPF dinámica (Refs. [2], [3]) en cambio encuentra 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 tiempo $t$ . Como se muestra en la sección Antecedentes, esto requiere la matriz de solapamiento $M_{ij}(t)$ entre estados evolucionados por Trotter y el solapamiento $L_i(t)$ con el estado exacto — ambos los estimamos usando backends de redes tensoriales (TeNPy) en qiskit_addon_mpf.

Para configurar el LSE dinámico necesitamos tres ingredientes:

Una fábrica de evolutor aproximado que el addon ejecutará para cada $k_j$ para producir $\rho_{k_j}(t)$ como un MPS/MPO. La construimos a partir de la estructura por capas del circuito de Trotter de orden $2$ (una capa por slice_by_depth), envuelta como un LayerwiseEvolver con parámetros de truncamiento de TeNPy.
Una fábrica de evolutor exacto que produce una referencia de alta precisión $\rho(t)$ . Usamos un circuito de Suzuki-Trotter de cuarto orden con paso temporal pequeño (dt=0.1, order=4) como sustituto de la evolución exacta.
Una fábrica de identidad y un MPS de estado inicial que inicializan la simulación TeNPy.

La celda a continuación construye la fábrica de evolutor aproximado.

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

aviso

Las opciones de LayerwiseEvolver que determinan los detalles de la simulación de la red tensorial deben elegirse con cuidado para evitar configurar un problema de optimización mal definido.

Aproximamos el estado exacto evolucionado en el tiempo con una fórmula de Suzuki-Trotter de cuarto orden usando un paso de tiempo pequeño dt=0.1. Los parámetros de truncamiento de TeNPy pueden afectar a la precisión, por lo que es importante explorar un rango de 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,
    },
)

Finalmente, definimos una identity_factory que produce el estado MPO inicial y preparamos el estado inicial de Néel como un MPS que coincide con la red usada por el modelo de Trotter por capas.

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

mps_initial_state = MPS_neel_state(models[0].lat)

Con las fábricas listas, ahora calculamos los coeficientes dinámicos en cada tiempo de evolución. Para cada $t$ , setup_dynamic_lse construye las matrices de solapamiento relevantes mediante TeNPy, y setup_frobenius_problem devuelve un cvxpy.Problem que minimiza el costo de la norma de Frobenius. El solucionador devuelve coeficientes $x_j(t)$ adaptados a ese tiempo; los recopilamos en mpf_dynamic_coeffs_list. Si el solucionador falla para un $t$ dado, recurrimos a coeficientes cero 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("")

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 MPF

Ahora evaluamos $\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 graficamos la serie temporal resultante frente a 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 of the previous code cell

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

Error de Trotter. Las fórmulas de producto individuales (marcadores grises) se desvían cada vez más de la curva exacta a medida que aumenta el tiempo. El circuito $k=1$ tiene la mayor desviación y es el más superficial, pero también ya está en el régimen donde $t/k \gtrsim 1$ , por lo que el término de error principal $1/k^{2}$ es grande. Las combinaciones MPF (marcadores coloreados) cancelan varios de estos términos principales de error de Trotter, por lo que siguen la curva exacta mucho más de cerca que cualquier circuito individual $k_j$ . La brecha restante refleja los términos de Trotter de orden superior que la MPF no cancela: una MPF estática de orden $2$ , $r=3$ solo elimina los dos primeros órdenes de error, y a $t/k_{\min}$ grande la cola no cancelada eventualmente domina — por lo que la MPF no garantiza que los circuitos muy superficiales permanezcan precisos en tiempos arbitrarios.
Error de muestreo. Las barras de error más anchas en las curvas MPF son una consecuencia directa de la combinación lineal: propagar errores estándar independientes por circuito $\sigma_{k_j}$ da una varianza total $\sigma_{\text{MPF}}^2 = \sum_j x_j^2 \, \sigma_{k_j}^2$ . Por lo tanto, cuanto mayor sea $\|x\|_2$ (y en la práctica $\|x\|_1$ , que es lo que controlamos), más disparos se requieren para alcanzar una incertidumbre objetivo dada. Este es el compromiso detrás de la opción del solucionador aproximado en la sección Antecedentes: limitamos $\|x\|_1$ para mantener manejable esta sobrecarga. Fundamentalmente, a diferencia del error de Trotter, el error de muestreo se reduce con $1/\sqrt{N_{\text{shots}}}$ , por lo que siempre se puede reducir invirtiendo más disparos.

En el ejemplo de hardware a gran escala a continuación, el ruido del hardware entra como una fuente de error adicional en cada $\langle A \rangle_{k_j}$ , que de manera similar se amplifica por los coeficientes MPF. Veremos cómo interactúa la mitigación de errores con las MPF en esa sección.

Ejemplo de hardware a gran escala

En esta sección escalamos el problema más allá de lo que es posible simular exactamente. Reproducimos algunos de los resultados mostrados en la Ref. [3], usando una cadena XXZ de 50 qubits en el tiempo $t = 3$ . Seguimos el mismo flujo de trabajo de cuatro pasos que el ejemplo a pequeña escala, ahora apuntando a hardware cuántico real con mitigación de errores. Como en la plantilla, cada paso está marcado en línea en el código, y un solo paso puede abarcar varias celdas cuando sus resultados intermedios valen la pena inspeccionar. El mapeo refleja el ejemplo a pequeña escala: definir un hamiltoniano, elegir parámetros de Trotter, calcular coeficientes MPF (estáticos y dinámicos), y construir circuitos. Las diferencias clave 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 de segundo orden simétrica 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 $t/k_{\min}=1$ , manteniendo los componentes superficiales dentro del régimen de convergencia de Trotter donde el modelo de error principal en el que se basa la MPF es válido.
Una ejecución de comparación de circuito único adicional con $k = 10$ pasos de Trotter, usada como línea base. Elegimos $k = 10$ porque su profundidad de dos qubits en hardware es mayor que la del componente MPF más profundo ( $k_{\max}=6$ ) más la sobrecarga de ejecutar múltiples circuitos MPF — suficientemente profundo para estar limitado por el ruido, que es el régimen en el que se espera que la combinación MPF supere a la línea base de circuito único. Es una comparación de "circuito único profundo" contra la combinación MPF, no un circuito que apunte al error de Trotter efectivo de la MPF (lo que requeriría muchos más pasos).

Ten en cuenta que, aunque todavía estamos en el Paso 1 aquí (mapeo y construcción de circuitos), también precalculamos los coeficientes dinámicos junto con los estáticos en esta celda. 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 mantener toda la configuración específica de 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)

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. Usamos el gestor de pases preestablecido de Qiskit con optimization_level=3, que selecciona automáticamente un buen conjunto de qubits físicos y enruta 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
]

<IBMBackend('ibm_fez')>

Ejecutar circuitos más profundos en hardware real requiere una mitigación de errores agresiva. Habilitamos el desacoplamiento dinámico, el twirling de compuertas y mediciones, la mitigación de errores de medición, y la extrapolación a ruido cero (ZNE). Ten en cuenta que los factores de ruido de ZNE que usamos aquí (1, 1.2, 1.4) son menores que en un escenario de circuito superficial, ya que los componentes MPF más profundos ya están cerca del umbral de ruido y grandes amplificaciones de ruido los llevarían más allá del punto en el que la extrapolación ZNE es confiable.

Enviamos los cuatro circuitos (tres componentes MPF en $k_j = [3, 4, 6]$ más la línea base $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 del resultado del trabajo, luego los combinamos con cada conjunto de coeficientes MPF exactamente como en el ejemplo a pequeña escala: $\langle A \rangle_{\text{MPF}} = \sum_j x_j \, \langle A \rangle_{k_j}$ , con 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)

[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,
)

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 of the previous code cell

Algunas observaciones sobre los resultados de hardware anteriores:

Ir más profundo no es gratis en hardware. Las líneas base de circuito único cuentan la historia directamente: el circuito $k = 6$ es esencialmente exacto ( $-0.256$ frente a la referencia $-0.244$ ), pero la línea base más profunda $k = 10$ es peor ( $-0.061$ , desviada en $\sim 0.18$ ), no mejor. Una vez que el error de Trotter ya es pequeño, agregar pasos principalmente profundiza el circuito y acumula más ruido de compuertas y decoherencia. Este es precisamente el régimen para el que están construidas las MPF: alcanzar la precisión de un circuito profundo usando solo componentes superficiales.
Una MPF de norma pequeña supera al circuito único profundo. La MPF estática-aproximada (limitada a $\|x\|_1 \approx 2$ ) llega a $-0.259$ , dentro de $\sim 0.015$ de la referencia y mucho más cerca que la línea base $k = 10$ . La MPF dinámica ( $-0.127$ ) también supera cómodamente esa línea base. Ambas combinan solo los circuitos superficiales $k_j = [3, 4, 6]$ , y sin embargo recuperan una respuesta que el circuito único profundo no pudo lograr.
La norma de los coeficientes importa más que la optimalidad matemática. La MPF estática-exacta tiene $\|x\|_1 = 4.66$ y es el peor estimador de todos ( $-0.567$ , desviada en más de $0.3$ ): la gran norma de coeficientes amplifica el ruido residual de compuertas, la decoherencia y el error de ZNE en cada $\langle A \rangle_{k_j}$ aproximadamente por el mismo factor, superando la cancelación de error de Trotter que consigue. Limitar la norma (el solucionador estático-aproximado, $\|x\|_1 \approx 2$ ) elimina esta sobrecarga y da la mejor estimación — aunque sus coeficientes ya no cancelan exactamente el error de Trotter principal.
Los circuitos superficiales individuales pueden seguir siendo competitivos. El único componente $k = 6$ ( $-0.256$ ) es en sí mismo esencialmente exacto aquí — en esta ejecución incluso está marginalmente más cerca que la MPF estática-aproximada. El problema es que no sabes de antemano cuál $k$ individual se encuentra en el punto óptimo de "convergido pero aún no limitado por el ruido", y la elección aparentemente segura de simplemente ir más profundo ( $k = 10$ ) para garantizar la convergencia de Trotter es exactamente la que falla. La MPF ofrece una combinación fundamentada de circuitos superficiales que no requiere adivinar la profundidad correcta.

La conclusión práctica es que en hardware, las MPF deben combinarse con una mitigación de errores sólida en cada $\langle A \rangle_{k_j}$ individual, la norma $L_1$ de los coeficientes debe mantenerse modesta (usa el solucionador aproximado, o la MPF dinámica), y los pasos de Trotter $k_j$ deben elegirse de modo que $t/k_{\min} \lesssim 1$ — aquí $k_{\min} = 3$ en $t = 3$ da $t/k_{\min} = 1$ , manteniendo los componentes dentro del régimen convergente donde el modelo de error principal en el que se basa la MPF estática es válido. Con esas elecciones, las MPF de norma pequeña aquí igualan a un circuito único convergido mientras que la línea base ingenua de "simplemente ir más profundo" no lo hace, recuperando la ventaja de profundidad frente a precisión mostrada en la Ref. [3]. Ten en cuenta también que las ejecuciones individuales tienen ruido — en un envío diferente del mismo trabajo (o un backend diferente), el orden exacto puede cambiar; las tendencias robustas son que las MPF con $\|x\|_1$ pequeño funcionan bien, la MPF estática-exacta con $\|x\|_1$ grande se amplifica por el ruido del hardware, y el circuito único demasiado profundo está limitado por el ruido.

Próximos pasos

Recomendaciones

Si te pareció interesante este trabajo, podrían interesarte los siguientes materiales:

Cómo elegir los pasos de Trotter para una MPF — orientación práctica sobre cómo seleccionar valores de $k_j$ para evitar inestabilidades
Cómo usar el modelo aproximado — ajuste de la restricción de norma $L_1$ y opciones del solucionador para la MPF estática aproximada
Referencia de API de qiskit-addon-mpf — documentación completa de los módulos estático, dinámico y de backend

Referencias

[1] Vazquez, A. C., Egger, D. J., Ochsner, D., & Woerner, S. Well-conditioned multi-product formulas for hardware-friendly Hamiltonian simulation. Quantum, 7, 1067 (2023)

[2] Zhuk, S., Robertson, N. F., & Bravyi, S. Trotter error bounds and dynamic multi-product formulas for Hamiltonian simulation. Physical Review Research, 6(3), 033309 (2024)

[3] Robertson, N. F., et al. Tensor network enhanced dynamic multiproduct formulas. arXiv:2407.17405 (2024)

Resultados de aprendizaje​

Requisitos previos​

Antecedentes​

¿Qué son las fórmulas multiproducto?​

¿Cómo se determinan los coeficientes?​

Cuándo usar MPF​

Cuándo las MPF podrían no ayudar​

Qué cubre este tutorial​

Requisitos​

Configuración​

Ejemplo de simulador a pequeña escala​

Paso 1: Mapear las entradas clásicas a un problema cuántico​

Construir circuitos de Trotter​

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

Paso 3: Ejecutar utilizando primitivas de Qiskit​

Paso 4: Post-procesar y devolver el resultado en el formato clásico deseado​

Coeficientes MPF estáticos​

Configurar el LSE​

Optimizar para xxx utilizando un modelo exacto​

Optimizar para xxx utilizando un modelo aproximado​

Coeficientes MPF dinámicos​

Combinar los valores esperados de Trotter con los coeficientes MPF​

Ejemplo de hardware a gran escala​

Próximos pasos​

Referencias​