MIPソルバーで求める”最適な”海外旅行計画 ~NY旅行の実例~

おことわり

これは真面目に書いた記事ですが、文献調査はあまりしてませんし数理最適化にも詳しくないので、内容の正確性は保証できません。ただ、タイトルのことをやってみたという記事です。なお、非常に初歩的な内容ですので特にコード等は含みません。

はじめに

最適な海外旅行計画を立てることは、多くのホモサピエンスにとっては困難であることが知られています。往復の航空便・宿泊先のホテル・訪問する場所などの旅行で決めるべき要素は数多く互いに依存しているため、ホモサピエンスが使える単純な手法では良い旅行計画が得られる保証はありません。なので、大体の場合はどこかで妥協せざるを得ないという悲しい状況があります。

そこで苦労せずとも良い旅行をするために、ヒューリスティクスを用いない解法を考えます（わざわざ旅行のためだけに解法は考えたくない）。ナイーブな発想では日帰り旅行では単にTSPTW（時間制約つき巡回セールスマン問題）として解けそうです。これをN泊に拡張し訪問先候補の選好と営業時間を導入して、海外旅行計画を求めてみて、実際に旅行に行きます。なお、N泊への拡張方法は、ポテンシャルをホテル出発時にリセットすることで実現しています。

問題設定

海外旅行計画は、次の手順で計画されるとします：

1. 旅行先（空港）の決定
2. 滞在先ホテルの決定
3. 訪問場所と訪問時刻の決定（狭義の旅行計画）

ここでMIPソルバーを用いる問題は、3の狭義の旅行計画です。仮に、2について、複数ホテルのうちどれが最も効率的に観光先を回れるかといった検討を行いたい場合は、それぞれのホテルについて、3を行ってその結果を比較することでおこないます。

また情報収集の方法としては、1,2はExpediaなどの比較サイト、3はInstagramなどのSNSやまとめサイトを想定しています。さらに3については、訪問先候補の選好度・営業時間・最小/最大滞在時間も適当に設定します。

解き方

単にMIPソルバーを叩くだけです。ただし、そのままでは非現実な解をはき出すので旅行者のフィードバックを何度か得てモデルを改良する必要があります。フィードバックの実例としては、一日に軽食が偏ってる/到着日の夜までバックパックを持ち歩きたくないから一度夕方にホテルに寄りたいなどです。これらをいい感じに定式化し、問題を繰り返し解き直す必要があります。

このフィードバックサイクルを10-20回ほどやらないと、単に最適化問題として解いた旅行計画は使いものになりませんでした。そのため往路の機内にラップトップを持ち込み、同行者と話し合いながら何度も計算を回すことをおすすめします。最終的にみんなが納得できればOKです。ただし、事前予約が必要なものについては予めある程度ちゃんと解いて予約しておきましょう。

移動方法はtransit（公共交通）かdrive（自動車）です。これらの移動時間はgoogle maps apiから取得します。またdriveはuberかlyftを使うことになります。これらはpickupまで時間がかかることがありますので、余裕を持って10-15分くらい時間時間を多めに見積もります。

具体的な制約の定式化については、特に触れませんが，いい感じの連立方程式/不等式を立てればだいたいOKです。目的関数は、選好や門限ペナルティなどいろいろ適当にフィードバックを反映しつつ設定しています（解法を考えない代わりに目的関数を頑張って考える）。

NY旅行の実例

去年と今年の2度、私的な旅行で上記の手法を実際に用いましたのでそれらについて記します。旅行先はロサンゼルス（LAX）とニューヨーク（JFK）で、どちらも3泊4日のスケジュールです（出発地は日本ではないため、3泊6日ではない）。

ここでは、後者の旅行計画を簡単に紹介します。なお実現可能性を疑うホモサピエンスも多いと思いますが、筆者と同行者が実際にこの計画の下で、大きな遅延なくすべての場所をちゃんと回って来ました。また、下記の旅行計画は我々の訪問先の候補と選好を組み込んだモデル内での最適解です。真の意味での最適解である保証はありませんのでご注意ください。ちなみに、この実例での問題規模は、変数・制約ともに数千程度です。

​​Day 1

