はじめに

量子力学は「ミクロな世界の不思議な法則」を説明する理論ですが、その性質は計算にも応用できます。
CIM（Coherent Ising Machine）はその一つで、レーザーの性質やゆらぎを使って組合せ最適化問題を解こうとするマシンです。

実際のCIMは光学実験室にある大がかりな装置ですが、その挙動を Python上でシミュレーションできるライブラリ として cim-optimizer があります。

この記事では、このライブラリを使って代表的な組合せ最適化問題である MAX-CUT問題を解くチュートリアルを、できるだけていねいに紹介します。
私自身まだ初学者ですが、同じように初めて触れる方と一緒に学びを楽しめたら嬉しいです。

環境構築

筆者は、ローカルとGoogle Colabの二つで試したのでそれぞれのケースを紹介します。

ローカル（Mac Mチップ）で試す場合

Anacondaを使ってローカルのMacのApple Silicon(Mチップ)で動作確認をした時のコードを紹介します。
なお、インストールにハマった場合は、cim-optimizer の公式のインストールベージも合わせて活用してください。

# 1) 仮想環境を作成（Python 3.10 例）
conda create -n cim-env python=3.10 -y
conda activate cim-env

# 2) 科学計算系の基本
conda install numpy matplotlib -y

# 3) PyTorch（MacはCUDAではなくMPS/CPU）
#    ※ Apple Silicon (M1/M2/M3) では pytorch-cuda は不要＆インストールできません
conda install pytorch torchvision torchaudio -c pytorch -y

# 4) cim-optimizer 本体
python -m pip install cim-optimizer

GPU(MPS)の動作確認

import torch
print("PyTorch:", torch.__version__)
print("MPS built?:", torch.backends.mps.is_built())
print("MPS available?:", torch.backends.mps.is_available())
print("CUDA available?:", torch.cuda.is_available())  # Macでは通常 False

Mac (Apple Silicon)
- MPS available? が True → Metal(GPU) が利用されています
- False → CPU 実行になります
Linux/Windows + NVIDIA GPU
- CUDA available? が True → CUDA(GPU) が利用されています
- False → CPU 実行になります

GPUを使いたい場合

基本的な動作確認や小さな問題なら CPU で十分です。

もし GPU を使いたい場合は、環境依存の注意点があるので
公式のインストールガイドに従って設定するのがおすすめです。

Google Colab で試す場合

すぐに試したい場合は Google Colab でも動作します。
GPU を使用する場合は、[ランタイム] → [ランタイムのタイプを変更] → [ハードウェアアクセラレータ] で GPU（例: T4 GPU）を選択してください。
CPU を使いたい場合は [ハードウェアアクセラレータ] で [CPU] を選びます。

Colab には NumPy や Matplotlib などの基本ライブラリはあらかじめ揃っているため、追加で必要なのは cim-optimizer だけです。

こちらの Colab ノートブックでも実行例を確認できます。

MAX-CUT問題の基礎

MAX-CUTとは？

MAX-CUT はグラフ理論で扱う組合せ最適化問題のひとつです。

グラフ理論では、点を頂点（vertex）、点と点を結ぶ線を辺（edge）と呼びます。向きのない線で構成されるものを無向グラフ（undirected graph）、矢印で向きを持つ場合は有向グラフ（directed graph）と呼びます。

MAX-CUT では、頂点を二つのグループに分けて、グループ間を結ぶ辺の重みの合計をできるだけ大きくすることを目指します。

この問題は NP 困難に分類され、規模が大きくなると厳密解を効率的に求めるのは難しいことが知られています。そのため近似解法や量子計算を応用した手法が研究されています。

数式での定義

MAX-CUT を数学で定式化すると次のようになります。

\text{Cut}(s) = \frac{1}{2} \sum_{(i,j)\in E} w_{ij}\,(1 - s_i s_j)

