こんにちは！量子計算の世界へようこそ。
これまでの連載で、量子ゲートの基礎やMax-Cut問題を学んできました。今回は、いよいよ組み合わせ最適化における「ラスボス」とも言える巡回セールスマン問題（TSP）に挑みます。
「4都市をどう回れば最短か？」という一見シンプルなパズルですが、量子コンピュータで解くには16個の量子ビットを巧みに操る必要があります。数式の展開から、Qiskitでの最新の実装まで、世界一わかりやすく解説します！
1. 変数設計：なぜ「16個」も量子ビットが必要なのか？
Max-Cut問題では「4つのノード＝4ビット」という単純な対応でした。しかし、TSPはそうはいきません。なぜなら、「いつ（Step）」「どこに（City）いるか」という2つの情報を同時に記録しなければならないからです。
- 都市： A, B, C, D（4つ）
- ステップ： 0, 1, 2, 3（4回）
これらを掛け合わせると、 4×4=16 マスの「スケジュール表」が出来上がります。この1マスずつに「0（行かない）」か「1（行く）」を割り当てるため、16量子ビットが必要になるのです。
Pythonコードでは、この「2次元の表」を「1次元のビット列」に変換する住所変換ツール（idx関数）を定義します。
# 都市番号(city)と時間(step)を、16個のビットの何番目かに変換する
def idx(city, step):
return city * n_steps + step
例えば、「都市B（1）にステップ3（2）で訪問する」なら、1×4+2 = 6 番目のビットを見れば良い、という仕組みです。
2. 数式の魔法：QUBOから物理（イジング）への変換
量子コンピュータは「0か1」という論理（QUBO）を直接は理解できません。計算できるのは「1か-1」という物理的なスピン（Z ゲート）の世界です。
そこで、すべての数式に x = (1-Z)/2を代入して展開（イジング化）します。
ここでは、ブログの核心である「なぜコードの係数が 0.5 や 0.25 になるのか」を解き明かします。
① 制約条件：1ステップに1都市ルール
「ある時間 tに、必ず1つの都市にしか行けない」というルールを数式にすると：
P = A(x_0 + x_1 + x_2 + x_3 – 1)^2
これに x=(1-Z)/2 を代入して地道に展開（※中略）すると、驚きの結果になります。
$$P = \dots = 2A – A \sum Z_i + 0.5 A \sum_{i<j} Z_i Z_j$$ (lateXですみません。。。)
ここがポイント！
- 1次の項（-A）: 単独の Z ゲートに掛ける係数。
- 2次の項（0.5A）: 2つの Zゲートの相互作用（ZZ）に掛ける係数。
コードを見ると、この理論値がそのまま実装されているのがわかります！
# 1次の項：-A
Hconst1 += -1.0 * Z_op(idx(city, step))
# 2次の項：0.5A
Hconst1 += 0.5 * ZZ_op(idx(city1, step), idx(city2, step))
② コスト条件：距離の最小化
「都市 i から j へ移動したときの距離 $W$」のペナルティは：
$$\text{Cost} = W \cdot x_{i, t} \cdot x_{j, t+1}$$ (lateXですみません。。。)
これも代入してバラすと…
$$\text{Cost} = \frac{W}{4} (1 – Z_1 – Z_2 + Z_1 Z_2)$$
すべての項に 0.25 が出てきましたね。
# weight / 4 (0.25) が至る所に！
Hcost3 += (weight * 0.25) * I_op()
Hcost3 += -(weight * 0.25) * Z_op(bit_i_t)
Hcost3 += -(weight * 0.25) * Z_op(bit_j_t_next)
Hcost3 += (weight * 0.25) * ZZ_op(bit_i_t, bit_j_t_next)
3. レベル3の進化：手動で量子回路を組み立てる
ハミルトニアン（エネルギーの設計図）ができたら、それを「量子回路」に翻訳します。QAOA（量子近似最適化アルゴリズム）の心臓部です。
ここで登場するのが CNOTサンドイッチのショートカット である rzz ゲートです。
# 1次の項（単一のZ） -> 単独の回転ゲート(rz)
for i, hi in nonzero_h:
qc.rz(2.0 * g * hi, i)
# 2次の項（ZZの相互作用） -> 相互作用ゲート(rzz)
for (i, j), Jij in nonzero_J:
qc.rzz(2.0 * g * Jij, i, j)
- なぜ2倍するの？: 数式の $e^{-i \gamma \dots}$ の形に合わせるため、係数を2倍（
2.0 * g * hi）にして角度として流し込みます。 - RXミキサー: 最後に
rxゲートを並べて「量子トンネル効果」を引き起こし、壁を突き抜けて最適解を探します。
4. 実践！C++シミュレータによる超高速計算
16量子ビットの探索範囲は2^16 = 65,536 通り。これを普通のPythonで計算すると日が暮れてしまいます。そこで、今回はIBMが提供する最新エンジン AerSimulator (C++製) と EstimatorV2 を使用します。
これは、本物の量子コンピュータ（実機）で動かすための「業界標準フォーマット」です。人間が書いた綺麗な回路を、マシンが理解できる形にコンパイル（翻訳）して実行します。
# C++シミュレータの準備
sim = AerSimulator()
pass_manager = generate_preset_pass_manager(3, sim)
isa_circuit = pass_manager.run(qc) # 実機向けに翻訳！
まとめ：このコードが「最強」である理由
このプログラムがすごいのは、Qiskitの自動機能に丸投げせず、人間が紙とペンで導き出した数式（0.5AやW/4）をそのままコードに打ち込んでいる点です。
- 物理のルールを文字（I, Z）で作る。
- 数式から係数を抜き出す。
- 抜き出した係数に合わせてゲート（rz, rzz）を置く。
このプロセスを完璧にマスターすれば、どんなに複雑なビジネス課題（物流最適化やシフト作成など）が来ても、自分の手で量子回路を設計できるようになります！
実行コードの全貌
import numpy as np
from qiskit.quantum_info import SparsePauliOp
from qiskit_algorithms import QAOA
from qiskit.primitives import StatevectorSampler as Sampler
from qiskit_algorithms.optimizers import COBYLA
n_cities = 4
n_steps = 4
n_bits = n_cities * n_steps
city_labels = [“A”, “B”, “C”, “D”]
distance = np.array([
[0, 2, 5, 8],
[2, 0, 4, 6],
[5, 4, 0, 6],
[8, 6, 6, 0]
])
def idx(city, step):
return city * n_steps + step
def I_op(n_bits=n_bits):
return SparsePauliOp.from_list([(“I” * n_bits, 1.0)])
def Z_op(bit_index, n_bits=n_bits):
pauli_str_list = [‘I’] * n_bits
pauli_str_list[bit_index] = ‘Z’
pauli_str_list.reverse() # Qiskitのエンディアン（右端が0）に合わせる
pauli_str = “”.join(pauli_str_list)
return SparsePauliOp.from_list([(pauli_str, 1.0)])
def ZZ_op(bit_i, bit_j, n_bits=n_bits):
pauli_str_list = [‘I’] * n_bits
pauli_str_list[bit_i] = ‘Z’
pauli_str_list[bit_j] = ‘Z’
pauli_str_list.reverse()
pauli_str = “”.join(pauli_str_list)
return SparsePauliOp.from_list([(pauli_str, 1.0)])
def zero_op(n_bits=n_bits):
return SparsePauliOp.from_list([(“I” * n_bits, 0.0)])
# Hconst1: 各ステップで訪問する都市は1つだけ
Hconst1 = zero_op()
for step in range(n_steps):
Hconst1 += 2.0 * I_op()
for city in range(n_cities):
Hconst1 += -1.0 * Z_op(idx(city, step))
for city1 in range(n_cities):
for city2 in range(city1 + 1, n_cities):
Hconst1 += 0.5 * ZZ_op(idx(city1, step), idx(city2, step))
# Hconst2: 各都市は1回だけ訪問される
Hconst2 = zero_op()
for city in range(n_cities):
Hconst2 += 2.0 * I_op()
for step in range(n_steps):
Hconst2 += -1.0 * Z_op(idx(city, step))
for step1 in range(n_steps):
for step2 in range(step1 + 1, n_steps):
Hconst2 += 0.5 * ZZ_op(idx(city, step1), idx(city, step2))
# Hcost3: 距離の最小化
Hcost3 = zero_op()
for t in range(n_steps):
t_next = (t + 1) % n_steps
for i in range(n_cities):
for j in range(n_cities):
if i != j:
weight = distance[i][j]
bit_i_t = idx(i, t)
bit_j_t_next = idx(j, t_next)
# weight * x_{i,t} * x_{j,t+1} の展開式
Hcost3 += (weight * 0.25) * I_op()
Hcost3 += -(weight * 0.25) * Z_op(bit_i_t)
Hcost3 += -(weight * 0.25) * Z_op(bit_j_t_next)
Hcost3 += (weight * 0.25) * ZZ_op(bit_i_t, bit_j_t_next)
H = Hconst1 + Hconst2 + 0.01 * Hcost3
# 重複したPauli文字列をまとめて計算を軽くする
H = H.simplify()
# ==========================================
# 3. 高速化のためのハミルトニアン事前分解
# ==========================================
print(“— ハミルトニアンを解析中 —“)
const = 0.0
h = np.zeros(n_bits, dtype=float)
J = defaultdict(float)
for P, c in zip(H.paulis, H.coeffs):
coeff = float(np.real(c))
if abs(coeff) < 1e-12:
continue
z_mask = P.z
zs = list(np.where(z_mask)[0])
if len(zs) == 0:
const += coeff
elif len(zs) == 1:
h[zs[0]] += coeff
elif len(zs) == 2:
i, j = zs
a, b = (i, j) if i < j else (j, i)
J[(a, b)] += coeff
print(f”解析完了: 定数項={const:.2f}, 1次項(Z)の数={np.count_nonzero(np.abs(h) > 1e-12)}, 2次項(ZZ)の数={len(J)}”)
# ==========================================
# 4. 手動量子回路（Ansatz）の超高速構築
# ==========================================
p = 10 # レイヤー数（10段で深い探索を行う）
theta = ParameterVector(“θ”, 2 * p)
gammas = theta[:p]
betas = theta[p:]
nonzero_h = [(i, float(h[i])) for i in range(n_bits) if abs(h[i]) > 1e-12]
nonzero_J = [((i, j), float(coeff)) for (i, j), coeff in J.items() if abs(coeff) > 1e-12]
qc = QuantumCircuit(n_bits)
qc.h(range(n_bits))
for layer in range(p):
g = gammas[layer]
b = betas[layer]
# コストレイヤー (gamma)
for i, hi in nonzero_h:
qc.rz(2.0 * g * hi, i)
for (i, j), Jij in nonzero_J:
qc.rzz(2.0 * g * Jij, i, j)
# ミキサーレイヤー (beta)
for q in range(n_bits):
qc.rx(2.0 * b, q)
# ==========================================
# 5. C++シミュレータ (AerSimulator) による最適化
# ==========================================
print(“\nQAOA計算をC++シミュレータで高速実行中…”)
sim = AerSimulator()
pass_manager = generate_preset_pass_manager(3, sim)
isa_circuit = pass_manager.run(qc)
H_isa = H
layout = getattr(isa_circuit, “layout”, None)
if layout is None:
layout = getattr(isa_circuit, “_layout”, None)
if layout is not None:
try:
H_isa = H.apply_layout(layout)
except Exception:
pass
est = EstimatorV2(options={“default_precision”: 0.0, “backend_options”: {“method”: “statevector”}})
isa_param_names = [p_.name for p_ in isa_circuit.parameters]
def energy(params):
by_name = {theta[i].name: float(params[i]) for i in range(2 * p)}
ordered = [by_name[name] for name in isa_param_names]
pub = (isa_circuit, H_isa, ordered)
job = est.run([pub])
return float(np.real(job.result()[0].data.evs))
optimizer = COBYLA(maxiter=100)
x0 = np.random.rand(2 * p)
result = optimizer.minimize(energy, x0=x0)
print(f”最小エネルギー: {result.fun:.4f}”)
# ==========================================
# 6. 正しいルートのデコードと出力
# ==========================================
def decode_tour(bit_int):
tour = []
for t in range(n_steps):
chosen = [c for c in range(n_cities) if ((bit_int >> idx(c, t)) & 1) == 1]
if len(chosen) != 1:
return None # ルール違反（1ステップに2都市など）
tour.append(chosen[0])
if sorted(tour) != list(range(n_cities)):
return None # ルール違反（同じ都市に2回行ったなど）
return tour
def tour_length(tour):
return float(sum(distance[tour[t], tour[(t + 1) % n_steps]] for t in range(n_steps)))
# 最適化されたパラメータを回路に代入して状態ベクトルを取得
bind = {theta[i]: float(result.x[i]) for i in range(2 * p)}
qc_opt = qc.assign_parameters(bind, inplace=False)
psi = Statevector.from_instruction(qc_opt)
probs = psi.probabilities()
# 制約を満たす（ルール違反をしていない）ルートの中で、最も確率が高いものを探す
best = None
for k in range(1 << n_bits):
tour = decode_tour(k)
if tour is None:
continue
P = float(probs[k])
if best is None or P > best[2]:
best = (tour, tour_length(tour), P, k)
if best is None:
print(“\n⚠️ 有効なルートが見つかりませんでした。ペナルティ係数の調整などが必要です。”)
else:
tour, L, P, k = best
tour_labels = [city_labels[c] for c in tour]
print(“\n✅ 最も確率の高い有効なルート:”)
print(” -> “.join(tour_labels + [tour_labels[0]]))
print(f”総移動距離: {L}”)
print(f”観測確率: {P:.6f}”)
いかがでしたか？16量子ビットのTSPを自力で解く快感、ぜひ皆さんも味わってみてください！
あくまで量子計算は確率なので何回か実行するかp = 50などにするなど工夫は必要ですが、最適解には到達します。
次もお楽しみに！