Risolvi il modello di Sherrington-Kirkpatrick con il Parity Twine Optimizer di ParityQC

Stima di utilizzo: 10 secondi su un processore Nighthawk r2. (NOTA: questa è solo una stima. Il tuo tempo di esecuzione potrebbe variare.)

Obiettivi di apprendimento

• Usa il Parity Twine Optimizer per risolvere il modello di Sherrington-Kirkpatrick.

• Scopri quali opzioni sono disponibili per il Parity Twine Optimizer e quali risultati vengono restituiti.

Contesto

Questo tutorial dimostra come risolvere il modello di Sherrington-Kirkpatrick usando il Parity Twine Optimizer di ParityQC.

Fornisce il codice per formulare il problema interamente in locale in una forma che possa interfacciarsi con il Parity Twine Optimizer.

Il modello di Sherrington-Kirkpatrick

Il modello di Sherrington-Kirkpatrick (SK) è un modello fondamentale nella meccanica statistica, in particolare nello studio dei vetri di spin. A differenza del modello di Ising standard, dove le interazioni sono solitamente limitate ai vicini più prossimi, il modello SK è un modello a raggio infinito, il che significa che ogni spin interagisce con ogni altro spin del sistema. Ciò porta a un paesaggio energetico altamente complesso e "accidentato", caratterizzato da molti minimi locali, che è il tratto distintivo del comportamento vetroso.

La caratteristica principale del modello SK risiede nella frustrazione. Nel modello, le intensità di interazione $J_{ij}$ tra gli spin sono distribuite casualmente tra valori positivi e negativi. Ciò porta a situazioni (per esempio, in disposizioni triangolari) in cui gli spin non possono essere disposti in modo da minimizzare simultaneamente tutte le interazioni. Nel modello SK, poiché ogni spin interagisce con ogni altro spin, questa frustrazione si aggrava globalmente, dando origine a una rete di vincoli in conflitto.

Formulazione matematica

Lo stato del sistema è definito da un insieme di $N$ spin di Ising, $s_i \in \{+1, -1 \}$. L'energia di una specifica configurazione è data dall'Hamiltoniana:

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

dove $J_{ij}$ è l'intensità di accoppiamento tra lo spin $i$ e lo spin $j$.

Nel modello SK, gli accoppiamenti $J_{ij}$ sono variabili casuali indipendenti e identicamente distribuite. Per garantire che l'energia rimanga estensiva (proporzionale a $N$) quando $N \rightarrow \infty$, la varianza degli accoppiamenti deve scalare con il numero di particelle:

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

Lo stato fondamentale per un dato $H$ è la specifica configurazione di spin ($s_1, s_2, ...s_n$), che minimizza l'energia. Trovare lo stato fondamentale è un problema di ottimizzazione NP-hard.

Requisiti

Prima di iniziare questo tutorial, assicurati che siano installati i seguenti pacchetti:

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

• Qiskit addon Optimization Mapper (pip install qiskit_addon_opt_mapper)

• NumPy (pip install numpy)

Ti serve anche l'autorizzazione per accedere alla funzione ParityQC Twine Optimizer. Per richiedere l'accesso, completa questo modulo.

Configurazione

(Questo codice presuppone che tu abbia già salvato il tuo account nel tuo ambiente locale.)

Innanzitutto, importa tutti i pacchetti richiesti per questo 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

Carica il Parity Twine Optimizer dal catalogo delle Qiskit Functions:

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

Passaggio 1: definisci il problema come funzione obiettivo

Invece di ottenere il problema SK da una libreria come facciamo per il problema Market Split, lo formuli direttamente.

La funzione generate_sk_problem formula il problema SK direttamente nel formato dizionario richiesto. L'unico input richiesto è n, il numero di spin nel modello.

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 funzione _generate_couplings_sk_model per calcolare i termini di accoppiamento casuali nel modello SK. Per un maggiore controllo sugli accoppiamenti, puoi usare argomenti opzionali, spiegati nella docstring della funzione.

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

Passaggio 2: risolvi il problema usando il Parity Twine Optimizer

Con le funzioni precedenti, puoi impostare il problema SK e trovare una soluzione usando il Twine Optimizer e un backend IBM Quantum® a tua scelta.

Per eseguire la funzione, scegli un backend adatto; per esempio, ibm_phoenix.

Facoltativamente, usa options per un controllo aggiuntivo sull'invio del job:

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

dove shots è un intero che specifica il numero di esecuzioni del circuito, postprocessing_level determina se il post-processing viene applicato al risultato, transpile_only sceglie se il problema viene solo transpilato in un circuito (e non risolto), e job_tags è l'etichetta per identificare il job sulla IBM Quantum Platform.

La dimensione del modello SK è definita da $N$, il numero di spin interagenti. Una volta scelto $N$, il codice precedente genera il problema per n_spins.

Esegui l'ottimizzatore:

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

Controlla lo stato del job:

# Monitor the job status
function_job.status()

Recupera i risultati:

result = function_job.result()

result

Il risultato è nella 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},
        },
    }
}

dove il dizionario solution corrisponde ai qubit definiti nel problema e fornisce i loro valori di spin ottimizzati per l'Hamiltoniana del modello SK. Questa specifica sequenza di spin nella soluzione ottimale rappresenta la configurazione che minimizza l'energia totale del sistema in base alle intensità di interazione casuali date. Nel modello SK, questa può essere considerata come lo stato a energia più bassa di un sistema magnetico disordinato.

metadata fornisce informazioni sulla transpilazione (conteggio/profondità delle gate a due qubit, gate usate, qubit attivi) e vari tempi di esecuzione.

Passi successivi

Consigli

• Richiedi l'accesso alla funzione completando questo modulo.
• Visita il riferimento API per questa Qiskit Function.
• Leggi la guida.
• Prova il tutorial per applicare il Parity Twine Optimizer al problema Market Split.
• Consulta il preprint ArXiv Connectivity-aware Synthesis of Quantum Algorithms, Drier et al. (2025).