ここで $s = (s_1, s_2, \dots, s_N)$ は各ノードのスピン値を並べたベクトルで、
それぞれの $s_i$ は $\pm 1$ をとります（2つのグループのどちらに属するかを表しています）。

$s_i s_j = +1$ → 同じグループ（寄与しない）
$s_i s_j = -1$ → 異なるグループ（寄与する）

MAX-CUT と Ising モデルの対応関係

Ising モデルでは次のように相互作用エネルギーが定義されます。

E(s) = \sum_{i<j} J_{ij}\, s_i s_j + \sum_i h_i s_i

ここで $J_{ij}$ はスピン間の結合、$h_i$ は外場を表します。

MAX-CUT を Ising モデルに対応づけるときには、外場項は現れないため
すべての $h_i$ を 0 とおくことができます。このとき、

E(s) = \sum_{i<j} J_{ij}\, s_i s_j

となります。さらに $J_{ij} = -w_{ij}$ と置くと

E(s) = \sum_{i<j} w_{ij}\, s_i s_j

エネルギー $E(s)$ を小さくするほど Cut の値が大きくなる関係となります。

小さな例で確認（N=4）

MAX-CUTの仕組みを、簡単な例で確認してみましょう。

下図を描画するためのコードはこちら

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from cim_optimizer.CIM_helper import brute_force  # 総当たり（Nが小さいとき用）


# --- 小さい重み付きグラフ（対称行列） ---
N = 4
W = np.array([
    [0, 2, 5, 1],
    [2, 0, 3, 6],
    [5, 3, 0, 4],
    [1, 6, 4, 0]
], dtype=float)

# NetworkX グラフへ（重み0は辺を張らない）
G = nx.Graph()
for i in range(N):
    for j in range(i+1, N):
        if W[i, j] != 0:
            G.add_edge(i, j, weight=W[i, j])

# --- Cut と Ising エネルギー ---
def cut_weight(spins, W):
    s = spins.reshape(-1, 1)
    M = (1 - s @ s.T) / 2
    return float(np.sum(np.triu(W * M, k=1)))

def ising_energy(spins, J, h=None):
    if h is None:
        h = np.zeros_like(spins)
    return float(-0.5 * np.sum(J * np.outer(spins, spins)) - np.dot(h, spins))

# --- Brute force で最適スピンを取得（MAX-CUT の写像：J = -W, h=0） ---
J = -W
spins_opt, E_opt = brute_force(J)
cw_opt = cut_weight(spins_opt, W)

# --- 可視化の下準備 ---
pos = nx.circular_layout(G)  # 均等（円形）配置

def split_edges(G, spins):
    cut = [(u, v) for u, v in G.edges() if spins[u] * spins[v] < 0]
    uncut = [(u, v) for u, v in G.edges() if spins[u] * spins[v] > 0]
    return cut, uncut

groupA = [i for i, s in enumerate(spins_opt) if s > 0]
groupB = [i for i, s in enumerate(spins_opt) if s < 0]
cut_edges, uncut_edges = split_edges(G, spins_opt)

# --- 2面並べ：左=入力 / 右=最適解 ---
fig, axes = plt.subplots(1, 2, figsize=(10, 5))


# まず全体を普通に描く
nx.draw_networkx(G, pos, with_labels=True, node_color="lightgray", ax=axes[0])

# --- 交差していない辺はまとめてラベル ---
other_edges = [(i, j) for i, j in G.edges() if (i, j) not in [(1,3), (0,2)] and (j,i) not in [(1,3), (0,2)]]
nx.draw_networkx_edge_labels(
    G, pos, ax=axes[0],
    edge_labels={(i, j): W[i,j] for i, j in other_edges},
    font_size=14, bbox=dict(facecolor="white", edgecolor="none", pad=0.2)
)

