はじめに

しばらくブログを書いてなかったのでAHC参加ネタを投稿。

先日Atcoder Heuristic Contest 030（AHC030）に参加し最終56位という結果だったが、提出した基本方針に二乗誤差を最小化する焼きなましが含まれていた。この方法では二乗誤差の計算を何度も行うことになるためその計算時間が肝になるのだが、線形代数ライブラリのEigenを使うと下手に差分計算するより愚直に全計算した方が速かったりする*1。とはいえ、少し工夫して差分計算にもEigenを活用すれば、状況によっては愚直計算よりも速くなる。

ということで、本記事では実際にAHC030で使ったEigenによる二乗誤差の差分計算の方法を説明する。このテクニックはAHC030に限らず二乗誤差を評価値として焼きなましする場合にも応用できる（かも？）。

AHC030の問題については以下を参照。 atcoder.jp

AHC030解法基本方針概要

解法基本方針についてはtwitterのリンクおよびそのスレッドを参照。

#AHC030 お疲れさまでした。暫定47位でした。

- ランダムにn回占った後、焼きなましでポリオミノの配置による各マスの合計値と占いで得た合計値の2乗誤差が最小になるポリオミノの配置を求めて出力
- 駄目だったらn'回更に占う

を正解するまで繰り返しました（epsとMが小さければ占いのみ行う）
— 甲斐性なし (@YamagenSakam) 2024年2月19日

二乗誤差最小化によるポリオミノの配置最適化

まずは評価値の計算を行列演算の形で表す。占いを $n$ 回行ったとして $i$ 回目の占いで得られる値を $y_i$ 、 $y_i$ を $i$ 番目の要素に持つ $n$ 要素のベクトルを $\mathbf{y}$ 、 $i$ 行目に $i$ 回目の占いにおいてどのマスを選んだか？という情報を持つ $n \times N^2$ の行列を $\mathbf{X}$ （各列がマスに対応し、選ばれたマスなら $1$ , そうでなければ $0$ ）とする*2。そして、 $\mathbf{w}$ を現状態のポリオミノの配置において、各マス何個のポリオミノが重なっているかを表す $N^2$ 要素のベクトルとし、これが最適化対象の変数である。これらを使って最適化問題として定式化すると、

$\displaystyle \begin{aligned} &{\rm arg}\min_{\mathbf{p}} \parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2 \\ &{\rm s.t. \ \ } \mathbf{w}はポリオミノの配置の組み合わせであり得る値 \end{aligned} \tag{1}$

となる。これを1つのポリオミノを選択して置く場所を変える焼きなましによって最適化を図る。そのため二乗誤差 $\parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2$ の計算を高速化できれば短い時間で焼きなましのループ回数が稼げ、より正解を出せる確率が上がる。

二乗誤差計算の高速化

Eigenの活用

C++においては行列計算を高速に実行できるライブラリEigenを活用すれば、愚直計算でもかなり速く計算できる。コードもこんな感じで1行で書ける。

double cost = (y - X * w).squaredNorm();

EigenではSIMD命令の活用、ループ変形等によりメモリアクセスの効率化などが行われているため演算回数だけ見れば遅くなりそうな処理も高速にできる場合がある。逆に差分計算はメモリアクセスが飛び飛びになりがちで、SIMD化もしにくいため処理時間上不利になることがある。

二乗誤差差分計算の行列表現

差分計算を行列演算で行うことを考える。ポリオミノの配置変更による $\mathbf{w}$ の変化を $\Delta \mathbf{w}$ とすると、二乗誤差の変化量 $D$ は以下のようになる。

$\displaystyle \begin{align} D &= \parallel \mathbf{y} - \mathbf{X} \mathbf{w} \parallel_2^2 - \parallel \mathbf{y} - \mathbf{X} \left(\mathbf{w} + \Delta \mathbf{w} \right) \parallel_2^2 \\ &=2 \left(\mathbf{y}^T \mathbf{X} \Delta \mathbf{w} - \mathbf{w}^T \mathbf{X}^T \mathbf{X} \Delta \mathbf{w} \right) - \Delta \mathbf{w}^T \mathbf{X}^T \mathbf{X} \Delta \mathbf{w} \end{align} \tag{2}$

この式の中で $\mathbf{X}^T \mathbf{X}$ は焼きなましの中で変化しないので、事前に計算しておくことができる。加えて、焼きなましの1ループでは1つのポリオミノしか動かさないので、 $\mathbf{w}$ の中で変化する要素はごく一部、つまり $\Delta \mathbf{w}$ の大半の要素は $0$ になる。このような多くの要素が $0$ のベクトルを扱うためEigenにはSparseVectorが用意されている。 $\Delta \mathbf{w}$ をSparseVector型として扱うことで、式（2）の計算がだいぶ速くなる。コードはこんな感じ。

Eigen::SparseVector<double> sdw(N * N);
vector<Eigen::Triplet<double>> trivec;
rep(i, N * N){
    if(abs(dw(i)) > 1e-5){
        sdw.insert(i) = dw(i);
    }
}
VectorXd XXdw = XX * sdw;
VectorXd Xdw = X * sdw;
double D = 2.0 * (y.dot(Xdw) - old_w.dot(XXdw)) - dw.dot(XXdw);