TERRYのブログ

最近は競技プログラミングにお熱。

Sponsored Link

トップ > AtCoder > HACK TO THE FUTURE 2025 (AHC040) 解説

HACK TO THE FUTURE 2025 (AHC040) 解説

AtCoder 競プロ マラソン

HTTF2025（AtCoder Heuristic Contest 040）お疲れ様でした。2,955,084,008,755点で6位です。前回に引き続き解説記事を書きます。

焼きなまし・ビームサーチ・モンテカルロ木探索（MCTS）・マルコフ連鎖モンテカルロ法（MCMC）・多変量正規分布・ベイズ推定・SIMD高速化・ガウス過程回帰といろんな要素が盛りだくさんです。それなりに上級者向け解説なので、分からないところは適宜調べて頂ければと思います。*1

ビジュアライザ (seed=0, score=753,883点)

問題概要

できるだけ小さい領域に長方形を詰め込んでね。

長方形の大きさは高橋くんが目分量で測るからバラツキがあるよ。

atcoder.jp

解法概要

推定パートと詰め込みパートに分けて解きました。

ラスト10ターン弱くらいまでは推定パートとして、箱を縦横にたくさん並べて長さの計測に専念しました。多変量正規分布のベイズ推定を行い、その分散・共分散ができるだけ小さくできるような配置を焼きなましで求めます。

その後、MCMCで N 個の長方形のセットを16セットサンプリングして、その16セットに対して並列にビームサーチとMCTSを行うことで詰め込み手順を求めました。16セット同時に詰め込み・評価を行うことで、「一部のケースだけ上手く行く」という詰め込み方が選ばれないようにし、ロバスト性を向上させる効果を狙っています。16セットにしたのはSIMDで処理するためで、これにより16セットを並列でシミュレーションすることができます。

詰め込みパートの概要

動画にするとこんな感じです。縦横に並べているのが推定パート、正方形に近い形に詰め込んでいるのが詰め込みパートです。

動画版ビジュアライザ（seed=0）

提出コードはこちらです。ライブラリ部分を除いても4500行以上ある超大作となっているので読むのはオススメしません……ｗ

atcoder.jp

推定パート解法

ほぼこれ↓です。

qiita.com

箱の縦横の長さを並べたベクトル \boldsymbol{X}=\lbrack w_0, h_0, w_1, h_1, \cdots, w_{N-1}, h_{N-1}\rbrack^T が多変量正規分布 \mathcal{N}(\boldsymbol{x_0}, \boldsymbol{\Sigma_0}) に従うので、この期待値 \boldsymbol{x_0} と分散共分散行列 \boldsymbol{\Sigma_0} をベイズ推定します。

事前分布

まず入力を受け取る前に、 w_i, h_i の生成方法が与えられているので事前情報として入れておきましょう。 w_i, h_i は w_i, h_i=\mathrm{rand}(L,U) により生成されるのですが、これを大胆にも正規分布 \mathcal{N}(\mu_0, \sigma_0^2) で近似してしまいます。

下記リンク先のようなPythonコードで実験をすると、概ね \mu_0=65000, \sigma_0=21280 くらいの値が得られるので、 \boldsymbol{x_0} を全要素が 65000 の 2N 次元ベクトル、 \boldsymbol{\Sigma_0} を対角成分が全て 21280^2 の 2N\times 2N 対角行列として初期化します。*2 *3

github.com

初期入力による更新

初期入力で各長方形の大きさの計測値 w_i', h_i' が与えられるので、これを用いて 2N 回ベイズ更新を行います。更新式はeivourさんのAHC003解説に書かれている通りで、記事内の y_t=w_i' 、\boldsymbol{c_t}=\lbrack 0, 0, \cdots, 1, \cdots, 0\rbrack^T として更新を行えばよいです。

配置の最適化と更新

標準偏差が \sigma の正規分布を k 個足し合わせるとその分散は k\sigma^2 となり、標準偏差は \sqrt{k}\ \sigma となりますが、高橋くんはどんなに箱を増やしても標準偏差 \sigma で計測してくれるので、たくさん並べた方がお得です。1回のクエリで縦と横の長さを教えてくれるので、縦横双方に長くなるよう下図のようにL字に並べることにしましょう。