# --- 特殊扱い：交差して見づらい 2 本だけ別にラベルをずらして描く ---
nx.draw_networkx_edge_labels(
    G, pos, ax=axes[0], edge_labels={(1,3): W[1,3]},
    font_size=14, bbox=dict(facecolor="white", edgecolor="none", pad=0.2),
    label_pos=0.55
)
nx.draw_networkx_edge_labels(
    G, pos, ax=axes[0], edge_labels={(0,2): W[0,2]},
    font_size=14, bbox=dict(facecolor="white", edgecolor="none", pad=0.2),
    label_pos=0.35
)


axes[0].set_title("Original weighted graph")
axes[0].axis("off")

# 右：ブートフォース最適解（色分け＋赤がカット）
nx.draw_networkx_nodes(G, pos, nodelist=groupA, node_color="tab:blue", ax=axes[1], label="Group A")
nx.draw_networkx_nodes(G, pos, nodelist=groupB, node_color="tab:orange", ax=axes[1], label="Group B")
nx.draw_networkx_edges(G, pos, edgelist=cut_edges,   edge_color="red",  width=2, ax=axes[1])
nx.draw_networkx_edges(G, pos, edgelist=uncut_edges, edge_color="gray", alpha=0.5, ax=axes[1])
nx.draw_networkx_labels(G, pos, ax=axes[1])


# --- 交差していない辺はまとめてラベル ---
other_edges = [(i, j) for i, j in G.edges() if (i, j) not in [(1,3), (0,2)] and (j,i) not in [(1,3), (0,2)]]
nx.draw_networkx_edge_labels(
    G, pos, ax=axes[1],
    edge_labels={(i, j): W[i,j] for i, j in other_edges},
    font_size=14, bbox=dict(facecolor="white", edgecolor="none", pad=0.2)
)

# --- 特殊扱い：交差して見づらい 2 本だけ別にラベルをずらして描く ---
nx.draw_networkx_edge_labels(
    G, pos, ax=axes[1], edge_labels={(1,3): W[1,3]},
    font_size=14, bbox=dict(facecolor="white", edgecolor="none", pad=0.2),
    label_pos=0.55
)
nx.draw_networkx_edge_labels(
    G, pos, ax=axes[1], edge_labels={(0,2): W[0,2]},
    font_size=14, bbox=dict(facecolor="white", edgecolor="none", pad=0.2),
    label_pos=0.35
)

axes[1].set_title(f"Brute force optimum\nCut={cw_opt:.1f},  E={E_opt:.1f}")
axes[1].legend()
axes[1].axis("off")

plt.tight_layout()
plt.savefig("max_cut_example.png", bbox_inches="tight")
plt.show()

print("=== Summary ===")
print("Optimum: Cut =", cw_opt, "  Ising Energy =", E_opt)

左図はノード（頂点）数 $N=4$ の重み付きグラフです。
点を2つのグループに色分けし、異なる色のペアが結ばれたときに「重み $w$」がカットに加算されます。
重みが大きい辺をできるだけグループ間に含めたい、というのが MAX-CUT の狙いです。

この例では、真ん中の辺（重み5と6）が特に大きいので、それらをまたぐようにグループ分けすると有利です。
実際に最適な分割は右図のようになります。

グループ分け: A = {1,2}（青）, B = {0,3}（橙）
このときカット（赤線）は次の4本です：

(1–0) = 2
(1–3) = 6
(2–0) = 5
(2–3) = 4

合計すると 2+6+5+4 = 17 → Cut = 17.0

Ising エネルギーで表すと？

同じ状態を「スピン」で表すと、青グループを $+1$、橙グループを $-1$ と置けます。
各辺ごとの寄与を計算すると：

(0–1): $2\times(-1)(+1) = -2$
(0–2): $5\times(-1)(+1) = -5$
(0–3): $1\times(-1)(-1) = +1$
(1–2): $3\times(+1)(+1) = +3$
(1–3): $6\times(+1)(-1) = -6$
(2–3): $4\times(+1)(-1) = -4$

合計すると -2 -5 +1 +3 -6 -4 = -13 → E = -13.0

合成MAX-CUTをCIMで解いてみる（N=20）