13:03 に JFK Terminal 4 (A) を出発
13:43 に Supermoon Bakehouse (1) に到着
14:13 に Supermoon Bakehouse (1) を出発
14:25 に Popbar (2) に到着
14:53 に Popbar (2) を出発
15:04 に チェルシー・マーケット (3) に到着
16:04 に チェルシー・マーケット (3) を出発
16:05 に ザ・ハイ・ライン (4) に到着
17:05 に ザ・ハイ・ライン (4) を出発
17:27 に ホテル (5) に到着
17:57 に ホテル (5) を出発
18:07 に タイムズスクエア (6) に到着
18:37 に タイムズスクエア (6) を出発
18:40 に マジェスティック劇場 (7) に到着
22:25 に マジェスティック劇場 (7) を出発
22:40 に ホテル (H) に到着

​​Day 2

9:03 に ホテル (H) を出発
9:21 に バッテリー・パーク (9) に到着
12:21 に バッテリー・パーク (9) を出発
12:29 に ブルックリン橋 (10) に到着
13:37 に ブルックリン橋 (10) を出発
13:43 に OddFellows Ice Cream Co. DUMBO (11) に到着
14:13 に OddFellows Ice Cream Co. DUMBO (11) を出発
14:28 に Morgenstern's Finest Ice Cream (12) に到着
14:48 に Morgenstern's Finest Ice Cream (12) を出発
15:00 に World Trade Center (13) に到着
17:00 に World Trade Center (13) を出発
17:11 に LROOM CAFE (14) に到着
17:41 に LROOM CAFE (14) を出発
17:51 に シェイク シャック (15) に到着
19:31 に シェイク シャック (15) を出発
19:43 に ホテル (H) に到着

​​Day 3

9:00 に ホテル (H) を出発
9:04 に エッサ ベーグル (17) に到着
9:35 に エッサ ベーグル (17) を出発
9:44に 5th Ave (18) に到着
13:31 に 5th Ave (18) を出発
13:44 に アメリカ自然史博物館 (19) に到着
15:59 に アメリカ自然史博物館 (19) を出発
16:02 に セントラル・パーク (20) に到着
17:06 に セントラル・パーク (20) を出発
17:19 に Cha Cha Matcha (21) に到着
17:34 に Cha Cha Matcha (21) を出発
17:40 に エンパイア・ステート・ビル (22) に到着
19:40 に エンパイア・ステート・ビル (22) を出発
19:52 に ホテル (H) に到着

​​Day 4

9:03 に ホテル (H) を出発
9:18 に Sarabeth's (24) に到着
10:40 に Sarabeth's (24) を出発
10:48 に ニューヨーク近代美術館 (25) に到着
13:03 に ニューヨーク近代美術館 (25) を出発
13:40 に JFK Terminal 4 (A) に到着

機内で1時間みっちりとフィードバックと解き直しを繰り返したので、（この計画をみても分からないと良し悪しが思いますが）個人的には完璧な計画が出来上がりました。しかしながら、この計画通りに旅行すると、濃すぎて疲れ、ホテルですぐに爆睡してしまいました。まあ、これが最適ということなのでしょう。本当に濃い旅行でした。

おわりに

MIPソルバーという、組合せ問題でのチートツールを用いることで満足度の高い旅行ができました。これならアウストラロピテクスでもできそうですので、誰にでもおすすめできます。

また、このような旅行に付き合ってくれた同行者には、深く感謝します。

なぜCatboostの推論は速いの？

前回の記事｢AutoML v.s. Catboost｣に出てくるCatboostは、XGBoostやLightGBMと比べて30-60倍も推論が速いという特徴があります。

推論時間は、kaggleなどのコンペでは推論に時間をかけられるのであまり気にしませんが、実サービスとなると重要ですよね。

推論時間の比較

以下のグラフは、3大GBDTライブラリでの推論時間比較です。Catboostがずば抜けて速いことがやかります。 そして学習時間の速さに定評があるLightGBMは、なんと最遅です。

この推論時間の速さは、CatboostがGBDTのベースとなる決定木に symmetric tree と呼ばれる、特殊な形の木を使っているからです。

ここで、symmetric treeとは以下の図の右側のような木です。左側は普通の決定木です。

