AHC025 参加記 - eijirouの競プロ参加記

AHC025で優勝しました。長期AHCの過去最高順位は10位だったということもあり、非常に嬉しく思います。

はじめに

解法だけでなく、考察や余談についても記載しました。
重要度が低いと思われるものは折りたたみにしました。必要に応じてクリック（タップ）してください。
ベイズ推定やMCMCについては私自身が詳しくないため、コンテスト後に少し勉強して理解したことを書きました。「お気持ち」だけ汲み取っていただければ幸いです。
アイテム集合 $A$ に対して、それらの重さ（の推定値）の和を $sum(A)$ と表記しました。行儀が悪い表現ですがお許しください。
環境によっては数式が読みにくいと思います。下の記事などを参考にすると読みやすくなるかもしれません。

syleir.hatenablog.com

問題

公式の問題文と入力生成方法を参考にしてください。

atcoder.jp

この記事では $N, D, Q, \lambda$ などを問題文と同様に定義して使用します。

処理の流れ

$(Q - N \log_2{N}) / N$ が大きいとき、クエリを消費しながらマージソートを行い、アイテムを重さ順にソートする。
$D$ 個の集合を初期化する。
1~500個ほどの推定器（各アイテムの重さを推定するもの）を初期化する。
山登り法で改善する。具体的には以下の処理を繰り返す。
1. 各推定器について、今までのクエリとの矛盾がなくなるように、反復法によって改善を行う。
2. 時間に余裕がある場合、推定器を大きく更新し、再び反復法を適用する。
3. 推定器の平均値を求める。
4. 推定器の平均値を元に、2つの集合の間でいくつかのアイテムを入れ替える遷移でよさそうなものを列挙する。
5. 推定器をサンプルとみなして各遷移のスコア改善量の期待値を求め、期待値が最も高いものを選択する。
6. 選ばれた遷移によってスコアが改善するか否かをクエリで判定し、改善する場合にはその遷移を採用する。
7. 残りの実行時間が十分短いとき、適当なクエリを投げてクエリを無駄に消費する。（TLE緊急回避モードです。基本的には実行されません。）
8. クエリを合計で $Q$ 回行った場合、山登り法を終了し、解を出力する。

各処理について説明します。

アイテムのソート

$(Q - N \log_2{N}) / N$ がおよそ3.0以上のとき、マージソートでアイテムをソートしました。

閾値の設定について

マージソートで $N \log_2{N}$ 回弱のクエリを消費します。

ソートした後、山登り法のために $N$ の定数倍ほどのクエリを残したいという理由から、このような閾値を設定しました。

実際には $D$ が大きいときに閾値を上げていますが、改善の度合いはごくわずかです。

ソートを行う理由

後で重さを推定するときに、各アイテム間の順序関係がわかっていると非常に推定の精度が上がります。

例えば 0-indexed で $i$ 番目に軽いアイテムの重さを $- \frac{\log{(1 - (i + 0.5) / N)}}{\lambda}$ で推定するということが考えられます。

式の導出

指数分布について検索すると、累積分布関数が

$F(x; \lambda) = 1 - e^{-\lambda x}$

であることがわかります。 $F(x; \lambda)$ の逆関数を求めます。

$y := F(x; \lambda)$

$e^{-\lambda x} = 1 - y$

$-\lambda x = \log{(1 - y)}$

$x = - \frac{\log{(1 - y)}}{\lambda}$

[0.0, 1.0] を $N$ 等分し、中央の値で代表すると $- \frac{\log{(1 - (i + 0.5) / N)}}{\lambda}$ が得られます。

$x = - \frac{\log{(1 - y)}}{\lambda}$ は一様分布から指数分布を生成することにも使用することができます。

seed = 1 の場合、次のような推定結果になりました。

グラフにおいて、アイテムの真の重さの昇順にソートされ、重さの和が1になるように正規化されています。以降の推定に関するグラフにも同様のことが行われています。

比較回数の最小化