小さい系 ($N=4$) ではブルートフォース（全探索）で最適解を直接求められました。
しかしノード数 $N$ が増えると、候補は $2^{N-1}$ 通りに指数関数的に増えてしまい、例えば $N=20$ で約50万通り、$N=30$ では5億通りを超えます。

そのため実際には「ヒューリスティック（近似解法）」を使って、短時間で良い解を探すのが一般的です。CIM もその一種です。
本来のCIM研究では「大規模なMAX-CUT」を物理デバイスで解くことが目的ですが、ここでは Pythonライブラリ cim-optimizer を使ってシミュレーションしてみましょう。

（公式サンプルは .npz ファイルを読み込む形ですが、pip 版では扱いが少し複雑になるため、本記事では簡単のためにランダムな対称行列を使った例を紹介します。）

サンプルコード

import numpy as np
import matplotlib.pyplot as plt
from cim_optimizer.solve_Ising import Ising
from cim_optimizer.CIM_helper import brute_force  # 小規模なら正解チェックに使える


# --- 合成の MAX-CUT 用データを作る ---
N = 20
rng = np.random.default_rng(0)
W = rng.standard_normal((N, N))
W = (W + W.T)/2
np.fill_diagonal(W, 0.0)

# MAX-CUT→Ising の定番変換：J にマイナスをつける
J = -W
h = np.zeros(N, dtype=float)

# （任意）小規模なら総当たりで下限エネルギーを確認
spins_g, E_g = brute_force(J)
print("Ground energy (brute force):", E_g)

# --- CIMで解く ---
# gamma は少し上げると収束しやすい例（公式ノートの流儀）
gamma = 0.010
res = Ising(J, h).solve(cac_gamma=gamma, hyperparameters_randomtune=False)

# バージョン差に強い取り出しヘルパー（←“さっきのエラー”対策）
def get_energy_spins(obj):
    r = getattr(obj, "result", obj)
    # エネルギー候補
    for k in ["lowest_energy", "best_energy", "energy", "final_energy"]:
        if hasattr(r, k):
            E = getattr(r, k); break
    else:
        E = None
    # スピン候補（環境で名前が違うことがある）
    for k in ["lowest_energy_spin_config",  # ← この名前の版がある
              "lowest_energy_spins",
              "best_spins",
              "lowest_energy_state",
              "best_state",
              "state",
              "spins"]:
        if hasattr(r, k):
            S = getattr(r, k); break
    else:
        S = None
    return E, S

E, S = get_energy_spins(res)
print("Best energy:", E)
print("Best spins: ", S)

# 進化の可視化（任意）
res.result.plot_spin_trajectories(plot_type="spins")
plt.show()

res.result.plot_spin_trajectories(plot_type="energy")
plt.show()

出力される結果

コンソール

コンソール結果

No External Field Detected
Target Ising Energy: -inf.
Best Ising Energy Found: -42.04461179559532.
Corresponding Spin Configuration: [ 1.  1. -1.  1. -1. -1. -1. -1. -1.  1.  1. -1.  1.  1. -1.  1.  1.  1.
 -1. -1.].
Time Elapsed: 0.08697199821472168.
Number of Runs Completed: 1.

コンソール出力の意味

No External Field Detected
外場 $h$ がゼロなので「外部からの影響はなし」という確認メッセージです。
Target Ising Energy: -inf.
目標エネルギーを指定していないので「未設定」を意味します。
Best Ising Energy Found: -42.04...
CIM が見つけた最小エネルギー値。これが最適解の指標になります。
Corresponding Spin Configuration: [ ... ]
そのときのスピン配置（±1 の並び）。+1 と -1 がグループ分けを表します。
Time Elapsed
計算にかかった時間。
Number of Runs Completed
実行した回数（ここでは1回）。

グラフの読み方

スピンの時間変化

