機械学習を用いた宇宙大規模構造の理論予言構築小林洋祐 (東大Kavli IPMU D2)

第6回High Performance Computing Physics勉強会
＠京都大学基礎物理学研究所, 2020. 1. 31

今回のトーク
このトークは京大基研の西道さん（2018年までIPMU所属）主導のプロジェクト”DarkQuest”の成果に基づいています。
[1] Nishimichi et al. 2019, “Dark Quest. I. Fast and Accurate Emulation of Halo Clustering Statistics
and Its Application to Galaxy Clustering” (arXiv:1811.09504)

使われている計算・機械学習手法の主要References
ラテン超方格法[2] Ba et al. 2015, “Optimal Sliced Latin Hypercube Design”, https://doi.org/10.1080/00401706.2014.957867

EMPCA[3] Bailey et al. 2012, “Principal Component Analysis with Noisy and/or Missing Data” (arXiv:1208.41122)

ガウス過程回帰[4] Rasmussen & Williams, “Gaussian Process and Machine Learning”, the MIT press, 2006

“

1. 宇宙大規模構造を用いた宇宙論
2. エミュレータの概念とその目的
3. エミュレータ構築の手法

4. まとめ宇宙の大規模構造
宇宙の ~1 Mpc (~ 100万光年) 以上のスケールに広がる物質や銀河の分布
宇宙初期に生成された量子ゆらぎが成長してできたと考えられている
大規模構造の形成は宇宙の大スケー
ルの枠組みを決める宇宙論パラメー
タに依存しており、大規模構造の情報から宇宙論を制限できる

Credit: Springel et al. (Virgo Consortium)
Simulation code: Gadget-2

大規模構造の情報①: 重力レンズ効果
2次元天球面上で見える銀河の像は、その前景にある重力源（主にダークマター）により歪められている

この歪み具合から、前景のダークマターハロー周りの物質の集まり具合が統計的にわかる

δg(x) ≡ng(x)
n̄g−1

ρm(x)
δm(x) ≡−1
ρ̄m

δ
(x)δ
(x)
≡ξ
(x−x)
gm 12
⟨ g 1 m 2⟩にセンシティブ

大規模構造の情報②: 銀河クラスタリング
クラスタリング＝「集まり具合」
大規模構造上に広がる銀河の数密度場の相関δg(x1)

δg(x2)

δg(x) ≡ng(x)
n̄g−1

δ
(x)δ
(x)
≡
ξ
(x−x)
g 12
⟨ g 1 g 2⟩宇宙論ではとくに2点相関が重要で、これを銀河の大域的分布から測定して解析する

Redshift z

Redshift z

ダークマターハローと銀河
ダークマター分布

銀河分布

比較的単純な物理（重力）

複雑な物理

大規模シミュレーションも大規模シミュレーションはできる

かなり厳しいWechsler & Tinker 2018

ダークマターハローと銀河ダークマター分布

銀河分布
本当に予言したいのはこっちだが比較的単純な物理（重力）

複雑な物理

大規模シミュレーションも大規模シミュレーションはできる

かなり厳しいWechsler & Tinker 2018

ダークマターハローと銀河ダークマター分布

銀河分布

比較的単純な物理（重力）

複雑な物理

大規模シミュレーションも大規模シミュレーションはできる

かなり厳しいここではこちらにフォーカス

Wechsler & Tinker 2018

宇宙大規模構造を用いた観測的宇宙論
宇宙大規模構造の観測（重力レンズ＆銀河クラスタリング）からどうやって宇宙論につなげるか？
基本的には、理論予言と観測データからベイズ推定をする

P(model | data) ∝ P(data | model)P(model)
事後分布

尤度関数

事前分布

Problems
ここでは専らパラメータ推定を想定

P(θ | data) ∝ P(data | θ)P(θ)→ MCMCサンプリングでパラメータ事後分布を推定する問題

これが目標

“

1. 宇宙大規模構造を用いた宇宙論
2. エミュレータの概念とその目的
3. エミュレータ構築の手法

4. まとめ宇宙大規模構造の理論予言
重力レンズ

〜

銀河クラスタリング 〜ξhm(r) = ⟨δh(x)δm(x + r)⟩
ξhh(r) = ⟨δh(x)δh(x + r)⟩

More+ 2015
重力レンズ
Problems

