特集「諸科学における統計数理モデリングの拡がりII」［原著論文］

統計数理（2023）
第 71 巻 第 2 号 159–187©2023 統計数理研究所

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出
竹内 努1,2 ・河野 海1 ・クレ スチェータ1,3,4 ・西澤 淳5,6,7 ・村上 広耶1 ・馬 海霞1 ・本武 陽一8（受付 2022 年 9 月 21 日；改訂 2023 年 3 月 16 日；採択 3 月 20 日）

要旨

宇宙空間において，銀河は銀河団，銀河群，フィラメント（1 次元的な銀河の連なり），ヴォイド（銀河の存在しない領域）といった構造を形作り不均一に分布している．これを宇宙の大規模構造と呼ぶ．宇宙論では通常の物質をバリオンと呼び，我々が直接観測するのはこのバリオンが形作る大規模構造である．
大規模構造を形成する主要因は重力不安定だが，これとは別種の構造の源となる機構が存在する．これがバリオン- 光子脱結合の時期に生じるバリオンの音波である．これはバリオン
音響振動（BAO）として知られ，銀河の空間分布に刻印されている．本研究では，銀河分布のBAO 信号を検証するために，位相的データ解析（topological data analysis: TDA）の方法で銀河の空間分布を解析した．TDA では，点群データを「穴」の集合として解釈し，その幾何学的構
造からパーシステントホモロジー群を構成することで，1 つ 1 つの穴構造の特徴を抽出する手法である．従って，それぞれの穴の情報を保持するパーシステントホモロジー群と穴の位置情報などを組み合わせることで，TDA によって点群データから定義される穴集合のそれぞれの穴のサイズ，位置などの特徴が定量評価できる．
我々はまず，バリオンの物理あり，なしの場合のシミュレーションデータセットを解析し，
TDA のパフォーマンスを調べた．この結果，TDA は，宇宙の大規模構造の中から BAO 信号を実際に検出できることがわかった．バリオンありのシミュレーションは BAO からの有意な信号を示したのに対し，バリオンなしデータはこの信号が検出されなかった．この結果をもと（SDSS DR14）の eBOSS（extendedに，次にスローンデジタルスカイサーベイデータリリース 14から z < 1.0 のクエーサー（活動銀河核の一種）を抽出
Baryon Oscillation Spectroscopic Survey）
してサンプルを構成し，TDA を適用した．この解析で，我々は r ∼ 150 [Mpc] のスケールで特徴的な穴（中空のシェル）を発見した．これは，銀河/クエーサー分布に刻印された BAO の信1

名古屋大学 素粒子宇宙物理学専攻：〒 464–8602 愛知県名古屋市千種区不老町2

統計数理研究所 客員：〒 190–8562 東京都立川市緑町 10–3

3

国立天文台 科学研究部：〒 181–8588 東京都三鷹市大沢 2–21–1

4

日本学術振興会特別研究員（PD）5岐阜聖徳学園大学 DX推進センター：〒 501–6194 岐阜県岐阜市柳津高桑西 1–1

6名古屋大学 高等研究院：〒 464–8601 愛知県名古屋市千種区不老町

7名古屋大学 素粒子宇宙起源研究所：〒 464–8602 愛知県名古屋市千種区不老町8

統計数理研究所（現 一橋大学 ソーシャル・データサイエンス研究科：〒 186–8601 東京都国立市中 2–1）

統計数理

160

第 71 巻 第 2 号 2023

号に正確に対応している．この解析は 2000 個のクエーサーからなるサブサンプルで実行した．
これは TDA が，たとえ非常に疎なサンプリングであっても中空の穴のような位相的構造を検出するのに有効であることを示している．
キーワード：位相的データ解析，パーシステントホモロジー，銀河分布，宇宙の大規
模構造，バリオン音響振動．
1.

はじめに

1.1 宇宙の物質分布の特徴づけと記述方法
各宇宙年齢における物質および銀河の分布には，宇宙の構造の形成と進化に関する基本的な情報が含まれている（e.g., Peebles, 1980; Efstathiou and Silk, 1983; Bernardeau et al., 2002）．
まず強調しておかなくてはならない事実として，宇宙論が理論的に予言するのは物質分布ある
いは銀河分布の統計的性質であり，時空間のどの位置にどの銀河が誕生するかなど具体的な天体についての予言力はない．よって，宇宙論の理論と観測の比較は統計的記述を通じて行うしかない．この理由から，物質分布のゆらぎの統計的特性を特徴づけるために，膨大な数の洗練された方法が提案されてきた（e.g., Peebles, 1980; Martinez and Saar, 2001; Bernardeau et al.,2002）．その中でも，n 点相関関数は，観測される銀河分布（e.g., Peebles, 1980; Martinez andの解析で最も一般的でよく研究されている方法である．宇
Saar, 2001; Bernardeau et al., 2002）
宙空間のある位置 x における銀河の数密度 n(x) を測定し，その統計平均 n と比較する．宇宙（アンサンブル：は我々の住んでいるもの 1 つしか観測できないが，仮想的に多数の宇宙の集合
ensemble）
を想定し，ある量 Q の多数の宇宙での平均を考えることを統計平均あるいはアンサ
ンブル平均と呼び，記号 Q で表す．これに対し，空間内で平均した値を体積平均とよび，記号 Q̄ で表す．宇宙論では，この 2 種類の平均が一致するという仮定（宇宙論的エルゴード仮説）に基づいて議論を展開する．
この準備のもと，d3 x1 , d3 x2 という 2 つの微小体積を考える．双方の微小体積中に銀河が見つかる同時確率を P(x1 , x2 )d3 x1 d3 x2 とする．銀河分布が一様，即ち無相関（クラスタリングがない）
場合，微小体積内に銀河が見つかる確率は
(1.1)

P(
x)d3 x = nd3 x ,なので，体積 d3 x1 内の銀河と体積 d3 x2 内に銀河が見つかる同時確率は
(1.2)

P(
x1 , x2 )d3 x1 d3 x2 = n2 d3 x1 d3 x2となる．銀河がクラスタリングしている場合，同時確率は
(1.3)

P(
x1 , x2 )d3 x1 d3 x2 = n(x1 )n(x2 )d3 x1 d3 x2 ≡ n2 [1 + ξ(x1 , x2 )] d3 x1 d3 x2となり，一様な場合に比べて n2 ξ(x1 , x2 )d3 x1 d3 x2 だけ高くなる．この ξ(x1 , x2 ) を銀河分布の 2 点相関関数と呼ぶ．同様にして，一般に n 点相関関数は n 個の微小体積内に銀河を見出す
同時確率を用いて定義される．
ここまでは銀河分布を点群として定式化したが，銀河分布を議論する際には数密度の分布を平滑化し，密度場として扱うことも多い．点群と密度場の数学的扱いは異なるが，宇宙物理ではこの違いは重視しないで同様に扱うことが多く，本稿でもそれに倣う．宇宙論的定式化のため，銀河分布の密度場 ρ(x) の密度ゆらぎ
(1.4)

δ(
x) ≡

ρ(x) − ρ̄ρ̄

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

161を考えるのが便利である．ここで ρ̄(x) は平均密度を表す．密度 ρ(x) は点群としての銀河の数
密度，あるいは平滑化した銀河分布の質量密度のどちらを表現する場合もありうる．密度ゆらぎ δ(x) を用いると，相関関数 ξ(x1 , x2 ) は
(1.5)