なぜsymmetric treeは速いか

｢同一の深さではすべての分岐条件が同じ｣ということが推論の並列化を可能としていて、これが推論時間の短縮を実現しています。

ちなみに、この木の構造が気持ち悪いという人はgrow_policyというパラメータを変えることで、変更できます（公式ドキュメントに詳細あり）。

---

おまけ：学習時間とGPU

CatboostはXGBoostと変わらない速度だからLightGBM使っとけ的な意見をよく耳にします。しかしそれは、CPUを用いた場合の話で、GPUを用いた場合は、学習時間でもCatboostが最速です。以下のグラフは、GPUでの学習時間を表しています。

並列化に向いたsymmetric treeと並列計算が得意なGPUの相性は抜群ということです。

皆さんも是非一度、deep learning以外でcolabのGPUを使い倒してみては如何でしょうか。

追記：GPU学習時間の検証

実際にcolabで、GPUを用いてLightGBMとCatboostの学習時間を計測しました。条件等は、下記のbinningの追記と同様です

公式の主張ほど学習時間の差はありませんが、それでもGPUがあればLightGBMより速いと言えそうです。

ただ、公式的にはデータセットがもっと大きいときに分散学習できるのが強みだと言っているので、そのような条件ではもっと差がつくのかもしれません。

レファレンス

図やグラフは、すべて以下の動画からの引用です。

CatBoost - the new generation of gradient boosting - Anna Veronika Dorogush - YouTube

この動画では、Yandexの開発者が直々にCatboostを初心者向けに解説してますので、このブログよりはるかに分かりやすいと思います。是非ご覧ください（ただし英語）。

また、各種時間の計測にかかる実験の条件等も説明されていますので、そちらに興味がある方もご覧ください。

---

追記：binningの時間も求めつつ実際に実験

一瞬だけbinningにかかる時間を考慮して実験してみてというツイートを観測したので，やってみた結果を書きます。

条件

• colabのGPUを使用
• sklearnのmake_classificationで特徴量次元128の二値分類データセット作成
• binningにかかる時間は，max_depth=1, n_rounds=1としたときの推論時間として近似的に求める*1
• 推論時間は，max_depth=6, n_rounds=1として学習させたモデルの推論時間。binningにかかる時間も含まれる

結果

左：binning時間・右：推論時間

• binningのみはLightGBMの方がCatboostよりも速い
• でも，推論全体にかかる時間はCatboostの方がLightGBMよりも速い
• ただし，symmetric treeがよく効くようにある程度max_depthが深くないとLightGBMが勝つ
• Catboostが主張してるほど推論は速くない（もっと適切な条件があるのかもしれないが）がLightGBMよりは速い

追追記：binning時間を近似なしで求める

ちゃんとLightGBMのDatasetをconstructしてbinningだけの時間を求めてみました（Catboostは遅延評価ではないっぽい）。すると，binningだけでもCatboostの方が速いという結果となりました。

exactなbinningにかかる時間

以下，考察です（コードを読んでないので間違ってるかもしれません）。

• Catboostは，推論時にもbinningを完全におこなう
• 一方LightGBMは，推論時には分岐条件に含まれるものだけbinningする等して効率化している
  • max_depth=1, n_rounds=1のときの推論時間（=binningの近似時間）の異様な速さに説明が付く

アプデ

本記事を受けて、中の人がコメントしてくれています。

We have fast inference because we are using special tree shape - oblivious trees, they have the same condition on one level, so it could be done with SIMD instructions. We are now making GPU model inference backend, wich is even faster on large models. Stay tuned!

— Stas Kirillov (@kizillof) 2019年6月5日

*1:正しい近似なのか自信ないので、おかしかったら突っ込んで下さい

AutoML v.s. Catboost (colab有)

前説

tjo.hatenablog.com

今話題のAutoMLは，有名な"データサイエンティスト"よりも優秀なようです。

しかし，この比較にはmamas16k氏に盛大に突っ込まれてるように

roundsを！！！チューニングしない！！！

— まます (@mamas16k) 2019年5月23日

という問題があります*1。これでは，天下のGoogle様が開発されたAutoMLが役不足ではないでしょうか。

