クソコード製造機

数理最適化とかPythonとか

最適化におけるスパース性を考慮したモデリングについて

最適化モデルを作成するときに、スパースな定式化を意識することにより、モデリングや前処理にかかる時間やメモリ使用量を大幅に改善できる場合がある。

これは、例えば、疎グラフにおいては、隣接行列でデータを持つよりも隣接リストなどでデータを持った方がオーダーレベルで時間やメモリ使用量が改善するのと同様である。
特に整数計画問題などとしてモデリングする場合は、変数の作成やメモリ確保の点で定数倍が大きいため、時には絶大な効果を発揮する。さらに、モデリングの仕方によっては宣言する必要のある制約が減少する可能性もあるため、モデリングや前処理がボトルネックとなっているような場合は、真っ先に検討する必要がある部分だと考える。

以下、gurobipyとpandasを例に、簡単な例を挙げてみる。

例えば、3頂点のグラフに対して、次のようなDataFrameで隣接行列が与えられている場合を考える。

import pandas as pd
df = pd.DataFrame([[0, 1, 1], [1, 0, 0], [0, 1, 0]])
print(df)
# 出力結果
#    0  1  2
# 0  0  1  1
# 1  1  0  0
# 2  0  1  0

このようなグラフ上に変数を定義したい場合に、変数 V \times Vで定義し、後からその変数を0で固定するというのが密な定式化である。

from itertools import product
import gurobipy as gp

model = gp.Model()
V = df.index
x = model.addVars(product(V, V), vtype="C", name="x")

for i, j in product(V, V):
    # 辺がない場合は0
    if df.loc[i, j] < 0.5:
        model.addConstr(x[i, j] == 0)

このような定式化は、変数の数がグラフのスパース性にかかわらず O(n^2)となるため、基本的に避けるべきであり、変数を E上に定義するべきである。
例えば次のようにstackメソッドを用いることで、 EをMultiIndexとして取り出すことができる。

ser = df.stack()
E = ser[ser > 0.5].index
print(E)
# 出力結果
# MultiIndex([(0, 1),
#             (0, 2),
#             (1, 0),
#             (2, 1)],
#            )

x = model.addVars(E, vtype="C", name="x")

さらに、gurobipy-pandasを使用すれば、変数の線形結合などを求めたいときに、pandasの集計関数を使えるようになり、便利になる。

import gurobipy_pandas as gppd

x = gppd.add_vars(model, E, name="x", vtype="C")
model.update()
print(x)
# 出力結果
# 0  1    <gurobi.Var x[0,1]>
#    2    <gurobi.Var x[0,2]>
# 1  0    <gurobi.Var x[1,0]>
# 2  1    <gurobi.Var x[2,1]>
# Name: x, dtype: object

Gurobi Python APIにおけるmodel.update()について

Gurobi Python APIのコードを見ていると、次のようなコードをよく見かける。

import gurobipy as gp

model = gp.Model()

# model.addVarで変数を作る処理を書く

model.update()

# model.addConstrで制約を書く

model.optimize()

このコード自体は問題なく動くが、model.update()を考えたい。実はこのコードはmodel.update()がなくても動作する。(最新バージョンの場合)

# 動く
import gurobipy as gp

model = gp.Model()

# model.addVarで変数を作る処理を書く

# model.addConstrで制約を書く

model.optimize()

つまり、なんらかの事情があって古いバージョン(v.6以前かな?)を使ったり、互換性を保ちたい場合を除いて、model.update()を呼び出すのは不要である。

そのうえ、不必要にmodel.update()を呼び出すことはパフォーマンス低下にも繋がってしまうため、避けるべきである。(遅延評価の良さを生かしきれない。)


それでも、限定的ではあるが、model.update()を呼びだす局面がある。例えば次のような状況である。

モデルのデバッグをするとき

制約のデバッグなどをするときに、作った変数などをprintしても役に立つ情報が得られない場合が多い。例えば次の状況である。

x = model.addVar(vtype="I", name="x")
y = model.addVar(vtype="I", name="y")
z = x + 2 * y
print(z)  # zの中身を確認したい。

遅延評価のため、この出力結果は<gurobi.LinExpr: <gurobi.Var *Awaiting Model Update*> + 2.0 <gurobi.Var *Awaiting Model Update*>>のようになってしまうが、print文の前にmodel.update()を呼びだしておくことで、<gurobi.LinExpr: x + 2.0 y>のような出力が得られるため、デバッグがしやすくなる。

変数のアトリビュートにアクセスしたいとき

これも用途はデバッグと近いが、変数のアトリビュートにアクセスするためには、model.update()を呼び出しておく必要がある。これはログ出力などをしたい場合にも注意が必要な点である。

