Resuelve el modelo de Sherrington-Kirkpatrick con el Parity Twine Optimizer de ParityQC

Estimación de uso: 10 segundos en un procesador Nighthawk r2. (NOTA: Esto es solo una estimación. Tu tiempo de ejecución puede variar.)

Resultados de aprendizaje

• Usar el Parity Twine Optimizer para resolver el modelo de Sherrington-Kirkpatrick.

• Aprender qué opciones están disponibles para el Parity Twine Optimizer y qué resultados se generan.

Contexto

Este tutorial muestra cómo resolver el modelo de Sherrington-Kirkpatrick usando el Parity Twine Optimizer de ParityQC.

Proporciona código para formular el problema de manera totalmente local, en un formato compatible con el Parity Twine Optimizer.

El modelo de Sherrington-Kirkpatrick

El modelo de Sherrington-Kirkpatrick (SK) es un modelo fundamental en mecánica estadística, específicamente dentro del estudio de los vidrios de espín. A diferencia del modelo de Ising estándar, donde las interacciones suelen limitarse a los vecinos más cercanos, el modelo SK es un modelo de alcance infinito, lo que significa que cada espín interactúa con todos los demás espines del sistema. Esto da lugar a un paisaje energético muy complejo y "accidentado", caracterizado por muchos mínimos locales, lo cual es el sello distintivo del comportamiento vítreo.

La característica central del modelo SK reside en la frustración. En el modelo, las intensidades de interacción $J_{ij}$ entre espines se distribuyen aleatoriamente entre valores positivos y negativos. Esto da lugar a situaciones (por ejemplo, en disposiciones triangulares) en las que los espines no se pueden organizar para minimizar todas las interacciones simultáneamente. En el modelo SK, dado que cada espín interactúa con todos los demás, esta frustración se agrava globalmente, dando lugar a una red de restricciones conflictivas.

Formulación matemática

El estado del sistema se define mediante un conjunto de $N$ espines de Ising, $s_i \in \{+1, -1 \}$. La energía de una configuración específica viene dada por el Hamiltoniano:

$H = - \sum_{1 \leq i \leq j \leq N} J_{ij} s_i s_j$

donde $J_{ij}$ es la fuerza de acoplamiento entre el espín $i$ y el espín $j$.

En el modelo SK, los acoplamientos $J_{ij}$ son variables aleatorias independientes e idénticamente distribuidas. Para garantizar que la energía siga siendo extensiva (proporcional a $N$) a medida que $N \rightarrow \infty$, la varianza de los acoplamientos debe escalar con el número de partículas:

$J_{ij} \sim \mathcal{N} \left( 0, \frac{J^2}{N} \right)$

El estado fundamental para un $H$ dado es la configuración específica de espines ($s_1, s_2, ...s_n$) que minimiza la energía. Encontrar el estado fundamental es un problema de optimización NP-difícil.

Requisitos

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

• Qiskit Functions Catalog IBM Client (pip install qiskit-ibm-catalog)

• Qiskit addon Optimization Mapper (pip install qiskit_addon_opt_mapper)

• NumPy (pip install numpy)

También necesitas permiso para acceder a la función ParityQC Twine Optimizer. Para solicitar acceso, completa este formulario.

Configuración

(Este código asume que ya has guardado tu cuenta en tu entorno local.)

Primero, importa todos los paquetes necesarios para este tutorial.

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

_needed = {"numpy": "numpy", "qiskit_ibm_catalog": "qiskit-ibm-catalog"}
_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 numpy as np

from qiskit_ibm_catalog import QiskitFunctionsCatalog

Carga el Parity Twine Optimizer desde el catálogo de Qiskit Functions:

catalog = QiskitFunctionsCatalog(channel="ibm_quantum_platform")
function = catalog.load("parityqc/parity-twine-optimizer")

Paso 1: Definir el problema como una función objetivo

En lugar de obtener el problema SK de una biblioteca como hacemos para el problema Market Split, lo formulas directamente.

La función generate_sk_problem formula el problema SK directamente en el formato de diccionario requerido. La única entrada necesaria es n, el número de espines en el modelo.