ξ(
x1 , x2 ) =x2 ) − ρ̄][ρ(x1 ) − ρ̄] [ρ(= δ(x1 )δ(
x2 )
ρ̄2

と表せる．なお，宇宙論では宇宙は一様等方であると仮定するので，対称性より
(1.6)

ξ(
x1 , x2 ) = ξ (|x1 − x2 |) ≡ ξ(r)として扱うことが一般的である．式（1.5）の形で表すと，相関関数が密度場の 2 次のキュミュラントになっていることが分かる．一般的に，ある空間上のベクトル x = (x1 , . . . , xn ) を添字ととよぶ．数学的には，ある確率場は無限次までのする確率変数 X(x) を確率場（random ﬁeld）
モーメントあるいはキュミュラントで一意に定まることから，キュミュラントと等価な n 点相関関数によって確率場が決まる．これは確率場の特性関数のキュミュラント展開定理を用いて示せる（e.g., Daley and Vere-Jones, 2003）．即ち，相関関数を測定することによって原理的には宇宙の密度場の統計的性質を特定できる．
しかし実際には，高次の n 点相関の推定は極めて難しく非現実的であるため，高次相関の情報を扱う別のアプローチが検討されている．銀河分布の 3 次元構造の発見（e.g., de Lapparent（トポロジー）を用いて大規模構造を定量化
et al., 1986）により，平滑化した密度場の幾何学
し，背後にある物理過程を特定しようというアイディアが提示された（e.g., Shandarin and．早くから用いられてき
Zeldovich, 1983; Gott et al., 1986, 1987, 1989; Boerner and Mo, 1989）
によるトポロジーの特徴付けである（Gott et al., 1986）．銀
たのがジーナス（種数，genus: G）（ν は定数）を閾値とすれば，νσ に対する上側レ
河分布の密度ゆらぎ δ(x) の標準偏差 σ の ν 倍（ここでは Eν と表記）とよ
ベル集合が各 ν について定義できる．これをエクスカーション集合
び，ρ(x) = νσ となる等密度面がその表面となる．密度場を考えているので，Eν は 3 次元の集合である．銀河分布のジーナスは，平滑化した銀河密度場のオイラー標数の −1/2，即ち

1
G=−
(1.7)
KdA ,
2
(K: 等密度面のガウス曲率，dA: 等密度面上の面積要素）と定義される．なお，ここでの定義
は宇宙論で用いられるもので，数学における通常の定義と 1 だけずれていることに注意されたい（松原, 1996）．エクスカーション集合 Eν の単位体積あたりのジーナスを ν の関数として表した量 G(ν) をジーナスカーブと呼ぶ．ジーナスカーブを銀河分布の集約統計量として用い，宇宙論的研究を行う方法がジーナス統計である．このアプローチは，大規模構造の重力による成長を密度分布の初期値からの微小なずれ（摂動）として扱うことで，大規模構造進化を記述する理論体系，いわゆる構造形成の摂動論的アプローチに関連して開発されてきた（e.g.,
Matsubara, 1994, 1996, 2003; Matsubara and Suto, 1996; Matsubara et al., 2022; Matsubara and．摂動論を用いると，ジーナスカーブ
Yokoyama, 1996; Pogosyan et al., 2009; Gay et al., 2012）
は大規模構造のスケールではモーメントの関数として解析的に書けるので（e.g., Matsubara,，これによって位相幾
1994, 1996; Matsubara and Suto, 1996; Matsubara and Yokoyama, 1996）
何的情報を定量的に表現できる．特にジーナスカーブはガウス場と他の確率場を区別できるため，大規模構造を解析するための一般的な道具として定着している．
（Minkowski
ジーナス統計の方法を含む，より一般的な体系がいわゆる “ミンコフスキー汎関数” である（e.g., Hikage et al., 2003; Matsubara, 2003; Sullivan et al., 2019）．ミンコ
functionals）
フスキー汎関数の大規模構造への応用でも式（1.4）で定義した平滑化した銀河密度場の密度ゆ

162

統計数理

第 71 巻 第 2 号 2023

らぎを考える．この場合ミンコフスキー汎関数 Mk は k = 0 については

1
(1.8)
M0 =
dV ,
V Eν
k = 1, 2, 3 については閾値 ν の等密度面上の曲率の積分によって定義されるので，原理的には
銀河密度場はミンコフスキー汎関数によって完全に特定できる（Schmalzing and Buchert, 1997;（これ以降，本論文を通して k = 0, 1, 2, 3 は次元の意味でのみ用
Schmalzing and Gorski, 1998）いる）．実際，この方法は宇宙論的解析に広く応用されている（密度場の解析: e.g., Mecke et al.,
1994; Kerscher et al., 1997; Beisbart et al., 2002; Hikage et al., 2003; Kerscher and Tikhonov,
2010; Wiegand et al., 2014; Wiegand and Eisenstein, 2017; 赤方偏移歪み: e.g., Appleby et al.,2018; 弱重力レンズ宇宙論: e.g., Spurio Mancini et al., 2018; Vicinanza et al., 2019; Parroni et al.,2020; 宇宙再電離: e.g., Yoshiura et al., 2017; Bag et al., 2018, 2019; Chen et al., 2019; 宇宙論/
重力理論の検証: e.g., Ling et al., 2015; Junaid and Pogosyan, 2015; Shiraishi et al., 2016; Fang
（CMB）
: e.g., Matsubara, 2010; Planck Collaboration
et al., 2017; および宇宙マイクロ波背景輻射
et al., 2016; Santos et al., 2016; Chingangbam et al., 2017; Ganesan and Chingangbam, 2017;．
Joby et al., 2019）
1.2 ミンコフスキー汎関数から位相的データ解析へ
このように，銀河密度場は密度場のトポロジーによって定量的評価ができる．上で述べた
ジーナス統計は，点群データの性質を特徴付ける．ジーナスのように，点群データ全体をその穴構造についての一つの特徴量で表現するのではなく，点群データに含まれる個別の穴構造についてのスケールのような特徴をそれぞれ保持できる表現がパーシステントホモロジー群で．
ある（TDA: e.g., Edelsbrunner et al., 2002; Zomorodian and Carlsson, 2004; Wasserman, 2018）
以下ではパーシステントホモロジー群に基づくデータ解析法である TDA についての簡易な説
明を行う．
データが持つ一つ一つの穴構造を評価できる TDA を用いることで，宇宙論におけるデータ
解析をより発展させることができる．本研究で対象となる点群データの穴構造は，銀河や任意の天体などで構成される N 体系の各点に配置された半径 r の n 次元球の集合によって形成される．宇宙論分野黎明期のこのアプローチによる研究は，基本的にはランダムパーコレーションの応用であった（Gott et al., 1986）．すなわち，各銀河の周りに球を置いて半径 r を増加させ，半径 r の関数として考えている体積内の連結成分の幾何学的構造を調べる方法である．これは TDA で用いられる考え方そのものである．
萌芽的な意味での TDA の概念の誕生は 90 年代に遡るが，実質的な誕生は 2000 年代に入ってからである（see e.g., Edelsbrunner et al., 2002; Edelsbrunner and Harer, 2010）．TDA では，
よくデータセットをパーシステントホモロジー群として表現することが行われる．ランダム
パーコレーションと同様，パーシステントホモロジー群では点群データの各点の周りに半径 rの（中空でない）球を考え，さまざまな次元の穴構造の生成と消滅を扱う．この過程を図 1 で概略的に説明する．2 次元平面内に 3 つの点からなる集合を考える．半径 r が小さい場合，これ．次に r がある値より大きくなると，は点の周りに 3 つの円盤が置かれた構造である（図 1: 左）
円盤が連結してサイクル構造が形成される（図 1: 中央），これが穴の誕生とみなされる．さら（図 1: 右）．に r が大きくなると，サイクル構造の内部が円盤で埋め尽くされて，穴は消滅する
ある穴が発生（誕生）する r を x 座標，消滅（死亡）する r を y 座標として，データに存在する全ての穴構造の生成と消滅をそれぞれプロットして構成されるのがパーシステンス図（persistenceと呼ばれるパーシステントホモロジー群の中心的な表現方法である．n 次元空間中のdiagram）

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

163

図 1．パーシステントホモロジー群の概念図（Kono et al., 2020）．位相的情報は，離散的な「穴」によって特徴づけられる．
データ点からの半径 r の n 次元の球で構成された

点群データに対してパーシステンス図は，0 次元の穴（連結成分），1 次元の穴（リング構造），2次元の穴（シェル構造）…といった複数の次元の穴に対して個別に定義できる．より厳密な定義と方法論は 2 章で示す．
天体物理学や宇宙論における近代的な TDA の応用例はここ数年で徐々に増加している．
Pranav et al.（2017, 2019a）, Feldbrugge et al.（2019）, Wilding et al.（2021）, そして Biagetti
et al.（2021）は一般の確率場への応用を念頭に置き，ガウスおよび非ガウス確率場の確率的ホモロジーに関する詳細な定式化を示した．この方法により，（Pranav et al., 2019b）は CMBに特異な信号を検出した．Cole and Shiu（2018）は，TDA を適用して CMB の非ガウス性を評価する方法を開発した．Xu et al.（2019）は，N 体シミュレーション（Libeskind et al., 2018）
から生成された大きな銀河データ集合に TDA を適用し，シミュレーションデータで 23 個のヴォイドを検出した．これは，TDA のヴォイド検出における有用性を示しており，従来の方法（e.g., CLASSIC: Kitaura and Angulo, 2012）の結果ともよく一致している．さらに彼らは，MassiveNu シミュレーション（Liu et al., 2018）と比較して宇宙論的パラメータを制限し，パラメータ間の縮退を解く方法について検証した．Elbers and van de Weygaert（2019）と Giri andは TDA による宇宙再電離期の電離バブルの構造と進化，およびサイズ分布とMellema（2021）
トポロジーの定量的推定の有効性について議論した．
1.3 バリオン音響振動（BAO）
宇宙の構造形成は質量において支配的なダークマターが駆動する．しかし，電磁波で観測できる宇宙，即ち目に見える宇宙の姿を形作るのは電磁波と相互作用する通常の物質のなす構造である．宇宙論では通常の物質をバリオンとよぶ慣習があるので，本論文でも以降通常物
質をバリオンと表現する．バリオンは，電磁相互作用（たとえば，輻射加熱/冷却，ガス圧，流体力学的過程）を介して複雑な進化をたどる．このような非自明な現象の典型例が，初期の宇（BAO）である（Peebles and Yu, 1970;
宙のバリオン- 輻射 2 流体で生成されたバリオン音響振動．ここでは，BAO と関連する観測量について簡単に説明する．
Sunyaev and Zeldovich, 1970）
初期宇宙のバリオン- 光子プラズマにおける点源による初期摂動を考える．プラズマでは，バ
リオンと光子がトムソン散乱（Thomson scattering）を解して結合し，単一の流体としてふるまう．光子は光速で動き，また非常に多数存在するため，バリオン- 輻射流体は密度に対してとてつもなく大きな圧力を持つ．圧力は周囲と均一になろうとし，その結果外側に向かって伝播す
る球形の音波を生じる．この時期の音速は，次式によって評価される．

164

(1.9)

統計数理

第 71 巻 第 2 号 2023

cs =  

c

3ρB
3 1+
4ρR

 .

ここで c は光速であり，ρ はエネルギー密度を，下付き文字 B と R はそれぞれバリオンと輻射は，バリオンと輻射の脱結合前の音速が，初期の光速の ∼ 57 % に及の成分を表する．式（1.9）
ぶことを示している．音響的地平線（sound horizon）の最終半径 rs は，次の式で与えられる．
 tdecrs =(1 + z)cs dt
0 ∞
cs

=
(1.10)dz .
ΩR,0 (1 + z)4 + ΩM,0 (1 + z)3 + ΩΛ,0zdec H0
ここで z は宇宙論的赤方偏移，t は宇宙年齢，H0 はハッブルパラメータ，ΩR,0 は輻射のエネ
ルギー密度パラメータ，ΩM,0 は物質のエネルギー密度パラメータ，ΩΛ,0 は宇宙項の密度パラメータ，下付き添字 dec は光子と物質の脱結合の時刻での値を，0 は現在時刻での値を意味す
る．それぞれの詳細は付録 A.3 節を参照のこと．我々が観測するのは，大規模構造に刻印された，原始揺らぎから生じた多くの音波の重ね合わせである．BAO は定義上バリオンに生じる現象だが，バリオンと暗黒物質の密度ゆらぎの構造は重力相互作用を通して作用し合うた: e.g., Ma andめ，最終状態では一致することが知られている（バリオンの追いつき（catch-up）．
Bertschinger, 1995; Eisenstein et al., 2007）
BAO のスケールは，共動座標（宇宙膨張による距離の拡大の影響を除いた座標系．付録 A を参照）
で一定である．このことから，BAO シグナルは基本的には銀河の 2 点相関関数で検出することができる．しかし，BAO スケールは銀河の典型的なスケールと比較して非常に大きい
ため，2 点相関に現れる BAO シグナルはかなり弱い．したがって，BAO シグナルを測定するには，広大な観測体積において高密度でサンプリングされた大規模銀河探査が必要とされる．
Sloan Digital Sky Survey (SDSS) は，全天の 3 分の 1 をカバーする最も広い可視測光・分光探査である． SDSS により，2 点相関関数で約 150 Mpc の BAO シグナルが検出された（Eisenstein．ここで単位 [pc] はパーセクを意味し，1 pc = 3.08 × 1018 cm である．Mpc, Gpc
et al., 2005）
はそれぞれメガパーセク，ギガパーセクと読み，106 pc, 109 pc を意味する． とはいえ，伝統的方法による BAO 解析は SDSS のような大規模広域探査データでのみようやく可能になる．
よって，データへの要求がより少なくてすむ BAO シグナルの検出法および定量化の開拓は，宇宙論モデルの制限において極めて重要な課題である．本研究ではまず TDA を導入し，銀河の空間分布データへの応用を紹介する．そして有意性検定の方法についても議論し，TDA がBAO スケールの評価に有用であることを示す．さらに，TDA の長所である逆解析についても議論する．これは，銀河分布上に刻まれたバリオン音響振動に対応するトポロジー構造の視覚化を可能にする初めての方法である．
本論文は次のように構成されている．2 章では，TDA の概念と方法，逆解析を紹介する．次に，3 章ではバリオン物理ありとなしの場合のシミュレーションデータを用い，BAO 研究の道具としての TDA の有効性を検証する．4 章では SDSS DR14（Pâris et al., 2018）の extended
Baryon Oscillation Spectroscopic Survey（eBOSS; Dawson et al., 2016）Large Scale Structure(LSS) catalog（Ata et al., 2018）に TDA を適用した結果を示す．5 章では，本論文で導入した方法論の結果と将来展望についてのまとめと議論を述べる．本稿の物理的内容は Kono et al.
（2020）の内容と議論を基礎にしている．
本 論 文 で は Planck 衛 星 に よ る CMB 観 測 で 推 定 さ れ た 宇 宙 論 パ ラ メ ー タ h =H0 /100 [km s−1 Mpc−1 ] = 0.6766, ΩΛ,0 = 0.6889, ΩM,0 = 0.3111, ln(1010 As ) = 3.047 を採位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

165

用する．この物理的意味については付録 A.3 節を参照していただきたい．
2.

位相的データ解析（TDA）

TDA では，N 次元ユークリッド空間 RN 中の有限個の点群 P の位相幾何学的性質を検証す
る（この論文全体を通して，R，N と N0 は，それぞれ実数，自然数，および非負自然数の集合
を表す）．点群 P が有限集合であるため，伝統的なトポロジーの議論を単純に適用しても P の特徴に関する重要な情報を抽出することはできない．TDA は，スケール（e.g., Edelsbrunner（persistent）性質を調べる
et al., 2002; Boissonnat et al., 2018）を変えてもある程度保持されることにより，点群データが持つ意味のある位相幾何学的特徴を抽出できる．このため，TDA で
基本的な役割を果たすこのような幾何学的モデルを導入する．この解析の主要な道具がパーシ
ステントホモロジー群である．
2.1 パーシステントホモロジー群
2.1.1 点データの「形」
複雑な構造を持つ大規模データの分析では，データの性質を代表する少数の特徴量によってデータを縮約して表現することが重要である．このため，次元縮小（dimensionality reduction）
が中心的役割を果たす．例えば，以下のようなユークリッド空間に設定された点群データを考
える．
(2.1)

P ≡ {xi ∈ RN | i = 1, . . . , m}.

例えばデータ P に周期的な構造がある場合はフーリエ変換が非常にうまく機能し，データ点がクラスタリングしている場合は相関関数による解析が有効である（e.g., Peebles, 1980）．
しかし，このようなアプローチは，非周期的，あるいは明確なクラスター構造を持たない複
雑なデータでは有効ではない．フーリエ変換や相関関数を用いる従来の解析方法では，このような点群の形状をうまく特徴付けることは難しい．代わりに，各データ点を中心とする球を置
き，その球の結合関係としてデータ点の接続関係を導入することでデータを特徴付ける方法が有効である．例えば接続関係を穴構造の数によって特徴づけると，それはドーナツとコーヒーカップを同一視するような同相変換に対して不変なホモロジー群として，点群データを大胆に縮約して表現したことに対応する．一方これでは，穴の大きさのようなスケールの情報が失われてしまう問題がある．そこで前述の，各データ点に置いた球の半径を大きくしていく際の穴構造の生成・消滅を抽出することが有用となる．これによって，穴の形状やサイズといった情報を保持した点群データの縮約表現が得られる．これがホモロジー群を拡張した “パーシステ
ントホモロジー群” である．
2.1.2 点群データにおける
「穴」
の代数表現
前述の通り，点群データの分析ではランダムパーコレーションのように，各点の周りに球を置いて半径を増加させ，半径の関数として考えている体積内の連結成分の穴構造を調べることが有用である．このような点群を元に構成された球の集合は，球という凸閉集合の集合と見なすことができる．凸閉集合の有限個の集まりで被覆された構造があった場合，その構造はトポロジカルな特徴を保存（ホモトピー同値）
しながら，代数操作が可能な単体複体による表現に変換できる．単体複体は図形の構造を表現したものであり，k 次元の三角形に対応する概念である単体の組み合わせとして構成される．
一般には単体複体は位相空間の構造を持たないが，ここではまず具体的なイメージを持ちやすいユークリッド空間を導入した単体複体について説明する．ユークリッド空間内において k

統計数理

166

第 71 巻 第 2 号 2023

単体は以下で定義される．
Definition 1.（ユークリッド空間内での単体（simplex）） Rn 内の k + 1 個の点 p0 , p1 , . . . , pk
で構成される k 個のベクトル p0p1 , p0p2 , . . . , p0pk が一次独立な関係にある場合，その k + 1 個の点によって k 次単体 |p0 , p1 , . . . , pk | は以下のように定義される．
 k
k


|p0 , p1 , . . . , pk | :=λi pi λi ≥ 0(i = 0, 1, . . . , k)， λi = 1
i=0

i=0

具体的には，0 次元の単体は点，1 次元の単体は三角形，2 次元の単体は四面体に対応する（図 2）．この単体をずれなく綺麗に貼り合わせたものが単体複体である．ずれなく貼り合わせるという操作に基づき，ユークリッド空間内における単体複体は以下のように定義される．
Definition 2.（ユークリッド空間内での単体複体（simplicial complex）） Rn 内の有限個の単体の集まり K が以下の条件を満たすとき，K を単体複体と呼ぶ．
K に属する単体 τ の面 σ が K に含まれる．
（1）
2 個の単体 τ, σ ∈ K の共通部分 τ ∩ σ は空集合でない限り，τ の面かつ σ の面である．
（2）（表面）の 4 枚の 3 角形という 1 次元低い単体
ここで単体の面とは，直観的には 4 面体なら境界の集合，線分なら両端の 2 点という 0 次元単体の集合を指す．図 3 にユークリッド空間内での単体複体の例を示す．
前述したように，単体複体そのものは位相空間によらず，以下のような組み合わせ論的構造として定義される．
Definition 3.（単体複体（simplicial complex）） 有限集合 V 上の単体複体 X とは，次の条件{a} ∈ X (∀ a ∈ V )，（1）
σ ∈ X のとき，τ ⊂ σ ならば τ ∈ X ．
（2）
を満たす V の部分集合の族である．
X における単体の最大次元は複体 X の次元と呼ばれ，dim X と書く．また，X の k 単体の集
合を Xk と表記する．
次に，点群データから単体複体を構成する手順を説明する．より具体的には，本研究における「穴」の抽出では，単体複体の一種である Čech 複体を構成する．
P 上の Čech 複体は次のように定義される．
Definition 4.（Čech 複体（Čech complex）） 点群 P ∈ RN を考える．それぞれの xi の周りに半径 r の球Br (xi ) = {x ∈ RN ; x − xi ≤ r} ,(2.2)

を置く．ここで x はユークリッドノルムである．Čech 複体 C(P, r) とは，P 上の単体複体
で，それを構成する単体が
k

(2.3)

{xi0 , . . . , xik } ∈ C(P, r) ⇔Br (xij ) = ∅j=0で定義されるものである．
すなわち点群 P の周りに置いた半径 r の球の関係（重なり）から単体を構成する．半径 r を変

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

167

図 2．単体（simplex）の例．図中の影をつけた部分は中身が詰まっていることを示している．

図 3．単体複体（simplicial complex）の例．図 2 と同様に，影をつけた部分は中空ではなく中身の詰まった図形であることを示す．

図 4．Čech 複体の例．左パネル: xi の周りの半径 r の球の和集合 B, 右パネル: 対応する
Čech 複体 C.

化させたときの Čech 複体の振舞を図 1 に示す．図 1 からもわかる通り，球の和集合と対応する Čech 複体は連続変形によって互いに移りあえる．これは，トポロジーの基本的定理の 1 つから導かれる（e.g., Kozlov, 2008）．つまり，点群の連結性をである脈体定理（Nerve Theorem）
モデル化した半径 r の球の和集合 kj=0 Br (xij ) と単体複体 C(P, r) は，連続変形によって互いに変形できることが保証される（ホモトピー同値）．したがって，単体複体には穴に関する情報が保存される．このようにして，点群データから代数操作が可能な「穴」の表現が実現される．
Čech 複体はデータ P と半径 r によって決まる．半径 r が小さいと Čech 複体は元のデータ
点群 P と同じ離散構造を持つ．一方 r を増大させるとデータ点は球によって連結し始め，r を十分大きく取ると最終的に全体が連結する（図 4）．つまり前述のように，r はデータ点間の連結性を扱うための “解像度” とみなすことができる．半径 r1 < r2 < · · · < rn に対応して Čech複体の増大列 C(P, r1 ) ⊂ C(P, r2 ) ⊂ · · · ⊂ C(P, rn ) が得られる．一般にこのような単体複体の増大列をフィルトレーションと呼ぶ．
実際の応用では，Čech 複体と共通の性質を持ちつつ計算機上での扱いに適した alpha 複体が使用されることが多い．Alpha 複体は，ヴォロノイ分割（Voronoi tessellation）を用いた Čech複体の改良版である．ここでは alpha 複体の厳密な導入はせず，図 5 で概念のみ示すにとどめる．以下の解析では alpha 複体を採用する．

統計数理

168

第 71 巻 第 2 号 2023

図 5．点データ集合の例（左上），その Čech 複体（右上），Delauney 複体（左下），alpha 複体（右下）（Kono et al., 2020）．

2.1.3 ホモロジー
次に，得られた単体複体から「穴」
構造を代数的に算出する方法を説明する．穴構造は，ホモロジー群として表現される．ホモロジー（homology）群は，穴に対応する幾何学的意味を持つ代
数的な対象である．ホモロジー群の代数的計算方法は複数あるが，その中の 1 つについて，以下に概説する．
単体複体 X に含まれる全ての頂点の集合を V = {1, . . . , m} とする．X に含まれる全ての k
単体の集合を Xk とする．Xk に含まれるある k - 単体 σk = |v0 , . . . , vk | の頂点に順序を設定する．つまり，頂点集合に順序 v0 < · · · < vk を設定する．このようにして向きを指定した k - 単体を σk  ≡ v0 , . . . , vk  とする．この向きを指定した k - 単体 σk  を基底として整数 Z を係数とするベクトル空間を導入する．


(2.4)
Ck (X) ≡aσk σk  aσk ∈ Z ,σk ∈Xk

ここで，k > dim X あるいは k < 0 ならば Ck (X) = 0 である．
単体複体の穴を定義するため，サイクル構造を導入する．サイクル構造とは，結合した単体を辿って 1 周回って元に戻ってくる構造である．k 次元の穴は，Xk に含まれる k 次元の単体を組み合わせてできるサイクル構造で，Xk+1 に含まれる k + 1 単体の境界となっていない構造のことを指す．従って，サイクル構造を抽出するには，Xk から境界のない k 単体の集合を見つければよい．このように，穴構造の抽出には単体の境界を算出することが必要となる．向
きが指定された k - 単体 σk  = v0 , . . . , vk  の境界は，以下で定義される境界作用素 ∂k を適用することで得られる．
Definition 5.（単体の境界作用素） 境界作用素 ∂k は
∂k v0 , . . . , vk  ≡

(2.5)
(2.6)

∂k


σk ∈Xkk(−1)i v0 , . . . , vˇi , . . . , vk  ,i=0aσk σk  =σk ∈Xkaσk ∂k σk 位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

169と定義される（ここで vˇi は i 番目の成分を取り除いたことを表す）．
前述したように，単体複体 X に含まれる k 次元の穴は，X に含まれる境界のない k - 単体の集
まりで，X に含まれる k + 1 - 単体の集まりの境界になっていないもののことをいう．サイクルに対応する X に含まれる境界のない k - 単体複体 c の集合は，
Ker ∂k := {c ∈ Ck (X) | ∂k (c) = 0 }(2.7)

で与えられる．また，X に含まれる k + 1 - 単体の境界の集合は，(2.8)

Im ∂k+1 := c ∈ Ck (X) c = ∂k+1 (c ), c ∈ Ck+1 (X)

で与えられる．従って，k 次元のサイクルで k + 1 次元の単体複体の境界でないものは，これらの商 Ker ∂k /Im ∂k+1 をとることで与えられる．
このように，単体複体 X に存在する k 次元の穴構造に対応する数学的意味を持つ，k 次ホモロジー群は以下のように定義される．
Definition 6.（ホモロジー（homology）） 商ベクトル空間(2.9)

Hk (X) = Ker ∂k /Im ∂k+1

は単体複体 X の k 次ホモロジー群である．なお k = 0 のホモロジー群 H0 と Hubble パラメー
ター H0 とを混同しないように注意されたい．
このように，データ P に内在する穴の情報は，その Čech 複体 C(P, r) のホモロジー群を計算することで抽出できる．しかし，この方法では穴のスケール情報を抽出することはできず，ま（persistence）にたある半径 r2 で穴が出現し，半径 r3 で消滅するといった穴のパーシステンス
ついての検証をすることもできない．こういった情報を抽出するためには，パーシステントホ
モロジー群を導入する必要がある．
2.1.4 パーシステントホモロジー群
穴の生成消滅の情報を持つパーシステントホモロジー群も，ホモロジー群と同様に代数的演算によって求められる．そのためにまず，半径 r を増大させ，新たな単体が加わった半径を記（0 ≤ t ≤ n）の単体が追加された半径 rt を用いた場合の単体複録する．また，ちょうど t 番目
体を X t := C(P, rt ) として記録する．このようにして，単体の発生回数 n の数の半径 rt と単体複体 X t のセットが得られる．次に，この rt と X t を用いて，半径 r の増大によって生成・消滅する穴構造を代数的に算出することを考える．X t ⊂ X t+1 であることから，穴の消滅は必ずその穴を境界とする高次の単体複体の誕生により生じる．従って，図形構造の変化を捉える
半径としては，単体の発生時刻を保持すれば十分である．そこで，単体が発生した時刻に，その単体を要素とする基底空間を設定する．例えば，rt で k 単体 σk  が誕生したとすると，それは以下のように表現される．
(2.10)eσk  =

(0,
1

0,
2

...,

σk ,
t

...,

0,
n−1

0)
n

単体複体 X t は，このような基底空間で表現された単体の和として代数的に表現することができる．
ホモロジー群の算出と同様に，上で与えた基底空間で表現された単体複体からサイクルや境界の発生・消滅を算出することで，穴構造の発生・消滅を抽出することを考える．その際に重要な点は，境界に対応する単体・単体複体の発生とサイクルのような図形構造に対応する単体

統計数理

170

第 71 巻 第 2 号 2023

複体の発生する半径が異なるという点である．単体の発生半径の情報しか持たない基底空間の元で，サイクル構造の存在を算出するために，基底に作用して時刻を発展（t を増大）させる作
用素 x を導入する．つまり，さきほどの単体 eσ̂k  に x を作用させると，以下のような変化を与える作用素である．
xeσk  = x(0,10, . . . , σk , . . . ,
0,
2
t
n−1

= (0,
1

0)
n

0, . . . , 0,2t0,
σk , . . . ,
t+1
n−1

0)
n

この作用素を用いた上で，境界の生成・消滅の情報を算出する境界作用素が以下で与えられる．
(2.11)∂k (eσk  ) =ki(xb(σk )−b(σk−1 ) )eσi

k−1 

i=0

, σk ∈ Xk ,i
ここで，b(·) は，その単体が生じた時刻 t を返す関数，σk−1
 は，k 単体 σk = v0 , . . . , vk  の面
v0 , . . . , vˇi , . . . , vk  を表すとした．パーシステントホモロジー群の境界作用素 ∂k に対して，ホモロジー群のときと同様に，サイクルに対応する Ker ∂k と k + 1 次元の単体複体の境界に対応する Im ∂k+1 の商を取ることで，生成時刻と生成から消滅までの区間情報を保持した穴構造が抽出される．このようにして，穴の生成・消滅の情報を保存するパーシステントホモロジー群が算出される．
一般的には，パーシステントホモロジー群 P Hk (X ) は X i のホモロジー群 Hk (X i ) の間の写像列として定義される．

Definition 7.（パーシステントホモロジー（persistent homology）群） k 次のパーシステン
トホモロジー群 P Hk (X ) は(2.12)P Hk (X ) : Hk (X 1 ) −→ · · · −→ Hk (X i ) −→ · · · −→ Hk (X n )として定義される．

ここで，ある単体 σ ∈ X n が誕生する「時刻」
を tσ ≡ min t| σ ∈ X t とした．
2.1.5 パーシステンス図
パーシステントホモロジー群 P Hk (X ) について，bi , di ∈ {1, . . . , n} (bi ≤ di ) および非負整数k が存在し，P Hk (X ) がホモロジー群の構造を保存しながら以下のように一意に分解できることが知られている．
(2.13)P Hk (X ) kI[bi , di ] .

i=1

ここで，区間 I[b, d] は k 次元の穴が Xb で誕生し，Xd で死滅することを表現し， はホモト
ピー同値を表すものとした．
式（2.13）で表される直和分解より，パーシステントホモロジー群は誕生と死滅の対 (bi , di )（パーシステント対）によって一意に定まるとわかる．パーシステント対の集合をプロットしたが定義できる．
パーシステンス図（persistence diagram）
Definition 8.（パーシステンス図（persistence diagram））(2.14)

Dk (X ) = {(bi , di ) | i = 1, . . . , k }を X の k 次パーシステンス図と呼ぶ．
においては bi < di であるため，パーシステンス図のすべてのパーシステント対（誕生 - 死滅対）

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

171

図 6．パーシステンス図の模式図．対角線近くに位置するデータ点は位相的にはノイズと見なされる一方，対角線から遠い点はデータ点群の持つ有意でロバストな構造を表す．

図 7．点データの密度がパーシステンス図上の構造の生成（birth）および消滅（death）に与える影響を表した概念図（Kono et al., 2020）．左から右の順に，点データ，データ上に（ball）を描いた図，半径 rdeath の球を描いた図，および対応する H1
半径 rbirth の球のパーシステンス図を示す．上段の 4 つの図が低密度データ点の場合，下段が高密度の場合である．

点は対角線よりも上に存在する．定義により，対角線の近くにある点の「寿命」は短く，誕生後すぐに死滅することが示されている．対照的に，対角線から遠く離れた点は長く生存できる穴（生成元）を表す．図 6 はこれを端的に示している．本研究は銀河分布の空間構造を議論するため，以下では誕生と死滅に対応する半径を b と d ではなく rbirth および rdeath と表記する．直感的には，rbirth は Hk を構成する点の平均間隔を表し，rdeath は構造のサイズを表している．
パーシステンス図の重要な特性の 1 つについて言及する．リング状の構造を持つ点群データを考える．図 7 ではこのような点群データにおける 2 つの場合のパーシステンス図が示されている．上のパネルはデータが低密度の場合のパーシステントホモロジー群と対応するパーシ
ステンス図，下のパネルは高密度の場合の同じ図を示している．低密度のデータは，実際には図 7 の高密度データの部分集合となっている．つまり，2 つのデータはリング構造に関して
基本的に同じ情報を共有している．ただし，高密度データの rbirth and rdeath の間隔は低密度

172

統計数理

第 71 巻 第 2 号 2023

データの間隔よりも小さい．これは，2 つのデータが同じ構造を持つ場合，データの密度が高い（サンプリングレートが高い）ほど検出された特徴が有意に長い寿命を持つ可能性が高くなることを示唆する．図 7 では低密度データの rbirth が大きく rdeath に近いことが見て取れる．これは rbirth が上記のデータ点間の平均距離の指標であることの直接的な帰結である．対照的に，
rdeath は注目する構造のサイズの指標であるためほとんど変化しない．この事実は，次章以降で示すシミュレーションおよび観測データの解析で検出された特徴の有意性についての議論と関連する．
2.1.6 ボトルネック距離
様々な応用において，異なる 2 つのパーシステンス図 D1 と D2 の差異を比較する方法が必要になる．パーシステンス図はベクトルとならないため，この比較には注意が必要である．これを可能にする方法の 1 つがボトルネック距離（bottleneck distance）である（e.g., Edelsbrunner．
and Harer, 2010; Kerber et al., 2017; Cohen-Steiner et al., 2007; Chazal et al., 2008）
Definition 9.（ボトルネック距離（bottleneck distance）） ボトルネック距離は(2.15)

db (D1 , D2 ) ≡ inf

supm x ∈|D |1x1 − m(x1 ) ∞1と定義される．ここで x ∞ はチェビシェフ距離（Chebyshev distance）
L∞ ≡ max{|x1 |, |x2 |}で，|D| はパーシステンス図 D に対角線 (x, x) 上を重なりを許して埋め尽くす無限個の生成元の集合 Δ を加えた集合 |D| ≡ D ∪ Δ, m は |D1 | と |D2 | の間の全単射である．対角線の多重集
合 Δ は，対応する相手のない生成元の対応相手として用意された，寿命が 0 の生成元の集合である．これによって，2 つのパーシステンス図の生成元の個数が異なるような場合でも直接比較が可能となる．
ボトルネック距離の概念を図 8 に示す．2 つのパーシステンス図が完全に同一である場合，
ボトルネック距離は最小値となる．ここでパーシステンス図の距離の安定性を定義する．
Definition 10.（パーシステンス図の距離の安定性） ハウスドルフ距離 DH を
(2.16)
DH = max sup (x, Y )，sup (y, X)
x∈X

y∈Y

とする．パーシステンス図の距離 d が安定とは，ある定数 L が存在し，任意の点集合 X = {xi },Y = {yj } に対し

図 8．2 つのパーシステンス図間のボトルネック距離を計算する方法の概念．他方のパーシス
テンス図に対応する点がない点は，対角線に存在するとみなして計算する．

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出(2.17)

173

d(D1 , D2 ) ≤ LDH (X, Y )が成り立つことをいう（Lipschitz 連続性）．
, パーシ
ボトルネック距離はこの意味での安定性が保証されており（Cohen-Steiner et al., 2007）
ステンス図の検証のために用いることができる．
2.2 パーシステンス図の逆解析
2 点相関関数による従来の BAO 解析では，データにおける宇宙論的な構造の位置や形状を特定することはできない．これは，この種の解析が相関の統計的な特徴のみを表現するため，
データの持つ平均的な特徴しか抽出できないからである．TDA では，得られたパーシステン
ス図の各生成元と対応する図形構造を特定することができる．このような逆解析は，解析結果の物理的議論を進める上で有用である．半径 r についてのフィルトレーションを通じて，対応する単体複体 σ(r) の位相的情報は保持される．この場合，逆解析の目的は特異的に長寿命な)穴（生成元）(rbirth,i , rdeath,i（つまりパーシステンス図に現れた有意な点に対応するデータ点群の特徴）に対応する単体複体 (σ(rbirth,i ), σ(rdeath,i )) を見つけることである．例えば，穴が消滅する際の穴を構成する単体複体やそのときの半径 rdeath,i から，その穴の特徴的な構造や位置，
スケールといった情報が推定できる．本研究では R TDA パッケージを用いて逆解析を行った．
R TDA では，逆解析だけでなく，異なるパーシステンス図の間の距離，パーシステンス図にある穴構造の存在についての検定などを実施することができる．R TDA はまた C++ で書かれている位相的データ解析アルゴリズム GUDHI（The GUDHI Project, 2015），Dionysus（Morozov,，PHAT（Bauer et al., 2012）のインターフェースも与える．
2007）
分析対象である 3 次元点群データから検出された 3 次元の穴（シェル）構造の逆解析のために，
R パッケージ SCHU（https://github.com/xinxuyale/SCHU; Xu et al., 2019）を用いた．SCHU
は Dionysus での位相的データ解析に最適化されている．これにより，パーシステンス図を元に検出された生成元に対応する特徴的な図形構造の位置と形状を特定することが可能になる．
逆解析法の研究は急速に発展しているため，他のさまざまな洗練された方法が提案されている（たとえば Obayashi et al., 2018 を参照されたい）．様々な方法の比較は本研究の範囲外であり，今後の研究で行う予定である．
3.

宇宙論的シミュレーションを用いた解析法の性能評価

3.1 データ
本研究で開発した方法の性能を検証するため，まず一連の暗黒物質のみの宇宙論的 N 体シ
ミュレーションを実行し，データを生成した．この計算では公開 N 体シミュレーションコー
ド Gadget-2（Springel, 2005）を用い，初期条件は Crocce et al.（2006）の方法を用いて生成した．初期条件は赤方偏移 z = 20 に設定し，シミュレーションボックスは通常の BAO スケール
∼ 150 Mpc よりも十分に大きい一辺 2 Gpc に設定した．BAO 検出のためには質量分解能は重要でないため，シミュレーションボックスには 2563 個の粒子を配置した．これは，粒子あたりの質量 6.4 × 1013 M （M は太陽質量）に対応する．このシミュレーションには冷たい暗黒
粒子しか含まれていないが，本研究で注目する時代は光子 - バリ
物質
（cold dark matter: CDM）
オン脱結合のかなり後であり，また見る必要があるのは重力だけなので，BAO の特徴は初期条件にのみ刻印しておけば一般的に十分である．初期ゆらぎのパワースペクトルは，線型ボルを用いて計算している．これをツマン方程式を数値的に解くコード CLASS（Lesgourgues, 2011）
元に z = 0 でのシミュレーションのスナップショットを取り，パワースペクトル解析を行うこ

174

統計数理

第 71 巻 第 2 号 2023

図 9．赤方偏移 z = 0 でのバリオンあり（右）およびなし（左）の場合のシミュレーションデー
タ．インサートは 100–250 Mpc のボックス領域の拡大図を表す．密度揺らぎは z 軸方
向に 100 Mpc の厚さの領域を射影しており，見やすくするために arcsinh をとった値を示している．

図 10．シミュレーションによる z = 0 でのパワースペクトル（Kono et al., 2020）．Eisenstein and Hu（1998）のバリオン振動のない線型パワースペクトルで規格化されている．黒丸はバリオン振動ありのシミュレーションデータ，点線は同じ宇宙論的パラ
メーターのセットを用い，ハローフィットモデル（Takahashi et al., 2012）によって計算した解析解である．比較のため，バリオン振動なしのシミュレーションデータを四角で示す．バリオンありと比べ，k < 0.2Mpc/h において振動がみられない．

とによってデータに BAO が確かに存在することが確認できる．
コントロールサンプルとして，BAO なしで初期パワースペクトルを置き換えたまったく同
じシミュレーションセットを生成した．これは初期パワースペクトルを計算する際，物質の総質量（CDM+ バリオン）を変更せずにバリオンの量を極限まで減らすことで実現している．すなわち，通常のモデルでは ΩB = 0.049 であるが，コントロールサンプルは ΩB = 0.002204 とした．以降便宜上両者を BAO あり/なし，やバリオンあり/なしと表現し区別することにする．
図 9 はバリオンがある場合とない場合のシミュレーションのスナップショットを示している．
どちらのシミュレーションセットにも同じランダムシードを使用しているため，目視によって違いを確認することは非常に難しい．一方で，図 10 は両者のシミュレーションセットのパ位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

175

図 11．シミュレーションデータのパーシステンス図．黒い点，赤い三角形，青いひし形はそ
れぞれ H0 , H1 , H2 を表す．黒い対角線は rbirth = rdeath である．左パネルはバリ
オンなし，右パネルはバリオンありシミュレーションの結果を示している．

ワースペクトルを示している．それぞれのパワースペクトルは見やすさのために，BAO の振動的な特徴を平滑化したパワースペクトル（Eisenstein and Hu, 1998）によって規格化されている．予想されるとおり，得られたコントロールサンプルのパワースペクトルは BAO の信号を示さない．本論文では，計算機システム上の制約のため 2563 粒子の親サンプルから 2000 個の粒子をランダムに抽出したサブサンプルを用いた．今後はコードを並列化すると同時に計算機をアップグレードすることで，はるかに大きなデータセットが処理可能になる予定である．
3.2 シミュレーションデータの解析結果
シミュレーションデータのパーシステンス図を検証しよう．検出された信号の統計的有意性を評価するには，信頼区間を計算する方法が必要になる．しかし，パーシステンス図の信頼区．その間推定は自明な議論ではなく，いくつかの方法が提案されている（e.g., Fasy et al., 2014）
ひとつとして，ここではボトルネック距離を使用する方法を採用する．本研究ではまず，ブー
トストラップリサンプリングでサンプルのセットを生成する．ここから対を抽出し，その 2 つのパーシステンス図の間の距離をボトルネック距離で測ることによって信頼区間を推定する．
シミュレーションで得られたパーシステンス図を図 11 に示す．左のパネルはバリオンなし，
右のパネルはバリオンありの結果を表す．黒い点，赤い三角形，青いひし形はそれぞれ 0 次(H0 )，1 次 (H1 )，2 次 (H2 ) のパーシステントホモロジー群に対応する．赤と青の点線は，それぞれ H1 と H2 の 90-% 信頼区間である．これらパーシステンス図の信頼区間はリサンプリング
回数 Nboot = 30 のブートストラップリサンプリングによる方法で計算したものである．この研究では，対応する信頼区間の上にある点を統計的に有意であるとみなす．
図 11 の右パネルで示されているように，バリオンありデータから得られた H2 には 4 個の（青いひし形）が検出された．一方バリオンなしデータでは有意な生成元，すなわち 3 次元の穴対応するスケールの有意な穴は検出されなかった．
バリオンありのサンプルで検出された有意な穴の rbirth と rdeath を表 1 に示す．検出された
中空の穴の平均半径は rdeath は 150.16 ± 8.46 [Mpc] であった．検出個数は 4 つだが，取得された r̄death は，BAO で期待される半径と一致する．図 12 に検出された生成元に対応するデータ
点群の構造の位置と形状を示す．この 3 次元構造は，中空の穴となっている．

176

統計数理

第 71 巻 第 2 号 2023

表 1．シミュレーションデータで検出された p < 0.2 の穴の rbirth と rdeath ．

図 12．逆解析の例．シミュレーションデータにおいて最も優位性の高い 4 個の H2 に対応する実空間での構造を示す．網掛け部分は H2 に対応する構造のすべての構成粒子の凸包であり，中空の構造となっている．

一方 H1 に属する生成元（赤い三角形）については，双方のデータで重要な特徴が得られた．
この生成元は 2 次元のループ構造に対応する．バリオンありのデータでは p < 0.2 の 17 個の有意な H1 が検出された．これら有意な生成元の平均 rdeath は 99.00 ± 2.26 [Mpc] である．バリオ
ンなしデータの H1 では 34 個の生成元が検出され，r̄death = 100.49 ± 3.24 [Mpc] となった．H1で検出された生成元に対応する図形構造の特徴的スケールの一致は，バリオンの存在によって影響されない 2 次元のループのような構造があることを示唆する．大規模構造のフィラメント
はこのような構造の候補の可能性がある．これとは対照的に，r̄birth はバリオンの有無で違いがあり，バリオンありで 46.13 ± 2.24 [Mpc], バリオンなしで r̄birth = 62.34 ± 2.84 [Mpc] となった．この解析から，ループを構成する銀河の間の距離はバリオンの存在によって減少すること
が分かる．
（黒い点）
についてはバリオンの有無でパーシステンス図上の分
そして，H0 に属する生成元
布に明らかな違いが見られた H0 は連結成分を現し，銀河団やフィラメントなどの 1 次元的な連続構造に対応する．バリオンなしのデータでは rbirth は対角線に沿って幅広く分布するのに対し，バリオンありでは rbirth は低い値に集中している．バリオンありの場合のような rbirthが局在したような分布は，連結構造がある特定の半径で誕生したことを意味する．このことは，銀河団やフィラメント構造を構成する銀河の近某銀河間距離が均等であり，従って特定の
rbirth で一気に結合して連結構造を形成することを意味する．一方で，rdeath の分布はバリオン
あり・なしに関わらず広く分布していた．このことは，銀河団やフィラメント間の距離はバリ位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

177

オンの影響を受けないことを示唆する．このように H0 のパーシステンス図の特徴は，バリオンの存在が銀河団やフィラメント構造内の近傍銀河間距離に何らかの制約を与える効果を持つ
ことを示唆する．更なる解析を行うことで，H0 と H1 のパーシステンス図上の特徴を物理的に解釈できることが期待される．
4.

Sloan Digital Sky Survey Data Release 14 による解析

4.1 SDSS 銀河探査データ
本研究では，Sloan Digital Sky Survey (SDSS) Data Release 14 (DR14) の extended Baryon
Oscillation Spectroscopic Survey (eBOSS)（Dawson et al., 2016）大規模構造カタログ（Ata et al.,2018）のクエーサーデータを用いた．これは SDSS-IV（Blanton et al., 2017）の一部である．観
測されたバンド（波長帯）は u, g, r, i, および z である．天球上の探査領域は 2044 deg2 , 赤方偏（NGC）で 1288 deg2 ，南銀河半球で移は 0.8 < z < 2.2 の範囲である．有効面積は北銀河半球
2995 deg となっている．クエーサーサンプルは Myers et al.（2015）
によって選択されたデータを使用した．このサンプルにおける BAO の存在は Ata et al.（2018）によって検証され，すでに確認されている．
サンプル全体には 147,000 個のクエーサーが含まれているが，本研究ではシミュレーション
データと同様 2000 個の銀河をランダム抽出して用いた．
4.2 SDSS データ解析結果
SDSS DR14 データのパーシステンス図を図 13 に示す．図 13 より，SDSS DR14 データでは H2に属する 4 つの有意な生成元すなわち「シェル」が検出された．平均の rdeath は 146.6±2.0 [Mpc]である（表 2）．有意な生成元の数は 4 個だが，r̄birth は BAO から予想される半径の期待値と一致している．H2 の逆解析の結果を図 14 に示す．
（ループ）として検出されている．これこの解析では，H1 に属する 19 個の生成元が有意な穴らの生成元の平均 rdeath は 101.82 ± 3.54 [Mpc] となり，シミュレーションデータで得られた
r̄birth と一致している．また観測データの rbirth は 57.92 ± 3.05 [Mpc] となった．そして H0 ホ
モロジーの分布は，興味深いことにバリオンなしシミュレーションの分布と類似している．これら低次元ホモロジーの解釈には今後のより詳細な研究が待たれる．

図 13．SDSS データのパーシステンス図．シンボルは図 11 と同様である．

統計数理

178

第 71 巻 第 2 号 2023

表 2．観測データ（SDSS DR14）で検出された p < 0.2 の rbirth および rdeath ．

図 14．実際の銀河探査（SDSS）データで検出された H2 に属する 4 つの有意な生成元に対応する実空間構造．網掛け部分の意味は図 12 と同様．

5.

まとめと議論本研究では，銀河分布の BAO 信号検出における位相的データ解析の有用性を検証した．本
研究の位相的データ解析で対象としたパーシステントホモロジー群は，様々なスケールに対応
した単体複体の集合列によって定義されたフィルトレーションを構成する．この性質はデータに内在する穴のスケールと形状を分析するのに特に有用である．
まず N = 2563 の N 体シミュレーションを用い，パーシステントホモロジー群に対するバリオンの影響を検証した．バリオンの物理を入れた場合とない場合のシミュレーションで，パーシステンス図は大幅に異なることが示された．バリオンありのデータからは 2 次のパーシステ（穴構造）が検出されたが，バリオンなしのデータからは検出さンス図から 4 つの有意な生成元れなかった．また H1 と H0 に属する生成元のパーシステンス図上での分布にも明らかな違いが見られた．検出された 1 次元の穴の r̄death は同程度の値を取るが，r̄birth は大きく異なる．同様に 0 次元の穴の r̄birth は，バリオンの有無で大きく異なる傾向を示した．
次に，実際の観測データとして SDSS のクエーサー探査から抽出した高々 2000 サブサンプ
ルで BAO 信号を検出することに成功した．得られた r̄death は H2 に対して 146.6 ± 2.0 [Mpc]であり，BAO スケールと一致する．このことは，銀河分布の大規模構造検出においてパーシ
ステントホモロジー群の計算コストが相関関数よりもはるかに安価であることを意味し，BAOの宇宙論的分析の新たな可能性を与える．
本研究の方法による BAO 検出は，伝統的な 2 点相関関数による解析と好対照をなす．2 点

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

179

相関関数の推定には，高密度にサンプリングされた銀河データが必要とされるが，位相的デー
タ解析は ∼ 2000 個の銀河のスパースサンプリングされたデータで BAO のような弱い信号を検出できる．これは，位相的データ解析を活用することで，スパースにサンプリングされた銀河探査も宇宙論的研究に使用できることを意味する．このことは次世代の宇宙論的銀河調査のデザインと戦略に新しい方向性を与えると期待される．
さらに，従来の解析と比較した位相的データ解析の著しい利点は，逆解析により BAO 構造を可視化できることである．パーシステンス図の逆分析が持つ機能は，従来の解析とは一線を画す明確な特長であり，強調に値する．たとえば本研究の直接的発展として，BAO 信号に寄与する実空間での銀河分布構造を簡単に特定できるため，伝統的方法に比べさまざまな宇宙論的議論がずっと簡単に可能になる．ただし，最も単純な位相的データ解析はノイズに対して頑健でないという弱点があり，今後のより詳細な応用ではこの点を克服する必要がある．ロバス
トな位相的データ解析の方法は現在も研究されており（たとえば Vishwanath et al., 2022）．将来の研究ではその成果を取り入れられると期待できる．
謝辞

本稿について大変貴重な改善点の指摘および助言を頂いた査読者のお 2 人，および担当編集
委員の島谷健一郎氏には心からの感謝を申し上げる．池田思朗，福水健次，栗木哲の各氏からは本研究に関する貴重な議論および有益な示唆を頂いた．ここに深く感謝する．本研究はJSPS 科学研究費助成事業（19H05076，20H01932，21H01128，21J23611），新学術領域研究「加，住友財団 2018 年度基礎科学研究助成金（180923），およ速宇宙」（15H05890 および 16H01096）
び統計数学研究所の共同利用「データサイエンス的方法による銀河進化研究の新展開」の補助を受けて行った．
付
A.

録

宇宙論の基礎A.1 Friedmann-Lemaître-Robertson-Walker 計量とスケール因子
一般相対性理論は，時空がその幾何学的構造によって特徴付けられることを示し，微分幾何学が使用するべき基本的な枠組みであることを示した．微分幾何において，多様体上の局所的性質を記述する基本的な概念が線素（line element）あるいは計量（metric）である．Friedmann（FLRW）
計量は，均一で等方的な時空の（局所的）幾何構造を表すたLemaître-Robertson-Walker
めに提唱された計量で，宇宙論のモデルとして広く用いられている．


dr2
2
2
2
2
dθ
ds2 = gμν dxμ dxν = −c2 dt2 + a2 (t)
+
r
+
sin
θdφ
(A.1)
1 − Kr2

ここで gμν が計量テンソル，t は時刻，r, θ, φ は極座標の動径距離および角度方向の座標を表す．また K はガウス曲率，a(t) はスケール因子である．スケール因子は宇宙膨張による空間のと規格スケールの変化を表現しており，宇宙論の慣習上 a0 ≡ a(t0 ) = 1（t0 : 現在の宇宙年齢）−2化して用いられる．ここでは K は（長さ） の次元を持つように定義を採用する．この場合 rは長さの次元を持ち，a(t) は無次元になる．
スケール因子を導入すると，宇宙の進化による天体間の距離の変化から宇宙膨張による「自明な」距離の変化を除いた正味の変化分を表現することができる．即ち，距離をあらわすベク
トルを r とすると

180

統計数理

(A.2)

第 71 巻 第 2 号 2023

r = a(t)
xと書ける．この x は宇宙膨張の影響を除いた距離の指標であり，共動座標（comoving coordinate）と呼ばれる．
A.2 赤方偏移と宇宙論パラメータ
過去のある時刻 t = tem で共動座標 (r, φ, θ) = (re , 0, 0) で放射された光が，時刻 t = t0 において原点 (r, φ, θ) = (0, 0, 0) に位置する観測者に到達する状況を考える．相対性理論では，光はに沿って移動する．ヌル測地線とは ds2 = 0 を満たす経路のことで，
ヌル測地線（null geodesic）
最小作用の原理で実現する軌道と考えてよい．式（A.1）を用いれば，距離と時間を
cdt
dr
= √
a(t)
1 − Kr2(A.3)のように関係付けられる．宇宙膨張の光への影響を議論するため，変数 r を
 rem
dr
√
(A.4)
χ≡
1 − Kr20
によって座標距離 χ に，時間 t を
η≡

(A.5)

 t

dt
a(t )

によって共形時間 η に変換すると，光の満たすべき方程式はc(ηobs − ηem ) = χ(A.6)となる．ここで ηobs は光が観測者に到達した時刻を表す．
時刻 η = η0 および η = η0 + δη0 で放出された光が，それぞれ η1 および η1 + δη1 で観測者にの右辺は η によらないことから到達したとする．式
（A.6）
(A.7)

δη0 = δη1

すなわち
δt0
δt1
=
a(t0 )a(t1 )(A.8)となる．時間間隔 δt0 および δt1 において，光の位相が保存することから，(A.9)

a(t0 )ν0 = a(t1 )ν1 ⇐⇒λ0
λ1=
a(t0 )a(t1 )

が得られる．ここで ν は光の振動数で，波長とは λν = c の関係がある．
膨張する宇宙では a(t0 ) > a(t1 ) なので，観測される時刻での波長は放射のときよりも長くなる，即ち λ0 > λ1 となる．この波長の伸びが宇宙膨張による赤方偏移と呼ばれている．赤方偏
移zは(A.10)z≡

λ0 − λ1a(t0 )=−1
λ1
a(t1 )と定義される．
ここで，上記のようにスケール因子 a(t) を a(t0 ) = 1 となるよう規格化する．このとき a(t)と z との関係は

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出z=(A.11)

181

1
−1,
a

あるいは
a(t) =(A.12)

1
1+zとなる．
A.3 フリードマン方程式と宇宙論パラメータ
ここまでの準備，およびアインシュタイン方程式（Einstein’s equation）を用いれば，膨張する宇宙を記述するための「運動方程式」
を導くことができる．アインシュタイン方程式とは，時空の歪みとエネルギー密度を関連付けた一般相対性理論の基本方程式であり，(A.13)

Rμν −

1 μν
8πG
Rg + Λg μν = 4 T μν
2cと書かれる．ここで Rμν はリッチテンソル（Ricci tensor），R ≡ Rα α はそのトレースでリッチ，G はニュートンの重力定数，Λ は宇宙定数である．詳細は相対性理論
スカラー（Ricci scalar）の文献に譲るが，リッチテンソルは計量 gμν とその微分量の関数であり，時空の歪みを表現する 2 階の共変テンソルである．計量が FLRW メトリックの場合，時空の運動方程式は
 2
ȧ
8πGρ
c2 Λ
c2 K
=
(A.14)
− 2 +
a
3
a
3


ä
4πG
3p
c2 Λ
(A.15)=−
ρ+ 2 +,a3c3となる．これはフリードマン方程式（Friedmann equations）と呼ばれている．ここで ρ は宇宙の（A.14）
はアインシュタイン方程式（式（A.13））の 00 成分，式（A.15）は密度，p は圧力を現す．式
トレースから求められる．
さらに，宇宙論パラメーターと呼ばれる以下の量を定義すると，宇宙膨張の物理量依存性がより具体的に記述できる．
• ハッブルパラメータ：(A.16)

H(t) ≡

ȧ(t)
.
a(t)

これは時間依存する量であるが，観測的には現在値 H0 が用いられることが多い．また無次元化した h = H0 /100 もよく用いる．
• 密度パラメータ：(A.17)

ΩM (t) ≡ρ(t)8πGρ(t)≡,ρc (t)
3H 2 (t)ここで ρc は臨界密度パラメータと呼ばれる量で，現在の宇宙では ρc = 1.88 × 10−29 h2 [g cm−3 ]と測定されている．
• 無次元宇宙項パラメータ：(A.18)• 曲率パラメータ：ΩΛ (t) ≡c2 Λ.
3H 2 (t)182(A.19)

統計数理

第 71 巻 第 2 号 2023

ΩK (t) ≡ −c2 K.
a(t)2 H 2 (t)式（A.16）（A.19）–
を用いると，式（A.14）
は(A.20)

ΩM (t) + ΩΛ (t) + ΩK (t) = 1と表せる．これは一般の宇宙年齢 t で成り立つので，現在の宇宙年齢 t0 では(A.21)

ΩM0 + ΩΛ0 + ΩK0 = 1となる．
導入で述べているように，本論文では Planck 衛星による CMB 観測で推定された宇宙論パラ
メータ h = H0 /100 [km s−1 Mpc−1 ] = 0.6766, ΩΛ,0 = 0.6889, ΩM,0 = 0.3111, ln(1010 As ) = 3.047を採用している．曲率パラメータ ΩK0 = 0 である．これらの値は，現在の宇宙は空間的に平坦で，宇宙定数（または暗黒エネルギー）により加速膨張している宇宙モデルを強く示唆する．
パラメータ As は宇宙初期のホライズンスケールの揺らぎの振幅を決めるもので，シミュレー
ションデータを生成する際にこの値を採用している．

参

考

文

献

Appleby, S., Chingangbam, P., Park, C., Yogendran, K. P. and Joby, P. K. (2018). Minkowski tensors
in three dimensions: Probing the anisotropy generated by redshift space distortion, AstrophysicalJournal, 863(2), 200.
Ata, M., et al. (2018). The clustering of the SDSS-IV extended Baryon Oscillation Spectroscopic Survey
DR14 quasar sample: First measurement of baryon acoustic oscillations between redshift 0.8 and
2.2, Monthly Notices of the Royal Astronomical Society, 473(4), 4773–4794.
Bag, S., Mondal, R., Sarkar, P., Bharadwaj, S. and Sahni, V. (2018). The shape and size distribution of
H II regions near the percolation transition, Monthly Notices of the Royal Astronomical Society,477(2), 1984–1992.
Bag, S., Mondal, R., Sarkar, P., Bharadwaj, S., Choudhury, T. R. and Sahni, V. (2019). Studying
the morphology of H I isodensity surfaces during reionization using Shapeﬁnders and percolation
analysis, Monthly Notices of the Royal Astronomical Society, 485(2), 2235–2251.
Bauer, U., Kerber, M. and Reininghaus, J. (2012). PHAT, a software library for persistent homology,
https://bitbucket.org/phat-code/phat/src/master/.
Beisbart, C., Dahlke, R., Mecke, K. and Wagner, H. (2002). Vector- and tensor-valued descriptors
for spatial patterns, Morphology of Condensed Matter, Lecture Notes in Physics, Volume 600,
238–260, Springer, New York.
Bernardeau, F., Colombi, S., Gaztañaga, E. and Scoccimarro, R. (2002). Large-scale structure of theuniverse and cosmological perturbation theory, Physics Reports, 367(1-3), 1–248.
Biagetti, M., Cole, A. and Shiu, G. (2021). The persistence of large scale structures I: Primordial
non-Gaussianity, Journal of Cosmology and Astroparticle Physics, 2021(4), 061.
Blanton, M. R., et al. (2017). Sloan Digital Sky Survey IV: Mapping the Milky Way, nearby galaxies,
and the distant universe, Astronomical Journal, 154(1), 28.
Boerner, G. and Mo, H. (1989). A percolation analysis of cluster superclustering, Astronomy andAstrophysics, 224(1-2), 1–9.
Boissonnat, J.-D., Chazal, F. and Yvinec, M. (2018). Geometric and Topological Inference, Cambridge
Texts in Applied Mathematics, Cambridge University Press, Cambridge.
Chazal, F., Cohen-Steiner, D., Glisse, M., Guibas, L. J. and Oudot, S. (2008). Proximity of persistence
modules and their diagrams, Research Report, RR-6568, INRIA, 29.

位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

183

Chen, Z., Xu, Y., Wang, Y. and Chen, X. (2019). Stages of reionization as revealed by the Minkowski
functionals, Astrophysical Journal, 885(1), 23.
Chingangbam, P., Ganesan, V., Yogendran, K. P. and Park, C. (2017). On Minkowski functionals of
CMB polarization, Physics Letters B, 771, 67–73.
Cohen-Steiner, D., Edelsbrunner, H. and Harer, J. (2007). Stability of persistence diagrams, Discrete& Computational Geometry, 37, 103–120.
Cole, A. and Shiu, G. (2018). Persistent homology and non-Gaussianity, Journal of Cosmology andAstroparticle Physics, 2018(3), 025.
Crocce, M., Pueblas, S. and Scoccimarro, R. (2006). Transients from initial conditions in cosmologicalsimulations, Monthly Notices of the Royal Astronomical Society, 373(1), 369–381.
Daley, D. J. and Vere-Jones, D. (2003). An Introduction to the Theory of Point Processes. Vol. I,
Probability and Its Applications, 2nd ed., Springer-Verlag, Berlin.
Dawson, K. S., et al. (2016). The SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Overview
and early data, Astronomical Journal, 151(2), 44.
de Lapparent, V., Geller, M. J. and Huchra, J. P. (1986). A slice of the universe, Astrophysical Journal,
302, L1–L4.
Edelsbrunner, H. and Harer, J. (2010). Computational Topology — An Introduction, American Mathematical Society, Providence, Rhode Island, http://www.ams.org/bookstore-getitem/item=MBK69.
Edelsbrunner, H., Letscher, D. and Zomorodian, A. (2002). Topological persistence and simpliﬁcation,
Discrete & Computational Geometry, 28(4), 511–533.
Efstathiou, G. and Silk, J. (1983). The formation of galaxies, Fundamentals of Cosmic Physics, 9,
1–138.
Eisenstein, D. J. and Hu, W. (1998). Baryonic features in the matter transfer function, AstrophysicalJournal, 496(2), 605–614.
Eisenstein, D. J., et al. (2005). Detection of the baryon acoustic peak in the large-scale correlation
function of SDSS luminous red galaxies, Astrophysical Journal, 633(2), 560–574.
Eisenstein, D. J., Seo, H.-J., Sirko, E. and Spergel, D. N. (2007). Improving cosmological distance measurements by reconstruction of the baryon acoustic peak, Astrophysical Journal, 664(2), 675–679.
Elbers, W. and van de Weygaert, R. (2019). Persistent topology of the reionization bubble network —I. Formalism and phenomenology, Monthly Notices of the Royal Astronomical Society, 486(2),1523–1538.
Fang, W., Li, B. and Zhao, G.-B. (2017). New probe of departures from general relativity using
Minkowski functionals, Physical Review Letters, 118(18), 181301.
Fasy, B. T., Lecci, F., Rinaldo, A., Wasserman, L., Balakrishnan, S. and Singh, A. (2014). Conﬁdence
sets for persistence diagrams, The Annals of Statistics, 42(6), 2301–2339.
Feldbrugge, J., van Engelen, M., van de Weygaert, R., Pranav, P. and Vegter, G. (2019). Stochastic homology of Gaussian vs. non-Gaussian random ﬁelds: Graphs towards Betti numbers andpersistence diagrams, Journal of Cosmology and Astroparticle Physics, 2019(9), 052.
Ganesan, V. and Chingangbam, P. (2017). Tensor Minkowski functionals: First application to the
CMB, Journal of Cosmology and Astroparticle Physics, 2017(6), 023.
Gay, C., Pichon, C. and Pogosyan, D. (2012). Non-Gaussian statistics of critical sets in 2D and 3D:
Peaks, voids, saddles, genus, and skeleton, Physical Review D, 85(2), 023011.
Giri, S. K. and Mellema, G. (2021). Measuring the topology of reionization with Betti numbers, Monthly
Notices of the Royal Astronomical Society, 505(2), 1863–1877.
Gott, I., Gott, J. R., III, Melott, A. L. and Dickinson, M. (1986). The sponge-like topology of large-scale
structure in the universe, Astrophysical Journal, 306, 341–357.
Gott, I., Gott, J. R., III, Weinberg, D. H. and Melott, A. L. (1987). A quantitative approach to the
topology of large-scale structure, Astrophysical Journal, 319, 1–8.

184

統計数理

第 71 巻 第 2 号 2023

Gott, I., Gott, J. R., III, et al. (1989). The topology of large-scale structure. III. Analysis of observations,Astrophysical Journal, 340, 625–646.
Hikage, C., et al. (2003). Minkowski functionals of SDSS galaxies I: Analysis of excursion sets, Publications of the Astronomical Society of Japan, 55, 911–931.
Joby, P. K., Chingangbam, P., Ghosh, T., Ganesan, V. and Ravikumar, C. D. (2019). Search for anomalous alignments of structures in Planck data using Minkowski tensors, Journal of Cosmology andAstroparticle Physics, 2019(1), 009.
Junaid, M. and Pogosyan, D. (2015). Geometrical measures of non-Gaussianity generated from single
ﬁeld inﬂationary models, Physical Review D, 92(4), 043505.
Kerber, M., Morozov, D. and Nigmetov, A. (2017). Geometry helps to compare persistence diagrams,
Journal of Experimental Algorithmics, 22, 1.4.
Kerscher, M. and Tikhonov, A. (2010). Morphology of the local volume, Astronomy and Astrophysics,
509, A57.
Kerscher, M., et al. (1997). Minkowski functionals of Abell/ACO clusters, Monthly Notices of the Royal
Astronomical Society, 284(1), 73–84.
Kitaura, F.-S. and Angulo, R. E. (2012). Linearization with cosmological perturbation theory, Monthly
Notices of the Royal Astronomical Society, 425(4), 2443–2454.
Kono, K. T., Takeuchi, T. T., Cooray, S., Nishizawa, A. J. and Murakami, K. (2020). A study on thebaryon acoustic oscillation with topological data analysis, arXiv:2006.02905.
Kozlov, D. (2008). Combinatorial Algebraic Topology, Springer, Berlin.
Lesgourgues, J. (2011). The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview,
arXiv:1104.2932
Libeskind, N. I., et al. (2018). Tracing the cosmic web, Monthly Notices of the Royal Astronomical
Society, 473(1), 1195–1217.
Ling, C., Wang, Q., Li, R., Li, B., Wang, J. and Gao, L. (2015). Distinguishing general relativity andf (R) gravity with the gravitational lensing Minkowski functionals, Physical Review D, 92(6),064024.
Liu, J., Bird, S., Zorrilla Matilla, J. M., Hill, J. C., Haiman, Z., Madhavacheril, M. S., Petri, A.
and Spergel, D. N. (2018). MassiveNuS: Cosmological massive neutrino simulations, Journal of
Cosmology and Astroparticle Physics, 2018(3), 049.
Ma, C.-P. and Bertschinger, E. (1995). Cosmological perturbation theory in the synchronous and conformal Newtonian gauges, Astrophysical Journal, 455, 7–52.
Martinez, V. and Saar, E. (2001). Statistics of the Galaxy Distribution, CRC Press, London.
Matsubara, T. (1994). Analytic expression of the genus in a weakly non-Gaussian ﬁeld induced bygravity, Astrophysical Journal, 434, L43–L46.
Matsubara, T. (1996). Statistics of isodensity contours in redshift space, Astrophysical Journal, 457,
13–17.
松原 隆彦 (1996). ジーナス統計を用いた宇宙の大構造の定量化, 日本物理學會誌, 51(4), 266–273.
Matsubara, T. (2003). Statistics of smoothed cosmic ﬁelds in perturbation theory. I. Formulation anduseful formulae in second-order perturbation theory, Astrophysical Journal, 584(1), 1–33.
Matsubara, T. (2010). Analytic Minkowski functionals of the cosmic microwave background: Secondorder non-Gaussianity with bispectrum and trispectrum, Physical Review D, 81(8), 083505.
Matsubara, T. and Suto, Y. (1996). Nonlinear evolution of genus in a primordial random Gaussian
density ﬁeld, Astrophysical Journal, 460, 51–58.
Matsubara, T. and Yokoyama, J. (1996). Genus statistics of the large-scale structure with non-Gaussian
density ﬁelds, Astrophysical Journal, 463, 409–419.
Matsubara, T., Hikage, C. and Kuriki, S. (2022). Minkowski functionals and the nonlinear perturbation
theory in the large-scale structure: Second-order eﬀects, Physical Review D, 105(2), 023527.
Mecke, K. R., Buchert, T. and Wagner, H. (1994). Robust morphological measures for large-scale位相的データ解析による銀河分布の定量化とバリオン音響振動抽出

185

structure in the universe, Astronomy and Astrophysics, 288, 697–704.
Morozov, D. (2007). Dionysus, a C++ library for computing persistent homolog, https://www.mrzv.org/
software/dionysus/.
Myers, A. D., et al. (2015). The SDSS-IV extended Baryon Oscillation Spectroscopic Survey: Quasar
target selection, Astrophysical Journal Supplement Series, 221(2), 27.
Obayashi, I., Hiraoka, Y. and Kimura, M. (2018). Persistence diagrams with linear machine learningmodels, Journal of Applied and Computational Topology, 1, 421–449.
Pâris, I., et al. (2018). The Sloan digital sky survey quasar catalog: Fourteenth data release, Astronomy
and Astrophysics, 613, A51.
Parroni, C., Cardone, V. F., Maoli, R. and Scaramella, R. (2020). Going deep with Minkowski functionals of convergence maps, Astronomy and Astrophysics, 633, A71.
Peebles, P. J. E. (1980). The Large-scale Structure of the Universe, Princeton University Press, New
Jersey.
Peebles, P. J. E. and Yu, J. T. (1970). Primeval adiabatic perturbation in an expanding universe,
Astrophysical Journal, 162, 815–836.
Planck Collaboration, Ade, P. A. R., et al. (2016). Planck 2015 results. XVII. Constraints on primordial
non-Gaussianity, Astronomy and Astrophysics, 594, A17.
Pogosyan, D., Gay, C. and Pichon, C. (2009). Invariant joint distribution of a stationary random ﬁeld
and its derivatives: Euler characteristic and critical point counts in 2 and 3D, Physical Review
D, 80(8), 081301.
Pranav, P., et al. (2017). The topology of the cosmic web in terms of persistent Betti numbers, Monthly
Notices of the Royal Astronomical Society, 465(4), 4281–4310.
Pranav, P., et al. (2019a). Topology and geometry of Gaussian random ﬁelds I: on Betti numbers, Euler
characteristic, and Minkowski functionals, Monthly Notices of the Royal Astronomical Society,
485(3), 4167–4208.
Pranav, P., et al. (2019b). Unexpected topology of the temperature ﬂuctuations in the cosmic microwave
background, Astronomy and Astrophysics, 627, A163.
Santos, L., Wang, K. and Zhao, W. (2016). Probing the statistical properties of CMB B-mode polarization through Minkowski functionals, Journal of Cosmology and Astroparticle Physics, 2016(7),029.
Schmalzing, J. and Buchert, T. (1997). Beyond genus statistics: A unifying approach to the morphology
of cosmic structure, Astrophysical Journal, 482(1), L1–L4.
Schmalzing, J. and Gorski, K. M. (1998). Minkowski functionals used in the morphological analysis
of cosmic microwave background anisotropy maps, Monthly Notices of the Royal Astronomical
Society, 297(2), 355–365.
Shandarin, S. F. and Zeldovich, I. B. (1983). Topology of the large-scale structure of the universe,
Comments on Astrophysics, 10, 33–45.
Shiraishi, M., Hikage, C., Namba, R., Namikawa, T. and Hazumi, M. (2016). Testing statistics of
the CMB B-mode polarization toward unambiguously establishing quantum ﬂuctuation of thevacuum, Physical Review D, 94(4), 043506.
Springel, V. (2005). The cosmological simulation code GADGET-2, Monthly Notices of the Royal
Astronomical Society, 364(4), 1105–1134.
Spurio Mancini, A., et al. (2018). 3D cosmic shear: Numerical challenges, 3D lensing random ﬁelds
generation, and Minkowski functionals for cosmological inference, Physical Review D, 98(10),103507.
Sullivan, J. M., Wiegand, A. and Eisenstein, D. J. (2019). The clustering of galaxies in the SDSS-III
Baryon Oscillation Spectroscopic Survey: Evolution of higher-order correlations demonstrated
with Minkowski functionals, Monthly Notices of the Royal Astronomical Society, 485(2), 1708–
1719.

186

統計数理

第 71 巻 第 2 号 2023

Sunyaev, R. A. and Zeldovich, Y. B. (1970). Small-scale ﬂuctuations of relic radiation, Astrophysics
and Space Science, 7(1), 3–19.
Takahashi, R., Sato, M., Nishimichi, T., Taruya, A. and Oguri, M. (2012). Revising the haloﬁt modelfor the nonlinear matter power spectrum, Astrophysical Journal, 761(2), 152.
The GUDHI Project (2015). GUDHI User and Reference Manual, http://gudhi.gforge.inria.fr/doc/latest/.
Vicinanza, M., Cardone, V. F., Maoli, R., Scaramella, R., Er, X. and Tereno, I. (2019). Minkowski
functionals of convergence maps and the lensing ﬁgure of merit, Physical Review D, 99(4), 043534.
Vishwanath, S., Sriperumbudur, B. K., Fukumizu, K. and Kuriki, S. (2022). Robust topological inference in the presence of outliers, arXiv:2206.01795.
Wasserman, L. (2018). Topological data analysis, Annual Review of Statistics and Its Application, 5(1),501–532.
Wiegand, A. and Eisenstein, D. J. (2017). The clustering of galaxies in the SDSS-III Baryon Oscillation
Spectroscopic Survey: Higher order correlations revealed by germ-grain Minkowski functionals,
Monthly Notices of the Royal Astronomical Society, 467(3), 3361–3378.
Wiegand, A., Buchert, T. and Ostermann, M. (2014). Direct Minkowski functional analysis of large
redshift surveys: A new high-speed code tested on the luminous red galaxy Sloan Digital Sky
Survey-DR7 catalogue, Monthly Notices of the Royal Astronomical Society, 443(1), 241–259.
Wilding, G., et al. (2021). Persistent homology of the cosmic web — I. Hierarchical topology in ΛCDM
cosmologies, Monthly Notices of the Royal Astronomical Society, 507(2), 2968–2990.
Xu, X., Cisewski-Kehe, J., Green, S. B. and Nagai, D. (2019). Finding cosmic voids and ﬁlament loops
using topological data analysis, Astronomy and Computing, 27, 34–57.
Yoshiura, S., Shimabukuro, H., Takahashi, K. and Matsubara, T. (2017). Studying topological structureof 21-cm line ﬂuctuations with 3D Minkowski functionals before reionization, Monthly Notices of
the Royal Astronomical Society, 465(1), 394–402.
Zomorodian, A. and Carlsson, G. (2004). Computing persistent homology, Proceedings of the Twentieth
Annual Symposium on Computational Geometry, SCG’04, 347–356, Association for Computing
Machinery, New York.

Proceedings of the Institute of Statistical Mathematics Vol. 71, No. 2, 159–187 (2023)

187

Quantification of Galaxy Distribution with Topological Data Analysis andDetection of the Baryon Acoustic Oscillation
Tsutomu T. Takeuchi1,2 , Kai T. Kono1 , Suchetha Cooray1,3,4 , Atsushi J. Nishizawa5,6,7 ,Koya Murakami1 , Hai-Xia Ma1 and Yoh-Ichi Mototake8
1 Division of Particle and Astrophysical Science, Nagoya University
2 The Institute of Statistical Mathematics
3 Division of Science, National Astronomical Observatory of Japan
4 Research Fellow of the Japan Society of the Promotion of Science (PD)
5 DX Center, Gifu Shotoku Gakuen University
6 Institute for Advanced Research, Nagoya University
7 Kobayashi Maskawa Institute, Nagoya University
8 The Institute of Statistical Mathematics; Now at Graduate School of Data Science, Hitotsubashi

University

Galaxies are distributed inhomogeneously in space, with clusters, groups, ﬁlaments
and voids. This is called the large-scale structure in the Universe. In astronomy, it
is customary to refer to the normal matter as baryons. What we directly observe isthe large-scale structure made of baryons. The large-scale structure has formed mainly
through the gravitational instability, but there is another origin of the structure exists,
which is the acoustic oscillation of baryons at the baryon-photon decoupling. This is
imprinted on the spatial distribution of galaxies in the Universe, known as the baryon
acoustic oscillation (BAO). In this work, we analyzed the spatial distribution of galaxies
with a method from the topological data analysis (TDA), in order to examine the BAO
signal in the galaxy distribution. The TDA provides a method to treat various types of“holes” in point set data, by constructing the persistent homology group from the geometric structure of data points and handling the topological information of the dataset.
We can obtain the information on the size, position, and statistical signiﬁcance of theholes in the data. A particularly strong point of the persistent homology is that it canclassify the holes by their spatial dimension. We ﬁrst analyzed the simulation datasetswith and without the baryon physics to examine the performance of the PH method.
We found that the persistent homology is indeed able to detect the BAO signal among
the large-scale structures in the Universe: simulation data with baryon physics present a
prominent signal from the BAO, while data without baryon physics does not show thissignal. Then, we applied the persistent homology to a quasar (a kind of active galaxies)
sample at z < 1.0 from extended Baryon Oscillation Spectroscopic Survey in Sloan Digital
Sky Survey Data Release 14 (SDSS DR14). We discovered a characteristic hole (a hollow
shell) at a scale r ∼ 150 [Mpc]. This corresponds exactly to the BAO signature imprinted
in the galaxy/quasar distribution. We performed this analysis on a small subsample of2000 quasars. This clearly demonstrates that the PH analysis is very eﬃcient in ﬁnding
this type of topological structures even if the sampling is very sparse.
Key words: Topological data analysis (TDA), persistent homology, galaxy distribution, the large-scale structure in the Universe, baryon acoustic oscillation.