チューニングをミスっているXgboost相手でデモンストレーションするのではなく，どうせならXgboost・LightGBMに続く最新の第3世代*2GBDT系のMLライブラリ Catboost を相手にやりましょう！

スコア比較（低い方が良い）

AutoML:0.8634979

TJO氏のブログからの引用です。

rmse(pred$shares, test$shares)
[1] 0.8634979

Xgboostでも勝てなかったこのAutoMLのスコアに，果たしてCatboostは敵うのでしょうか・・・

Catboost:0.8625799

あ，勝っちゃいましたね。

pred = model.predict(test_X)
np.sqrt(mean_squared_error(pred, test_y))
> 0.8625798651724433

でもこれ，前処理としてダミー変数をカテゴリカル変数に戻して*3，後はCatBoostRegressorに突っ込んだだけです。 ハイパラもoptunaやベイズなんちゃらは使わずに，完全感覚チューニングでlr=0.01, depth=8, l2_leaf_reg=30と適当です。

詳しくは以下のcolabを見て下さい。

colab.research.google.com

また，学習には20分弱３０分で，roundsが50006000くらいです。この圧倒的速さを最近流行の棒グラフで表現しました*4。

AutoMLとCatboostの学習時間の比較

加えて先のAutoMLは12時間・一桁万円課金の結果らしいので，速度・精度・カネの３観点を総合してCatboostの圧倒的勝利といったところではないでしょうか。（今回のデータセットにかぎりますが）

---

追記：random seedを固定して学習しなおしたところ学習時間・roundsが増加しましたので訂正しました。またグラフも訂正しています。

---

追追記：DS界の井上尚弥ことonodera氏がワンパンしました。さすがGMです。

炎上マーケに釣られてみた。ノールックで0.861は出るので0.8705はちょっとアレhttps://t.co/KEwTF7OD4J

— onodera (@0verfit) 2019年5月24日

---

追追追記

この記事の動機と成果を、一部の人向けに書いておきます。

某JO氏はノーダメージ説を観測しましたが，プロダクトのポジショニング変更を余儀なくされているのでG社マーケ的にはダメージなはずです。

「専門家いらずで上出来」→「だれでもそこそこの精度」→「検証・改善にMLエンジニアやDSが若干名必要」

※若干名いればAutoMLの精度上回るから意義なくね？ pic.twitter.com/bW6eY8smpR

— atksh (@atksh1) 2019年5月24日

私はただ「専門家いらずで上出来」という主張を見て言い過ぎなのではと思い検証して記事を書いただけです。

— atksh (@atksh1) 2019年5月24日

実は中にいても完全な詳細は知らされていないので、今後どうなるかはPMへのリクエスト次第という部分もあります（今回バズったのはPMも気付いたらしく問い合わせのメールが飛んできました笑）

— TJO (@TJO_datasci) 2019年5月24日

中の人にご迷惑をおかけしたみたいですね。

*1:指摘に気づいてやり直したみたいですが，なぜかスコアは悪化しています。

*2:勝手に名づけた

*3:Catboostのドキュメントがonehotにはしないでねって言ってるのでそれだけやりました。

*4:厳密には、私がおこなったチューニング時間を足す必要がありますが。

PyTorchのMulti-GPUでメモリ使用量の偏りと改善策

よくあるGPUメモリの偏り

PyTorchでマルチGPUをすると， こんな感じで １つのGPUのメモリだけ他に比べて多い （２倍以上）場合って OOMでバッチサイズが増やせずに悲しいですよね。

解決法

以下の方法で，8 GPUでOOMしないバッチサイズが16->32になりました。

◯forwardにlossの計算までを入れる

i.e.,

class Model(nn.Module):
  def __init__(self, d_in, n_class):
    super().__init__()
    self. n_class = n_class
    self.fc = nn.Linear(d_in, n_class)

  def forward(self, x, labels=None):
    x = self.fc(x)
    if labels:
      ### ここが重要！！！ ###
      return x, F.cross_entropy(x.view(-1, self.n_class), labels.view(-1))
    else:
      return x

model = DataParallel(Model(100, 10))
for x, labels in dataloader:
  pred, loss = model(x, labels)
  optimizer.zero_grad()
  loss.backward()
  optimizer.step()

