ベイズ最適化: 適応的実験計画の基礎と実践 輪読会 3.5 トンプソン抽出

https://www.kindaikagaku.co.jp/book_list/detail/9784764906631/ の輪読会で発表した資料です。

Transcript

1. ベイズ最適化 適応的実験計画の基礎と実践 輪読会 第３回 3.5 トンプソン抽出 tatsukawa 1/20
2. 多腕バンディット問題の選択アルゴリズム アイデア： 1. 事後分布のサンプリング 2. を次の点に設定 これまでの獲得関数と異なり、関数のサンプリングに起因するランダム性を含む これが原因で素直に勾配法などで獲得関数の最小値が求められない トンプソン抽出 (Thompson Sampling) f ∼ p(f ∣D ​ ) n (1) x = ∗ arg min ​ f (x) x∈X トンプソン抽出の獲得関数は α ​ (x) = TS −f (x) (n) で定義される。ただし、 である。 定義 3.5.1 (トンプソン抽出) f ∼ (n) GP(μ ​ , κ ​ ) n n (2) (1): はトライアル までのデータ履歴 (2): は を与えたときのガウス過程の事後分布を表す D ​ = n {(x ​ , y ​ )} ​ i i i=1 n n GP D ​ n 2/20
3. サンプリングした関数毎に最小点が変わる（ランダム性がある）ため、 で表すことが可能 右辺は目的関数の最小値の従う事後分布からサンプリングすることを式で表したもので、次節のエン トロピー探索獲得関数において重要な役割を持つ 獲得関数の解釈 : エントロピー探索との関係 「獲得関数の最大化 によって得られる 」 = 「 」 (1) x ​ n+1 x ​ ∼ n+1 p ​ (x ​ ∣D ​ ) min ∗ n (2) p ​ min (3) (1) : 獲得関数はサンプリングした目的関数に負の符号をつけているので目的関数の最小化をしようとすると獲得関数の最大化になる。ややこしい。 (2) : 結局のところ、 1. 関数のサンプリング 2. 最小化 の 2 プロセスをひとまとめにして書いたのが右辺というだけ (3) : 予測エントロピーの獲得関数において新しい観測点が得られたときのエントロピーの減少度合いを測るときのエントロピーの確率分布に を採用 p ​ min 3/20
4. 1. 素朴な方法 計算が単純で実装が楽 ガウス過程の学習の計算量が 事後分布のサンプリングの計算量が 勾配法が使えない 2. 疎スペクトル近似に基づく方法 （ と分解して効率よく計算する） ガウス過程の学習の計算量が 事後分布のサンプリングの計算量が 勾配法が使える 平行移動不変カーネルにしか使えない 獲得関数の実装方法 O(n ) 3 O(L ) 3 k(x, y) = ϕ(x) ϕ(y) ⊤ O(nm ), m < 2 n (1) O(mL) (2) (1): が より小さすぎると外層領域での振る舞いが悪化する variance starvation が起こる (2): ほとんどの場合で不変カーネルが使われるのでそこまで問題ではないかも m n 4/20
5. アイデア：入力空間 を 個の点集合（テキストだと 個）で近似してトンプソン抽出 ​ ​ X f (n) (y ​ , … , y ​ ) 1 L ≃ {x } ​ i i=1 L ≃ ​ : x ​ ↦ y ​ f ~(n) i i ∼ N (μ ​ , κ ​ + σ I) n n 2 獲得関数も と近似する 獲得関数の実装方法： 1. 素朴な方法 X L Δ α (x) ≃ TS − ​ (x ​ ) f ~(n) i ​ ​ ​ α ​ (x) x∈X arg max TS ≃ ​ − ​ (x ​ ) x ​ ∈{x ​ ,…,x ​ } i 1 L arg max ( f ~(n) i ) = ​ (x ​ ) x ​ ∈{x ​ ,…,x ​ } i 1 L arg min f ~(n) i 計算量は （ 共分散行列の逆行列でのコレスキー分解） トンプソン抽出の素朴な実装方法 O(L ) 3 ∵ 5/20
6. ボホナーの定理を活用した効率的なサンプリング方法 不変カーネル であることが前提 のままだと少々扱いづらいので、確率密度関数 として扱いたい → そのために の分母を求めたい → 実はこれは なので簡単に求まる とすると、 と書ける はy軸に対称 （ ) から言える） p (−ω) = s σ ​ e k(x)dx = 0 −2 ∫ −2πix (−ω) ⊤ σ ​ e k(y)dy = 0 −2 ∫ −2πiy ω ⊤ p ​ (ω), y := s −x 獲得関数の実装方法： 2. 疎スペクトル近似に基づく方法 k(x, y) = k(x − y) 連続関数 が半正定値である ある非負測度 が存在して、 が成立 定理 3.5.1 (ボホナーの定理) k : R → d R ⇔ s(⋅) k(x) = e s(ω)dω ∫ 2πix ω ⊤ (1) s(ω) p ​ (ω) s p ​ (ω) = s ​ s(ω )dω ∫ ′ ′ s(ω) k(0) = e s(ω)dω = ∫ 2πi0ω s(ω)dω ∫ k(0) = σ ​ 0 2 k(x) = σ ​ e p ​ (ω)dω = 0 2 ∫ 2πix ω ⊤ s σ ​ E[e ] 0 2 2πix ω ⊤ p ​ (ω) s k(x) = k(−x (1): のフーリエ変換が であるということ : s k k = F [s] 6/20
7. 個のサンプル を使ってカーネル関数を近似 ​ ​ k(x, y) = σ ​ E ​ [e ] 0 2 p ​ s 2πi(x−y) ω ⊤ = σ ​ E ​ e e 0 2 p ​ s [ 2πix ω ⊤ ( 2πiy ω ⊤ ) ∗ ] ≃ ​ ​ e e + e e 2m σ ​ 0 2 j=1 ∑ m [ 2πix ω ​ ⊤ j ( 2πiy ω ​ ⊤ j ) ∗ { −2πix ω ​ ⊤ j ( −2πiy ω ​ ⊤ j ) ∗ } ∗ ] = ​ ​ Re e e m σ ​ 0 2 j=1 ∑ m [ 2πix ω ​ ⊤ j ( 2πiy ω ​ ⊤ j ) ∗ ] = ​ ​ cos (2π(x − y) ω ​ ) m σ ​ 0 2 j=1 ∑ m ⊤ j (モンテカルロ法による近似 ) (1) 不変カーネルの式変形 m ω ​ ∼ j p ​ (ω) s (1): 前ページの が 軸対称であることを利用して と を MC のサンプル平均に使っている。また、このおかげでその後は実部だけ見れば良い p80 の (3.28) 式 の上から２行目の は の誤植 は複素共役で、 は虚数 （ とすると、 ） p ​ (ω) s y ω ​ j −ω ​ j ω ​ J ω ​ j ∗ i (e ) = iθ ∗ cos θ − i sin θ = e−iθ z + z = ∗ 2Re[z] ∵ z = a + bi z + z = ∗ a + bi + a − bi = 2a = 2Re[z] 7/20
8. 特徴量ベクトル の積でカーネル関数を表現した方が都合が良いため、もう少し式変形を進める に対し次の式が成立 : これを用いるとカーネル関数は ​ ​ k(x, y) ≃ ​ ​ cos (2π(x − y) ω ​ ) m σ ​ 0 2 j=1 ∑ m ⊤ j = ​ ​ E ​ 2 cos (2π(x ω ​ ) + b) cos (2π(y ω ​ ) + b) m σ ​ 0 2 j=1 ∑ m b [ ⊤ j ⊤ j ] ≃ ​ ​ 2 cos (2π(x ω ​ ) + b ​ ) cos (2π(y ω ​ ) + b ​ ) m σ ​ 0 2 j=1 ∑ m ⊤ j j ⊤ j j = ϕ(x) ϕ(y). ⊤ 不変カーネルの式変形 ϕ(x) b ∼ U[0, 2π] (1) E ​ [2 cos (2πx ω ​ + b) cos (2πy ω ​ + b)] b ⊤ j ⊤ j (1): テキストだと からサンプリングしているが、自分が見た関連論文の殆どが で取るようにしていたため、そちらに合わせた。そのためテキス トでは の係数 に がかかっていたが、このスライドでは としている。 U [0, 1] U [0, 2π] b 2π 1 8/20
9. 二乗指数カーネル の は ​ ​ ​ s(ω) = F[k(x)] = ​ k(x)e dx ∫ Rd −2πix ω ⊤ = σ ​ ​ e e dx 0 2 ∫ Rd −x Λx/2 ⊤ −2πix ω ⊤ = σ ​ e ​ e dx 0 2 −2π ω Λω 2 ⊤ ∫ Rd − ​ ∥x−2πiΛ ω∥ ​ 2 1 −1 2 2 = σ ​ (2π) ∣Λ∣ e . 0 2 n/2 −1/2 −2π ω Λω 2 ⊤ (平方完成) なる確率密度関数を求める だと思うと、係数は なの で、 で全体を割れば OK. 例：二乗指数カーネルのフーリエ変換 k(x, y) = k(x − y) = σ exp(−(x − 0 2 y) Λ(x − ⊤ y)/2) s(ω) s(ω) ∝ p ​ (ω) s C = −1 4π Λ 2 −1 ​ = ​ (2π) ∣C∣ n 1 ​ = (2π) (1/(4π )) n/2 2 n/2 ∣Λ∣−1/2 (2π) ∣Λ∣ n/2 −1/2 σ ​ 0 2 テキストの (3.31) の係数は多分誤植 と思っていたら元論文の Sparse Spectrum Gaussian Process Regression でも同じ間違いをしていた。 9/20
10. 疎スペクトル近似を用いて事後平均関数 と事後共分散関数 を書き直す ​ ​ μ ​ (x) n κ ​ (x, x ) n ′ = k ​ (x) (K ​ + σ I ​ ) y ​ n ⊤ n 2 n −1 n ≃ ϕ(x) Φ ​ (Φ ​ Φ ​ + σ I ​ ) y ​ ⊤ n ⊤ n n ⊤ 2 n −1 n = k(x, x ) − k ​ (x) (K ​ + σ I ​ ) Φ ​ ϕ(x ) ′ n ⊤ n 2 n −1 n ′ ≃ ϕ(x) ϕ(x ) − ϕ(x) Φ ​ (Φ ​ Φ ​ + σ I ​ ) Φ ​ ϕ(x ) ⊤ ′ ⊤ n ⊤ n n ⊤ 2 n −1 n ′ = ϕ(x) {I ​ − Φ ​ (Φ ​ Φ ​ + σ I ​ ) Φ ​ }ϕ(x ) ⊤ m n ⊤ n n ⊤ 2 n −1 n ′ 目的関数の事後分布 ​ ​ k(x, y) ϕ(x) = ϕ(x) ϕ(y) ⊤ = [ ​ cos {(2πx ω ​ ) + b ​ }] ​ ∈ R , Φ ​ = [ϕ(x ​ )] ​ ∈ R ​ m 2σ ​ 0 2 ⊤ j j j m n i i n×m 疎スペクトル近似 μ ​ n κ ​ n 10/20
11. 事後共分散関数がまだ複雑なので次の定理を用いて単純にする の に Woodbuly の式を適用 とすると、 (A + U CV ) = −1 I ​ − m I Φ ​ (σ I ​ + m n ⊤ 2 n Φ ​ Φ ​ )Φ ​ = n n ⊤ n Ψ ​ n −1 と書けるので、 . を使って平均関数をもう少し整理することが可能 ​ ​ Ψ ​ Φ ​ n −1 n ⊤ μ ​ (x) n = ⋯ = σ Φ ​ (Φ ​ Φ ​ + σ I ​ ) 2 n n n ⊤ 2 n −1 ≃ ϕ(x) Φ ​ (Φ ​ Φ ​ + σ I ​ ) y ​ = ϕ(x) Ψ ​ Φ ​ y ​ ⊤ n ⊤ n n ⊤ 2 n −1 n ⊤ n −1 n ⊤ n 目的関数の事後分布 に対し以下が成立 (A + U CV ) = −1 A − −1 A U (C + −1 −1 V A U ) V A −1 −1 −1 定理 3.5.2 (Woodbury の恒等式) A ∈ R , C ∈ m×m R , U ∈ n×n R , V ∈ m×n Rn×m κ ​ (x, x ) = n ′ ϕ(x) { ​ }ϕ(x ) ⊤ Ψ ​ と定義しておく n −1 ​ I ​ − Φ ​ (Φ ​ Φ ​ + σ I ​ ) Φ ​ m n ⊤ n n ⊤ 2 n −1 n ′ Ψ ​ n −1 A = I ​ , U = m Φ ​ , C = n ⊤ σ I ​ , V = −2 n Φ ​ n Ψ ​ = n A + U CV = I ​ + m σ Φ ​ Φ ​ −2 n ⊤ n Ψ 11/20
12. 計算量の改善 ( ) 計算量が支配的な部分は で、これは の計算量は まとめ：疎スペクトル近似による事後分布の近似 ​ ​ μ ​ (x) n κ ​ (x, x ) n ′ Ψ ​ n Φ ​ n ϕ(x) ≃ σ ϕ(x) Ψ ​ Φ ​ y ​ −2 ⊤ n −1 n ⊤ n ≃ ϕ(x) Ψ ​ ϕ(x ) ⊤ n −1 ′ = I ​ + σ Φ ​ Φ ​ ∈ R m −2 n ⊤ n m×m = [ϕ(x ​ )] ​ ∈ R i i n×m = [ ​ cos (2πx ω ​ ) + b ​ ] ​ ∈ R ​ m 2σ ​ 0 2 ⊤ j j j m 疎スペクトル近似による GP の事後確率 ∵ m < n Φ ​ Φ ​ n ⊤ n O(nm ) 2 Ψ−1 O(m ) 3 12/20
13. 当たり前だが だと疎スペクトル分解の嬉しさは殆どない（ が支配的） 各点の分散も計算したが殆ど小さくて見えなくなった これが実装ミスなのかこういうものなのかよくわかっていない 実装： を変えたときの平均関数 m m = 1000 ∵ O(m ) 3 13/20
14. 元論文 : "The approximation is equivalent to replacing the original spectrum by a set of Dirac deltas of amplitude distributed according to . Thus, we “sparsify” the spectrum of the GP" 元のスペクトル密度関数 を に従って分布する振幅 のデルタ関数の集合で置換 ↑ これをスパース化と呼んでいる 余談：疎スペクトル近似の疎は何を意味しているか s(ω) σ ​ 0 2 p ​ (ω) s (1) s(ω) p ​ (ω) s σ ​ 0 2 (1): 記号はこのスライドに合わせて変更した 14/20
15. 今までは事後分布の計算量が改善されたことを言っただけ トンプソン抽出の計算量改善については何も触れていない 実は疎スペクトル近似を使うことで、 のサンプリングも効率的にできる 1. 素朴な方法と同じように入力空間 を 個の点集合で近似 2. 個の入出力 は疎スペクトル近似した事後分布 に従うはず ​ ​ ​ ​ μ ~ n ​ κ ~ n Φ ​ L = σ ϕ(x ​ )Ψ ​ Φ ​ y ​ ​ = σ Φ ​ Ψ ​ Φ ​ y ​ [ −2 i n −1 n ⊤ n ] i=1 L −2 L n −1 n ⊤ n = [ϕ(x ​ ) Ψ ​ ϕ(x ​ )] ​ j = Φ ​ Ψ ​ Φ ​ i ⊤ n −1 j ij L n −1 L ⊤ = [ϕ(x ​ )] ​ ∈ R i i=1 L L×m 3. で表される が に従うことを利用して をサンプリング とすると、 で .. なので を からサンプリングした後 で OK. 3.5.4 疎スペクトル近似を用いた計算方法 f ∼ (n) GP(μ ​ , κ ​ ) n n X L L x ​ , y ​ i i N ( ​ ​ , ​ ) μ ~ n κ ~ n y = Ax + b, x ∼ N (μ, K) y N (Aμ + b, AKA ) ⊤ y ​ = i f (x ​ ) (n) i A = Φ ​ , b = L 0, μ = σ Ψ ​ Φ ​ y ​ , K = −2 n −1 n ⊤ n Ψ ​ n −1 Aμ + b = σ Φ ​ Ψ ​ Φ ​ y ​ , AKA = −2 L n −1 n ⊤ n ⊤ Φ ​ Ψ ​ Φ ​ L n −1 L ⊤ y = Ax ∼ N ( ​ ​ , ​ ) μ ~ n κ ​ n ~ θ N (μ, K) y ​ = i ϕ(x ​ ) θ, i = i ⊤ 1, … , L 15/20
16. 全体の計算量 のサンプリング : , （ と の行列積があるため） の計算 : 一般的には で設定 について線形まで計算量を削減できたが、 が 次元で各次元で 分割することを考えると となるため、 が大きい場合（e.g. ）ではまだ不十分 を有限個の点集合で近似せずに求めたい 3.5.4 疎スペクトル近似を用いた計算方法 トンプソン抽出： を計算してその最小値を求めることが目的 ​ ​ X θ ∈ Rm y Φ ​ L ≃ {x ​ } ​ i i=1 L ∼ N (σ Ψ ​ Φ ​ y ​ , Ψ ​ ) −2 n −1 n ⊤ n n −1 = Φ ​ θ L = [ϕ(x ​ )] ​ ∈ R i i=1 L L×m 疎スペクトル近似でトンプソン抽出 y O(nm + 2 mL) θ O(nm ) 2 m × m m × n y = Φ ​ θ L O(mL) m < n ≪ L L X d n L = nd d d > 10 → X 16/20
17. 別に有限個の点集合を取ってこなくても最小点を求められる。 を一つサンプリングしておけば、 ​ ​ f (n) α ​ (x) TS ≃ ϕ(x) θ ⊤ ↓ ≃ −ϕ(x) θ ⊤ という形で獲得関数が（近似ではあるが）解析的に求まるので、後は勾配法で最小化すれば良い。 なので 獲得関数の勾配 は、 ​ ​ ​ ∂x(k) ∂α ​ TS = − ​ ​ θ ​ ​ cos 2π ​ x ω ​ + b ​ ∂x(k) ∂ j=1 ∑ m j ​ m 2σ ​ 0 2 {( l=1 ∑ d (l) j (l)) j } = 2π ​ ​ θ ​ ω ​ sin {(2πx ω) + b ​ } ​ m 2σ ​ 0 2 j=1 ∑ m j j (k) ⊤ j 勾配法の反復回数を とすると計算量は であり、 なら大幅に高速化 3.5.5 疎スペクトル近似を用いた効率的な計算方法 θ ϕ(x) = ​ cos {(2πx ω ​ ) + b ​ } ​ m 2σ ​ 0 2 ⊤ j j ​ ∂x(k) ∂α ​ TS T O(nm + 2 T m) T ≪ L 17/20
18. 実装：近似した獲得関数のプロット 18/20
19. 多腕バンディットだと実用的に高い性能を発揮するアルゴリズムとして知られている ベイズ最適化だとそれほど使われていない 低次元では比較的上手くいく 高次元だと性能が悪化 トンプソン抽出獲得関数が探索と活用の探索に寄りすぎているから、という指摘 [B. Shahriari, et al. (2016)] 他のアプローチとして、 からサンプリングした値を初期値として、尤度最大化によってガウス過程 を最適化する手法も提案されている 過学習のリスクはあるが、目的関数に合わせて が変化するため、比較的小さな で目的関数 を近似できることがあることが報告されている [A. Yang, et al. (2020)] 3.5.6 トンプソン抽出獲得関数の性質 p ​ s ω ​ j m (1) (1) : 多分テキストの参考文献は間違ってそう 19/20
20. 1. 素朴な方法 計算が単純で実装が楽 ガウス過程の学習の計算量が 事後分布のサンプリングの計算量が 勾配法が使えない 2. 疎スペクトル近似に基づく方法 ガウス過程の学習の計算量が これが一番重要 事後分布のサンプリングの計算量が 勾配法が使える 平行移動不変カーネルにしか使えない （再掲）獲得関数の実装方法 O(n ) 3 O(L ) 3 O(nm ), m < 2 n (1) ← (3) O(mL) (2) (1): が より小さすぎると外層領域での振る舞いが悪化する variance starvation が起こる (2): ほとんどの場合で不変カーネルが使われるのでそこまで問題ではないかも (3): 近年のベイズ最適化の手法でよく使われている (Random Features 等) m n 20/20