なるべく比較回数が少ないソートを選ぶことにより、山登り法でより多くのクエリを使用することができます。

私は比較回数が少なそうなマージソートを使用しました。

ソートの種類について

ソートの性能を評価する指標として次のようなものがあります。

今回のコンテストで重要な指標
- 比較回数
- ソートを途中で止めたときの情報
一般的な指標
- 実行速度
- 安定性

代表的なソートアルゴリズム

クイックソート

一般的に速いとされていますが、安定性はありません。（安定性をもたせる実装だと遅くなると思います。）

理想的な場合の比較回数はあまり多くないですが、実際には理想的に動くとは限らないため、比較回数はやや多くなります。また、ピボットの選択が乱択でない場合、計算量が $O(N^{2})$ になるケースを意図的に生成できます。

ソートを途中で止めた場合、いくつかのグループに分かれ、グループ間の大小関係がわかる状態になります。

軽いアイテムのソートを優先したり、逆に重いアイテムのソートを優先したりすることが可能です。

ヒープソート

小さいもの（あるいは大きいもの）から得ることができます。

ソートを途中で止めた場合、軽いアイテム（あるいは重いアイテム）の順序関係を正しく把握できます。

マージソート

安定性があり、比較回数は少ないです。

トップダウン方式とボトムアップ方式があります。

トップダウン方式は全体を半分ずつに区切る方式で、再帰関数を用いることで実装しやすくなります。

ボトムアップ方式は隣り合うグループをマージしていくような方式（Segment Tree をイメージするとわかりやすいかもしれません）で、通常は非再帰で実装します。

ボトムアップ方式はダブルバッファリングとの相性がよく、トップダウン方式よりも若干速い場合があります。

一方でトップダウン方式のほうが比較回数は少ないです。ボトムアップ方式だと $N = 33$ のときに、最後に 32 のグループと 1 のグループに分かれ、単純にマージすると遅い可能性が高いです。

私が今回のコンテストで使用したのはトップダウン方式です。

途中で止めると、いくつかのグループに分かれ、グループ内でソートされた状態になります。クイックソートと対照的ですね。

merge-insertion sort

比較回数の最悪ケースが少ないらしいです。よく知りません。

ソートを途中で止める場合に、軽いアイテムと重いアイテムのどちらの情報が重要か

よくわからないというのが私の結論です。

直感的には重いアイテムのほうが影響が大きいため、重いアイテムの情報が重要そうに思えます。

一方で、山登り法で頻繁に動かすのは軽いアイテムであり、軽いアイテムの大小関係だけ把握するという考え方もあります。

ソートするかしないかの二択が離散的すぎて気持ち悪いという感覚があり、途中までソートすることを試したのですが、結局のところうまく使えませんでした。

ソート方法のまとめ

今回はマージソートか merge-insertion sort を使うのがいいと思います。

集合の初期化

ソートを行ったかどうかにかかわらず、アイテム $i$ をグループ $i \mod D$ に割り当てました

マージソートを行った場合、各アイテムの重さに関してある程度正しい推定ができ、初期解のコストを下げることが可能です。しかし、初期解をよくすると最終的なコストが増加しやすかったため、初期解は常に $i \mod D$ で設定しました。

~~局所探索の初期解はランダム性が高いほうがよいのでしょうか。私には理由がわかりませんでした。~~

（2023/10/30 追記: 初期解がよいとアイテムの入れ替え回数が減り、同一集合内のアイテム同士の比較が行われず、推定の精度が下がるからだと思います。）

推定

反復法

まず、今までのクエリ結果と矛盾しない（このことを制約と呼ぶことにします）重さを求めるための反復法について説明します。

推定の初期値は指数分布で生成し、ソートした場合には推定値もソートします。

以下に反復の流れを示します。

次の処理を最大30回行う。
1. 今までの全てのクエリについて次の処理を行う。
  1. 集合 $l, r$ の和をそれぞれ求める。
  2. 求めた和とクエリの結果が矛盾している場合、重さの和が入れ替わるように $l, r$ の各要素の重さをスケーリングする（修正量をアイテムの重さに比例させる）。
  3. 重さが1.0未満のものは1.0にし、重さが 10⁵ N / D 以上のものは 10⁵ N / D にする。
