1．デ ー タ の
形
膨大な実験画像データ，Protein Data Bank，大規模 シミュレーションで再現される原子配置データ等々，現 在の諸科学・産業の現場ではデータがあふれている．こ のようなデータ群から考察対象に対して新たな発見や理 解を深めることを可能にするデータ科学が昨今人気であ る．このようなデータ駆動型の問題では「データが多い」 ことよりもむしろ「データが複雑」であることが問題で ある場合が多い．その場合一見複雑な構造に見えるデー タを何らかの形で変換し低次元化できれば，よりデータ 構造の本質に迫ることが可能になるであろう．複雑な データに対して適切な記述子を開発することは現在の データ科学における重要な課題となっている． 例えば N 次元ユークリッド空間の有限個の点の集ま り P＝{ xi∈RN｜i ＝1, …, m} （1） からなるデータを考える．簡単な状況として P が周期的 に配置された点の集まりの場合は，フーリエ解析が代表 的な記述子を与える．また P がいくつかのクラスタの集 まりで構成されているときは，各クラスタを記述子とし データを表現することも可能であろう．周期系 P では点 達の周期的なつながり，クラスタ系では点達の近さに関 するつながりといったように，点達のつながり具合を特 徴付けることでデータの構造は表現できる．しかしなが ら，後述するガラスの原子配置などからも明らかなよう に，現在のデータ科学が直面している諸問題に対してこ れらの古典的な解析手法は力不足であり，より詳細な構 造を表現できる記述子の開発が望まれている． 例えば図 1 左で与えられる平面上のデータを考えてみ る．この場合クラスタとしての特徴付けではなく，「穴」 の存在としてデータのつながりを特徴付けるほうが自然 に見える．このような図形のつながり具合をおおざっぱ に調べることは，現代数学の一分野であるトポロジーが 得意とするところである．特に，穴の存在に着目した特 徴抽出は，従来ホモロジーと呼ばれる概念がその役割を 担ってきた．実際，図 1 のようなデータに対しては，各 点を適当に膨らませた球で置き換えることで，一次ホモ ロジーを用いてデータに内在する穴を検出することがで きる＊1．さらにホモロジーの概念を拡張したパーシステ ントホモロジーを用いることで，穴の存在だけでなくそ の形や大きさなどのより詳細な幾何情報も得ることがで き，「データの形」を調べることが可能になる． 本稿では大学 1 年生程度の基本的な線形代数の知識の
パーシステントホモロジーと材料構造解析
Persistent Homology and Structural Analysis in Materials Science
平岡 裕章
京都大学高等研究院高等研究センター・ヒト生物学高等研究拠点（WPI-ASHBi），Yasuaki Hiraoka 理化学研究所革新知能統合研究センター
KUIAS, Kyoto University. / AIP, RIKEN. [email protected]
大林 一平
理化学研究所革新知能統合研究センターIppei Obayashi AIP, RIKEN.
赤木 和人
東北大学材料科学高等研究所Kazuto Akagi AIMR, Tohoku University.
Keywords:
persistent homology, topological data analysis, materials informatics. 「マテリアルズインフォマティクス」＊1 クラスタ解析は 0 次ホモロジーに対応する．
図 1 「穴」の構造をもつデータ（左）とそれを球体で膨らませた モデル（中央，右）．
みを仮定して，パーシステントホモロジーについての解 説を試みる．その後，著者らが近年取り組んでいるシリ カガラスへの応用，材料構造の中距離秩序の抽出，およ び機械学習との融合についても紹介する．なお，ここで 取り扱う内容のより詳細については文献 [Edelsbrunner 10, 平岡 13, 平岡 17] を参照されたい．また，我々のグ ループではここで紹介するパーシステントホモロジー の 計算ソフトウェア「HomCloud」も開発している [HomCloud]． このソフトウェアはトポロジーの予備知識を仮定しな い素人でも，パーシステントホモロジーを用いて材料の 構造解析までを楽しめる．本記事と合わせて HomCloud を実際に計算させることで，より理解が深まると思われ る．なお，本解説は [平岡 17, 平岡 19] に加筆修正を施 したものである．
2．
幾 何 モ デ ル
点列データ（1）に対して適当なデータ処理を施しそ の形を調べていこう．その際，（i）穴の抽出が可能，（ii） 計算機での処理が可能，の 2 点を満たす数学的概念に Čech 複体と呼ばれる単体複体（多面体）がある． ここで有限集合 V 上の単体複体 X とは，V の部分集 合の集まりで，（i）すべての点 a ∈ V について {a}∈ X， （ii）σ∈ X, τ⊂σ⇒τ∈ X，を満たすもので与えられる． 直感的には頂点，辺，三角形，高次元の三角形を張り合 わせてできる多面体を想像したらよい． Xの要素σで点の個数が｜σ｜＝k＋1 個となるものを k単体といい，k を単体σの次元という（｜σ｜はσの点 の個数）．単体複体 X の次元 dim X を，X に含まれる単 体の次元の最大値で定め，また記号 Xkで X の k 単体の 集まりを表すことにする． このとき Čech 複体C（P, r）とは P 上の単体複体であり， その単体は {xi0, …, xik}∈ C（P, r）⇐⇒ k a＝0B（xr ia）≠Ø で割り当てられる．ここで B（xr i）＝{ x∈RN｜x−xi r} は点 xiを中心とする半径 r の球を，記号 Ø は空集合を 表す．つまり半径 r の球の交わりを調べることで単体を 構成していくのである． 図 2 は，パラメータ r を変化させた際の，Čech 複体 の変化を表示している．例えば一番左の解像度では，球 の交わりがないので，辺や三角形が現れない 3 頂点のみ で構成される単体複体が Čech 複体として得られる．真 ん中の解像度では，二つの球の交わりはあるのですべて の辺が現れているが，三つの球の交わりはないので三角 形は現れていない．ここで，各半径 r ごとに上段の球の 和集合と下段の Čech 複体は，互いに連続変形で移り合 えることに気がつく．これは脈体定理と呼ばれるトポロ ジーの基本定理からの帰結であり， m i＝1B（xr i）と C（P, r） は切り張りしない連続変形で移り合える＊2ことが保証さ れている．よって，特に穴の情報は保存される． このように Čech 複体はデータ P とパラメータ r によ り定まる．パラメータ r が十分小さいとき，Čech 複体 は P と同じ離散的なデータ構造を与える．一方，r が徐々 に大きくなると，データ間のつながりが生まれ始め，十 分大きくするとすべてが完全につながった状態になる． このように，r はデータのつながり具合を操る解像度と みなすことができる． また，パラメータの増大列 r1＜…＜rnに応じて，Čech 複体の増大列 C（P, r1）⊂…⊂C（P, rn）が得られることに も注意する．このような単体複体の増大列のことをフィ ルトレーションと呼ぶ． 実際の数値計算では，ここで紹介した Čech 複体と同 じ性質をもち，かつ計算効率の良いアルファ複体がよく 用いられる [Edelsbrunner 10]．その説明にはいくつか 準備が必要なためここでは割愛するが，ソフトウェアと しては CGAL [CGAL] が公開されている．また本稿では 式（1）の形の入力データのみを扱うが，その他の重要 な幾何モデルとして方体複体があげられる．方体複体に 対しても次章以降の話題は適用でき，それらは画像解析 などへの応用が強く期待されている．方体複体の詳細に ついては文献 [Kaczynski 04] をあげておく．3．ホ モ ロ ジ ー
ホモロジーとは，単体複体の中に存在する穴を代数的 に抽出する道具である．ここでは，説明の簡略化と次章 以降の話題を考慮して，ベクトル空間としてホモロジー を導入する． 単体複体 X の頂点集合を V＝{1, …, n} としよう．こ こで単体σの頂点を順序付きで並べたものを〈σ〉＝〈v0… vk〉と書くことにする．つまり v0＜…＜vkであり，不 ＊2 ホモトピー同値． 図 2 上段は球の増大列． 下段は対応する Čech 複体の増大列等号は V＝{1, …, n} 上に定まる通常のものを用いる． まず単体の次元 k ごとに，k 単体を基底としてもつ体 K 上のベクトル空間＊3 C（X）＝k σ ∈ Xk a〈σ σ〉｜aσ∈K , 0≤k≤dimX を導入する．ここで k＜0, k＞dim X では C（X）＝0 とk しておく．ホモロジーの考え方では，穴を抽出するため にまずは境界を特徴付ける．そこで，境界作用素と呼ば れる線形写像∂k：C（X）Ck k−1（X）を k 〈v0…vk〉＝ k i＝0 （－1）〈vi 0…vi…vk〉, k a〈σσ〉＝ aσ 〈kσ〉 ∂ ∂ ∂ で導入する．ここで〈v0…v̂ …vi k〉は viを除くことを意 味する．また 2 番目の式は 1 番目の式の線形拡張であ る．定義に従って計算をすると，∂k◦ ∂k＋1＝0 を確認でき， これより im∂k＋1⊂ ker ∂kが従う＊4．このとき，単体複 体 X の k 次ホモロジーは，商ベクトル空間
H（X）＝ker ∂k k/im∂k＋1＝{[c]｜c ∈ ker ∂k}
で定められる．ここで [c]＝c＋im ∂k＋1 は同値類を表す 記号である．また H（X）のベクトル空間としての次元k をβ（X）と書き k 次ベッチ数と呼ぶ．次の例で確かめk るが，H（X）は X に存在する k 次元の穴をベクトル空k 間として表現し，またベッチ数β（X）は k 次元の穴のk 数を表すことになる． Example 3.1 頂点集合 V＝{1, 2, 3} 上の単体複体 X＝{{1, 2}, {1, 3}, {2, 3}, {1}, {2}, {3}} のホモロジーを計算してみる．境界作用素は ∂k＝0（k≠1） と 〈12〉〈13〉〈23〉 1＝ －1 －1 0 〈1〉 1 0 －1 〈2〉 0 1 1 〈3〉 ∂ で与えられる．よって，一次ホモロジーは H（X）＝ker ∂1 1＝span{〈12〉−〈13〉＋〈23〉} K となり，一次元の穴（リング）〈12〉−〈13〉＋〈23〉から 生成される．同様に計算をしてみると，0 次ホモロジー も一次元ベクトル空間 H0（X） K となり，0 次元の穴 （連結成分）が一つ存在することを示している． 包含関係にある単体複体の組 X ⊂ Y には，それぞれ ホモロジーが定まるが，さらに線形写像 H（X）∋[c] [c]∈ Hk （Y）k が誘導される．ここで定義域の[c]はH（X）での同値類を，k 値域の [c] は H（Y）での同値類をそれぞれ示している．k Example 3.2 X を Example 3.1 と同じ単体複体とし， Y＝X ∪ {1, 2, 3} とする．単体複体 Y では，im∂2＝ker∂1＝ span{〈12〉−〈13〉＋〈23〉}より，H（Y）＝ 0 である．よっ1 て，誘導される線形写像 H（X）H1 （Y）は零写像とな1 る．これは，包含写像 X Yで穴が消滅していること を示している．一方，Y＝X ∪{4}∪{1,4} としておくと， 誘導される線形写像 H（X）H1 （Y1 ）は同型写像となる． これは包含写像 X Yでの穴の存続を意味する． Čech 複体 C（P, r）のホモロジーを計算することで，デー タ P に内在する穴の情報が取り出せることになる．例え ば図 1 のデータについては，左，中央，右の三つの半径 パラメータ r1＜ r2＜ r3ごとに，一次ホモロジーは H（C（P, r1 1））＝0, H（C（P, r1 i）） K（i＝2, 3） で与えられる．ここで，r2では右下のほうに存在する小 さな穴，r3ではデータ全体で構成されている大きな穴が， 一次元のベクトル空間として表現されている．しかしな がら，これらの一次元ベクトル空間を個別に調べるだけ では，それぞれの穴の大きさや，r2での穴が消滅し r3で 新たな穴が発生しているといった存続性を調べることは できない． このような穴の性質を調べるには，Example 3.2 の ように，包含関係が誘導する線形写像 H（C（P, r1 2）） H（C（P, r1 3））を調べればよい．この場合の線形写像は零 写像で与えられ，穴の消滅と発生がきちんと記録されて いることに注意されたい． このように，ホモロジーと包含写像が誘導する線形写 像をうまく組み合わせることで，データに存在する穴と その幾何学的な情報を抽出することが可能になる．次に 説明するパーシステントホモロジーは，この考え方をさ らに発展させた道具である．
4．パーシステントホモロジー
前章のホモロジーは単体複体に対して定められたが， パーシステントホモロジーは単体複体のフィルトレー ション X：X1⊂…⊂ Xt⊂…に対して定義される．ここ でフィルトレーション X は，あるパラメータ n 以降で は Xt＝X（t n）であると仮定する．また，単体n σ∈ Xnの発生時刻を t σ＝min{t｜σ∈ Xt}で定める．このと きパーシステントホモロジー H（X）は，各 Xk iでのホモ ロジー Hk（Xi）とその間に誘導される線形写像の列 H（X）k ：H（Xk 1）…H（Xk i）…H（Xk n） （2） ＊3 K は実数体 R や有理数体 Q と思ってよい． ＊4 im，ker は線形写像の像と核を表す．として定められる．線形写像が包含関係から誘導されて いることから，直感的には k 次元の穴がフィルトレー ション内でどの程度存続しているかが記録されているで あろう．実際，パーシステントホモロジー H（X）に対k して bi, di∈ {1, …, n}（bi di）および非負整数 lkが一意 に定まり H（X）k lk i＝1I[bi, di] （3） の形に分解できることが知られている．ここで区間 I[b, d] は 0 および一次元ベクトル空間 K の n 個の列とそれら の間の線形写像として I[b, d]：0…K…K…0 で与えられ，左端の K は b 番目，右端の K は d 番目に 位置している．ここで中央に連続して現れる KK は恒 等写像で与えられる．そのトポロジカルな意味は，Xbで 発生し Xdまで存続しその後消滅する k 次元の穴を表し ている． 一意な分解（3）より，発生と消滅の組（bi, di）を指 定することでパーシステントホモロジーは同型の差を除 いて決まる．そこでこの組を平面上の多重集合 D（X）＝{（bk i, di）｜i＝1, …, lk} として表示したものを，X の k 次パーシステント図と呼 ぶ． 定義より，パーシステント図の点はすべて対角線より 上側に位置する．対角線付近の点（bi, di）は存続時間 li＝di−biが短く，一方で，対角線から離れた点は存続 時間が長い穴に対応する．このように，パーシステント 図を描かせることで，視覚的に穴の存続性を理解するこ とができる．例えば，図 1 のデータから Čech 複体のフィ ルトレーションを構成し，その一次パーシステント図を 描かせると，中央の図に存在する右下の小さい穴は対角 線付近に現れ，一方で右図にある大きな穴は対角線から 離れたところに位置することになる（図 3 参照）．パー システント図に現れる発生・消滅対（b, d）（もしくは（3） に現れる I[b, d]）をパーシステント図の生成元ともいう． パーシステントホモロジーにはいくつかの導入の仕方 があるが，前章と同様に境界作用素を用いた定義や，行 列の基本変形を用いた表示（3）の導出については，例 えば著者の [平岡 13] もしくは [平岡 15] を参照されたい．
5．シリカガラスの
構造解析
本章では著者らが最近出版した論文 [Hiraoka 16] で 議論した，パーシステントホモロジーを用いたガラスの 構造解析に関する成果を部分的に紹介する． よく知られているように，ガラス状態の原子配置は一 見無秩序に配置されているように見える．しかしながら 動径分布関数や散乱パターンを見ると，液体状態とは近 距離や中距離オーダで違いが見られる．このようなガラ スの幾何構造を明確に特徴付ける新たな記述子の開発は 科学・産業の両側面から強く期待されている． そこでガラスの原子配置を空間点データ（1）とみな しパーシステント図を描かせてみよう．図 4 は分子動力 学シミュレーションで計算したシリカの結晶，ガラス， 液体状態の原子配置に対して一次パーシステント図を計 算したものである（シミュレーションを含む計算の詳細 は原論文 [Hiraoka 16] を参照されたい）．一次のパーシ ステント図ということでリング構造に着目していること に注意しておく． まず結晶状態（左）では実現できるリングの形が限ら れていることから，パーシステント図はいくつかの孤立 図 3 図 1 の一次パーシステント図 図 4 シリカの一次パーシステント図（左：結晶，中央：ガラス，右：液体）． 色付けは対数重複度を用いている．中央図にはそれぞれの領域で特徴的な原子リング構造を示している（大きい球： 酸素原子，小さい球：ケイ素原子） 図 4 シリカの一次パーシステント図（左：結晶，中央：ガラス，右：液体）． 色付けは対数重複度を用いている．中央図にはそれぞれの領域で特徴的な原子リング構造を示している（赤球：酸素 原子，青球：ケイ素原子）島から構成されている（完全結晶の場合は 0 次元的な点 の集まりに収束する）．また液体状態ではリング形状を ある程度自由に変形できることが反映されて，パーシス テント図が二次元的に散布されている．一方で，ガラス のパーシステント図にのみ曲線（CP, CT, CO）が現れて いることに注意する．いささか乱暴に要約すると，パー システント図のサポートの次元が 0, 1, 2 に応じて結晶， ガラス，液体の 3 状態が識別されていることになる．つ まりパーシステント図のサポート次元はガラス転移を幾 何学的に説明する「秩序変数」として振る舞っているよ うにも見える＊5． さてここまでは与えられた原子配置に対してパーシ ステント図を求めたが，得られたパーシステント図の生 成元を指定した際，その生成元を構成している原子の集 まりを再現することが可能である．これはパーシステン トホモロジーの生成元についての逆問題とみなすことが でき，これによりパーシステント図の中で興味ある特徴 （例えば曲線など）が得られた際に，具体的にその特徴 を生み出す原子配置を特定することが可能になる．ここ で生成元は商ベクトル空間の元を用いて記述されている ことから，一般には指定された生成元を与える点データ は複数存在するが，適切な最適化問題の設定を導入する ことで逆を一意に定めることが可能となる．詳細は論文 [Escolar 15]を参照されたい． 次にガラスのパーシステント図に現れる曲線の意味を 考えてみる．まずこれらの曲線は，シミュレーションの パラメータを保ったまま時間発展させた場合，その形状 は不変に保たれる．つまり熱的な揺らぎのもとでパーシ ステント図の生成元は曲線の接線方向には動けるが，曲 線の法線方向には動くことができないことが観測され る．ここでパーシステント図の生成元は原子配置のリン グ構造から定められていることを思い出そう．よってこ こで観察される法線方向に動けないという制約条件は， 熱揺らぎのもとでの原子のリング構造の幾何学的な制約 にほかならない．また同様に接線方向はリング形状の許 容される変形を示していることになる．つまりガラスの 硬さに関係する制約条件がパーシステント図に埋め込ま れているとも思える．実際に，論文 [Hiraoka 16] では 加圧シミュレーションのもとでガラスのパーシステント 図に現れている曲線が不変であることが示されている． ここではシリカガラスの構造解析のみを取り扱った が，その他の材料科学への適用例として論文 [Ichinomiya 17, Kimura 18, Saadatfar 17]をあげておく．
6．
中距離秩序の抽出
原子間の距離を用いて二体相関を表す動径分布関数 は，一次元ながら遠方までの秩序情報をわかりやすく示 す．一方，三次元の原子座標から得られるパーシステン ト図は，前述のように 0 次（つながり），一次（リング）， 二次（空洞）の構造的特徴を表すため 2 体以上の情報（形 の情報）が含まれると期待されるが，実際の解析で使用 するにあたっては若干の注意が必要である．動径分布関 数とパーシステント図との対応を見ておくことは，複雑 な系を調べるうえで有用である． パーシステントホモロジーを効率良く計算するため にアルファ複体が用いられることはすでに述べたとおり であるが，その準備としてのボロノイ分割により，特に 粗密の差が小さい系において得られる情報が近距離のも のに偏る副作用をもたらす．面心立方（fcc）構造と六 方最密充填（hcp）構造を例に考えよう．これらはどち らも三角格子を積層して得られる構造であるが積み上げ 方に違いがあり，fcc 構造は ABCABC... の，hcp 構造は ABABAB...の周期性をもつ（図 5（a），（b））．hcp 構造 の方が対称性が低いため，原子間距離を 2.0 とした両構 造に対する動径分布関数は，余分なピークの有無に着目 して両者を明確に区別できる（図 5（c），（d））． 一方，fcc 構造と hcp 構造に対する 0 次，一次，二次 のパーシステント図（PD0, PD1, PD2）は，いずれも両 者で違いを生じない（図 6（a），（b）の 1 段目）．これは， 最近接距離がすべて等しく（PD0），単一サイズの正三 角形のみが含まれ（PD1），充填率が同じであるため正 四面体や正八面体の比率も同じである（PD2）ことから の帰結であり，粉体の結晶化プロセスとその密度依存性 を明らかにした論文 [Saadatfar 17] ではこの性質を積極 ＊5 この考察を正当化することは興味深い課題であると思われる． fcc hcp (a) (b) 1.0 1.5half distance2.0 2.5 fcc 1.0 1.5half distance2.0 2.5 hcp (c) (d) 図 5 （a）面心立方（fcc）構造と（b）六方最密充填（hcp） 構造の模型．薄網で示した正八面体と正四面体 はいずれの構造にも含まれる．（c）fcc 構造と， （d）hcp 構造の動径分布関数．パーシステント 図との比較のため 2 体距離の半分で示している 図 5 （a）面心立方（fcc）構造と（b）六方最密充填（hcp） 構造の模型． 黄色で示した正八面体と正四面体はいずれの構造に も含まれる．（c）fcc 構造と（d）hcp 構造の動径分 布関数．パーシステント図との比較のため 2 体距離 の半分で示している的に生かしている．ではパーシステント図に基づいて両 者を区別することはできないのだろうか． PD0は二つの点がつながるイベントを記録したもので あるため，何もしなければ最近接距離の r＝1.0 のみに 値をもつ．そこで原子をランダムに間引いて原子間の平 均距離を大きくし，遠距離の情報を「回復する」ことを 考える．図 6 は約 5 000 原子の fcc 構造，hcp 構造から ランダムに原子を間引いた系についてのパーシステント 図の変化を示したものである．一つのパーシステント図 に含まれる情報量の減少を補うため，複数の構造を生成 して積算した．50％間引きの PD0 では動径分布関数の 2番目のピーク位置と同じ r＝1.4 付近に新たなピークが 出現するが，両構造はまだ識別できない．一方 PD1 の 情報量は飛躍的に増えており，hcp 構造では動径分布関 数の r＝1.6 付近のピークに対応して birth time が 1.6 の生成元が現れている．fcc 構造と hcp 構造の違いは明 らかである．death time（縦方向）の広がりは，辺の数 が多く等方的な多角形の出現を意味する．PD2 でも生成 元の増加が見られるが，この多様な多角形が生み出す空 洞構造はさらに複雑になるため，fcc と hcp の識別性は PD1に劣る．さらに多くの原子を間引けば遠方まで見え るようになり，PD0 に基づいた fcc 構造と hcp 構造の分 類が可能になる（図 6（a）の 3 段目）が，PD1 や PD2 から局所的な形の情報を取り出すことはむしろ困難で あった．2 体相関以上の情報を活用して解析の幅を広げ るためには，系に応じた間引きを適用する必要がある． 動径分布関数とトポロジカル解析を相補的に使うと， 適切な間引き率を見つけやすくなる．図 7 は高温条件下 でアモルファス炭素からグラファイト（黒鉛）が生成す る様子を第一原理分子動力学法によって加速シミュレー ションしたもので，周期境界条件を課された三次元的な ネットワーク構造から層状構造への変化が見て取れる （図 7（a），（b））．しかし動径分布関数の変化は小さく， 何らかの秩序が発達した様子がうかがえるが，それ以上 の解析は難しい（図 7（e））．これに対して，時間発展 に伴う PD1 の変化を差分表示すると，動径分布関数の 第 1 ピーク位置に対応する birth time が 0.075 nm 付近 1.0 1.5 2.0 2.5 birth 1.0 1.5 2.0 2.5 deat h 1.0 1.5 2.0 2.5 birth 1.0 1.5 2.0 2.5 death 1.0 1.5 2.0 2.5 birth 1.0 1.5 2.0 2.5 deat h 1.0 1.5 2.0 2.5 birth 1.0 1.5 2.0 2.5 death 1.0 1.5 2.0 2.5 birth 1.0 1.5 2.0 2.5 deat h 1.0 1.5 2.0 2.5 birth 1.0 1.5 2.0 2.5 deat h 1.0 1.5 death2.0 2.5 1.0 1.5 death2.0 2.5 1.0 1.5 death2.0 2.5 fcc/hcp 100% PD1 (a) (b) 1.0 1.5death2.0 2.5 1.0 1.5death2.0 2.5 1.0 1.5 death2.0 2.5 fcc/hcp 100% PD2 hcp 50% PD1 hcp 50%PD2 fcc 50% PD1 fcc 50%PD2 fcc 100% PD0 hcp 100%PD0 fcc 50% PD0 hcp 50%PD0 fcc 25% PD0 hcp 25%PD0 図 6 （a）各間引き率での 0 次のパーシステント図． （b）一次（左）および二次（右）のパーシステント図 (a) (b) (e) (f) 10 20 30 40 50 Time [ps] 6-membered rings inter-layer order 1 2 3 4 0.05 0.10 0.15 0.20 0.25 [nm] 0 ps 50 ps RDF PD1 (100%) PD1 (50%) (d) (c) 0 ps 50 ps 図 7 アモルファス炭素の分子動力学シミュレーションにおける （a）初期と（b）後期のスナップショット． （c）間引きなし，および（d）半数間引きでの一次の PD の 差分．（e）初期と後期の動径分布関数（RDF）．パーシステ ント図（PD）との比較のため横軸は 2 体距離の半分で表示． （f）六員環と層状秩序に対応する生成元の数の時間変化
（炭素間結合距離の半分）のところに大きな変化が見ら れる（図 7（c））．生成元が増えている濃い領域はよくま とまっており，リング構造の形や大きさが収束した様子 がわかる．次に，間引き率を変えながら炭素原子をラン ダムに取り除いていくと，約 50％のあたりで動径分布 関数の第 2，第 3 ピークのあたりに対応する PD1 の変 化を捉えることができた（図 7（d））．間引きをするこ とで，近距離の情報量は減っているが，代わりに 0.2 ∼ 0.3 nm程度の birth time をもつ生成元の減少（薄い網） と 0.35 nm 程度の birth time をもつ生成元の増加（濃 い網）が見られ，逆問題解析によってそれが層状構造の 成長に由来することが確かめられた．この解析を小さな 時間刻みで行うことで，炭素原子の六員環が増えてから 黒鉛の層状の秩序が発達する様子が定量的に捉えられる （図 7（f））．動径分布関数を参照しつつ間引き率を変化 させて PD の「焦点距離」を調節することで，中距離秩 序の定量的な解析が可能になったといえよう．
7．パーシステント
図の機械学習
上述のようにパーシステント図を直接目で比較して データの特徴的パターンを発見できるというのは幸運な 場合で，実際にはそれではよくわからない場合も多い． そこでどうするかというと，多くのパーシステント図を 用意してそれを統計的に処理することで何らかの特徴的 パターンを発見するのである．本章では [Obayashi 18] で導入された手法について解説する． ここで紹介する機械学習の手法の全体像を図 8 に示 す．一般的な機械学習の手法では，個々の入力データを データの特徴を抽出したようなベクトルに変換し，その ベクトルに対して解析をする（このベクトルは特徴量と 呼ばれる）．パーシステント図をうまくベクトルに変換 することで，入力データの特徴的な幾何パターンを取り 出すことができる．パーシステント図の変換手法はさま ざまな手法が提案されているが，ここでは Persistence Image（PI）と呼ばれるシンプルな手法を利用する．また， 機械学習の手法としては単純な線形モデルを考える．PI と線形モデルの組合せは，そのシンプルさゆえに図 8 の ような逆解析を利用した手法と相性が良く，学習結果を 理解するのに有効である． まず PI について解説する．PI は論文 [Adams 17] で 導入された手法で，基本的なアイディアは「パーシステ ント図のヒストグラムの各ビンの頻度数をベクトルの各 要素とみなす」というものである．ただ，このような安 直な方法だと入力への微小なノイズに対してベクトルが 大きく変化してしまうという欠陥がある．またパーシス テント図の上の各点の重要度（対角線から離れるほど重 要度が高い）が反映されないという欠点もある．これら の問題を解決するため，PI では次のような平面上の分 布関数ρをヒストグラムの代わりに考える． w（b, d）＝arctan（C（d－b）p） ρ（x, y）＝ （bi, di）∈ D（X）k w（bi, di）Nai, b（x, y）i Nai, b（x, y）＝exp －（i x－bi）2＋（y－di）2 2σ2 ⎟⎟⎠ ⎞ ⎜⎜ ⎝ ⎛ （4） ここでσ, C, pはパラメータで，データに合わせて適当 に事前に決める必要がある．これは D（X）の各点に wk の重みの二次元ガウス分布を置いて足し合わせたもので ある．w は対角線との距離によって決まり，対角線に 近くなるほど点の重要度が低くなることを反映してい る．またガウス分布を使うことで実用上はパーシステン ト図の二次元平面をグリッドに切ってこの関数を有限次 元に近似して利用する．すなわちパーシステント図の注 目する領域をグリッドで区切って，各グリッドにおける ρ（x, y）の代表値（グリッドの中心の値やρ（x, y）のグ リッド内での平均値など）を求め，その代表値を一列に 並べたものをベクトルとみなすものである．コンピュー タ上で計算するにはパーシステント図のヒストグラムに 図 8 パーシステントホモロジーと機械学習を組み合わせた手法の概要w（b, d）に対応する重みをかけ，それに二次元のガウシ アンフィルタをかけて得られた「ぼかしたヒストグラム」 の各ビンの値を一列に並べたものをベクトルとみなせば よい．このようにすることでノイズの問題や重要度の問 題に対応する． 次に機械学習について準備する．機械学習，特にここ で用いる教師あり機械学習の目的は （v1, t1）, …,（vn, tn）∈ RN×R という n 個のデータからその背後に隠された v と t の間 の関係 t＝f（v）を発見することである．線形モデルは f に f （v）＝ a·v＋b＋（noise） （5） という形を仮定し，a と b をデータからうまく調整する （すなわち「学習」する）ことで適切な f を見いだす． この二つのアイディアを組み合わせるために，v ∈ RN として PI で得られたベクトルを利用する．パーシステ ント図として抽出された幾何的特徴をベクトルに符号化 して機械学習と組み合わせることで，幾何構造と t との 関連を見いだすのである．式（4）と（5）を組み合わせ， （5）に現れる a をパーシステント図の平面上の符号付き 分布とみなすことで，t とパーシステント図の関係は近 似的に t a（x, y）w（bi, di）Nai, b（x, y）dxdyi ＋b＋（noise） ＝ （bi, di）∈ D（X）k （6） と表されることがわかる． ここで a（x, y）は，PI の離散化のルールを逆に考えて ベクトル a を平面上の分布とみなしたものである．この 式の a（x, y）w（bi, di）Nai, b（x, y）dxdyi の部分はパーシステント図の各点（bi, di）の t の値への 貢献度を表しており，その正負や大小を見ることで t の 値に大きく影響する（bi, di）を特定することができる． bi, diの各値はその幾何構造の特徴的空間スケールに対 応しているため，その値から t と関連の深い重要な空間 スケールを特定したことになる．そして，この重要度が 高い（bi, di）を上述の逆問題の手法で入力データに戻 すことで t の値に重要な形を具体的に知ることもできる のである．材料の問題にこの手法を適用する場合には， 注目している物性値を t として使う．論文 [Kimura 18] では焼結鉱の X 線 CT 画像にこの解析手法を適用し，v として CT 画像中の酸化鉄やカルシウムフェライトと いった各成分の分布形状から計算されたパーシステント 図をベクトル化したもの，t として画像から取り出した クラック（ひび割れ）の大きさを使うことで，クラック と相関の強い形を特定することに成功した．この形は焼 結鉱のき裂発生サイトの特定に有用であると考えられ る． 実際の問題に適用する際には，式（6）のような線形 の関係が成立していると期待できない場合もある．そう いった場合にもロジスティック回帰のような一般化線形 モデルを利用することで問題を解決できる可能性があ る．結局（6）のように各（bi, di）の重要度が計算でき ることが大事なので，そのような性質をもった非線形変 換を PI で得られるベクトルに適用するという手段も考 えられる．また主成分分析のような教師なし機械学習の 利用も可能であり，[Kimura 18] では主成分分析を使っ て焼結鉱の還元反応による変化の特徴付けが行われた． これらの手法や機械学習をする場合の注意点などについ ては，[Obayashi 18] を参照してほしい．そこでは正則 化項による過学習の抑制，l1正則化項を使うことによる スパースな学習，ロジスティック回帰，パーシステント 図の情報に＋αの情報を付け加えることによる性能の向 上，といった問題を議論している．また，パーシステン ト図に対するカーネル法の開発については論文 [Kusano 16]に述べられている．
8．お
わ り に
本解説では，離散点データの形を扱うパーシステント ホモロジーの基本的な性質を概観し，幾何学的特徴に注 目することでこれまで難しかった複雑な系の記述に道が ひらけることを示した．従来の手法と相補的に利用する ことで，より多くの情報を取り出すこともできる．また， 入力データが内在する秩序への気付きを得る順問題だけ でなく，それがどこに起因するのかを調べる逆問題が扱 えることは本手法の大きな強みである．これを生かしつ つ機械学習などと組み合わせることで，その適用範囲は さらに広がるものと期待される． 謝 辞 本 解 説 記 事 は，JST/CREST「 ソ フ ト マ タ ー 記 述 言 語 の 創 造 に 向 け た 位 相 的 デ ー タ 解 析 理 論 の 構 築 （15656429）」，JST/ イノベーションハブ構築支援事業 「情報統合型物質・材料開発イニシアティブ（MI2I）」， JST/戦略的イノベーション創造プログラム（SIP）「革 新的構造材料 D72 ユニット」における研究開発成果に 基づくものである．◇ 参 考 文 献 ◇
[Adams 17] Adams, H., Emerson, T., Kirby, M., Neville, R., Peterson, C., Shipman, P., Chepushtanova, S., Hanson, E., Motta, F. and Ziegelmeier, L.: Persistence images: A stable vector representation of persistent homology, J. Mach. Learn.
Res., Vol. 18, pp. 218-252（2017）
www.cgal.org/
[ E d e l s b r u n n e r 10] E d e l s b r u n n e r, H . a n d H a r e r, J. :
Computational Topology: An Introduction, AMS（2010） [Escolar 15] Escolar, E. G. and Hiraoka, Y.: Optimization in the
real world: Towards solving real-world optimization problems,
Mathematics for Industry, Vol. 13, pp. 79-96, Springer（2015） [平岡 13] 平岡裕章：タンパク質構造とトポロジー : パーシステン
トホモロジー群入門，共立出版（2013）
[平岡 15] 平岡裕章：データに潜む幾何構造∼パーシステントホモ ロジー∼，数理科学，Vol. 53, pp. 48-53，サイエンス社（2015） [Hiraoka 16] Hiraoka, Y., Nakamura, T., Hirata, A., Escolar, E. G., Matsue, K. and Nishiura, Y.: Hierarchical structures of amorphous solids characterized by persistent homology, Proc.
Natl. Acad. Sci. USA, Vol.113, pp. 7035-7040（2016） [平岡 17] 平岡裕章，西浦廉政：ランダムの中に見る秩序─パーシ ステントホモロジーとその応用─，日本物理学会誌，Vol. 72, No. 9, pp. 632-640（2017） [平岡 19] 平岡裕章，大林一平：パーシステントホモロジーの基礎 と材料工学への適用例，日本金属学会会報「まてりあ」，Vol. 58, No. 1, pp. 17-22（2019）
[HomCloud] HomCloud: https://www.wpi-aimr.tohoku. ac.jp/hiraokalabo/homcloud/
[Ichinomiya 17] Ichinomiya, T., Obayashi, I. and Hiraoka, Y.: Persistent homology analysis of craze formation, Phys. Rev. E, Vol. 95, 012504（2017）
[Kaczynski 04] Kaczynski, T., Mischaikow, K. and Mrozek, M.:
Computational Homology, Springer（2004）
[Kimura 18] Kimura, M., Obayashi, I., Takeichi, Y., Murao, R. and Hiraoka, Y.: Non-empirical identification of trigger sites in heterogeneous processes using persistent homology, Sci. Rep., Vol. 8, p. 3553（2018）
[Kusano 16] Kusano, G., Hiraoka, Y. and Fukumizu, K.: Proc.
33rd Int. Conf. on Machine Learning, PMLR, Vol. 48, pp.
2004-2013（2016）
[Obayashi 18] Obayashi, I., Hiraoka, Y. and Kimura, M.: Persistence diagrams with linear machine learning models, J.
Appl. Comput. Topology, Vol. 1, pp. 421-449（2018）
[Saadatfar 17] Saadatfar, M, Takeuchi, H., Robins, V., Francois, N. and Hiraoka, Y.: Pore configuration landscape of granular crystallization, Nat. Commun., Vol. 8, 15082（2017）
2019年 3 月 13 日 受理
