はじめに

AtCoder上で開催されたフューチャー社主催のヒューリスティックコンテストHACK TO THE FUTURE2022（HTTF）に参加した*1。問題等詳細は以下を参照。

この問題は大きく分けて各作業者のスキル推定と作業スケジューリングという2つの要素があるが、その中でもスキル推定の方に以前から気になってたリーマン多様体上の最適化なるものが使えそうだったので試してみたというのが今回の話。実のところ自分も多様体は何ぞやとか数学的なことはよく分かっていないが、今回行った単位超球面上の最急降下法であればそこまで難しくないし、以下の資料を参考に比較的簡単に実装できる。

http://www.orsj.or.jp/archive2/or60-9/or60_9_549.pdf

結果としては割とうまく推定できるが普通の最急降下法でやった場合とスコア的にはさほど差がなく、最終的にはスコアが高く出た普通の最急降下法を採用した（誤差範囲の気もする）。

以下、具体的行ったことを書いていく。

スキル推定問題の定式化

問題全体は上記のリンク先を見ていただくとして、ここではスキル推定問題について定式化をしておく。スキルは作業者間で独立なので、ここでは1人の作業者のスキルを推定する問題を定式化する。作業者は $K$ 種類のスキルを有しており、その $k$ 番目のスキルを $s_k \geq 0$ とする*2。また、タスクの要求スキルも同様に $K$ 種類あり、タスク $i$ の要求スキルを $d_{ik} \geq 0$ とする。そして、実際にその作業者がタスク $i$ を実行したときにかかった時間は $t_i$ であり、この時間は $\sum_{k = 1}^{K}\max(0, d_{ik} - s_k) = 0$ なら $t_i = 1$ 、 $\sum_{k = 1}^{K}\max(0, d_{ik} - s_k) > 0$ ならこの値に $-3 \sim 3$ の乱数を加えたもので得られる（この場合も $0$ 以下になるなら $t_i = 1$ とする）。 $n$ 個タスクが作業者に割り振られ、 $\left\{(\mathbf{d}_1, t_1), \cdots, (\mathbf{d}_n, t_n) \right\}$ が得られているので、この集合から作業者のスキル $\mathbf{s}$ を推定したい。

ということで、このスキル推定問題を2乗誤差最小化問題として以下のように定式化する。

$\displaystyle {\rm minimize} \ \ E(\mathbf{s}) = {\rm minimize} \ \ \frac{1}{2n} \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - |s_k|) \right\}^2 \tag{1}$

ここで制約付き最適化になるとややこしいので、 $s_k \geq 0$ という制約はつけず負値を取ることを許して、式の中で絶対値を取るようにしている。

通常の最急降下法でスキル推定

最終的に採用したのはこっちの方。式（1）の中にmaxやら絶対値やらが入っているので一部の人から怒られそうだが、気にせず $\mathbf{s}$ で微分して勾配を求める*3。

$\displaystyle \nabla E(\mathbf{s}) = \frac{1}{n} \mathbf{h} \odot \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - |s_k|) \right\} \odot \mathbf{g}_i \tag{2}$

ただし、 $\mathbf{h}$ と $\mathbf{g}_i$ は $K$ 次元のベクトルで、その $k$ 番目の要素をそれぞれ $h_k, g_{ik}$ とすると