x = model.addVar(vtype="I", name="x", ub=30)
model.update()  # 無いと以下がエラーになる。
print(x.vtype)
print(x.varname)
print(x.ub)

Gurobiの階層型多目的最適化の制御

Gurobiで多目的最適化をする際に、それぞれの目的関数における振る舞いを制御したいときがある。

今回は例として、目的関数が3つの多目的最適化について、優先度を設定して順々に最適化することを考える。

import gurobipy as gp

model = gp.Model()

# モデルをがんばって定義する

# obj0, obj1, obj2にはそれぞれ目的関数を表すLinExprオブジェクトを入れる

model.setObjectiveN(obj0, 0, 2)
model.setObjectiveN(obj1, 1, 1)
model.setObjectiveN(obj2, 2, 0)

こうすることで、obj0, obj1, obj2の順で最適化が行われる。第3引数が優先度であり、この値が大きい目的関数が優先して最適化される。第2引数はindexであり、後ほど目的関数を指定するときに使う。

今回は次のようにソルバーを制御することを考える。

  • 各目的関数の時間制限として、10分を設定する。
  • 各目的関数において、分枝限定法の前にヒューリスティクスにおける解の探索(NoRelHeur)を行う。これはobj0では300s、他2つでは200s行うことにする。
  • 相対許容誤差(MIPGap)を設定し、誤差以内になったら求解を終了する。許容誤差は、obj0とobj2では1%、obj1では5%とする。

このような場合には、次のように設定すれば良い。

# 全体共通の設定
model.setParam("TimeLimit", 2000)
model.setParam("NoRelHeurTime", 200)
model.setParam("MIPGap", 0.01)

# 各目的関数における設定
env0 = model.getMultiobjEnv(0)
env0.setParam("TimeLimit", 600)
env0.setParam("NoRelHeurTime", 300)

env1 = model.getMultiobjEnv(1)
env1.setParam("TimeLimit", 600)
env1.setParam("MIPGap", 0.05)

env2 = model.getMultiobjEnv(2)
env2.setParam("TimeLimit", 600)

基本的には、目的関数個別に設定されていないパラメータについては、全体のパラメータが参照される。

注意が必要なのがTimeLimit(WorkLimitを使う場合も)であり、全体のモデルに適用した時間制限は、すべての目的関数を合わせた制限時間となる。たとえば、ここで

model.setParam("TimeLimit", 600)

と設定してしまうと、600s経過した時点ですべての最適化が打ち切られてしまうため、obj1の最適化が一切行われない可能性がある。また、TimeLimitを600sに設定していても、実際は600sを少しオーバーすることがあるので、1800sから少し余裕を持たせている。

今回は例示のために、設定情報をマジックナンバーとしてコードに埋め込んでいるが、実際のシステムでは、yamlファイルなど別ファイルに書き出したものを読み取る形にするのが、保守性的には良いと思われる。

Pythonの標準ライブラリでトポロジカルソートをする

Pythonでグラフの操作をしたい場合は、たいていnetworkxをインストールして使えば事足りる。(フローや最短路、最小全域木など)

ただし、なぜかトポロジカルソートは以下の標準ライブラリとして使用可能である。

from graphlib import TopologicalSorter

使い勝手も非常に独特である。TopologicalSorter.add("a", "b", "c")とすると、b → aとc → aのエッジがはられる。(ノードがなければ追加される。)
ts.static_order()で、トポロジカルソートの結果のジェネレータが得られる。

ts = TopologicalSorter()
ts.add("a", "b", "c")
ts.add("b", "c")
print(list(ts.static_order()))  # => ["c", "b", "a"]

辞書でグラフを構築することもできるが、この場合でも、それぞれのノードへの入力辺を持つノードを指定することになる。

graph = {"a": {"b", "c"}, "b": {"c"}}
ts = TopologicalSorter(graph)
print(list(ts.static_order()))  # => ["c", "b", "a"]

おそらくこのライブラリは、タスク関係の管理に重点を置いているのだと思われる。
実際、実行済みタスクを登録して、現在実行可能なタスクを取得するようなメソッドも実装されている。

なお、graphlibという大層な名前になっているが、現状実装されているのはTopologicalSorterのみであるため、基本的にはnetworkxを使うことになるだろう。
今後拡充されていく予定などがあるのだろうか。

GurobiのLazy Constraintで制約を逐次追加する

組合せ最適化問題について、整数計画の定式化をしていると、表現したい条件についての制約が非常に多くなってしまうことがある。

このような場合には、必要な制約のみを逐次追加していく切除平面法が用いられることがある。

本記事では、それをGurobiのLazy Constraintという機能を用いて実装していく。

