Lock in $30 Savings on PRO—Offer Ends Soon! ⏳

ヒューリスティックコンテストで機械学習しよう

nagiss
September 25, 2023

 ヒューリスティックコンテストで機械学習しよう

2023/8/10のDeNA/Go AI技術共有会の発表資料です。

nagiss

September 25, 2023
Tweet

More Decks by nagiss

Other Decks in Technology

Transcript

  1. AI 6 ▪ お題に対して、なるべく良いスコアを得るプログラムを書 いて提出する ▪ プログラムの実行時間制限は1テストケースあたり数秒 ▪ 長期開催(1〜2週間)と短期開催(4〜8時間)がある ▪

    今回のターゲットは長期 ▪ 長期開催はこれまでに 13 回 ▪ 余談 ▪ 「ヒューリスティックコンテスト」という言葉はAtCoderが作っ た (はず) ▪ それ以前は「マラソン」しかなかった (Topcoder由来) AtCoder Heuristic Contest (AHC) とは ※ レーティング付与対象のみカウント
  2. AI 7 ▪ RECRUIT 日本橋ハーフマラソン 2023夏(AtCoder Heuristic Contest 022) ▪

    明日から開催!!!!!!!乗り遅れるな!!!!!! ▪ 2023-08-11(金) 12:00 ~ 2023-08-20(日) 19:00 おすすめコンテスト https://atcoder.jp/contests/ahc022
  3. AI 10 ▪ AHC003 ▪ MCMCによるパラメータ推定 (1位) ▪ 多変量ガウス分布でモデル化して推定 (2位)

    ▪ 誤差逆伝播によるパラメータ推定 (3位) ▪ 線形回帰, UCB1 など (多くの上位解法) ▪ HACK TO THE FUTURE 2022 予選 ▪ MCMCによる推定 (1位) ▪ HACK TO THE FUTURE 2023 予選 (AHC016) ▪ Transformer Encoderによる予測 (1位) ▪ RECRUIT 日本橋ハーフマラソン 2023冬(AHC018) ▪ ガウス過程回帰で推定 (7位) ▪ いつも ▪ Optunaでパラメータチューニング ▪ 複数の解法を用意して、入力の変数の値からどの解法が最も適しているか予測 統計的手法が活躍した回はいっぱいある
  4. AI 12 ▪ 目的 ▪ 道案内アプリをリリースするので、目的地までの なるべく短い経路を出力するプログラムを書いてください ▪ 条件 ▪

    各道の長さは分かっていません ▪ 経路を提示するたびに到着までにかかった時間が得られ、 次回以降に道の長さを推測するヒントとして使えます ▪ 1テストケースにつき1000回経路を出力 AHC003 - 問題概要 (1/2) https://atcoder.jp/contests/ahc003/tasks/ahc003_a
  5. AI 14 ▪ MCMCによるパラメータ推定 (1位) ▪ 正規分布での近似から脱却したのはこの解法だけではなかろうか ▪ 多変量ガウス分布でモデル化して推定 (2位)

    ▪ 誤差逆伝播によるパラメータ推定 (3位) ▪ 線形回帰, UCB1 など (多くの上位解法) AHC003 - 上位解法 1位 https://qiita.com/contramundum/items/b945400b81536df42d1a 2位 https://speakerdeck.com/yos1up/ahc003-2wei-jie-fa
  6. AI 15 ▪ 目的 ▪ 1000個のタスクをチームメンバーに割り振って、なるべく早く 全て完了させてください ▪ 条件 ▪

    タスクとメンバーは、10~20種類の技能レベルを持ちます ▪ 技能が不足するほど、タスクの所要時間が長くなります ▪ タスクの各技能レベルは与えられますが、メンバーの各技能レベ ルは与えられないので、かかった時間から推定してください ▪ タスクには依存関係があります HTTF2022予選 - 問題概要 https://atcoder.jp/contests/future-contest-2022-qual/tasks/future_contest_2022_qual_a
  7. AI 16 ▪ MCMCによる推定 (1位など) ▪ AHC003のことがあったので1位以外もMCMCしようとしていた 人はちらほらいた ▪ AHC003と比べると推定の正確さは重要でなかった模様

    ▪ タスクの依存関係からを効率良く割り振りを決める部分の方が 重要だった HTTF2022予選 - 上位解法 https://qiita.com/contramundum/items/52609b5a4c943bc6a275
  8. AI 19 ▪ 目的 ▪ なるべく体力を消費しないで全ての家を水源と接続してください ▪ 条件 ▪ 200x200のグリッドに家と水源が

    配置されています。 ▪ 全てのマスに岩盤が存在し、好きな 順で好きな場所を掘削できます。 ▪ 岩盤には頑丈さが設定されており、 頑丈なほど掘削に体力を要します。 ▪ 頑丈さは掘削するまでわかりません。 AHC018 - 問題概要 画像引用元: https://atcoder.jp/contests/ahc018/tasks/ahc018_a
  9. AI 20 ▪ ガウス過程回帰で推定 (7位など) ▪ すごい ▪ 試し掘り等の戦略の方が重要だった ▪

    推定に時間をかけるとその分本質に時間をかけられなくなる罠 AHC018 - 上位解法 https://docs.google.com/presentation/d/1JEcyHLw8XrDqL4FHUGYIVQC63KSZ2eaHRjO0E2y1WeU
  10. AI 23 ▪ 1. リッジ回帰 + データ更新 ▪ 2. ガウス過程回帰

    + データ更新 ▪ 3. MCMC (メトロポリス法) ▪ 4. sklearnのランダムフォレスト埋め込み よく現れる問題の実装例を説明するよ
  11. AI 24 ▪ 問題文 ▪ 説明変数の行列 X と目的変数のベクトル y が与えられます。

    ▪ リッジ回帰してください。 ▪ ただし、データは1行ずつ渡され、データが渡されるごとに回帰 係数を求め直す必要があります。 リッジ回帰 - 問題の定義 これが役立つ問題 … AHC003 など ※ リッジ回帰の説明は省略します
  12. AI 25 ▪ リッジ回帰の係数は           で求められる ▪ 計算量は ▪ データが渡されるごとにこの計算を

    やり直す? ▪ もっと良い方法があります リッジ回帰 - 問題の定義 N … データ数 M … 説明変数の数 X … 説明変数のN×M行列 y … 目的変数のN次元ベクトル λ … リッジ回帰の正則化の係数 I … 単位行列 w … 回帰係数 ▪ 問題文 ▪ 説明変数の行列 X と目的変数のベクトル y が与えられます。 ▪ リッジ回帰してください。 ▪ ただし、データは1行ずつ渡され、データが渡されるごとに回帰 係数を求め直す必要があります。
  13. AI 26 ▪ 逆行列を安易に計算するべきではない ▪ 本当に必要なのは   ではなく大体 ▪ 行列分解を使うと効率的に計算できる ▪ 今回は修正コレスキー分解

    (LDL 分解) を使う ▪ 修正コレスキー分解 ▪ 対称正定値行列  を       と分解する ▪ 回帰問題は対称正定値行列が 出て来がち ▪      を求める計算量が逆行列経由と比較して ⅙ 程度 リッジ回帰 - 解法の前に L … 下三角行列、対角成分は全て1 D … 対角行列 ※正定値じゃないと分解できないわけではない
  14. AI 27 ▪ 分解の実装はこれだけ ▪ 対角成分に D の対角成分を入れている ▪ Python

    なら Scipy にも実装されている ▪ C++ なら Eigen に実装されている リッジ回帰 - 修正コレスキー分解 def ldl(A): n = len(A) diag = A.diagonal() A *= np.tri(n) for i in range(n): A[i:, i] -= A[i:, :i] @ (diag[:i] * A[i, :i]) A[i + 1 :, i] *= 1.0 / diag[i] ← 参照なので、A が変更されると diag の中も変わる AtCoderの次回のLanguage UpdateでC++にEigenが入る
  15. AI 28 ▪ 分解された L, D から以下のように求められる ▪ 計算量は O(n2)

    リッジ回帰 - 修正コレスキー分解を使って   を求める def ldl_solve(LD, b): n = len(LD) x = b.copy() for i in range(n): x[i] -= np.dot(LD[i, :i], x[:i]) x /= LD.diagonal() for i in range(n - 1, -1, -1): x[i] -= np.dot(LD[i + 1 :, i], x[i + 1 :]) return res
  16. AI この形の行列は対称正定値 29 ▪ 修正コレスキー分解を使うと、リッジ回帰           は以下のように実装できる ▪ データの追加をどのように扱うか? リッジ回帰

    - 実装 class Ridge: def __init__(self, X, y, lambda_): """O(NM^2 + M^3)""" self.n_data, self.n_features = X.shape self.LD = X.T @ X # O(NM^2) np.fill_diagonal(self.LD, self.LD.diagonal() + lambda_) if self.n_data != 0: ldl(self.LD) # O(M^3) self.b = X.T @ y # O(NM) def get_weight(self): """O(M^2)""" return ldl_solve(self.LD, self.b)
  17. AI 30 ▪ データ 1 行 (xT, y) が追加されると、           は、

              となる ▪ 右側はどうにでもなるが、問題は左側 ▪ 単純に行列分解や逆行列を再計算すると O(M3) の計算量 ▪ 実は、A の行列分解が計算済みなら、A+xxT の行列分解を O(M2) で計算できる ▪ ランク 1 更新と呼ばれる操作 ▪ 対称行列でなければ xyT でも良い リッジ回帰 - データの追加
  18. AI 31 ▪ 修正コレスキー分解の場合の実装 ▪ リッジ回帰で、データを追加する実装 リッジ回帰 - 修正コレスキー分解のランク 1

    更新 def ldl_rank_one_update(LD, x): a = 1.0 for k in range(len(x)): xk = x[k] dk = LD[k, k] LD[k, k] += a * xk * xk x[k + 1 :] -= xk * LD[k + 1 :, k] LD[k + 1 :, k] += xk * a / LD[k, k] * x[k + 1 :] a *= dk / LD[k, k] class Ridge: def add_data(self, x, y): self.b += x * y ldl_rank_one_update(self.LD, x) Q. なんでこうなるの? A. わからん
  19. AI 32 ▪ この問題を O(NM2) で 解くことができた リッジ回帰 - まとめ

    ▪ 問題文 ▪ 説明変数の行列 X と目的変数のベクトル y が与えられます。 ▪ リッジ回帰してください。 ▪ ただし、データは1行ずつ渡され、データが渡されるごとに回帰 係数を求め直す必要があります。 N … データ数 M … 説明変数の数 勾配法の方がもっと効率が良い可能性もあります (cf. AHC003 の 5 位解法)
  20. AI 33 ガウス過程回帰 - 問題の定義 ▪ 問題文 ▪ 説明変数の行列 X

    と目的変数のベクトル y が与えられます。 ▪ ガウス過程回帰してください。 ▪ ただし、データは1行ずつ渡され、データが渡されるごとに回帰 係数等を求め直す必要があります。 これが役立つ問題 … AHC018 AHC003 の 2 位解法もやることは近いはず
  21. AI 34 ▪ 大雑把に言えば (大雑把すぎてほぼ嘘) ▪ いくつかのデータ点の 間をそれっぽく補間する ▪ ここまでは

    カーネル回帰も同じ ▪ 補間の信頼度を 計算してくれる ガウス過程回帰 - そもそもガウス過程ってなんだよ (1/2) 画像引用元: https://scikit-learn.org/stable/auto_examples/gaussian_process/plot_gpr_noisy_targets.html
  22. AI 35 ▪ もうちょっと詳しく ▪   の事前分布として多変量ガウス分布+ノイズを仮定 ▪   と

      の共分散は以下で表されるとする ▪   ▪ x 同士が近いほど共分散が大きい ▪ これは RBF カーネルと呼ばれるものだが、 他のカーネル関数を使っても良い ▪ 上記の仮定の元で、未知のデータ点に対して y の分布を求める ▪ これはガウス分布になる ガウス過程回帰 - そもそもガウス過程ってなんだよ (2/2) X … 説明変数の行列 y … 目的変数のベクトル x i T … Xのi行目 γ … ハイパーパラメータ 平均は0を仮定するが、それ以外の場合 最初に引いて最後に足せば良い 未知のデータが複数ある場合はそれらの多変量ガウス分布を 求められるが、今回その共分散は求めないことにする
  23. AI 36 ▪ 3 変量正規分布 N(μ,Σ) から y 0 ,

    y 1 , y 2 がサンプルされる ▪ μ=[0,0,0], Σ=[[1, 0.8, 0.84], [0.8, 1, 0.8], [0.84, 0.8, 1]] ▪ これ以上の情報が無ければ y 0 , y 1 , y 2 は大体 0 と予想できる ▪ Q. y 0 =-1, y 1 =0 と判明した時、 y 2 の事後分布は? ガウス過程回帰 - 例 多変量正規分布の周辺分布 0.8 0.8 0.4096
  24. AI 37 ▪ 3 変量正規分布 N(μ,Σ) から y 0 ,

    y 1 , y 2 がサンプルされる ▪ μ=[0,0,0], Σ=[[1, 0.8, 0.84], [0.8, 1, 0.8], [0.84, 0.8, 1]] ▪ これ以上の情報が無ければ y 0 , y 1 , y 2 は大体 0 と予想できる ▪ Q. y 0 =-1, y 1 =0 と判明した時、 y 2 の事後分布は? ▪ A. N(0.64, 0.462) ガウス過程回帰 - 例
  25. AI 38 ▪ 以下の式を計算すれば良い ▪ 平均 ▪ 分散 ガウス過程回帰 -

    どうやって計算するの? X … 説明変数の行列 y … 目的変数のベクトル x i T … Xのi行目 x*T … 予測したいデータ K … i行j列の要素はk(x i T, x j T) k* … i番目の要素はk(x i T, x*) k** … k(x*, x*) λ … ノイズの分散 I … 単位行列 カーネル関数が適切であれば対称正定値
  26. AI 39 ▪ 修正コレスキー分解の行の追加 ▪ O(N2) ガウス過程回帰 - 実装 (1/3)

    def ldl_add_row(LD, x): n = len(LD) assert len(x) == n + 1 for i in range(n): x[i] -= np.dot(LD[i, :i], x[:i]) x[:n] /= LD.diagonal() x[n] -= np.dot(np.square(x[:n]), LD.diagonal()) なんでこれで動くのかはよくわからない 英語版 Wikipedia のコレスキー分解のページにあっ た式を弄ったらできた
  27. AI 40 ▪ 初期化とカーネル関数 ガウス過程回帰 - 実装 (2/3) class GaussianProcess:

    def __init__(self, X, y, n_max_data, gamma, lambda_): n_data, self.n_features = X.shape self.n_data = n_data self.gamma = gamma self.lambda_ = lambda_ self.X = np.empty((n_max_data, self.n_features)) self.X[:n_data] = X self.y = np.empty(n_max_data) self.y[:n_data] = y self.LD = np.empty((n_max_data, n_max_data)) self.LD[:n_data, :n_data] = self.kernel_func(X, X) np.fill_diagonal(self.LD[:n_data, :n_data], self.LD[:n_data, :n_data].diagonal() + lambda_) ldl(self.LD[:n_data, :n_data]) def kernel_func(self, A, B=None): if B is None: return np.ones(len(A)) D = A[:, None, :] - B[None, :, :] np.square(D, out=D) K = D.sum(2) K *= -self.gamma np.exp(K, out=K) return K
  28. AI 41 ▪ 予測とデータ追加 ガウス過程回帰 - 実装 (3/3) class GaussianProcess:

    def predict(self, x): k_star = self.kernel_func(x[None], self.X[: self.n_data]).ravel() K_inv_k_star = ldl_solve(self.LD[: self.n_data, : self.n_data], k_star) mean = np.dot(K_inv_k_star, self.y[: self.n_data]) k_star_star = self.kernel_func(x[None]).reshape(()) var = k_star_star - np.dot(K_inv_k_star, k_star) return mean, var def add_data(self, x, y): n_data = self.n_data self.X[n_data] = x self.y[n_data] = y k_star = self.kernel_func(x[None], self.X[:n_data]).ravel() self.LD[n_data, :n_data] = k_star k_star_star = self.kernel_func(x[None]).reshape(()) + self.lambda_ self.LD[n_data, n_data] = k_star_star ldl_add_row(self.LD[:n_data, :n_data], self.LD[n_data, : n_data + 1]) self.n_data += 1
  29. AI 42 ▪ 以下の計算量で処理できた ▪ データ追加 O(N2) ▪ 未知データに対する予測 O(N2)

    ガウス過程回帰 - まとめ N … 既知データの数 ▪ 問題文 ▪ 説明変数の行列 X と目的変数のベクトル y が与えられます。 ▪ ガウス過程回帰してください。 ▪ ただし、データは1行ずつ渡され、データが渡されるごとに回帰 係数等を求め直す必要があります。 ※ N’ 件の未知データに対して平均だけを予測する場合、   計算量を O(N’N + N2) にできる。この計算はカーネル回帰に等しい。
  30. AI 43 ▪ これ結構扱いづらい ▪ ハイパーパラメータに敏感 ▪ → L-BFGS-B、共役勾配法、MCMC 等により

      尤度を最大化するハイパーパラメータを探索 ▪ → カーネル関数を学習 ▪ 計算量が大きい ▪ → 補助変数法などによる近似 ▪ むずい ▪ めんどい ガウス過程回帰 - 所感 このあたりよくわかってない 理想(sklearnにある例。sklearnのガウス過程 は自動でハイパーパラメータも最適化される) ハイパーパラメータが適当だと こうなる
  31. AI 44 MCMC - 問題の定義 ▪ 問題文 ▪ [-2.0,2.0]の連続一様分布からサンプルされた隠しパラメータ  があります。

    ▪ 標準正規分布によるノイズが加わった として、-0.5と0.1の2つ の観測値を得ました。 ▪  はどれくらいだと推測できますか? ▪ また、      はどれくらいだと推測できますか? 事後期待値 EAP これが役立つ問題 … AHC003, HTTF2022予選
  32. AI 45 ▪ は-0.2くらいな気がするけど      の期待値は想像し辛い ▪ 方針:   の事後分布    

    から  を大量にサンプルすることで解く MCMC - 問題の定義 ▪ 問題文 ▪ [-2.0,2.0]の連続一様分布からサンプルされた隠しパラメータ  があります。 ▪ 標準正規分布によるノイズが加わった として、-0.5と0.1の2つ の観測値を得ました。 ▪  はどれくらいだと推測できますか? ▪ また、      はどれくらいだと推測できますか? 事後期待値 EAP
  33. AI 46 ▪ 方針:   の事後分布     から  を大量にサンプルする ことで解く

    ▪ Q1.     ってどうやって計算するのさ ▪ A1. ベイズの定理を使おう! MCMC - ベイズの定理 比例 θ … 隠しパラメータ x … 観測 p(θ|x) … 事後分布 (xが観測され た時のθの分布) p(θ) … 事前分布 (xを観測する前 からわかっていたθの分布) p(x|θ) … 尤度 (そのθから、xがど れくらいの確率で観測されるか) 事後分布 尤度 事前分布 ※ この式で求められるのは「事後分布に比例する関数」であって   正確な事後分布ではないが、メトロポリス法を使うにはこれで十分 定数倍は無視していいとか、まるで競プロerみたいですね
  34. AI def log_likelihood(theta): x1, x2 = -0.5, 0.1 p1 =

    -(x1 - theta) * (x1 - theta) / 2.0 p2 = -(x2 - theta) * (x2 - theta) / 2.0 return p1 + p2 def log_prior(theta): if -2.0 <= theta <= 2.0: return 0.0 else: return -math.inf def log_posterior(theta): return log_likelihood(theta) + log_prior(theta) 47 MCMC - ベイズの定理 - 実装 比例 事後分布 尤度 事前分布 ※ 色々と都合が良いので対数で実装している ※ 所々で定数倍を省略して簡略化している 同時確率(積)の対数なので和 そのθからxがサンプルされる確率密度(の定数倍の対数) 対数なので和
  35. AI 48 ▪ 方針:   の事後分布     から  を大量にサンプルする ことで解く

    ▪ Q2.     からどうやってサンプルするのさ ▪ A2. メトロポリス法を使おう! MCMC - メトロポリス法 MCMC の一種 ▪ 言葉で書くよりコード見た方が早い def metropolis(): sampled_theta = [] current_theta = 0.0 for _ in range(1000000): next_theta = current_theta + random.uniform(-0.5, 0.5) if math.log(random.random()) < log_posterior(next_theta) - log_posterior(current_theta): current_theta = next_theta sampled_theta.append(current_theta) 適当に初期化 current_thetaからnext_thetaが選ばれる確率が、 next_thetaからcurrent_thataが選ばれる確率と 等しくなるように近傍を選ぶ (重要) 遷移が棄却されてもリストに追加する
  36. AI 49 ▪ 方針:   の事後分布     から  を大量にサンプルする ことで解く

    ▪ サンプルできたのであとは統計量を計算するだけ ▪ それぞれ-0.19くらいと0.19くらいになる MCMC - 実装 mean_theta = sum(sampled_theta) / len(sampled_theta) mean_ramp_theta = sum(max(theta, 0) for theta in sampled_theta) / len(sampled_theta) の期待値 の期待値
  37. AI 50 ▪ メトロポリス法でサンプルすることで  が -0.19 くらい、    が 0.19 くらいと推測できた MCMC

    - まとめ ▪ 問題文 ▪ [-2.0,2.0]の連続一様分布からサンプルされた隠しパラメータ  があります。 ▪ 標準正規分布によるノイズが加わった として、-0.5と0.1の2つ の観測値を得ました。 ▪  はどれくらいだと推測できますか? ▪ また、      はどれくらいだと推測できますか?
  38. AI 53 ▪ こんな感じで実装すれば、 ランダムフォレスト - C++側の実装 // 回帰木 template

    <int max_node_count = 32> struct DecisionTree { array<int, max_node_count> children_left, children_right, feature; array<double, max_node_count> threshold, value; template <typename Array> double Predict(const Array& x) const { auto node = 0; while (children_left[node] != -1) node = (x[feature[node]] <= threshold[node] ? children_left : children_right)[node]; return value[node]; } }; // ランダムフォレスト回帰 template <int n_trees, int max_tree_size> struct RandomForest { array<DecisionTree<max_tree_size>, n_trees> trees; template <typename Array> double Predict(const Array& x) const { auto res = 0.0; for (const auto& tree : trees) res += tree.Predict(x); return res / n_trees; } };
  39. AI 54 ▪ こんな感じに埋め込める ▪ 結構かんたん ランダムフォレスト - Python側の実装 def

    dump_forest(reg): print(f"RandomForest<{reg.n_estimators}, 32>" "{") for tree in reg.estimators_: tree = tree.tree_ print(" DecisionTree<32>{") print(" {" + ",".join(map(str, tree.children_left)) + "},") print(" {" + ",".join(map(str, tree.children_right)) + "},") print(" {" + ",".join(map(str, tree.feature)) + "},") print(" {" + ",".join(map(str, tree.threshold)) + "},") print(" {" + ",".join(map(str, tree.value.ravel())) + "},") print(" },") print("}") reg = RandomForestRegressor() reg.fit(X, y) dump_forest(reg)
  40. AI 55 ▪ 1. リッジ回帰 + データ更新 ▪ 2. ガウス過程回帰

    + データ更新 ▪ 3. MCMC (メトロポリス法) ▪ 4. sklearnのランダムフォレスト埋め込み よく現れる問題の実装例を説明したよ
  41. AI 57 ▪ MC Digital プログラミングコンテスト2022 (AtCoder Heuristic Contest 008)

    ▪ グリッド上で逃げるペットを追い込むゲーム ▪ estie プログラミングコンテスト2022 (AtCoder Heuristic Contest 014) ▪ グリッド上に四角形をたくさん作るゲーム これ強化学習強いんじゃないの?って問題が時々出る 画像引用元: https://img.atcoder.jp/ahc008/f828b9475ffb41d54f05619db6ccbd4f.html?lang=ja&show=example https://atcoder.jp/contests/ahc014/tasks/ahc014_a 内の「例を見る」 参考: https://twitter.com/search?q=%23AHC008%20%23visualizer&src=typed_query&f=live https://twitter.com/search?q=%23AHC014%20%23visualizer&src=typed_query&f=live