$\displaystyle h_k = \left\{ \begin{array}{ll} 1 & (s_k > 0) \\ 0 & (s_k = 0) \\ -1 & (s_k < 0) \end{array} \right., \ \ \ g_{ik} = \left\{ \begin{array}{ll} 1 & (d_{ik} - |s_k| > 0) \\ 0 & (d_{ik} - |s_k| \leq 0) \end{array} \right.$

である。また、 $\odot$ はアダマール積を表す。スキル $\mathbf{s}$ はこの勾配を用いて、

$\displaystyle \mathbf{s}^{(l + 1)} = \mathbf{s}^{(l)} - \alpha \ \nabla E(\mathbf{s}^{(l)}) \tag{3}$

という更新式を適当な回数まわして推定した。本来なら収束判定を勾配のノルムとかで行うべきだが、手を抜いてループ回数を固定している。また、ステップ幅 $\alpha$ についても本来なら線形探索などで適切な値を求めるべきだが、これまた手を抜いて適当な値に固定している。

超球面上の最急降下法でスキル推定

ここからが本題。すぐ後に示すように式（1）を少し変形するとノルム $1$ 制約の最適化問題が出てくるので、そこに超球面上の最急降下法を使った。

スキルの成分分解

問題文をよく見ると真のスキルの生成方法が書かれている。これを見ると、 $\mathbf{s} = \frac{{\rm randdouble}(20, 60) | \mathbf{s}' |}{\parallel \mathbf{s}' \parallel_2} \ \ (\mathbf{s}' \sim \mathcal{N}(\mathbf{0}, \mathbf{I}))$ という式で生成されている。従って、ノルムが $1$ のベクトル $\mathbf{w}$ と実数 $q \ \ (20 \leq q \leq 60)$ を用いて、 $\mathbf{s} = q|\mathbf{w}|$ と書くことができる。これを使うと式（1）は

$\displaystyle \begin{align} {\rm minimize} \ \ E_{d}(\mathbf{s}) = {\rm minimize} &\ \ \frac{1}{2n} \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - q|w_k|) \right\}^2\\ {\rm s.t.} &\ \ \sum_{k=1}^K w_k^2 = 1 \\ & 20 \leq q \leq 60 \end{align} \tag{4}$

という $q$ と $\mathbf{w}$ に関する最適化に変換される。この式（4）の最適化のため、 $\mathbf{w}$ を固定して $q$ を最適化する処理と、 $q$ を固定して $\mathbf{w}$ を最適化する処理を交互に繰り返す方法を取った。

$q$ の最適化

$\mathbf{w}$ を固定すると、式（4）は $q$ について凸な関数になる。また、 $q$ は1次元の実数。ということで三分探索で最適解を求めた。ただ、タイムリミットの関係上三分探索のループは5回にしてるので、かなり粗い探索になっている。

$\mathbf{w}$ の最適化

$\mathbf{w}$ はノルムが $1$ という制約条件があり、これは $\mathbf{w}$ が単位超球面上に存在する制約だといえる。このような制約条件つき最適化問題に対しては様々なアルゴリズムが提案されているが*4、今回は冒頭で挙げた資料の通り超球面上での制約条件なしの最適化問題とみなす手法を実装した。

まずは一応式（4）の $\mathbf{w}$ に関する通常の勾配を書いておくと、

$\displaystyle \begin{align} \nabla E_{d}(\mathbf{w}) = \frac{q}{n} \mathbf{h}' \odot \sum_{i = 1}^n \left\{(t_i - 1) - \sum_{k = 1}^K \max(0, d_{ik} - q|w_k|) \right\} \odot \mathbf{g}_i' \end{align} \tag{5}$

ただし、

$\displaystyle h_k' = \left\{ \begin{array}{ll} 1 & (w_k > 0) \\ 0 & (w_k = 0) \\ -1 & (w_k < 0) \end{array} \right., \ \ \ g_{ik}' = \left\{ \begin{array}{ll} 1 & (d_{ik} - q|w_k| > 0) \\ 0 & (d_{ik} - q|w_k| \leq 0) \end{array} \right.$

となる。この勾配方向に式（3）のように動かしても、超球面のことは全く考慮されていないので明後日の方向に行ってしまう可能性がある。ということで、

なるべく超球面から離れないよう現在点における接平面に沿った方向に動かす
それでも超球面からはみ出るので超球面上に戻す（レトラクション）

ということを繰り返すのが当該のアルゴリズム。以下もう少し具体的に書く。

動かす方向

今超球面上の点 $\mathbf{w}^{(l)}$ にあって、そこから動かす方向を決める。 $\mathbf{w}^{(l)}$ における接平面の法線ベクトルは $\mathbf{w}^{(l)}$ である。この法線ベクトルが張る空間の直交補空間は接平面に平行な部分空間（接ベクトル空間）であり、 $\nabla E_{d}(\mathbf{w})$ のその部分空間へ直交射影したベクトルを $\mathbf{d}$ とすると、

$\displaystyle \mathbf{d} = \nabla E_d(\mathbf{w}^{(l)}) - \langle \mathbf{w}^{(l)}, \nabla E_d(\mathbf{w}^{(l)}) \rangle \mathbf{w}^{(l)} \tag{6}$

となる。これを用いて、

$\displaystyle \mathbf{w'}^{(l)} = \mathbf{w}^{(l)} - \alpha \mathbf{d} \tag{7}$

にまずは動く。こうすれば、とりあえずは接平面に沿った移動になるので、超球面から離れないという観点で言うと少なくとも $-\nabla E_d(\mathbf{w}^{(l)})$ 方向に移動するよりはマシになるはず、というのが直感的な理解。なお、 $\alpha$ はステップ幅でこちらも本来は（曲線）探索するべきだが実装では手を抜いて固定値にした。

超球面へのレトラクション

式（7）で得た $\mathbf{w'}^{(l)}$ は、それでもやはり超球面からははみ出ている。従ってこれを超球面上に戻す必要がある。この操作をレトラクションと呼ぶ。これは単純に $\mathbf{w'}^{(l)}$ のノルムで割ればよい。よって、

$\displaystyle \mathbf{w}^{(l + 1)} = \frac{\mathbf{w'}^{(l)}}{\parallel \mathbf{w'}^{(l)} \parallel_2} \tag{8}$

として $\mathbf{w}$ を更新できる。

結果

seed=0の時のビジュアライザ動画を貼りたかったけどはてなブログに動画がアップロードできなかったので、最後の状態を画面キャプチャした画像を貼る。

通常の最急降下法
成分分解+超球面上の最急降下法

うーん、どっちがいいとも悪いとも言いにくい。結局のところそこまで変わらないので、処理時間的にも短い通常の最急降下法を最終的に採用したわけだが。気持ちとしては後者の方がテクくて好きだったので何とかして採用したかった。。。

大半の人はスキル推定に焼きなましやMCMCを使っているので、実際のところそちらの方がいいのかもしれない。

あと、これを書いてて思ったが、振られたタスクの要求スキルが小さいと各スキルは過小評価されるので、最初のうちは満遍なくいろいろなタスクを振ってスキル推定の精度を向上させた方がよかったかもしれない。

ちなみに

作業スケジューリングについては、重み付き二分マッチング問題を誰かがタスクを終える度に実行して作業の割り当てを行った。この際、辺の重みは予測される作業時間に加えて、なるべくそのタスクを終えることで次に実行できるタスクが増えるように調整している。スキル推定はそこまで重要じゃない気もするので、作業スケジューリングの方にもっと時間を割いて考えるべきだったかもしれない。

所謂ヒューリスティックコンテスト形式のコンテストはAtCoder Heuristic Contestが始まってからまともに参加するようになったが、今回のように連続最適化の手法も使えるコンテストがたまにあって中々楽しい。

*1:最終結果はスコア5,664,033で129位。

*2:スキルとして取りうる値は整数だが、ここでは実数値として考え、後から四捨五入する方法を取った。

*3:一応、最初は絶対値ではなく2乗で0以上を保証するなど微分可能な式を立ててやってたが結果的には今の定式化の方がよかった。今や深層学習なんかでReLU等をガンガン微分してるのでまぁ。。。

*4:以前書いた内点法とかyamagensakam.hatenablog.com