Apache2.0 ライセンスとなった数理最適化ソルバー SCIP と Cbc を手持ちの問題で比較する

2022/12/21に公開

今年の11月に開催された数理最適化ソルバーの SCIP の ワークショップ において、 SCIPのライセンスが Apache 2.0 になることが発表されました。 SCIP は元々 ZIB という独自のライセンスで、学術利用でのみフリーで利用可能でしたが、今回のライセンス改定により、無償での商用利用も可能となります。アナウンスのあった11月の時点では、あくまで次期バージョンからの改訂という話だったのですが、先日その 8.0.3 がリリースされ、名実ともにApache 2.0 ライセンスとなったようです。

https://www.scipopt.org/index.php#license

SCIP の他にもオープンソースな数理最適化ソルバーはいくつかあります[1]が、商用のもの[2]と比較すると、お世辞にも性能が良いとは言えない状況でした。そのため、解ける問題のサイズも限られたものになってしまい、数理最適化が普及しにくい要因の一つとなっていました。今回のライセンス改定により、オープンソースの範囲内で解ける問題サイズが増えそうです。

そこで、本記事では、最近自分が定式化して遊んだ問題を題材にして、どれくらい性能が上がるのかを見てみたいと思います。

SCIP のベンチマーク性能

独自の問題を解いて比較する前に、標準的なベンチマーク問題で複数のソルバーの性能を比較した結果を見てみましょう。
Cbc と比較して SCIP は頑張っていますが、Gurobi とは大きな差があることがわかります。今回目指せるのもこの辺りになりそうです。

https://mattmilten.github.io/mittelmann-plots/
http://plato.asu.edu/ftp/milp.html

SCIP のインストール

SCIP のバイナリは、ダウンロードページ から入手できます。
ただし、以下の記述がある通り、環境が Mac の場合は libgfortran などを用意する必要があります。

For some of the mac packages you need the libgfortran and libtbb libraries, and bison. You can install them for example via homebrew: brew install bison, brew install gcc and brew install tbb before installing the scipoptsuite.

libgfortran は gcc パッケージに含まれていますが、対応している gcc のバージョンは 7 か 10 と少し古いものになっていて、M1 Mac で brew を使ってインストールしようとすると、エラーとなってしまいます。

そこで、今回はソースコードからインストールすることにしました。といっても、インストール自体はとても簡単で、以下のコードで一発でした

# Boostがないとおこられる
brew install boost

# SOPLEXのインストール
git clone https://github.com/scipopt/soplex.git
cd soplex
mkdir build
cd build
cmake ..
make
sudo make install

# SCIPのソースコードをclone(ここでは master ブランチを利用)
cd ../..
git clone https://github.com/scipopt/scip.git

# ビルドディレクトリ作成
cd scip
mkdir build
cd build

# ビルド. DAUTOBUILD=ON で足りない独自ライブラリを探してダウンロードしてくれる
cmake .. -DAUTOBUILD=ON
make

# システムにインストール
sudo make install

インストール後に別ターミナルを立ち上げて scip と打ち込み、プロンプトが表示されることが確認できれば、インストールは完了です。

比較1: 披露宴座席配置問題

数理最適化 アドベントカレンダー 2022 9日目の記事で、披露宴の座席配置を数理最適化を用いて決定した話が紹介されていました。その記事の中で高々20人であっても数理最適化ソルバーを使って厳密解を見つけるのは難しい点に言及されていました。

https://zenn.dev/akira_t/articles/seat-opt-via-gw

その話を受けて、数理最適化ソルバーを使ってどれくらいの問題を解くことができるか、現実的にはどういう戦法を取るのかという記事も公開されました。

https://qiita.com/iwanaga-jiro/items/361976589916fdcd1a23
https://qiita.com/iwanaga-jiro/items/96a9898d4790fcd14d3a

実は自分も「自分だったらどうするかなぁ」ということを考えて定式化していたので、この問題を取り上げたいと思います。

解きたい問題

座席間の近さ C_{a,b} と披露宴のゲスト同士の親しさ S_{i,j} が与えられた時に、知り合い同士をできるだけ近い座席に配置したい

定式化

定式化は上記の記事のものとよく似ていますが[3]、論理積と \min が等しくなる性質を使っています。

x y x\land y \min(x,y)
0 0 0 0
0 1 0 0
1 0 0 0
1 1 1 1

定数

  • S_{i,j}: ijの知り合い度合い
  • C_{a,b}: 座席abの席の近さ

変数

  • x_{i,a} \in \{0, 1\}
    • ゲストiが座席aにアサインされたときに1となるフラグ
  • z_{i,j,a,b} \in \{0, 1\},\ \mathrm{for}\ i < j
    • ゲストiが座席a、ゲストjが座席bにアサインされたときにのみ1となる補助変数
    • 言い換えると z_{i,j,a,b} = \min(x_{a,i},x_{j,b})

目的関数

ゲスト i が座席 a、ゲスト j が座席 b に座ったときの好ましさS_{i,j}C_{a,b} と表せるとして、以下のようにします

  • \mathrm{maximize} \sum_{i,j,a,b}z_{i,j,a,b}S_{i,j}C_{a,b}

制約条件

  • \sum_{a} x_{i,a} = 1 (各ゲストに1つ座席が与えられる)
  • \sum_{i} x_{i,a} = 1 (全ての座席がうまる)
  • z_{i,j,a,b} \leq x_{i, a} (x_{i,a}0 のとき、 z_{i,j,a,b}0)
  • z_{i,j,a,b} \leq x_{j, b} (x_{j,b}0 のとき、 z_{i,j,a,b}0)