最初は全スピンが 0 付近にあり（グループ未定）、時間が進むと ±1 に分かれて安定します。これは「どのノードがどちらのグループに入るか」が決まっていく過程を表します。途中で ±1 を少し超えるのは、CIMの連続的な振幅進化による自然なオーバーシュートです。

エネルギーの時間変化

エネルギーは段階的に下がっていき、最後に一定値で安定します。この段階的な下降は「より良いグループ分けに乗り換えた瞬間」を意味し、最終的に基底エネルギーに一致したことで、CIM が正しい最適解を見つけたことが確認できます。
今回、最終的なエネルギー値（約 -42）は brute forceで求めた基底エネルギーと一致しており、CIM が正しく MAX-CUT の解を見つけられたことが確認できます。

コードの注意点

フィールド名の違いに注意
環境やバージョンによって、
- lowest_energy_spins
- lowest_energy_spin_config
  どちらかしか存在しない場合があります。
  → 事前に print(dir(res.result)) で確認するか、前述の「ヘルパー関数」で吸収すると安全です。
グラフ保存時の余白
res.result.plot_spin_trajectories(plot_type="spins") をプロットして plt.show() すると、x 軸ラベルが下で切れてしまうことがあります。
→ plt.savefig("cim_spins.png", bbox_inches="tight") を使えば、余白込みできれいに保存されます。

オプションを試す

少し使い方に慣れてきたら、solve() のオプションを増やして実行してみましょう。
たとえば num_runs（試行回数）や num_parallel_runs（並列実行数）を指定すると、結果の安定性や計算時間が変わります。

また、出力に含まれる Time Elapsed（経過時間） を見れば、CIM がどれくらいの速さで解に到達したかも確認できます。
問題のサイズを大きくしたときに「時間がどう伸びるのか」を比べてみると面白いと思います。

並列実行（CPU/MPS でもOK）

res = Ising(J, h).solve(
    cac_gamma=gamma,
    num_runs=5,              # 何回トライするか
    num_parallel_runs=5,     # 同時並列（CPU/MPSでも加速することがある）
    hyperparameters_randomtune=False
)

目標エネルギーを入れて達成判定（ベンチしたいとき）

res = Ising(J, h).solve(
    cac_gamma=gamma,
    target_energy=float(E_g),  # 既知の下限を入れると「達成/未達」が分かる
    hyperparameters_randomtune=False
)

ハマりどころ（Q&A集）

（Mac）CUDA を入れようとしてコケる
→ Apple Silicon では CUDA なし。pytorch-cuda は不要。conda install pytorch ... -c pytorch でOK。
pip3 install ... で “externally-managed-environment”
→ その pip3 は Homebrew の Python を見てる。必ず仮想環境を有効化してから python -m pip install ... を使う。
lowest_energy_spins が無いと言われる
→ lowest_energy_spin_config を使う版があります（この記事ではヘルパー関数を用意して回避）。
Jupyter 以外で %matplotlib inline を書いてエラー
→ それはノートブック専用のマジック。スクリプトなら不要。
公式 .npz が読み込めない
→ pip には普通入っていません。リポジトリを clone して instances/MC_Instances_NPZ を使うか、上のように合成 Jで試す。
GPU を本当に使ってるか不安
→ Mac は torch.backends.mps.is_available() が True なら MPS（Metal）を使う。False でも CPU で動けばOK（まずは動作優先で大丈夫）。

まとめ

この記事では、CIM（Coherent Ising Machine）の Python 実装ライブラリ cim-optimizer を使って、代表的な組合せ最適化問題である MAX-CUT を解く流れを紹介しました。

まずは MAX-CUT と Ising モデルの対応づけを確認し、そのうえで CIM の仕組みを踏まえ、実際に cim-optimizer を用いて解く手順を体験しました。

本記事が、cim-optimizer を動かしてみたい方の参考になれば幸いです。

参考

CIMの最適化の原理については、こちらの記事でも簡単に紹介しています。
興味のある方は、あわせてご覧ください。