2. 更新が行われなかった場合、反復法を終了する。

クエリの結果が等号だった場合、重さの和が等しくなるようにスケーリングしました。稀にしか起こらないと思います。

重さのスワップを行った経緯

元々は最急降下法によって重さを推定しようと考えていました。

不等式が満たされない場合に、重さの和が等しくなるように推定値を更新するときの修正量を「学習率が1.0である」と定義することにします。すなわち、学習率を $\alpha$ とすると、 $sum(l)$ が $sum(l) + \alpha \frac{sum(r) - sum(l)}{2}$ になるように更新します。

学習率の値をいくらぐらいに設定するのがよいか考えます。

学習率が1.0未満の場合、満たされない制約が1つである場合などに更新を無限に繰り返すことになるため、不適切であると考えられます。（機械学習の文脈では基本的に学習率を1.0未満にしますが、今回は学習率の定義が異なることやデータが正確であることに注意してください。）

学習率が2.0より大きい場合、重さが負になる場合が考えられ、発散する場合もあるかもしれません。不適切です。

学習率は1.0以上2.0以下がよさそうです。パラメータ調整を行ったところ、学習率は2.0になりました。学習率2.0というのは重さの和を入れ替えるということになります。学習率を2.0にした場合、学習率を1.9にした場合よりもかなり性能がよく、学習率は敏感なパラメータだったと言えます。

連立一次方程式を解く際に、修正量を拡大して収束を加速させるという手法として、SOR法というものがあるらしいです。今回は制約が不等式のためSOR法ではありませんが、修正量を大きめに設定するという考え方は似ているかもしれません。ちなみにSOR法（かそれに近い何か）はトヨタコン決勝の懇親会でsaharanさんがAHC022（リクルートハーフマラソン）の文脈で教えてくださいました。ありがとうございます。

重さの和をスワップすることの利点

$sum(l) + sum(r)$ が不変のため、アイテムの重さの和が保たれます。重さの和が保たれない場合、クエリに含まれないアイテムの値を相対的に変更することになり、あまりよくない感じがします。

クエリの比較対象が要素1個同士の場合、単純に2つのアイテムの重さを入れ替えるだけなので、初期値が指数分布であれば指数分布が保たれることになります。分布が保たれるというのは非常によい性質です。実際には比較対象が1個同士とは限らないため、分布は少しずつ崩れると思います。

ベイズ推定の考え方

ソートを行った場合に、全ての制約を満たす重さとして次の2つの解を見つけたとします。

事前分布が指数分布であることを踏まえると、右側の推定のほうが確率が高そうです。

制約を満たす解の中でも、指数分布に近いものを見つけたいというモチベーションが生まれます。

理想的には解析的に事後分布を求めたいですが、おそらくそのような計算は困難です。そのためサンプリングを行って分布を近似することを考えます。

MCMCの考え方

サンプリングを行う度に指数分布から値を生成し（ソートした場合はソートも行い）、全ての制約を満たすもののみをサンプルとして扱うという方針が考えられますが、制約が増えるにつれて受理率が非常に低くなることが予想されます。

そこで、既に得られたサンプル（特に直前のサンプル）から新たなサンプルを生成して使用することにします。既存のサンプルから新たなサンプルを生成する方法として多くのバリエーションが考えられています。

サンプリング

並列に $\max{\left( 1, \frac{2 \times 10^{7} \times D}{N Q^{2}} \right)}$ 個ほどのサンプル（推定器）を用意します。

推定器の個数について

実行時間が2秒に収まる範囲で、なるべく推定器の個数を増やしたいです。

プログラムのボトルネックは、反復法で集合の和を求め、推定値を更新するところです。

集合の平均要素数はおよそ $N / D$ です。

クエリの数は $Q$ 回です。