上記に加えて、 x_{i,a}=1 かつ x_{j, b}=1 の時に、 z_{i,j,a,b} = 1 という条件が必要そうに見えますが、目的関数を見ると z_{i,j,a,b} は勝手に大きくなろうとしてくれるので不要です。

実装

model = pulp.LpProblem('wedding_seats_assignment', sense=pulp.const.LpMaximize)


x = np.array(
    pulp.LpVariable.matrix('x', (guest_names, seats), cat=pulp.LpBinary),
)

z = np.array(
    pulp.LpVariable.matrix(
        'z', (guest_names, guest_names, seats, seats), cat=pulp.LpBinary)
)
for i in range(x.shape[0]):
    for j in range(x.shape[0]):
        if i >= j:
            z[i, j] = 0

z = np.tril(z.T).T
np.fill_diagonal(z, 0)

# 目的関数の設定
z_s = np.tensordot(guest_distances, z, axes=([0, 1], [0, 1]))
z_s_c = np.tensordot(z_s, seat_distances, axes=([0, 1], [0, 1])).sum()

model += z_s_c

# 制約条件の追加
for i in range(x.shape[0]):
    model += pulp.lpSum(x[i, :]) == 1

for a in range(x.shape[1]):
    model += pulp.lpSum(x[:, a]) == 1

for i in range(x.shape[0]):
    for j in range(x.shape[0]):
        if i >= j:
            continue
        for a in range(x.shape[1]):
            for b in range(x.shape[1]):
                model += z[i, j, a, b] <= x[i, a]
                model += z[i, j, a, b] <= x[j, b]

model.solve(SOLVER)

比較結果

前述の通りゲストが20人いると単純に解くことはできないので、ゲスト数 N を変えて実験しました

N=10 N=11 N=12 N=13
Cbc 157.53sec 233.13sec 711.03sec > 6000sec
SCIP 10.23sec 12.82sec 27.22sec 166.24sec

ざっくりと5倍から10倍10倍以上高速化していそうです[3:1]。また、この手の問題はサイズが大きくなると急激に時間がかかるようになるものですが、その「急激に時間がかかるようになる」サイズも SCIP のほうが大きそうです。
なお Cbcでは N=13 でまったく歯が立たなくなっていますが、ログを見ると 384sec 程度で最適解は見つけられているが最適解である証明に時間がかかっているようです。なので、最適解(である保証)をあきらめればもっと短時間で解が得られます。

比較2: 最長しりとり問題

今年のアドベントカレンダーでさだまさしさんの曲タイトルでしりとりをしました。

https://zenn.dev/ohtaman/articles/sada_shiritori

上記の記事のポテンシャル最適化のコードをSCIPで動かすとどうなるかをみてみようと思います

解きたい問題

さださんの楽曲(ネット上でタイトルを確認できた 572曲)のタイトルでしりとりをした場合に最大でどれくらい続けられるかを求める

定式化

以下のように定式化しました。詳細は記事をご確認ください。

定数

  • N\in \mathbb{N}: タイトル数

変数

  • x_{ij} \in \{0, 1\}: タイトル i からタイトル j につなげる時は 1、つなげない場合は 0 となる変数
  • z_{i} \in \{0,\dots,N\}: タイトル i の訪問順

目的関数

曲と曲がつながった数、つまり x_{ij} = 1 の数を増やしたいので、以下のようにします。

  • \mathrm{maximize} \sum_{ij}{x_{ij}}

制約条件

  • \sum_{J} x_{JI} \leq 1 \quad \mathrm{for}\ \forall I (各タイトルは最大1回だけ利用可能)
  • \sum_{J}x_{IJ} \leq 1 \quad \mathrm{for}\ \forall I (あるタイトルにつづくタイトルは1つだけ)
  • z_{I} + 1 - N(1 - x_{Ij}) \leq z_{j}\ \forall I,\ j (部分巡回路除去)
  • x_{i,j} = 0 \ \forall (i, j) \notin C (しりとりにならない経路は利用できない)

比較結果

こちらもタイトル数を変えて実行時間を測ってみました

N=100 N=200 N=500 N=572
Cbc 0.38sec 3.83sec 54.34sec 46.06sec
SCIP 0.21sec 2.48sec 142.59sec 366.43sec

面白いことに、 N が大きいとき、今度は Cbc の方が早いという結果になりました。問題の構造によって得意不得意がありそうです。

まとめ

Apache 2.0 ライセンスになった SCIP のインストール方法の紹介と Cbc を、過去自分が書いたコードベースでの比較結果をまとめました

  1. SCIP のインストールはとても簡単でお手軽。PuLP であれば、そのままソルバーを切り替えられるので、お手軽に高速化できます
  2. ベンチマークを見る限り、Cbc よりも SCIP の方が高速です。ただし、すべての問題において一律に高速というわけではなないこともわかりました。PuLP を使っている限りではソルバーの切り替えは簡単なので、問題毎に比較して良さそうなものを選択するとよさそうです
脚注
  1. 最も有名なのは Cbc. 最近は HiGHs も話題にあがります. HiGHs については、ベンチマーク上 SCIP と同等か高速なようですので、時間があれば比較・紹介記事を書きたいと思います ↩︎

  2. たとえば Gurobi. 最近は COPT も頭角を表してきているという噂です. ↩︎

  3. この結果はどちらも並列化しておらず、シングルスレッドで動作しています。並列処理することで結果が変わる可能性があります ↩︎ ↩︎

Discussion