criteriaを使ってloss計算すると，メインGPUのみでlossの計算が おこなわれることが，先のメモリ使用量の偏りを引き起こしているらしいです*1。

よく見てみたら，BERT等のPyTorch実装も同様のforward実装になってました。最近のベストプラクティスなのかもしれない？

*1:https://discuss.pytorch.org/t/dataparallel-imbalanced-memory-usage/22551/26

PyMC3でWBIC

​​◯タイトル通り，PyMC3でWBICを求めてみました。 なお，WAICはpymc3.stats.waicで求められるので*1，やっていません。

◯元ネタは，以下の記事です。 RのstanでやられていたのをPythonのPyMC3に移植してみました。

statmodeling.hatenablog.com

◯私が頑張った*2ポイントとしては，以下の３点です。（詳しくはコードを見て下さい）

• 逆温度をもつ分布が標準実装されてないので，Custom Distributionで作った
• observed=Yと一括で指定することはせずに，データ数だけ分布を用意してobserved=Y[i]としたか
• WBICにデータごとの対数尤度が必要になるので，pm.Deterministicでobs[i].logpt / betaして無理矢理取ってくる

◯実際に用いたコードは文末に添付しますので参考にしてください。実行すると

model1｜WBIC = 201.884
model2｜WBIC = 194.659

のような結果が得られ，ちゃんと混合正規分布のmodel2が選べています。

◯結果の信憑性については，先のstanでやられていた方の結果と整合しました。 が，間違っているところやベストプラクティスがあれば是非コメントしてください。

PyMC3で混合正規分布の推定。正規分布と混合正規分布の２つを用いて，WBICでモデル選択

※モデル内のy_predは，事後分布を可視化する用に突っ込んでいます。WBICの計算には不要です。

---

追記

◯このツイートでアクセス数が急増しました。 普段フォローしている方のRTもあり，びっくりしています。しかし，ちょっと嬉しいですね。言及をいただいたこと，この場でお礼申し上げます。

"PyMC3でWBIC - 簡単なレポート" https://t.co/Cp4fVGSEim

— Kentaro Matsuura (@hankagosa) 2019年2月6日

*1:公式Example｜https://docs.pymc.io/notebooks/model_comparison.html

*2:どう実装するかわからず，無能すぎて2日も悩んでしまいました。

Model Based k-Meansで直線への割付を推定

やったことを，とりあえずグラフを載せます。

↑のように，2つの直線にノイズがのったデータから， 元の直線の推定とそれに基づいた分類をModel Based k-Meansでやってみました。

（元の直線を表示するの忘れましたが，求めています。）

Model Based k-Meansって何？

The Top Ten Algorithms in Data Miningによると

Another broad generalization is to view the “means” as probabilistic models instead of points in R^d. Here, in the assignment step, each data point is assigned to the model most likely to have generated it. In the “relocation” step, the model parameters are updated to best fit the assigned datasets.

とあります。 ざっくりというと「通常のk-Meansの割付ステップで最も 近いクラスタを選ぶ代わりに，その近さを 確率で測ってクラスタを割付け， 更新ステップでは，クラスタ 重心の更新の代わりに，モデルの パラメタ更新を行う」ということです。

アルゴリズム

コードは載せるまでもないと思いますので，アニメーションで示します。

「初期化 → 割付（一番誤差が少ない直線へ）→ 更新（切片と傾きの推定）→ 割付 → ･･･ → 更新 → 割付」 の初めの5ステップ分です。

また，割付については，確率の代わりに誤差を用いています。 つまり一番近い（yの差が小さい）直線に割り付けられるようにしています。

拡張可能性

1. 割付をソフトにしてデータのすべてを重み付けで学習
2. モデルを非線形*1にして表現力↑
3. ソフトな割付を別のモデルで学習して，予測のアンサンブルの重みに

3つ目に付いてですが，kaggleなどでアンサンブルをよくやりますが，通常そのウエイトをサンプルによらずに固定するところ， これが動的に変えられるようになると，精度あがったりするのかなと思います。

つまり，各データごとに適切なモデルを選ぶことができる可能性*2があるということです。

*1:GBDTやNNなど

*2:もちろん過学習の問題も大きくなりますが