制約の数は平均して $Q / 2$ 回程度です。

以上より、推定器の個数を $\frac{N Q^{2}}{D}$ に反比例させるとよさそうだということがわかります。

実行時間に余裕があるとき、直前のサンプルから新たなサンプルを生成します。

具体的には次のような手順です。

指数分布から $N$ 個の値を生成する。
アイテムの大小関係が元のサンプルと等しくなるように生成した値を並べ替える。
反復法を適用し、制約を満たす解へ収束させることを目指す。

ソートを行っていた場合、アイテムの大小関係はサンプルによらないため、元のサンプルを無視して新たにサンプルを生成することと等価になります。

なお、実装上は反復法が収束したか否かはほとんど区別しませんでした。処理を収束したものに限定するよりもサンプル数を減らさないことのほうが重要だったようです。

ちなみに私よりも高得点だったwriterのwataさんは、ギブスサンプリングでサンプルを更新したそうです。ここで差をつけられた可能性もありますね。

サンプリング結果

ソートしない場合に、いくつかのアイテムの推定を大きく間違えていました。

山登り法

初期解を山登り法で改善します。

2つの集合の間でいくつかのアイテムを入れ替え、コストが減少したなら採用し、そうでないなら元に戻します。

コストが減少したか否かを判定する方法について述べます。

元の集合を $A \cup B, C \cup D$ とし、新たな集合を $A \cup D, C \cup B$ とします。 $A, B, C, D$ の間に共通するアイテムはないものとします。

$sum(A \cup B) \lt sum(C \cup D)$ の場合、 $sum(A) \lt sum(C)$ かつ $sum(B) \lt sum(D)$ ならばコストが減少します。 $sum(A \cup B) \gt sum(C \cup D)$ の場合は不等号が逆になります。

最初の条件は必要か

コストが減少する必要十分条件は

$sum(A) \lt sum(C)$ かつ $sum(B) \lt sum(D)$

または

$sum(A) \gt sum(C)$ かつ $sum(B) \gt sum(D)$

です。2回のクエリで判定することができ、 $sum(A \cup B) \lt sum(C \cup D)$ で場合分けする方法よりも効率的に思えるかもしれません。

しかし、実際には場合分けを行ったほうがよいと思います。元のグループの大小関係はわかっている場合が多く、そのような場合は $sum(A \cup B)$ と $sum(C \cup D)$ の比較を省略できます。 $sum(A \cup B) \lt sum(C \cup D)$ で $sum(A) \gt sum(C)$ の場合に、 $sum(B)$ と $sum(D)$ の比較を省略でき、比較が1回で終わることがあります。

最大3回の比較を要しますが、集合間の大小関係を把握できるメリットもあるため、集合同士の比較も行ったほうがよいと考えられます。

不採用のときに比較を1回で済ませるために、 $sum(A) \lt sum(C)$ 、 $sum(B) \lt sum(D)$ のうち、成り立ちにくそうなほうから比較を行いました。

クエリの省略

結果がわかる場合はクエリを省略します。省略するパターンを以下に示します。

片方の集合が空のとき
既に同じクエリを投げているとき
集合の比較で、推移律から結果がわかるとき

集合の大小関係が $sum(A) \lt sum(L) \lt sum(B) \lt sum(R) \lt sum(C)$ となっていて、集合 $L, R$ を変更してコストが下がったとします。このとき、 $sum(L), sum(R)$ に関して以下の関係が保存されます。

$sum(A) \lt sum(L)$
$sum(A) \lt sum(R)$
$sum(L) \lt sum(C)$
$sum(R) \lt sum(C)$

コンテスト上位の方の中には、より複雑な比較の省略を行った人もいたようでした。

入れ替え候補の列挙

2つのグループを選択し、重さの推定値が与えられるものとすると、partition problem に帰着されます。partition problem はNP完全で、擬多項式時間アルゴリズムやヒューリスティックな手法が存在します。