def generate_sk_problem(
    n: int,
    coupling_mean: float = 0.0,
    coupling_std: float = 1.0,
    local_fields_mean: float = 0.0,
    local_fields_std: float = 0.0,
    edge_density: float = 1.0,
    ensure_extensivity: bool = False,
    seed: int | None = None,
) -> dict:
    """Generate the Sherrington-Kirkpatrick (SK) model with varying
     edge density.

    Samples couplings and local fields via ``generate_couplings_sk_model``
    and assembles the corresponding Ising Hamiltonian

        H = -∑_{i<j} J_ij z_i z_j - ∑_i h_i z_i,

    where z_i ∈ {-1, +1}.

    Args:
        n: Number of spins (>= 2).
        coupling_mean: Mean coupling before optional SK scaling.
        coupling_std: Coupling std before optional SK scaling.
        local_fields_mean: Mean longitudinal field.
        local_fields_std: Std of the longitudinal fields.
        edge_density: Fraction of non-zero couplings, in ``[2/n, 1]``.
        ensure_extensivity: Whether to apply the SK 1/n scaling.
        seed: random number generator seed.

    Returns:
        An encoding the SK Hamiltonian.

    Raises:
        ValueError: If ``n < 2``, ``coupling_std < 0``, ``local_fields_std < 0``,
            or ``edge_density`` is outside ``[2/n, 1]``.
    """
    couplings, local_fields = _generate_couplings_sk_model(
        n=n,
        coupling_mean=coupling_mean,
        coupling_std=coupling_std,
        local_fields_mean=local_fields_mean,
        local_fields_std=local_fields_std,
        edge_density=edge_density,
        ensure_extensivity=ensure_extensivity,
        seed=seed,
    )

    # Handle quadratic terms: coupling[i, j] * zj[i] * zj[j]
    # Only iterate over the upper triangle (i < j)
    sk_problem = {
        str((i, j)): float(couplings[i, j])
        for i in range(n)
        for j in range(i + 1, n)
        if couplings[i, j] != 0
    }

    # Handle linear terms: local_fields[i] * zj[i]
    sk_problem.update(
        {
            str((i,)): float(local_fields[i])
            for i in range(n)
            if local_fields[i] != 0
        }
    )

    return sk_problem

Usa la función _generate_couplings_sk_model para calcular los términos de acoplamiento aleatorios en el modelo SK. Para tener más control sobre los acoplamientos, puedes usar argumentos opcionales, que se explican en el docstring de la función.

def _generate_couplings_sk_model(
    n: int,
    coupling_mean: float = 0.0,
    coupling_std: float = 1.0,
    local_fields_mean: float = 0.0,
    local_fields_std: float = 0.0,
    edge_density: float = 1.0,
    ensure_extensivity: bool = False,
    seed: int | None = None,
) -> tuple[np.ndarray, np.ndarray]:
    """Generate random couplings and local fields for an Ising / SK model.

    Couplings are Gaussian. With ``ensure_extensivity=True`` they follow the
    Sherrington-Kirkpatrick scaling ``J_ij ~ N(coupling_mean/n, coupling_std^2/n)``
    (extensive energy, O(n)); otherwise ``J_ij ~ N(coupling_mean, coupling_std^2)``
    (energy O(n^2)). Fields are ``h_i ~ N(local_fields_mean, local_fields_std^2)``.

    ``edge_density`` sets the fraction of the ``n*(n-1)/2`` possible couplings that
    are non-zero (1 = fully dense). The kept edges always include a random spanning
    tree, so the interaction graph is guaranteed connected. This requires at least
    ``n-1`` edges, so ``edge_density`` must be at least ``2/n``.

    Args:
        n: Number of spins (>= 2).
        coupling_mean: Mean coupling before optional SK scaling.
        coupling_std: Coupling std before optional SK scaling.
        local_fields_mean: Mean longitudinal field.
        local_fields_std: Std of the longitudinal fields.
        edge_density: Fraction of non-zero couplings, in ``[2/n, 1]``.
        ensure_extensivity: Whether to apply the SK 1/n scaling.
        seed: random number generator seed.

    Returns:
        Tuple ``(couplings, fields)``: a symmetric ``(n, n)`` matrix with zero
        diagonal, and an ``(n,)`` field vector.

    Raises:
        ValueError: If ``n < 2``, ``coupling_std < 0``, ``local_fields_std < 0``,
            or ``edge_density`` is outside ``[2/n, 1]``.
    """
    if n < 2:
        raise ValueError(f"n must be >= 2, got {n}")
    if coupling_std < 0 or local_fields_std < 0:
        raise ValueError(
            "coupling_std and local_fields_std must be non-negative"
        )

    # A connected graph on n nodes needs at least n-1 of the n*(n-1)/2 possible
    # edges, so edge_density has a hard lower bound of 2/n.
    min_edge_density = 2.0 / n
    if not min_edge_density <= edge_density <= 1.0:
        raise ValueError(
            f"edge_density must be in [{min_edge_density:.4g}, 1] for n={n} "
            f"(at least n-1 edges are needed to keep the graph connected), "
            f"got {edge_density}"
        )

    rng = np.random.default_rng(seed)

    j_loc, j_scale = (
        (coupling_mean / n, coupling_std / np.sqrt(n))
        if ensure_extensivity
        else (coupling_mean, coupling_std)
    )

    upper_idx = np.triu_indices(n, k=1)
    n_edges = len(upper_idx[0])

    # Select which edges are present.
    if edge_density < 1.0:
        n_keep = int(round(edge_density * n_edges))
        # Map each (i, j) node pair to its position in the flat upper-triangle list.
        pair_to_flat = {
            (int(i), int(j)): idx
            for idx, (i, j) in enumerate(
                zip(upper_idx[0], upper_idx[1], strict=False)
            )
        }

        # Random spanning tree: node perm[k] links to a random earlier node.
        perm = rng.permutation(n)
        keep = np.zeros(n_edges, dtype=bool)
        for k in range(1, n):
            child, parent = perm[k], perm[rng.integers(0, k)]
            i, j = min(child, parent), max(child, parent)
            keep[pair_to_flat[(int(i), int(j))]] = True

        # Fill the remaining budget with random non-tree edges.
        remaining = n_keep - (n - 1)
        if remaining > 0:
            keep[
                rng.choice(
                    np.flatnonzero(~keep), size=remaining, replace=False
                )
            ] = True
    else:
        keep = np.ones(n_edges, dtype=bool)

    n_present = int(keep.sum())
    if j_scale == 0.0:
        vals = np.full(n_present, j_loc)
    else:
        vals = rng.normal(loc=j_loc, scale=j_scale, size=n_present)

    couplings = np.zeros((n, n))
    couplings[upper_idx[0][keep], upper_idx[1][keep]] = vals
    couplings += couplings.T  # symmetrize; diagonal stays zero

    fields = (
        np.full(n, local_fields_mean)
        if local_fields_std == 0.0
        else rng.normal(loc=local_fields_mean, scale=local_fields_std, size=n)
    )

    return couplings, fields

