数理最適化におけるChatGPTの活用 〜 その1：最小全域木問題＋OSS最適化ソルバー編

数理最適化では、従来より専門家による定式化やツール（モデリング言語・ソルバー）の利用、場合によってはアルゴリズムの構築が必要でした。特に、よくあるプログラミングのアルゴリズムと異なり、外部のツールを利用するため、若干特殊な知識が必要になる背景があります。

一方で、ChatGPTを利用してみたところ、従来人の手で行っていた部分が大幅に短縮されました。まだ、完全にChatGPTのみで動作する状態ではありませんが、数理最適化での活用の方法について、本記事にまとめます。

最小全域木問題に対する、PuLPを用いた実装の出力

最小全域木問題は、比較的シンプル・有名な数理最適化問題です。
有名な問題に対して、ChatGPTがPuLPによる実装を出力できるのか、検証を行いました。

対象とする問題

ここでは、下記のような閉路除去定式化を考えます。
式としての意味は、こちらを参照してください。: https://scmopt.github.io/opt100/05mst.html

\begin{equation*}
\begin{aligned}
& \underset{x}{\text{minimize}}
& & \sum_{e\in{E}}c_ex_e \\
& \text{subject to}
& & \sum_{e\in{E}}x_e = n-1, \\
&&& \sum_{e\in{E(S)}} x_e\leq{|S| - 1} , \quad \forall{S} \subset{V}, \\
&&& x_e\in\{0, 1\}, \quad \forall{e} \in {E}
\end{aligned}
\end{equation*}

入出力の結果

ChatGPTに対しては、シンプルに次のような文章を入力します。

最小木問題に対する閉路除去定式化を、PuLPで書いてください。

ソースコードの出力は次のように得られました。

import pulp

# ノード数
n = 6

# 辺の重み
c = [
    [0, 1, 3, 4, 0, 0],
    [1, 0, 2, 0, 5, 0],
    [3, 2, 0, 0, 6, 0],
    [4, 0, 0, 0, 0, 7],
    [0, 5, 6, 0, 0, 8],
    [0, 0, 0, 7, 8, 0]
]

# 変数を定義
x = pulp.LpVariable.dicts('x', [(i, j) for i in range(n) for j in range(i + 1, n)], cat=pulp.LpBinary)

# 最小全域木問題の目的関数
objective = pulp.lpSum([c[i][j] * x[(i, j)] for i in range(n) for j in range(i + 1, n)])

# 最小全域木問題の制約式
constraints = []
for i in range(n):
    # 次数が2以下のノードは木の端点であるため、必ず選ばれる
    if sum(x[(i, j)] for j in range(i + 1, n)) <= 1:
        constraints.append(x[(i, i + 1)] == 1)
    elif sum(x[(i, j)] for j in range(i + 1, n)) >= 2:
        constraints.append(x[(i, i + 1)] == 0)

# 閉路除去制約を追加
for i in range(n):
    for j in range(i + 1, n):
        if c[i][j] > 0:
            constraints.append(x[(i, j)] + x[(j, i)] <= 1)

# 問題を定義
problem = pulp.LpProblem('Minimum Spanning Tree with Cycle Removal', pulp.LpMinimize)

# 目的関数を最小化するように設定
problem += objective

# 制約式を追加
for constraint in constraints:
    problem += constraint

# PuLPを使用して問題を解く
status = problem.solve()

# 解の出力
print('Status:', pulp.LpStatus[status])
print('Objective:', objective.value())
for i in range(n):
    for j in range(i + 1, n):
        if x[(i, j)].value() == 1:
            print('Edge:', i, '-', j, 'Weight:', c[i][j])

結果の解釈

6点からなる全結合無向グラフに対して、cという重みを表す隣接行列が定義されています。
無向グラフであるため、対称行列になっています。

決定変数の定義であるx = pulp.LpVariable.dicts('x', [(i, j) for i in range(n) for j in range(i + 1, n)], cat=pulp.LpBinary)についても、問題なく実装できています。両端にあるノードのIDによって、エッジを定めるような定義であり、もともとの定義では$x_e$となされていた部分です。