問題としては、無向グラフ G = (V, E)と辺の重み d_eについての最小全域木問題に対する、カットセット不等式を用いた定式化を考える。
また、問題例としては、重みをランダムに定めた格子グラフを生成する。
その生成コードは次の通り。

import networkx as nx
import random
import gurobipy as gp
from gurobipy import GRB

n, m = 5, 5
lb, ub = 1, 100

random.seed(42)
G = nx.grid_2d_graph(n, m)
V = n * m
nodes = list(G.nodes)

d = {}
for e in G.edges:
    d[e] = random.randint(lb, ub)

カットセット不等式を用いた定式化は、次のようになる。


 \displaystyle
\begin{aligned}
&\mathrm{minimize} &&\sum_{e \in E} d_e x_e &\\
&\mathrm{s. t.} & &\sum_{e \in E} x_e = |V| - 1,& \\
&&& \sum_{e \in \delta(S)} x_e \geq 1  &\{\} \neq S \subsetneq V \\
&&& x_e \in \{0, 1\}, & e \in E
\end{aligned}

これをGurobiで素直に実装すると次のようになる。

model = gp.Model()
x = model.addVars(G.edges, vtype=GRB.BINARY)

model.setObjective(x.prod(d))

model.addConstr(gp.quicksum(x) == V - 1)

for s in range(1, 1 << V - 1):
    S = [nodes[i] for i in range(V) if s & (1 << i)]
    cut_x = [x[(u, v)] for u, v in G.edges if (u in S and v not in S) or (v in S and u not in S)]
    x_cutsum = gp.quicksum(cut_x)
    model.addConstr(x_cutsum >= 1)

model.optimize()

カットセット不等式は頂点数の指数本になるため、頂点数が20~30程度で限界がくる。

よって、カットセット不等式を逐次追加していくことを考える。
方針は次のようにする。

  1. ある頂点 s \in Vを定める。
  2. 分枝限定法で暫定解が得られたときに、 S \subset Vを、得られた暫定解によって sから到達可能な頂点集合とする。
  3.  V = Sならば、何もしない。そうでなければ、 Sに対応するカットセット不等式を追加する。

この実装は次のようになる。深さ優先探索で、到達可能な頂点集合を求めている。
model.cbLazyで制約を追加するためには、事前にパラメータ「LazyConstraints」を1に設定する必要があることに注意。

model = gp.Model()
model.setParam("LazyConstraints", 1)
x = model.addVars(G.edges, vtype=GRB.BINARY)

model.setObjective(x.prod(d))

model.addConstr(gp.quicksum(x) == V - 1)

def callback(model, where):
    if where == GRB.Callback.MIPSOL:
        x_cur = model.cbGetSolution(x)
        S = [(0, 0)]
        que = [(0, 0)]
        while len(que) > 0:
            v = que.pop()
            for u in G[v]:
                if u in S:
                    continue
                if ((u, v) in x.keys() and x_cur[(u, v)] >= 0.5) \
                        or ((v, u) in x.keys() and x_cur[(v, u)] >= 0.5):
                    S.append(u)
                    que.append(u)

        if len(S) < V:
            cut_x = [x[(u, v)] for u, v in G.edges if (u in S and v not in S) or (v in S and u not in S)]
            x_cutsum = gp.quicksum(cut_x)
            model.cbLazy(x_cutsum >= 1)


model.optimize(callback)

頂点数を64(8×8)に設定したときのGurobiのログの抜粋が以下のようになる。

ログを見ると、255本のカットセット不等式が追加されていることがわかる。

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0 1703.50000    0    4          - 1703.50000      -     -    0s
     0     0 1715.50000    0    4          - 1715.50000      -     -    0s
     0     2 1726.50000    0    8          - 1726.50000      -     -    0s
* 2910   462              13    1861.0000000 1831.00000  1.61%   2.3    0s
* 2923   462               9    1842.0000000 1831.00000  0.60%   2.3    0s
* 3068   464              11    1834.0000000 1833.00000  0.05%   2.3    0s
* 3541     0              12    1833.0000000 1833.00000  0.00%   2.4    0s

Cutting planes:
  Lazy constraints: 255

カットセット不等式は 2^{64} - 2本となるため、本質的に必要な制約が少ないことがわかる。
また、クラスカル法でもとめた重みと一致していることも確認できるため、最小全域木として必要な制約は追加できていることもわかる。

ただし、これでもMiller-Tucker-Zemlin制約を用いた定式化のような本質的に制約の本数が少ない定式化に比べると劣ってしまうため、定式化が難しい場合の最終手段ととらえるべきだろう。

CSVの文字コード問題