L字に配置する例

さて、この配置方法ではランダムに箱を並べてもそこそこよい推定ができるのですが、せっかくなら最も効率が良い配置でクエリを投げたいです。今考えている多変量正規分布のバラツキの大きさが小さくなると嬉しいのですが、この分布のバラツキの大きさは分散共分散行列 \boldsymbol{\Sigma_t} の固有値の大きさと捉えることができます。ここで、行列の対角成分の和をトレースといい、トレースは行列の固有値の和と等しくなるため、分散共分散行列のトレースが小さくなるようなクエリを投げると良さそうです。*4

トレースの減少量の式は既にAHC003の解説記事に書かれているので、これを評価関数として配置方法を焼きなまします。この式は愚直に計算すると \Theta(N^2) かかってしまいますが、適切に実装すると O(N) で差分更新が可能です。

ここで注意点として、毎回箱0を左上に置いてしまうと「箱0ありの計測値」しか取れず、「箱0なしの計測値」が分からないので、箱0の長さの推定精度が悪化します。それを避けるため、下図のようにT字型の配置も許容して焼きなましを行いました。*5

T字に配置する例

なお、L字部分の箱同士の干渉には注意が必要です。シミュレーションすればだいたいのケースで大丈夫なのですが、時々下図のように想定と異なる干渉の仕方をしているケースがあります。下図の例では事前情報（緑）から h_0+h_2\gt h_1 と推定され、 \boldsymbol{c_t}=\lbrack 0,1,0,0,0,1,\cdots\rbrack^T とするのが正しいように見えますが、実際の干渉状況を見ると h_0+h_2\lt h_1 であり、 \boldsymbol{c_t}=\lbrack 0,0,0,1,0,0,\cdots\rbrack^T とすべきで、前者を使うと誤った更新となってしまいます。この部分は後述パートのMCMCで吸収します。

L字部分の干渉が想定と異なる例（事前情報からは箱0+箱2>箱1と推定されていたが、実際は箱0+箱2<箱1のため誤った推定となる）

以上を実装すると、それなりに精度よく長方形の大きさを推定することができます。seed=0では入力で与えられる標準偏差が \sigma=5574 、焼きなましなしの完全ランダムだと \sigma=1920\ (\pm 240) 程度（括弧内は推定した標準偏差 2N 個の標準偏差）ですが、焼きなましを入れると \sigma=1400\ (\pm40) 程度まで下がるので、この後の配置パートがかなりやりやすくなります。

詰め込みパート解法

ターンごとに以下の1.～3.を繰り返し、その結果を出力することを繰り返しました。

  1. MCMCで N 個の箱を16セットサンプリング
  2. ビームサーチで中盤（N-15 個くらいまで）までの詰め込み方法を求める
  3. ビームサーチで求めた詰め込みを固定して、最後の仕上げ部分をMCTSで求める

詰め込みパートの概要（再掲）

MCMC

MCMC（マルコフ連鎖モンテカルロ法）を知らない人もいると思いますので最初に簡単に説明すると、ある確率分布から効率良くランダムにサンプリングするための手法です。実装内容は「対数尤度をスコアとした温度 1 固定の焼きなまし」と思ってもらえればだいたい大丈夫です。*6 詳しく知りたい方は以下の書籍がオススメです。

ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門 (KS理工学専門書)

ゼロからできるMCMC マルコフ連鎖モンテカルロ法の実践的入門 (KS理工学専門書)

Amazon

推定パートで多変量正規分布モデルを使って推定を行ってきたので、多変量正規分布からボックス=ミュラー法とコレスキー分解を使って直接サンプリングしてもよいです。ただ、多変量正規分布という仮定では詰め込みパートに入った後の計測結果を上手く使うことができず、特に T の小さいケースでかなりの情報を捨ててしまうことが懸念されるため、MCMCを用いて改めてサンプリングを行いました。

MCMC自体は何の変哲もないメトロポリス法です。箱 i の大きさを変化させたとき、t ターン目の詰め込みをシミュレーションすることで詰め込み結果の W_t, H_t が分かり、これと t ターン目の計測結果 W_t', H_t' を比較することで尤度が分かるので、これを使ってメトロポリステストを行います。

