ParityQC Parity Twine Optimizer を使用してSherrington-Kirkpatrickモデルを解く

使用量の目安: Nighthawk r2 プロセッサーで10秒。(注: これはあくまで目安であり、実際の実行時間は変動する場合があります。)

学習成果

Parity Twine Optimizer を使用してSherrington-Kirkpatrickモデルを解く。
Parity Twine Optimizer で利用可能なオプションと、出力される結果について学ぶ。

背景

このチュートリアルでは、ParityQC Parity Twine Optimizer を使用してSherrington-Kirkpatrickモデルを解く方法を示します。

このチュートリアルでは、Parity Twine Optimizer とインターフェースできる形式で、問題を完全にローカルで定式化するコードを提供します。

Sherrington-Kirkpatrickモデル

Sherrington-Kirkpatrick(SK)モデルは、統計力学、特にスピングラスの研究における基礎的なモデルです。相互作用が通常最近接に限定される標準的なイジングモデルとは異なり、SKモデルは無限レンジモデルであり、システム内のすべてのスピンが他のすべてのスピンと相互作用することを意味します。これにより、多数の局所最小値を特徴とする非常に複雑で「起伏の激しい」エネルギー地形が生じ、これはガラス的挙動の特徴です。

SKモデルの中核となる特徴はフラストレーションにあります。このモデルでは、スピン間の相互作用強度 $J_{ij}$ が正と負の値の間でランダムに分布します。これにより、 (例えば三角形配置において)すべての相互作用を同時に最小化するようにスピンを配置できない状況が生じます。SKモデルでは、すべてのスピンが他のすべてのスピンと相互作用するため、このフラストレーションは系全体で積み重なり、互いに矛盾する制約が複雑に絡み合った状態になります。

数学的定式化

システムの状態は、 $N$ 個のイジングスピン $s_i \in \{+1, -1 \}$ の集合によって定義されます。特定の配置のエネルギーは、次のハミルトニアンによって与えられます。

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

ここで $J_{ij}$ はスピン $i$ とスピン $j$ 間の結合強度です。

SKモデルでは、結合 $J_{ij}$ は独立同分布の確率変数です。 $N \rightarrow \infty$ のときにエネルギーが示量的(つまり $N$ に比例する)であり続けるようにするためには、結合の分散は粒子数に応じてスケールする必要があります。

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

与えられた $H$ に対する基底状態とは、エネルギーを最小化するスピンの特定の配置( $s_1, s_2, ...s_n$ )です。基底状態を見つけることはNP困難な最適化問題です。

要件

このチュートリアルを開始する前に、以下がインストールされていることを確認してください。

Qiskit Functions Catalog IBM Client (pip install qiskit-ibm-catalog)
Qiskit addon Optimization Mapper (pip install qiskit_addon_opt_mapper)
NumPy (pip install numpy)

また、ParityQC Twine Optimizer 関数にアクセスするための権限が必要です。アクセスをリクエストするには、このフォームに記入してください。

セットアップ

(このコードは、アカウントを保存する作業がすでにローカル環境で完了していることを前提としています。)

まず、このチュートリアルに必要なパッケージをすべてインポートします。

# 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

Qiskit Functions Catalog から Parity Twine Optimizer を読み込みます。

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

ステップ1: 問題を目的関数として定義する

マーケットスプリット問題の場合と異なり、SK問題はライブラリから取得するのではなく、直接定式化します。

generate_sk_problem 関数は、SK問題を必要な辞書形式で直接定式化します。必要な入力は、モデル内のスピン数を表す n のみです。

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

_generate_couplings_sk_model 関数を使用して、SKモデル内のランダムな結合項を計算します。結合をより細かく制御するには、関数のdocstringで説明されているオプション引数を使用できます。

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

ステップ2: Parity Twine Optimizer を使用して問題を解く

上記の関数を使用することで、SK問題を設定し、Twine Optimizer と選択したIBM Quantum®バックエンドを使用して解を求めることができます。

関数を実行するには、ibm_phoenix などの適切なバックエンドを選択します。

必要に応じて、options を使ってジョブ送信をさらに細かく制御できます。

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

ここで、shots は回路の実行回数を指定する整数、postprocessing_level は結果に後処理を適用するかどうかを決定し、 transpile_only は問題を回路にトランスパイルのみ行うか(解かないか)を選択し、job_tags はIBM Quantum Platform上でジョブを識別するためのラベルです。

SKモデルのサイズは、相互作用するスピンの数である $N$ によって定義されます。 $N$ を選択すると、上記のコードは n_spins に対する問題を生成します。

オプティマイザーを実行します。

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

ジョブのステータスを確認します。

# Monitor the job status
function_job.status()

結果を取得します。

result = function_job.result()

result

結果は次の形式になります。

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

ここで、solution 辞書は問題内で定義された量子ビットに対応し、SKモデルのハミルトニアンに対する最適化されたスピン値を示します。この最適解におけるこの特定のスピンの並びは、与えられたランダムな相互作用強度に基づいてシステム全体のエネルギーを最小化する配置を表します。 SKモデルでは、これは無秩序な磁性システムの最低エネルギー状態と見なすことができます。

metadata はトランスパイルに関する情報(2量子ビットゲートの数/深さ、使用されたゲート、アクティブな量子ビット)や各種の実行時間を示します。

学習成果​

背景​

Sherrington-Kirkpatrickモデル​

数学的定式化​

要件​

セットアップ​

ステップ1: 問題を目的関数として定義する​

ステップ2: Parity Twine Optimizer を使用して問題を解く​

次のステップ​