私は半分全列挙をベースとして partition problem を解きました。理想的には厳密解が得られますが、集合の要素数が多いときは長い時間がかかります。そこで、各集合の入れ替え可能なアイテムを、重さが小さいものを中心に 2 ~ 10 個に限定しました。

「重さが小さいものを中心に」の補足

より正確には、元の集合における和の差の推定値の 1/4 に近い重さのアイテムを優先して選びました。集合の差はすぐに小さくなるため、重さが小さいものを優先しているのとあまり変わりません。

優先度には乱数をかけました。

具体的な優先度の付け方を以下に示します。

アイテムの重さの推定値を $w$ とする。
優先する重さを $\Delta := \frac{|sum(l) - sum(r)|}{4}$ とする。
乱数の大きさを調節する変数（温度）を $t:= \frac{N Q'}{23000 \times D}$ とする。ただし $Q'$ は残りのクエリ数を表すものとする。
優先度を $rand\_double(0, t) + | \log_2{\frac{w}{\Delta}} |$ で定義し、優先度が小さいものを優先して選ぶ。

入れ替え可能なアイテム数について

$Q$ が大きいほど入れ替え可能なアイテム数を増やしました。

$Q$ が小さいときは推定の精度が悪かったからか、変更を小さくしたほうがスコアが上がりました。

半分全列挙の実装上の工夫

全列挙パートでは要素を1つずつ追加します。その際、マージ手法を使うことによりソート関数を呼ばずにソートできます。また、復元用に全ての配列は削除せずに保持します。

2つの全列挙された配列はソートされているため、尺取法の要領で最適な値を計算できます。要するに $O(2^{n / 2} n)$ ではなく、 $O(2^{n / 2})$ で計算できます。

解の復元は、全列挙パートで追加したときと逆の順序で行います。和が前の配列にあればその要素を使う必要がなく、そうでないときは使う必要があります。配列がソート済みであることに注意すると、値が含まれるか否かは二分探索でわかるため、復元は $O(n^{2})$ で行えます。

数値誤差に対処するのは面倒なため、重さは全て整数として扱いました。

関連する AtCoder Beginner Contest の問題 atcoder.jp

私の提出コード atcoder.jp

3-partition problem

2-partition problem へ帰着すると解きやすい問題になりますが、本来解くべき問題は D-partition problem です。

2-partition problem しか解かない場合、山登り法の局所解に陥るケースがありそうだったため、3-partition などのより多様な近傍を使用することをコンテスト中に考えました。しかし、私は 3-partition の近傍でスコアを上げることができませんでした。

うまくいかなかった理由として以下のことが考えられます。

そもそも 2-partition の近傍で十分よい解が得られる。
半分全列挙をベースとする方法で 3-partition problem を解くことに無理があった。Largest Differencing Method などの別の手法を使うべきだった。（Largest Differencing Method はコンテスト後に知りました。）

入れ替え候補の選択

以上のことを踏まえて、次のような流れで入れ替え候補の選択を行います。

サンプルの平均値を計算する。
2つの集合の組を全列挙し、次の処理を行う。ただし、時間がないときは集合の和の差が大きいペアについてのみ処理を行う。
1. サンプルの平均値を推定値とし、半分全列挙ベースの手法でよさそうな入れ替えを求める。
2. 入れ替えることによる改善の期待値をサンプルから（擬似的に）求める。
期待値が最も高い入れ替えを求める。

全てのサンプルが制約を満たすとき、サンプルの平均値も制約を満たします。これは制約に線形性があることから言えます。

線形性の補足

具体例で説明します。

$x + y \lt z$

という制約があり、

$(x, y, z) = (x_1, y_1, z_1), (x_2, y_2, z_2)$

が制約を満たすとします。

このとき、

$x_1 + y_1 \lt z_1$

の両辺を2で割ると

$\frac{x_1}{2} + \frac{y_1}{2} \lt \frac{z_1}{2}$

が得られます。よって

$(x, y, z) = (\frac{x_1}{2}, \frac{y_1}{2}, \frac{z_1}{2})$

も制約を満たします。より一般的に、正の定数をかけても制約を満たすということが言えます。

さらに

$x_1 + y_1 \lt z_1$ と $x_2 + y_2 \lt z_2$ の両辺を足すと

$(x_1 + x_2) + (y_1 + y_2) \lt z_1 + z_2$

が得られます。よって

$(x, y, z) = (x_1 + x_2, y_1 + y_2, z_1 + z_2)$

も制約を満たします。

解に正の定数をかけたり解同士の加算を行ったりしたものも解になることから、平均や加重平均をとってもよいということがわかります。

また、この記事においてグラフを描画するときに重さの和が1になるように正規化している理由は線形性があるからです。

期待値の求め方

基本的には、全てのサンプルについて $\max{(0, 改善量)}$ を計算し、その平均値をとればよいです。

私の実装では、サンプルの平均もサンプルとして扱っています。サンプル数を $n$ として、サンプルの平均は $\sqrt{n} / 2$ 個分のサンプルとして期待値を計算しました。加重平均の重みが普通のサンプルの $\sqrt{n} / 2$ 倍になっているということです。

サンプルの平均もサンプルとして扱うことで、サンプル数が少ないときに正則化効果が生まれたと推測しています。

その他

集合の表現方法

3つの方法を紹介します。

配列

アイテムの番号を整数として配列に保持します。最も単純な表現方法だと思います。

C++だと std::vector<int> などです。

ソートしない場合、要素の追加は $O(1)$ ですが、要素の削除や含まれるかの判定が $O(n)$ で遅いです。

Index Set

tomerunさんのブログを参考にしてください。

topcoder-tomerun.hatenablog.jp

今回の問題では $N \leq 100$ であるため、平衡二分木などよりも高速で便利だと思います。

bitset

C++における std::bitset<100> です。

次のような利点があります。

コピーが高速である。
集合の演算をビット演算で簡潔に記述でき、処理が高速である。

使い分け

配列と bitset を使用しました。イテレーションを回すことが多い場所では配列を使用し、そうでないところでは bitset を使用しました。

bitset を配列に変換する方法

C++20 の std::countr_zero などを使用することで、要素数の線形時間で変換できます。定数倍が軽くはなく、 $N$ 回ループを回したときと速度があまり変わらないかもしれません。

vector<int> bitset_to_vector(const bitset<100>& l) {
    ull low = (l & bitset<100>(ULLONG_MAX)).to_ullong();
    ull high = (l >> 64).to_ullong();

    vector<int> ret;
    ret.reserve(l.count());

    while (low) {
        int v = countr_zero(low);
        ret.push_back(v);
        low ^= (1ull << v);
    }
    while (high) {
        int v = countr_zero(high);
        ret.push_back(v + 64);
        high ^= (1ull << v);
    }
    return ret;
}

順位表の使い方

ACを一部のテストケースに制限し、自分が苦手なタイプのケースを特定していました。

assert(Q < 8 * N);

コンテスト終盤は $Q / N$ が小さいケースに弱いことがわかっていたため、 $Q / N$ が小さいケースを中心に改善点を探りました。

結局のところ、最後まで $Q / N$ が小さいケースに弱かったため、サンプリングをうまく使えていなかったのかもしれません。

スコアの評価方法

各テストケースのコストの対数を求め、その和を最小化しました。

パラメータ調整

大雑把なパラメータ調整は実装したときに行い、細かいパラメータ調整はコンテスト最終日に行いました。

$N, D, Q$ によって決定されるパラメータ調整は、関数の探索になるのでかなり面倒でした。

パラメータ調整は全て手動で行いました。Optuna のような自動最適化ツールは使っていません。

高速化の寄与

高速化によってどれだけコストが下がるかは、実行時間を長くすることにより、高速化を行う前に推測できます。

私の方針だと高速化は非常に重要で、参照渡しにし忘れていた部分を直しただけで0.8%ほど順位表のスコアが上がったことがありました。

提出コード