1ターン目のMCMC法の初期サンプルとしては多変量正規分布モデルからボックス=ミュラー法とコレスキー分解を使ってサンプリングした値を、2ターン目以降の初期サンプルは前ターンの最終的なサンプリング結果を用いました。なお、MCMC法では箱が16セット得られればよいので必要以上にサンプリングする必要はなく、自己相関がなくなる程度にサンプリングを進めれば十分です。

①多変量正規分布からサンプリングを行った手法 ②MCMCを行うが詰め込みパートの情報は捨てる手法 ③MCMCを行い、詰め込みパートの情報を使って都度尤度の計算式をアップデートする手法 を比較すると、①と③は1000ケースで0.1%ほど、②と③は1000ケースで0.05%ほどスコアに差が付きました。①と③は前述のL字部分の干渉が想定と異なるケースも考慮できた点、②と③は詰め込みパートの情報も有効活用できた点で差が付いたものと思われます。*7

ビームサーチ

詰め込みの中盤まではビームサーチを行いました。よくあるBottom-Left法ベースの解法です。*8

Bottom-Left法の説明については以下文献などをご参照ください。

Bottom-Left法ベースの詰め込み

詰め込み方針

下から・右からの双方から自由に詰め込めるスタイルだと上手く行かなかった*9 ので、予め領域の幅を \alpha\sqrt{\sum_iw_ih_i} （\alpha は 1.05 ～ 1.15 くらいでチューニング）だけ確保し、下からの投入だけを考えることにしました。

ビームサーチの状態遷移先については、有望そうな候補として「箱の左側面が壁か他の箱にピッタリくっついているもの」に絞って遷移しました。下図の例だと、操作 (U, 0) や (U, 2) はいかにも中途半端な位置に配置されていて、ここから良い解が得られる見込みは低そうです。そこでこのような遷移はNGとしました。

遷移先の候補

各遷移先について、\Theta(N) かけて長方形の配置される y 座標を求めて、そのときの評価値を計算しました。

評価関数

評価関数は「詰め込み領域の高さ」の重み付き和としました。*10 箱 i まで詰め込んだときの領域の高さを H として、これを16セット分昇順にソートしたものを H_0, H_1, \cdots, H_{15} とします。その上で適当な定数 \beta\ (0.5\le\beta\le1) を用いて、 H_0+\beta H_1+\beta^2H_2+\cdots+\beta^{15}H_{15} を評価値としています。平均ではなく重み付き和としているのは、解を複数回提出可能なことから平均値を改善するよりも上振れ狙いをした方が良い結果が得られると考えたためです。

ここまでの計算量は、

  • ビーム幅を B 、長方形のセット数を P\ (=16) とおく
  • ビームのターン数は長方形の数と同じ N
  • 各ターンの遷移先候補は N 個程度
  • 各遷移先について、
    • 配置先の y 座標を P セット求めるのに O(PN)
    • 領域の高さ H のソートに O(P\log P)

であることから、O(BPN^2(N+\log P)) となります。

候補数の削減

ここまでは全ての遷移先候補の操作について処理を行っていましたが、上記の操作 (U, 0) や (U, 2) のように「長方形 i に沿って投入したのに長方形 i より手前で止まってしまう」という操作は無駄です。そのような操作が今後有効な操作になることはないため、有効な操作の集合保持しておきそのような操作についてのみ遷移を試みることにします。

下部に露出している長方形に沿って投入する操作のみが有効な操作となりうることを考えると、そのような操作は平均で \sqrt{N} 個程度となることが期待されます。このとき、

  • ビーム幅を B 、長方形のセット数を P\ (=16) とおく
  • ビームのターン数は長方形の数と同じ N
  • 各ターンの遷移先候補は \sqrt{N} 個程度
  • 各遷移先について、
    • 有効かどうかの判定に O(PN)
    • 配置先の y 座標を P セット求めるのに O(PN)
    • 領域の高さ H のソートに O(P\log P)