銀河クラスタリング宇宙大規模構造の理論予言
解析的なアプローチは難しい -> ゆらぎの摂動展開が小さいスケールで破綻するから2点相関などの宇宙論統計量は空間スケール r [Mpc] の関数なので、摂動論が有効な大スケールだけでなく、小スケールの情報も解析に取り入れたい
このとき、
1. ゆらぎの非線形効果が効く広範なスケールにわたって精密に
Problems
2. MCMCサンプリングによる宇宙論パラメータ推定ができるくらい高速に重力レンズやクラスタリングの統計量を予言したい！
→ 宇宙論的N体シミュレーションなら、ダークマター分布を数値的に生成可能で、解析的な予言の難しい非線形効果まで正確に取り入れられる

N体シミュレーションに基づく“エミュレータ”構築“エミュレータ”とは？
工学分野では「代理モデル (Surrogate model)」とも言う→計算量のかかるシミュレーションをやらずに、好きなパラメータを入力してシミュレーションと同じように物理量の値を出力できるようにしたもの

地球統計学では「クリギング (Kriging)」という名前
Problems→特にエミュレータ構築の手法であるガウス過程回帰を指してクリギングと言うらしい

これを宇宙論の問題に適用した例がこのトークの内容“いろいろな宇宙論パラメータで大規模構造形成シミュレーションを行い、宇宙論に依存するダークマターハローの統計量を測ってそのエミュレータを作る”

Dark Quest (Nishimichi+ 2019)
エミュレータ（呼称：”Dark Emulator”）構築のためのN体シミュレーション群

2

Ωch2
ΩΛΩbh

101個の宇宙論パラメータ点でダークマター分布の重力進化を解き、宇宙大規模構造を数値的に生成する

現状、kmax ∼ 0.15 hMpc−1

までのデータでMCMC解析できたダークマター分布からハローを見つけ、カタログにする

ハローカタログとダークマター分布のデータから

wde

ns

Cosmological parameter spaceAs

興味ある統計量（2点相関関数など）を測定得られたデータを宇宙論パラメータ空間上で補間、シミュレーションをやっていない宇宙論でも欲しい値を出力できるようにするDark Emulator

ξhm(r) = ⟨δh(x)δm(x + r)⟩

ξhh(r) = ⟨δh(x)δh(x + r)⟩

Dark Emulatorの学習精度

ξhm(r) = ⟨δh(x)δm(x + r)⟩r

ξhh(r) = ⟨δh(x)δh(x + r)⟩r

“

1. 宇宙大規模構造を用いた宇宙論
2. エミュレータの概念とその目的
3. エミュレータ構築の手法

4. まとめシミュレーション実験計画法
機械学習にかけるシミュレーションデータを生成するとき、宇宙論パラメータの値をどうとるか？
- 単純に「パラメータ1つにつきm個」とすると、6パラメータではm^6個の組み合わせ- あまりに膨大な数のシミュレーションはできないが、考えているパラメータ領域は満遍なくサンプルしたい

ラテン超方格法 (Latin hypercube design)「パラメータ1つにつきm個」かつ「6パラメータでもm個」
特色Space-filling property : パラメータ空間から一様にサンプルできる
Good projection property：低次元空間へ射影したときにムダが少ないラテン超方格法
p次元パラメータ空間でn回のシミュレーションをするとして、まず、各々のパラメータで、n個全て異なる値にとる → good projection propertyn = 10, p = 2の場合

Crombecq+ 2009

学習・検証データを分けるラテン超方格法
機械学習では、学習データと検証データを分ける必要がある。このとき(1) 学習データも検証データも、それぞれ考えているパラメータ領域からの一様なサンプリングになっていて欲しい(2) 加えて、学習データと検証データを合わせたデータ集合がまたパラメータ領域からの一様なサンプリングになっていて欲しい→ “スライス’’という概念
1つのラテン超方格が複数の互いに素なラテン超方格を含むようなものを構成できる

“Sliced” Latin Hypercube Design
ラテン超方格は n × p 行列で表せる。下図はパラメータが2つ (p = 2) の場合3つの小さいラテン超方格（スライス）から1つの大きなラテン超方格を作る操作

slice 1slice 2slice 3

一番良いSliced Latin Hypercube Designはどれか？
12シミュレーション、2パラメータ、3スライスの場合、下のどれもありうるが。。。

“Maximin” criteria
Sliced Latin hypercube構成のoptimalityをどんな尺度で決めるか？
→ ナイーブには、データ点の最小距離（のようなもの）が最大（＝Maximin）になるように、つまりを最小化するようにパラメータ点を入れ替えていく
n: データ数

d: データi,j間の距離加えて、各スライスの中でも一様なサンプリングになっていてほしいので、を最小化すべき目的関数にする
t: スライスの数

Nishimichi+ 2019

