Résoudre le modèle de Sherrington-Kirkpatrick avec le ParityQC Parity Twine Optimizer

Estimation d'utilisation : 10 secondes sur un processeur Nighthawk r2. (REMARQUE : il ne s'agit que d'une estimation. Ton temps d'exécution peut varier.)

Objectifs d'apprentissage

• Utiliser le Parity Twine Optimizer pour résoudre le modèle de Sherrington-Kirkpatrick.

• Découvrir quelles options sont disponibles pour le Parity Twine Optimizer et quels résultats sont produits en sortie.

Contexte

Ce tutoriel montre comment résoudre le modèle de Sherrington-Kirkpatrick à l'aide du ParityQC Parity Twine Optimizer.

Il fournit du code permettant de formuler le problème entièrement en local, dans une forme compatible avec le Parity Twine Optimizer.

Le modèle de Sherrington-Kirkpatrick

Le modèle de Sherrington-Kirkpatrick (SK) est un modèle fondamental en mécanique statistique, en particulier dans l'étude des verres de spin. Contrairement au modèle d'Ising standard, où les interactions sont généralement limitées aux plus proches voisins, le modèle SK est un modèle à portée infinie, ce qui signifie que chaque spin interagit avec tous les autres spins du système. Cela conduit à un paysage énergétique très complexe et « accidenté », caractérisé par de nombreux minima locaux, ce qui est la marque caractéristique du comportement vitreux.

La caractéristique centrale du modèle SK réside dans la frustration. Dans le modèle, les forces d'interaction $J_{ij}$ entre les spins sont réparties aléatoirement entre des valeurs positives et négatives. Cela conduit à des situations (par exemple, dans des arrangements triangulaires) où les spins ne peuvent pas être disposés de manière à minimiser toutes les interactions simultanément. Dans le modèle SK, comme chaque spin interagit avec tous les autres spins, cette frustration s'amplifie à l'échelle du système entier, donnant lieu à un réseau de contraintes contradictoires.

Formulation mathématique

L'état du système est défini par un ensemble de $N$ spins d'Ising, $s_i \in \{+1, -1 \}$. L'énergie d'une configuration donnée est donnée par l'hamiltonien :

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

où $J_{ij}$ est la force de couplage entre le spin $i$ et le spin $j$.

Dans le modèle SK, les couplages $J_{ij}$ sont des variables aléatoires indépendantes et identiquement distribuées. Pour garantir que l'énergie reste extensive (proportionnelle à $N$) lorsque $N \rightarrow \infty$, la variance des couplages doit dépendre du nombre de particules :

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

L'état fondamental pour un $H$ donné est la configuration spécifique des spins ($s_1, s_2, ...s_n$) qui minimise l'énergie. Trouver l'état fondamental est un problème d'optimisation NP-difficile.

Prérequis

Avant de commencer ce tutoriel, assure-toi que les éléments suivants sont installés :

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

• Qiskit addon Optimization Mapper (pip install qiskit_addon_opt_mapper)

• NumPy (pip install numpy)

Tu as également besoin d'une autorisation pour accéder à la fonction ParityQC Twine Optimizer. Pour demander l'accès, remplis ce formulaire.

Configuration

(Ce code suppose que tu as déjà enregistré ton compte dans ton environnement local.)

Commence par importer tous les paquets requis pour ce tutoriel.

# 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

Charge le Parity Twine Optimizer depuis le catalogue Qiskit Functions :

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

Étape 1 : définir le problème comme une fonction objectif

Au lieu d'obtenir le problème SK à partir d'une bibliothèque comme nous le faisons pour le problème Market Split, tu le formules directement.

La fonction generate_sk_problem formule le problème SK directement dans le format de dictionnaire requis. La seule entrée requise est n, le nombre de spins dans le modèle.

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

Utilise la fonction _generate_couplings_sk_model pour calculer les termes de couplage aléatoires du modèle SK. Pour un contrôle plus fin des couplages, tu peux utiliser des arguments facultatifs, expliqués dans la docstring de la fonction.

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

Étape 2 : résoudre le problème avec le Parity Twine Optimizer

Avec les fonctions ci-dessus, tu peux configurer le problème SK et trouver une solution en utilisant le Twine Optimizer et un backend IBM Quantum® de ton choix.

Pour exécuter la fonction, choisis un backend adapté ; par exemple, ibm_phoenix.

Tu peux aussi utiliser options pour mieux contrôler la soumission :

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

où shots est un entier qui spécifie le nombre d'exécutions du circuit, postprocessing_level détermine si un post-traitement est appliqué au résultat, transpile_only indique si le problème est seulement transpilé en circuit (sans être résolu), et job_tags est l'étiquette utilisée pour identifier le job sur IBM Quantum Platform.

La taille du modèle SK est définie par $N$, le nombre de spins en interaction. Une fois que tu choisis $N$, le code ci-dessus génère le problème pour n_spins.

Exécute l'optimiseur :

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

Vérifie le statut du job :

# Monitor the job status
function_job.status()

Récupère les résultats :

result = function_job.result()

result

Le résultat est de la forme :

{
    '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},
        },
    }
}

où le dictionnaire solution correspond aux qubits définis dans le problème et donne leurs valeurs de spin optimisées pour l'hamiltonien du modèle SK. Cette séquence spécifique de spins dans la solution optimale représente la configuration qui minimise l'énergie totale du système en fonction des forces d'interaction aléatoires données. Dans le modèle SK, cela peut être considéré comme l'état de plus basse énergie d'un système magnétique désordonné.

metadata fournit des informations sur la transpilation (nombre/profondeur des portes à deux qubits, portes utilisées, qubits actifs) et divers temps d'exécution.

Étapes suivantes

Recommandations

• Demande l'accès à la fonction en remplissant ce formulaire.
• Consulte la référence API pour cette Qiskit Function.
• Lis le guide.
• Essaie le tutoriel pour appliquer le Parity Twine Optimizer au problème Market Split.
• Consulte le préprint ArXiv Connectivity-aware Synthesis of Quantum Algorithms, Drier et al. (2025).