であることから、平均計算量は O(BPN\sqrt{N}(N+\log P)) に落ちます。今回の問題では N が最大で 100 程度となることを考えると、 \sqrt N 倍高速化は最大でざっくり 10 倍程度の高速化となり、結構嬉しいです。

MCTS

ビームサーチの評価関数で高さしか考慮していないため、特に終盤の詰め込み方がイマイチになりがちでした。そのため、終盤をMCTSで探索することでごまかします。考え方としては「ビームサーチの評価関数をランダムプレイアウトにする」のと感覚が近いかもしれません。*11

ビームサーチは最後の仕上げ部分が弱く、だからといって最初からMCTSをすると候補数が多すぎて最終手までノードの展開が行えないため、中盤までビームサーチ・終盤をMCTSとすることで特性を補い合うことを狙っています。

長方形詰め込み問題におけるMCTSのアイデアは下記のbowwowforeachさんの記事を参考にしました。

bowwowforeach.hatenablog.com

候補手のノード

だいたい3回くらい訪れたノードは展開します。

各ターンで取れる行動（配置位置×回転有無）のうち、置いた箱の y 座標が小さいものから順に5個程度を候補としてノードに持たせ、 UCB1-Tuned に従ってノードの選択を行いました。

プレイアウト

プレイアウトでは各ターンで回転をランダムに決め打ち、置いた箱の y 座標が最も小さい位置に置く貪欲をしました。

プレイアウト中は実行した操作手順を保持しておき、ベストスコアを更新したら暫定解を更新します。これを時間いっぱい繰り返し、最終的に得られた解を出力しました。

SIMD高速化

この解法の非自明おもしろパートです。

まずSIMDとはSingle Instruction, Multiple Dataの略で、CPUの1命令で複数のデータをまとめて処理することができます。SIMDについて詳しくはリンク先記事をご参照ください。

AVX2では、CPU内の256bitレジスタを使って複数のデータに対してまとめて演算ができます。データを16bit整数で持つと、下図のような16bit整数16個の足し算が普通の足し算1回と遜色ない速度で完了します。すごいですね。

SIMDによる16要素同時足し算

SIMDは画像データなどの大量のデータを高速に処理するために使われることが多いかと思いますが、このような少量かつ固定長のデータに対して使うケースはあまり見ない気がするので、それなりにおもしろポイントかなと思います。

データの持ち方

AVX2では256bitレジスタを用いて計算するのですが、計算対象のデータの型によって同時に演算できるデータの数が変わってきます。16bit整数なら 256\div16=16 要素を同時に処理することができますが、これが32bit整数だと 256\div32=8 要素しか同時に処理することができません。そのため、できるだけ16bit整数で処理を行いたいです。

長方形の大きさはそのままだと符号なし16bit整数に収まらないため、長方形の大きさを一律で 128 で割って符号なし16bit整数に収めました。その分精度は犠牲になりますが、そもそも扱う長方形の大きさは不確かな推定値であるため、少しズレたくらいで特に問題はありません。

また、各データは AoS (Array of Structures) ではなく SoA (Strucute of Arrays) の形で持つようにします。

ビットマスクを使って並列化

SIMD処理では、if文の代わりにビットマスクを使って条件分岐を行うということをよくやります。

詳細は下図をご覧頂ければと思いますが、このように処理することでif文を使わず記述することができ、高速に16並列で処理することができます。

長方形16セットに対して操作をシミュレート

配置位置を求める処理①

配置位置を求める処理②

以上により、平均計算量は O(BN\sqrt{N}(N+P\log P)) に落ちました。

速さの代償として、MCMC・ビームサーチ・MCTSの全てで以下のようなコードが延々と記述されることになります。*12 提出コードで使っているintrinsics関数（ _mm で始まるSIMD系関数）の数を数えたら268個ありました。*13 さすがにここまで多いと結構大変です。

github.com

// 長方形がどこに置かれるかを調べる
let mut x0 = _mm256_setzero_si256();

