Skip to main content
IBM Quantum Platform

Resuelve el modelo de Sherrington-Kirkpatrick con el optimizador «Parity Twine» de ParityQC

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


Resultados del aprendizaje

  • Utiliza el optimizador Parity Twine para resolver el modelo de Sherrington-Kirkpatrick.
  • Descubre qué opciones están disponibles en el Parity Twine Optimizer y qué resultados se obtienen.

En segundo plano

Este tutorial muestra cómo resolver el modelo de Sherrington-Kirkpatrick utilizando el optimizador «Parity Twine» de ParityQC.

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

El modelo de Sherrington-Kirkpatrick

El modelo de Sherrington-Kirkpatrick (SK) es un modelo fundamental de la mecánica estadística, concretamente en el ámbito del estudio de los vidrios de espín. A diferencia del modelo de Ising estándar, en el que 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 numerosos mínimos locales, lo cual es el rasgo distintivo del comportamiento vítreo.

La característica fundamental del modelo SK radica en la frustración. En el modelo, las intensidades de interacción JijJ_ij entre los espines se distribuyen aleatoriamente entre valores positivos y negativos. Esto da lugar a situaciones (por ejemplo, en disposiciones triangulares) en las que no es posible organizar los espines de forma que se minimicen todas las interacciones al mismo tiempo. En el modelo SK, dado que cada espín interactúa con todos los demás espines, esta frustración se agrava a nivel global, lo que da lugar a una red de restricciones contradictorias.

Formulación matemática

El estado del sistema viene definido por un conjunto de NN espines de Ising, si{+1,1}s_i \in \{+1, -1 \} . La energía de una configuración concreta viene dada por el hamiltoniano:

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

donde JijJ_{ij} es la intensidad del acoplamiento entre el espín ii y el espín jj.

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

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

El estado fundamental de un sistema dado HH es la configuración específica de los espines ( s1,s2,...sns_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 que tengas instalados los siguientes elementos:

  • Qiskit Functions Catalog IBM Cliente (pip install qiskit-ibm-catalog)
  • Complemento de Qiskit «Optimization Mapper» (pip install qiskit_addon_opt_mapper)
  • NumPy (pip install numpy)

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


Configuración

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

En primer lugar, importa todos los paquetes necesarios para este tutorial.

import numpy as np

from qiskit_ibm_catalog import QiskitFunctionsCatalog

Carga el «Parity Twine Optimizer» del catálogo « 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 con el problema «Market Split», lo formulas directamente.

La función generate_sk_problem formula el problema SK directamente en el formato de diccionario requerido. El único dato nque hay que introducir es el número de giros del 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 :func:`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:
        A ``ProblemRepresentation`` 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

Utiliza la función _generate_couplings_sk_model para calcular los términos de acoplamiento aleatorios en el modelo SK. Para tener un mayor control sobre los acoplamientos, puedes utilizar argumentos opcionales, que se explican en la cadena de documentación 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: Resuelve el problema utilizando el «Parity Twine Optimizer»

Con las funciones anteriores, puedes plantear el problema SK y encontrar una solución utilizando el optimizador de Twine y el backend de IBM Quantum® que elijas.

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

Si lo deseas, puedes utilizar las opciones para tener un mayor control sobre el envío:

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

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

El tamaño del modelo SK viene definido por NN, el número de espines que interactúan. Una vez que elijas « NN », 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 trabajo:

# Monitor the job status
function_job.status()

Obtener resultados:

result = function_job.result()

result

El resultado tiene la siguiente 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 se corresponde con los qubits definidos en el problema y proporciona sus valores de espín optimizados para el hamiltoniano del modelo SK. Esta secuencia concreta de giros 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 menor energía de un sistema magnético desordenado.

metadata ofrece información sobre la transpilación (número de puertas de dos qubits/profundidad, puertas utilizadas, qubits activos) y diversos tiempos de ejecución.


Próximos pasos

Recomendaciones
¿Le ha resultado útil esta página?
Informe de un error, de una errata o solicite contenido en GitHub.