Dark Questでは、宇宙論パラメー
タ20点ずつが1つのスライス全部で5つのスライスがある
Dark Quest Simulation- Gadget-2 2048^3 particlesガウス過程回帰を用いたエミュレータ構築
ダークマターハローの各種の統計量に対して、
Input: 宇宙論パラメータ
Output: その統計量のデータベクトルという対応を学習させる（教師あり学習）
問題点
ガウス過程回帰では、普通は1つのガウス過程で1つのスカラー値関数を回帰する
シミュレーションから測ったデータベクトルはかなり多くの成分を持っているので、その成分一つ一つに対してガウス過程を構成していくのは現実的でない→データベクトルの成分数（次元）を落として、その低次元データを回帰する

前処理：主成分分析(PCA)
高次元なデータベクトルを線形変換によって低次元ベクトルで要約する
このとき、低次元ベクトルから元のデータベクトルを最も忠実に復元するような変換を選ぶことができる
X: (データの次元 × データの個数)
P: (データの次元 × 主成分の個数)
C: (主成分の個数 × データの個数)→元のベクトルと復元されたベクトルとの間でを最小にするPとCをとる

EMPCA
データにノイズや欠損がある場合、それを考慮する重みをつけて主成分分析をしたい
PとCを交互に最適化していくことで、データのノイズ共分散を含んだ χを最小化できる
参考：EM (Expectation Minimization)ア
ルゴリズムの一例
Cの最適化: E-step
Pの最適化: M-step2

ガウス過程回帰①
多次元パラメータ空間上の関数の非線形回帰

xn ↦ f(xn)

入力ベクトル

出力値

ガウス過程とは・・・- 確率変数の集合であって、その任意の有限個が多変数ガウス分布に従うというもの

入力空間の各点 x に対応する出力値 f(x) を1つ1つ確率変数であると思って、その統計的性質が平均:-> 簡単のため、通常は線形変換によってゼロにとる共分散:であるようなガウス過程になっていると「見なす」

ガウス過程
グラフィカルモデルによるガウス過程 (Rasmussen & Williams 2006)
ノイズが乗った出力（e.g. シミュレーションデータ）

ノイズなしの出力値（予言したい量）

入力ベクトル（e.g. 宇宙論パラメータ）

ガウス過程回帰② 事前分布
実際にデータでガウス過程回帰をするには：
まず入力xに対して出力値 f(x) が対応するデータが n 点あるとして、この出力の集合

f := {f(xi) | i = 1,...,n} がある共分散関数 k(x, x′) をもつガウス過程であるとみなす。
このとき、n × n 共分散行列が得られる

この時点ではまだきちんとした予言はできない→事前分布

事前分布からサンプルされた1変数関数

ガウス過程回帰③ 事後分布
ガウス過程が回帰に有用なのは、学習データをもとに事後分布を作れるから学習データ

事前分布

未知のデータ

事後分布: 学習データからの条件付き確率として得られる

事後分布からサンプルされた1変数関数

ガウス過程回帰④ ノイズありの場合
データのもつノイズを考慮すると、データの出力値は、実際にはだと思える
このとき、n × n 共分散行列は学習データ

事前分布

新しい予言値

事後分布

事前分布の共分散

この分だけ共分散が狭まる

ガウス過程回帰⑤ 学習
ガウス過程回帰のとき、共分散関数（あるいは、カーネルとも）には具体的な関数形を仮定する
e.g. ガウス関数:

ハイパーパラメータ
ここで導入したパラメータ θ0 や Θ を通じて共分散関数の形に自由度を持たせられる

ガウス過程回帰⑤ 学習
上のガウス型の共分散関数で、ハイパーパラメータℓを動かす

ガウス過程回帰⑤ 学習
ガウス過程回帰のとき、共分散関数（あるいは、カーネルとも）には具体的な関数形を仮定する
e.g. ガウス関数:

ハイパーパラメータ
ここで導入したパラメータ θ0 や Θ を通じて共分散関数の形に自由度を持たせられる

ハイパーパラメータを最適化すること＝ガウス過程による「学習」
具体的には、marginal likelihoodを最大化するハイパーパラメータを見つける

Dark Emulator

ξhm(r) = ⟨δh(x)δm(x + r)⟩

ξhh(r) = ⟨δh(x)δh(x + r)⟩

まとめ
宇宙大規模構造の統計量を精密かつ高速に予言するコード（Dark Emulator）を、N体シミュレーション群と機械学習手法を用いて構築
ハロー数密度場と物質密度場との相関、ハロー数密度場の2点相関などを、（解析的模型を使ったモデリングが難しい）小さいスケールまで精密に予言可能
具体的なテクニック：
1. ラテン超方格を使って効率的なシミュレーション計画
2. 主成分分析を使ってシミュレーションデータの次元削減
3. ガウス過程回帰を使って宇宙論パラメータ空間上で学習データを補間