目的関数についても、PuLPの文法に沿う記述ができています。objective = pulp.lpSum([c[i][j] * x[(i, j)] for i in range(n) for j in range(i + 1, n)])はfor文直書きではなく、リスト内包表記とpulp.lpSumとを組み合わせて書かれています。少しこなれた印象のある書き方ですね。

問題は制約条件です。制約条件は、任意の部分集合を表す$S\subset{V}$が難しさの原因で、人が実装するときにも悩むポイントではないでしょうか。実際、上記のモデルを動作させても、ここが誤りであるため、動作が詰まってしまいます。

正しくは、制約条件部分を次のように修正する必要がありますね（下記は筆者による人力の実装）。目的関数の定義や、考え方のフレームとなる実装・添字のとり方は提示されるため、（慣れた人ならば）容易に修正することができます。

# エッジの数=ノードの数-1
# Σ_{e} x_e = n-1
constraints = []
left = 0
for i in range(n):
    for j in range(i+1, n):
        left += x[(i, j)]
constraints.append(left == n-1)

# 任意の誘導部分グラフに対し、その誘導部分グラフ内の枝の総和=誘導部分グラフのサイズ-1
all_subsets = []
for l in range(2, n):
    for subset in itertools.combinations([i for i in range(n)], l):
        all_subsets.append(subset)

for subset in all_subsets:
    left = 0
    for idx, s in enumerate(subset):
        if idx == n:
            continue
        for t in subset[idx+1:]:
            left += x[s, t]
    constraints.append(left <= len(subset)-1)

感想

数理最適化だけではなく、競技プログラミングなどの領域でもそうですが、ChatGPTが完全に単独で作っても、現状ではまだエラーが発生します。一方で、考え方のフレームが瞬時に提示されるメリットは非常に大きく、PuLPの文法や変数定義の仕様をほぼ意識せずにコードが完成します。ChatGPTが出力したコードを鵜呑みにすると、仮に動いたとしても思わぬ制約の見落としなどありそうですが、工数削減に対して有効という感想です。

networkXを用いた可視化コードの出力

合わせて、最小木問題に対するnetworkXによる可視化コードの作成を依頼します。

入出力の結果

上記に対して更に、最小木問題の結果が可視化されるコードを書いてください。

次のような出力が得られました。

# 可視化
G = nx.Graph()
G.add_nodes_from(range(n))
for i in range(n):
    for j in range(i + 1, n):
        if x[(i, j)] == 1:
            G.add_edge(i, j, weight=c[i][j], color='r')  # 赤色で最小木の部分
        else:
            G.add_edge(i, j, weight=c[i][j], color='k')  # 黒色で最小木の部分でない部分
pos = nx.circular_layout(G)
edge_colors = nx.get_edge_attributes(G, 'color').values()
nx.draw_networkx_nodes(G, pos, node_size=700)
nx.draw_networkx_edges(G, pos, edge_color=edge_colors)
nx.draw_networkx_edge_labels(G, pos, edge_labels={(i, j): c[i][j] for i, j in G.edges()})
nx.draw_networkx_labels(G, pos, font_size=20, font_family='Yu Mincho')
plt.axis('off')
plt.show()

結果の解釈

これは、ほぼ完璧に近いコードなのですが、実は一点だけ誤りがあります。決定変数の中身の値を取得するときは、if x[(i, j)] == 1:ではなく、if x[(i, j)].value() == 1:であるため、修正が必要です。なお、この部分の修正も、エラーの原因をデバッグしながら、人手で実施しました。

感想

networkXのAPIは、Pandasやmatplotlibと記憶が混ざってしまって、書くことが面倒くさいです。
しかし、煩雑な部分も見やすいように調整しながら、かなり正確にグラフを作っていることがわかります。
networkXはPuLPよりも有名なライブラリであるため、その分学習が正確である影響もあるのだろうか？という印象もあります。

最適化のモデルでは数値で出力されるため、その可視化を行うことが面倒くさいことが多いのですが、かなり実用的な印象です。
定式化側と異なり、ある程度可視化は定石があるはずであり、PuLPとnetworkXの繋ぎの部分のような領域だけ、人が介在すると良さそうです。