for (p_x1, p_y0, p_y1) in izip!(placements_x1, placements_y0, placements_y1) {
    // やっていることはジャッジコードと同じ
    //
    // let x = if y0.max(p_y0) < y1.min(p_y1) {
    //     p_x1
    // } else {
    //     0
    // };
    //
    // x0 = x0.max(x);
    let max_y0 = _mm256_max_epu16(y0, p_y0.load());
    let min_y1 = _mm256_min_epu16(y1, p_y1.load());

    // AVX2には符号なし16bit整数の比較命令がないので符号ありで代用……
    let gt = _mm256_cmpgt_epi16(min_y1, max_y0);
    let x = _mm256_and_si256(p_x1.load(), gt);
    x0 = _mm256_max_epu16(x0, x);
}

let x1 = _mm256_add_epi16(x0, rect_w);

ちなみにコンパイル後はこんな感じのアセンブリになります。 ymm といっぱい出ているのは256bitレジスタを活用できている証拠です。*14

.LBB140_9:
 vpmaxuw ymm1, ymm5, ymmword, ptr, [r10, +, r9, -, 32]
 vpminuw ymm2, ymm7, ymmword, ptr, [rbx, +, r9, -, 32]
 vpcmpgtw ymm1, ymm2, ymm1
 vpand   ymm1, ymm1, ymmword, ptr, [r11, +, r9, -, 32]
 vpmaxuw ymm2, ymm5, ymmword, ptr, [r10, +, r9]
 vpminuw ymm3, ymm7, ymmword, ptr, [rbx, +, r9]
 vpmaxuw ymm1, ymm6, ymm1
 vpcmpgtw ymm2, ymm3, ymm2
 vpand   ymm2, ymm2, ymmword, ptr, [r11, +, r9]
 add     r12, 2
 vpmaxuw ymm6, ymm1, ymm2
 add     r9, 64
 cmp     r13, r12
 jne     .LBB140_9

バイトニックソート

P 要素のソートを普通に行うと O(P \log P) ですが、バイトニックソートというアルゴリズムを使うと比較・スワップを並列に行うことができ、 P 要素のソートが O( (\log P)^2) ステップで完了します。具体的な実装はGitHubに置きました。

github.com

実装すると、ランダム16要素に対するソートが標準ライブラリの5倍くらいの速度で完了するようになりました。*15

これにより、平均計算量は O(BN\sqrt{N}(N+(\log P)^2)) に落ちました。*16 最終的にseed=1でもジャッジサーバー上でビーム幅2000程度は出せるようになっています。ソート部分がどのくらいボトルネックになっているかは謎ですが、速いと気持ちが良いのでよしとします。

パラメータチューニング

以下の手順に従ってパラメータチューニングを行いました。入力パラメータ \boldsymbol{p}=\lbrack N, T, \sigma\rbrack^T に対して最適なパラメータベクトル \boldsymbol{q} を求める関数 \boldsymbol{q} = f(\boldsymbol{p}) をガウス過程回帰で求めるのが目的です。Google Cloudで88コアのC3インスタンスを2つ使用してチューニングを行っています。*17

  1. (N, T, \sigma) をランダムサンプリング
  2. その設定でテストケースを生成
  3. 生成した入力に対してoptunaで最適なパラメータを求める
  4. 1.～3.を繰り返し、データが十分集まったらガウス過程回帰を用いて関数 f を求める
  5. 実行時は入力パラメータ \boldsymbol{p}=\lbrack N, T, \sigma\rbrack^T に対して \boldsymbol{q} = f(\boldsymbol{p}) によりパラメータを決定する

できなかったこと

箱同士の大小関係を考慮した推定

推定パートでは箱をL字に並べて x,y 方向長さを計測していましたが、下図のように配置することで箱の大小関係も調べることができます。下図の例の場合、箱16を箱4に沿って右から投入しており、箱16が箱12に引っかかれば箱12の方が大きく、引っかからなければ箱4の方が大きいことが分かります。

単純にL字に並べるよりも多くの情報が得られるはずですが、のちのMCMCパートで適切なサンプリングを行うことが難しく断念しました。*18

x, y方向の長さに加えて箱4と箱12の大小関係も観測

AVX512

AVX512を使えば32セット並列にシミュレーションが可能ですが、Rustはnightly版コンパイラでないとAVX512に対応していないため断念しました。*19

やらない方がよかったこと

デスノート

コンテスト中の貴重な4時間を溶かしました。夜神月許せねえ……。

