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) 모델은 통계 역학, 특히 스핀 유리(spin glass) 연구에서 기초가 되는 모델입니다. 상호작용이 보통 최근접 이웃으로 제한되는 표준 이징 모델과 달리, SK 모델은 무한 범위 모델로, 시스템 내 모든 스핀이 다른 모든 스핀과 상호작용함을 의미합니다. 이는 유리 상태(glassy) 거동의 특징인 많은 국소 최솟값(local minima)을 가진 매우 복잡하고 "울퉁불퉁한" 에너지 지형(energy landscape)으로 이어집니다.

SK 모델의 핵심 특징은 좌절(frustration)에 있습니다. 이 모델에서 스핀 간 상호작용 강도 $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$ 일 때 에너지가 크기 변수(extensive), 즉 $N$ 에 비례하는 양으로 유지되도록 하려면 결합의 분산이 입자 수에 따라 스케일링되어야 합니다.

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

주어진 $H$ 에 대한 바닥 상태(ground state)는 에너지를 최소화하는 스핀의 특정 구성( $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 카탈로그에서 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 모델의 무작위 결합 항을 계산합니다. 결합에 대한 더 세밀한 제어가 필요하면 함수의 독스트링에 설명된 선택적 인수를 사용할 수 있습니다.

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를 사용해 문제 해결하기​

다음 단계​