数理最適化を含むデータ分析するシステムを作る際に、いつも悩むのが入出力データの形式だ。
選択肢としては、次の3つがあるだろう。

ユーザーが扱うことをだけを考えると、最もなじみがあるExcelファイルが無難だろう。
ただし、Excelファイルはバイナリファイルであるため、Git管理したときに差分が表示されない。基本的にGit管理するファイルにバイナリファイルは含めたくない。
CIで自動テストをすることを考えると、テストデータをGit管理する必要がある。コードレビューするためにブランチを取ってきてExcelで開くのは非常にめんどくさい。

そうなるとCSVファイル(文字コード:UTF-8)が第2候補に挙がってくる。
CSVファイルはテキストファイルであるため、Git管理との相性もよく、プログラムとの親和性も高い。

ただし、ユーザーがデータをいじろうとすると問題が発生する。基本的にユーザーのPCにメモ帳以外のエディターはインストールされていないので、CSVファイルを渡すとExcelで開こうとするのだが、なんとExcelはBOMのないCSVファイルを受け取るとデフォルトでCP932でエンコードしてくる。しかもPowerQueryなどでない限り文字コードを指定することもできない。(やはりExcelはクソ)
ユーザーにメモ帳でデータをいじらせるのは流石に酷である。

最終手段として、もともとCP932でエンコーディングされたCSVファイルを使うという発想になってくる。こうなると、プロジェクト内で入出力データだけがCP932になるという非常に奇怪な状態になる。文字コードデファクトスタンダードUTF-8とされている現状でCP932のファイルをわざわざ使うのもはばかられる。
Excelでの編集の問題も解決しておらず、CSVで開き、編集して保存すると、Excelが数値や時刻だと認識したものが勝手にフォーマットが変換される。(やはりExcelは滅ぼさなければならない)
これでもユーザーがメモ帳でCSVデータをいじるという苦行が発生する可能性がある。

そうなると結局Excelファイルで、となって堂々巡りになってしまう。この問題は常に頭痛の種だ。

結局、WindowsがさっさとCP932を廃止して、ExcelUTF-8をデフォルトにしてくれれば問題の大部分が解決する。それでも、もし将来UTF-8は古いと言われるようになってしまったら、また同じ問題が起こるのだろうか。

線形順序付け問題と貪欲法

線形順序付け問題は、行列 A \in \mathbb{R}^{n \times n}_{\geq 0}が与えられたときに、置換 \sigma \in S_nについて、目的関数


 \displaystyle
 f(\sigma) = \sum_{i < j} A_{\sigma(i)\sigma(j)}
を最大化する問題として表される。もし i > jについて、 A_{\sigma(i)\sigma(j)} = 0となるような \sigmaが取れる場合は、DAG上のトポロジカルソートになる。

線形順序付け問題はNP困難であるため、厳密解を求めたければ、整数計画として定式化してソルバーに投げることになるだろう。

厳密解がいらない場合では、ヒューリスティクスも数多く提案されている。ここでは単純な解法として、貪欲法に基づく1/2-近似アルゴリズムを述べておく。

  1.  S = \{1, 2, ..., n\}, \sigma(i) = 0\, (i = 1, 2, ..., n)と初期化しておく。
  2.  i = 1, 2, ..., nに対して以下を繰り返す。
    1.  \sum_{j \in S \setminus \{k\}} \left( A_{kj} - A_{jk} \right)が最大となるような kをとる。
    2.  \sigma(i) = k, S = S \setminus \{k\}とする。
  3.  \sigmaを出力して終了。

各ステップにおいて、 \sum_{j \in S \setminus \{k\}} \left( A_{kj} - A_{jk} \right) \geq 0なる kが少なくとも1つは存在することが背理法によって示せる。(左辺を k \in Sに渡って和を取ると0になる。)

よって、任意の iについて、


 \displaystyle
 \sum_{j = i + 1}^{n} A_{\sigma(i)\sigma(j)} \geq \sum_{j = i + 1}^{n} A_{\sigma(j)\sigma(i)}
が成立し、 iについて両辺和をとると、

 \displaystyle
 f(\sigma) = \sum_{i < j} A_{\sigma(i)\sigma(j)} \geq \sum_{j < i} A_{\sigma(i)\sigma(j)}
が成立することがわかる。

明らかな最適値の上界として、


 \displaystyle
 \sum_{i \neq j} A_{ij} = \sum_{i < j} A_{\sigma(i)\sigma(j)} + \sum_{j < i} A_{\sigma(i)\sigma(j)}
が与えられることから、1/2-近似アルゴリズムであることがわかる。

計算量は \mathcal{O}(n^2)になる。

より精度がよいアルゴリズムヒューリスティクスなどはちゃんとしたサーベイを見てほしい。