DEATH NOTE デスノート(1) (ジャンプ・コミックス)

DEATH NOTE デスノート(1) (ジャンプ・コミックス)

Amazon

コンテスト中のメモ

コンテスト中の考察メモはこちらです。かなり雑に書いていてまとまっていないので、暇な人だけご覧ください。

github.com

2024年まとめ

余談ですが、GP30 2024年Heuristic部門で1位になりました！年の中盤はちょっと微妙な成績だったのですが、序盤と終盤は自分にしては出来すぎな成績だったように思います。

GP30 2024年Heuristic部門

今年のAHCの成績はこんな感じでした。1桁か破滅かの極端な取り組み方しかしてないですね。まあ行けると思った回はちゃんと成果を出せているので良かったかなと思います。

2024年AHC成績

感想

推定回は大好きなのでとても楽しく取り組めました。特に今回は自分の持っている知識を全部使ったといっても良いくらい色々な手法を1つの解法に詰め込めたので、かなり満足感が高いです。SIMDの実装はパズルみたいでついついのめり込んでしまいますね。

コンテスト中はてっきり推定精度で負けているのかと思っていましたが、単純に詰め込みの最適化性能で負けていましたね。帯状に置く手法はあまり充填率を上げられないと思い込んでいましたが、想像以上によい詰め込みを実現できると聞いてびっくりしました。ビームサーチ中に下界が分かるのも強いなあと感じています。

2025年からはレーティング減衰が適用されるようですし、今後も継続して良い順位を取っていきたいですね。AWTF本戦出場も確定なので、そちらも頑張ります！

*1:初心者向けに書こうとすると時間がいくらあっても足りなさそうでした。ごめんなさい……。

*2:事前分布のσが大きいので、分散が比較的大きいseed=1でもσ=7992がσ=7482に下がる程度ですが、入れて損することはないでしょう。

*3:入力生成方法を盛大に誤読しており、Lがテストケースごとではなく長方形ごとに生成されると勘違いしていました……。テストケースごとにLが分かるならもう少し良い事前分布を設定したり、複数のモデルを用意したりといった工夫ができたかもしれません。

*4:という理解です。トレースという概念も分散共分散行列の固有値・固有ベクトルが表すものも今回初めて知ったので間違っていたらご指摘ください……。

*5:T字型にするほか、箱を削除する近傍も入れましたがどこまで効いているかは謎です。

*6:これは一番シンプルなメトロポリス法の話で、ギブスサンプリングやHMC法などはもう少し凝ったことをやります。

*7:0.05%スコアが下がると7位、0.1%スコアが下がると8位に落ちるのでそれなりに大きいです。

*8:今回の問題についてはBottom-Left法はあまり筋が良くなさそうでした。公式解説放送の解法等もご覧ください。

*9:なるべく正方形になるように詰め込むのが良いのですが、自由に詰め込めるスタイルだと細長い長方形になることがありました。

*10:できれば無駄な領域などを評価値に入れたかったのですが、高速に求める方法が思い付かず断念しました……。

*11:ちなみにMCTSをまともに使うのは3年半前のCodinGame（CodinGame Spring Challenge 2021参加記 - TERRYのブログ）以来だったので、実装方法の全てを忘れていました……。

*12:基本16並列処理ですが、MCMCの尤度計算など一部では32bit浮動小数点数を取り扱う関係で8並列x2回処理となります。

*13:コピペしている箇所もあるので実際はもっと少ないとは思います。

*14:余談ですが、コンパイラくんが自動でループアンローリングをして2個ずつ処理してくれているようです。えらい！

*15:Google Cloud C3インスタンス上での雑検証で155ns→31ns。

*16:Pが計算量のボトルネックになり得る箇所はソート以外にもスコアの重み付き計算などがありますが、2^i要素ずらして足し合わせることを繰り返すことでP要素の和の計算を log P 回の加算で求められることに注意してください。

*17:88コアインスタンスを延べ20～30時間くらい実行し、出費は4500円くらいでした。

*18:一度サンプルが迷子になってしまうと、尤もらしい状態に戻ってくるのが困難になるためと思われます。

*19:一応アセンブリを直接書くという手段もなくはないですが……。