Paso 2: Resolver el problema usando el Parity Twine Optimizer

Con las funciones anteriores, puedes configurar el problema SK y encontrar una solución usando el Twine Optimizer y un backend de IBM Quantum® elegido.

Para ejecutar la función, elige un backend adecuado; por ejemplo, ibm_phoenix.

Opcionalmente, usa options para tener control adicional sobre el envío:

options = {
    "shots": 100000,
    "postprocessing_level": 1,
    "transpile_only": False,
    "job_tags": ["sk"],
}

donde shots es un entero que especifica el número de ejecuciones del circuito, postprocessing_level determina si se aplica posprocesamiento al resultado, transpile_only elige si el problema solo se transpila a un circuito (y no se resuelve), y job_tags es la etiqueta para identificar el job en IBM Quantum Platform.

El tamaño del modelo SK se define mediante $N$, el número de espines que interactúan. Una vez que eliges $N$, el código anterior genera el problema para n_spins.

Ejecuta el optimizador:

n_spins = 50
sk_problem = generate_sk_problem(n_spins)

function_job = function.run(
    problem=sk_problem,
    variable_type="spin",
    backend_name="ibm_phoenix",
    options=options,
)
print(f"Job ID: {function_job.job_id}")

Comprueba el estado del job:

# Monitor the job status
function_job.status()

Recupera los resultados:

result = function_job.result()

result

El resultado tiene la forma:

{
    'solution': {'0': 1, '1': 1, '10': 1, '11': 1, ... },
    'objective_value':  -240.5425312543882,
    'solution_bitstring': '00001101110100100111001110101001101111111011001110',
    'metadata': {
        'circuit_metrics': {
            'depth': 523,
            'gate_count': 10118,
            'two_qubit_gate_depth': 196,
            'two_qubit_gate_count': 2499,
            'num_qubits': 50,
            'operations': {'sx': 3353, 'rz': 3320, 'cz': 2499, 'delay': 894, 'measure': 50, 'x': 2},
        },
        'solver_info': {
            'variable_mapping': {'0': 0, '1': 1, '10': 2, '11': 3, ... },
            'bitstring_distributions': {
                'before_postprocessing': {'011101110010110111001110011000': 1, ... },
                'after_postprocessing': {'011011110000110101001111011000': 1, ... }
            },
            'best_parameters': {
                'beta': [-0.4602084830507902],
                'gamma': [1.8500357096574955]
            }
        },
        'resource_usage': {
            'RUNNING: MAPPING': {'CPU_TIME': 290.272},
            'RUNNING: OPTIMIZING_FOR_HARDWARE': {'CPU_TIME': 0.494},
            'RUNNING: WAITING_FOR_QPU': {'CPU_TIME': 8.775},
            'RUNNING: EXECUTING_QPU': {'QPU_TIME': 31.0},
            'RUNNING: POST_PROCESSING': {'CPU_TIME': 162.96},
        },
    }
}

donde el diccionario solution corresponde a los qubits definidos en el problema y proporciona sus valores de espín optimizados para el Hamiltoniano del modelo SK. Esta secuencia específica de espines en la solución óptima representa la configuración que minimiza la energía total del sistema en función de las intensidades de interacción aleatorias dadas. En el modelo SK, esto puede considerarse como el estado de energía más bajo de un sistema magnético desordenado.

metadata proporciona información sobre la transpilación (recuentos/profundidad de puertas de dos qubits, puertas utilizadas, qubits activos) y varios tiempos de ejecución.

Próximos pasos

Recomendaciones

• Solicita acceso a la función completando este formulario.
• Visita la referencia de la API para esta Qiskit Function.
• Lee la guía.
• Prueba el tutorial para aplicar el Parity Twine Optimizer al problema Market Split.
• Revisa el preprint de ArXiv Connectivity-aware Synthesis of Quantum Algorithms, Drier et al. (2025).