位相的データ解析入門 ― パーシステント ホモロジーの基礎と応用 ―

2.1. 鎖複体とホモロジー

加群の準同型 \(f: A \rightarrow B\) の核 (kernel) とは、\(f(x) = 0\) を満たす \(x \in A\) 全体の集合です。 準同型 \(f\) の核を \(\operatorname{Ker}(f)\) と書きます。 これは \(A\) の部分加群です。 また、\(f(x) \,(x \in A)\) 全体の集合を像 (image) と呼び、\(\operatorname{Im}(f)\) と書きます。 これは \(B\) の部分加群です。
　次のような加群 \(C_n\) と準同型 \(\partial_n\) の列を考えましょう。 $$ \cdots \overset{\partial_{n+1}}{\rightarrow} C_n \overset{\partial_n}{\rightarrow} C_{n-1} \overset{\partial_{n-1}}{\rightarrow} \cdots $$ この列が鎖複体 (chain complex) であるとは、\(\partial_n \circ \partial_{n+1} = 0 \, (^\forall n)\) を満たすことを言います。 この時、\(\operatorname{Im}(\partial_{n+1}) \subset \operatorname{Ker}(\partial_n)\) となります。 すると、商加群 \(H_n = \operatorname{Ker}(\partial_n) / \operatorname{Im}(\partial_{n+1})\) を定義することができます。 これを鎖複体の \(n\) 次ホモロジー群 (homology group) と呼びます。 上記の鎖複体を \(C\) と書くならば、\(C\) の \(n\) 次ホモロジー群は \(H_n(C)\) などと書くことも多いです。

2.2. ホモロジーの意味

　ホモロジーの直感的な意味について説明しましょう。 結論から言えば、ホモロジーは空間の「穴」に対応しています。
　先ほどの鎖複体について、各 \(C_n\) は \(n\) 次元単体 (つまり、多面体の最小単位) の組み合わせの集まりです。 \(C_n\) の基底が単体に対応し、各ベクトルは複数の単体に重みをつけて組み合わせたものと考えられます。 単体を組み合わせたものですから、それ自体が「多面体のようなもの」と考えても良いかもしれません。 各準同型 \(\partial_n\) は、単体の境界を与えます。 すると、\(\operatorname{Im}(\partial_n)\) は「多面体のようなもの」の境界からなる集合と考えられます。 一方、\(\operatorname{Ker}(\partial_n)\) は、境界を持たない「多面体のようなもの」の集合と考えられます。 例えば、三角錐の境界を考えると、これは境界を持たない多面体ですね。 一般に境界の境界は空集合ですから、\(\partial_n \circ \partial_{n+1} = 0\) を満たす必要があるのです。 境界を持たない「多面体のようなもの」を、サイクル (cycle) と呼びます。
　以上より、ホモロジー群 \(H_n = \operatorname{Ker}(\partial_n) / \operatorname{Im}(\partial_{n+1})\) とは、 別の何かの境界にはならないサイクルを集めたものと考えられます。 各サイクルは、空間の内部に埋め込まれた、それ自体が「穴」を持った図形と考えられます。 それがもしも別の何かの境界になるのであれば、その「穴」は空間の内部で「埋めることができる穴」と考えられます。 一方、別の何かの境界にならないとは、その「穴」は「埋めることができない穴」です。 すなわち、ホモロジーは、空間そのものに空いた「穴」を検知しているわけです。

　鎖複体の具体例としては、単体鎖複体 (simplicial chain complex) や特異鎖複体 (singular chain complex) などがあります。 単体複体 (simplicial complex) とは簡単に言えば多面体のことであり、多面体を構成する各単体の集まりを基底として \(C_n\) を生成することができます。 これが単体鎖複体です。 一方、任意の位相空間に対して定義できるのが特異鎖複体です。 特異鎖複体の定義はここでは述べませんが、単体複体を位相空間とみなせば、単体鎖複体のホモロジーと特異鎖複体のホモロジーは同型です。
　また、これらの鎖複体については、非負整数 \(n\) に対してのみ \(n\) 次ホモロジー群が定義されています。 \(0\) 次ホモロジーについて見てみましょう。 \(0\) 次元単体は頂点の集まりです。 頂点は境界を持ちませんから、すべてがサイクルです。 \(1\) 次元単体は線分ですから、線分の 2 つの端点を同一視することが、\(1\) 次元単体の境界で割ることに対応します。 つまり、\(0\) 次ホモロジーの基底は、(弧状) 連結成分に対応しています。

3.1. パーシステント ホモロジーの定義

　半順序集合 \(P\) で添え字づけられた加群のフィルトレーション (filtration) とは、 各要素 \(p \in P\) に対して加群 \(A(p)\) が対応する族 \( \{ A(p) \}\) であって、 \(p \le q\) ならば \(A(p) \subset A(q)\) となるものを言います。
　同様に、鎖複体のフィルトレーションとは、鎖複体の族 \( \{ C(p) \}\) であって、 \(p \le q\) ならば \(C(p) \subset C(q)\) となるものと定義できます。

　一般に、半順序集合 \(P\) で添え字づけられたパーシステンス加群 (persistence module) とは、 各要素 \(p \in P\) に対して加群 \(A(p)\) が対応し、 関係 \(p \le q\) に準同型 \(A(p) \rightarrow A(q)\) が対応して以下を満たすものです。

• 関係 \(p \le p\) には恒等写像 \(id: A(p) \rightarrow A(p)\) が対応する。
• \(p \le q \le r\) のとき、関係 \(p \le r\) に対応する準同型は、 \(p \le q\) と \(q \le r\) それぞれに対応する準同型の合成に一致する。

　鎖複体のフィルトレーション \( \{ C(p) \}_p\) に対し、各鎖複体ごとにホモロジー群を取れば、パーシステンス加群 \( \{ H_n(C(p)) \}_p\) が得られますね。 このパーシステンス加群こそがパーシステント ホモロジー (persistent homology) です。

3.2. 具体例: チェック複体

　半順序集合 \(P\) で添え字づけられた単体複体のフィルトレーション (filtration) とは、 各要素 \(p \in P\) に対して単体複体 \(X_p\) が対応する族 \( \{ X_p \}\) であって、 \(p \le q\) ならば \(X_p \subset X_q\) となるものを言います。

　例えば、距離空間 \((X, d)\) 上の \(N\) 個の点 \(x_1, x_2, \cdots, x_N\) について考えましょう。 各点 \(x_i\) を中心とする半径 \(r\) の閉球 \(B(x_i, r)\) は、以下のように定義されます。 $$ B(x_i, r) = \{ y \in X \,|\, d(x_i, y) \le r \} $$ ここから、半径 \(r\) を固定して、単体複体 (多面体) を構成してみましょう。 まず、\(0\) 次単体 (頂点) は各点 \(x_1, x_2, \cdots, x_N\) になります。 さらに、\(B(x_i, r) \cap B(x_j, r) \neq \emptyset\) のとき、2 つの頂点 \(x_i, x_j\) の間に \(1\) 次単体 (辺) が存在するものとします。 同様に、 \(\cap_{1 \le k \le p} B(x_{i_k}, r) \neq \emptyset\) のとき、頂点の族 \(\{ x_{i_k} \}_k\) の間に \(p\) 次単体が存在するものとします。 このようにして得られた単体複体を、チェック複体 (Čech complex) と呼びます。

3.3. 具体例: 高さ関数と劣位集合

　半順序集合 \(P\) で添え字づけられた位相空間のフィルトレーション (filtration) とは、 各要素 \(p \in P\) に対して位相空間 \(X_p\) が対応する族 \( \{ X_p \}\) であって、 \(p \le q\) ならば \(X_p \subset X_q\) となるものを言います。

　例えば、位相空間 \(X\) が、アフィン空間 \(\mathbb{R}^n\) の中に埋め込まれていたとしましょう。 \(n\) 次元方向の "高さ" を取り出す関数 \(h(x_1, \cdots, x_n) = x_n\) について、 劣位集合 (sublevel set) \(X_t = h^{-1}((-\infty, t])\) を考えます。 すると、実数体 \(\mathbb{R}\) で添え字づけられた空間のフィルトレーション \(\{ X_t \}_t\) が得られます。 ここで、高さ関数 \(h\) は、実際には連続関数であればなんでもよいことに気づくでしょう。 \(X\) がコンパクトであれば、フィルトレーションはある時刻 \(T\) で止まることが保証されます。 つまり、ある実数 \(T \in \mathbb{R}\) に対し、\(T \le t\) ならば \(X_T = X_t\) となるということです。
　専門的な内容になりますが、特に、コンパクト多様体 \(M\) 上のモース (Morse) 関数 \(h\) を考え、同様に\(M_t = h^{-1}((-\infty, t])\) と置きましょう。 すると、以下のようなフィルトレーションが得られます。 $$ M_{t_0} \subset M_{t_1} \subset M_{t_2} \subset \cdots \subset M_{t_N} $$ ここで、\(M_{t_{i-1}} \subset M_{t_i}\) の間にただ 1 つのモース特異点を超えることとすれば、 上記のフィルトレーションはモース理論による多様体 \(M\) のハンドル分解の構成そのものであることに気づくでしょう。

　位相空間のフィルトレーションが存在すれば、特異鎖複体などを対応させることにより、鎖複体のフィルトレーションを得ることができます。 これにより、やはりパーシステント ホモロジーを考えることができます。

3.4. バーコードとパーシステンス図

　全順序集合 \(P\) の区間 (interval) とは、部分集合 \(I \subset P\) であって、次を満たすものを言います。

• 凸性 (convexitiy): \(p \le q \le r\) について、\(p, r \in I\) ならば、\(q \in I\) である。

　全順序集合 \(P\) で添え字づけられた、パーシステンス加群 \( \{ A(p) \}_p\) に対して、以下の定理が成り立つことが知られています。

定理 ([5] を参照) 各 \(A(p)\) が有限次元ならば、 \( \{ A(p) \}_p\) は区間加群の直和で書ける。

3.5. 古典論とのつながり

　古典的なパーシステント ホモロジーは、自然数全体の集合 \(\mathbb{N}\) で添え字づけられた鎖複体のフィルトレーションに対して定義されていました。 古典的な定義をご存じでない方は本節を飛ばして構いません。 (古典的な定義についての詳細は [6] などが参考になります。) 本節では、これまで概説した定義と、古典的な定義との整合性を確かめてみましょう。

　まずは一般論から始めます。 モノイド \(M\) は、その演算により \(M\) 自身への作用を持ちます。 この作用により、モノイド作用の圏 \(\mathscr{C}_M\) が定まります。 つまり、圏 \(\mathscr{C}_M\) は、\(M\) の各要素を対象とし、射 \(s: x \rightarrow y \) は \(s \cdot x = y\) を満たす要素 \(s\) として定義されています。
　もちろん、圏についても、半順序集合と同様に、鎖複体のフィルトレーションを定義することができます。 つまり、圏 \(\mathscr{C}\) で添え字づけられた鎖複体のフィルトレーションとは、\(\mathscr{C}\) から鎖複体の圏への関手 (で、各準同型が単射のもの) のことです。 以降はモノイド作用の圏 \(\mathscr{C}_M\) で添え字づけられた鎖複体のフィルトレーションについて考えましょう。
　体 \(K\) に対し、モノイド \(M\) が生成するモノイド環 \(K[M]\) を考えます。 次数付き \(K[M]\) 加群とは、\(K[M]\) 加群 \(A\) であって、以下の条件を満たす分解 \(\{ A_x \}\) を備えたものです。

• \(\{ A_x \}\) は、モノイド \(M\) で添え字づけられた、\(K\) ベクトル空間としての直和分解 \(A = \oplus_{x \in M} A_x\) である。
• モノイド \(M\) の各要素 \(s\) に対し、\(s \cdot A_x \subset A_{s \cdot x}\) である。

　具体例として、自然数全体の集合 \(\mathbb{N}\) を加法によりモノイドとみなしましょう。 モノイド作用の圏 \(\mathscr{C}_\mathbb{N}\) は、通常の順序を備えた \(\mathbb{N}\) 自身と同型であることは簡単に確かめられます。 また、モノイド環 \(K[\mathbb{N}]\) は多項式環 \(K[X]\) に同型であることも簡単に確かめられます。 パーシステント ホモロジーは区間加群の直和に分解できるのでした。 \(\mathbb{N}\) の区間は閉区間 \([b, d]\) か、半直線 \([b, \infty)\) のいずれかでしかありません。 閉区間 \([b, d]\) の区間加群に対応する次数付き \(K[X]\) 加群は \((X^b)/(X^d)\) です。 ただし、\((X^b)\) とは、単項式 \(X^b\) が生成する、多項式環 \(K[X]\) のイデアルです。 これは、\(K[X]\) 加群としては \(K[X]\) 自身に同型ですが、次数付き \(K[X]\) 加群としては同型でないことに注意しましょう。 また、半直線 \([b, \infty)\) の区間加群に対応する次数付き \(K[X]\) 加群は \((X^b)\) です。 以上は、[7] において示された次数付き \(K[X]\) 加群の構造定理により、パーシステント ホモロジー群が \((X^b)/(X^d)\) と \((X^b)\) の直和で書けることに対応しています。

4.1. 実装例: 散布図のデータ解析

　ユークリッド空間 \(\mathbb{R}^n\) の中から、有限個の点の集まり \(x_1, x_2, \cdots, x_N\) が選ばれていたとします。 各 \(x_i\) が \(n\) 個の数値からなる何らかのデータを表すものと思えば、それらがユークリッド空間上にプロットされている状況と考えられますね。 ユークリッド空間は、通常のノルム距離を考えれば、距離空間です。 よって、第 3.2 節のようにチェック複体を考えることができます。
　また、ユークリッド空間上ではアルファ複体 (alpha complex) を定義することができ、これはチェック複体と同じホモロジーを定めます。 ここでは定義を述べませんが、構成のアイデアはチェック複体とほぼ同じです。 しかし、アルファ複体の場合、単体複体の次元がユークリッド空間の次元 \(n\) で上から抑えられるため、チェック複体よりも高速に計算することができます。 (チェック複体の場合、データの数 \(N\) を用いて、最大で \(N-1\) 次元の単体が現れることに注意しましょう。)

　実際に、Python から GUDHI ライブラリを使用した実装例を見てみましょう。

        
            import gudhi
            import matplotlib.pyplot as plt
            import numpy as np

            # サンプル サイズを設定します。
            sample_size = 300  # サンプルとしてランダムに取る点の数です。

            # 平面上のランダムな点をとります。
            x = np.random.rand(sample_size)  # x 座標です。
            y = np.random.rand(sample_size)  # y 座標です。

            # 散布図を描画します。
            plt.scatter(x, y, s=1)
            plt.show()

            # アルファ複体を生成します。
            alpha_complex = gudhi.AlphaComplex(points=[[x[i], y[i]] for i in range(sample_size)])

            # アルファ複体から単体木を生成します。
            simplex_tree = alpha_complex.create_simplex_tree()

            # パーシステント図を取得します。
            diagram = simplex_tree.persistence()

            # パーシステント図を表示します。
            gudhi.plot_persistence_diagram(diagram)
            plt.show()
        
    

　次に、ユークリッド空間の代わりに、距離空間 \((X, d)\) の中から、有限個の点の集まり \(x_1, x_2, \cdots, x_N\) が選ばれていたとします。 この時、第 3.2 節で述べた通り、チェック複体を定義することができます。 しかし、チェック複体の計算は非常に重く、多くの場合、実用的ではありません。 GUDHI ライブラリでは、比較的計算の軽いヴィートリス-リップス複体 (Vietoris-Rips complex) が提供されています。 その定義はチェック複体に似ていますが、以下の点が異なります。 チェック複体の場合、\(\cap_{1 \le k \le p} B(x_{i_k}, r) \neq \emptyset\) のとき、頂点の族 \(\{ x_{i_k} \}_k\) の間に \(p\) 次単体が存在するものとしました。 ヴィートリス-リップス複体では、\(B(x_{i_k}, r) \cap B(x_{i_l}, r) \neq \emptyset (^\forall k, l)\) のとき、頂点の族 \(\{ x_{i_k} \}_k\) の間に \(p\) 次単体が存在するものとします。 それ以外はチェック複体と全く同じです。 しかし、ヴィートリス-リップス複体のホモロジーは、チェック複体のホモロジーと同型であるとは限りません。

　ヴィートリス-リップス複体の実装例も見てみましょう。 ヴィートリス-リップス複体は、距離行列 \(K = (d(x_i, x_j))\) などから生成することもできます。

        
            import gudhi
            import matplotlib.pyplot as plt
            import numpy as np

            # サンプル サイズを設定します。
            sample_size = 300  # サンプルとしてランダムに取る点の数です。

            # 平面上のランダムな点をとります。
            x = np.random.rand(sample_size)  # x 座標です。
            y = np.random.rand(sample_size)  # y 座標です。

            # 散布図を描画します。
            plt.scatter(x, y, s=1)
            plt.show()

            # 距離行列を生成します。
            scatter = np.array([[x[i], y[i]] for i in range(sample_size)])
            K = np.exp(-scatter @ scatter.T)

            # ヴィートリス-リップス複体を生成します。
            rips_complex = gudhi.RipsComplex(distance_matrix=K)

            # ヴィートリス-リップス複体から単体木を生成します。
            simplex_tree = rips_complex.create_simplex_tree()

            # パーシステント図を取得します。
            diagram = simplex_tree.persistence()

            # パーシステント図を表示します。
            gudhi.plot_persistence_diagram(diagram)
            plt.show()
        
    

4.2. 実装例: 形状解析

　位相的データ解析は、データの位相に注目した解析手法です。 位相に注目するとは、距離や角度、曲率、体積のような特徴量を忘れるということです。 一見、形状解析との相性は悪そうに思えます。 しかし、パーシステント ホモロジーは、位相的情報であるホモロジーの他にフィルトレーションの添え字を持っています。 この添え字方向の特徴量については、実は十分豊富な情報が保存されているのです。 むしろ、特徴量を「抽出する」という点においては、不要な情報を上手く忘れられる位相的手法こそが有効に働くこともあります。

　例えば、次の 2 変数関数 \(f\) について考えてみましょう。 $$ f(x, y) = (\sin(2x) + 2\cos(3y))(2\sin(5x) + \cos(7y)) $$ \(z = f(x, y)\) のグラフを 3 次元空間上にプロットすると、以下のようにデコボコした曲面が得られます。

　実際に計算してみましょう。 今回の場合、単体木の代わりに立方体複体 (cubical complex) クラスを使用するのが便利です。

        
            import math

            import gudhi
            import matplotlib.pyplot as plt
            import numpy as np

            # 描画範囲と精度を設定します。
            minimum = 0  # 描画範囲の最小値です。
            maximum = np.pi  # 描画範囲の最大値です。
            grid_num = 200  # 軸ごとに等間隔に取る点の数です。


            # 解析対象のグラフを与える関数を定義します。
            def f(x, y):
                return (math.sin(2 * x) + 2 * math.cos(3 * y)) * (
                    2 * math.sin(5 * x) + math.cos(7 * y)
                )


            # 等間隔に xy 座標を取得します。
            x = np.linspace(minimum, maximum, grid_num)  # x 座標です。
            y = np.linspace(minimum, maximum, grid_num)  # y 座標です。

            # メッシュグリッドを作成します。
            X, Y = np.meshgrid(x, y)

            # 関数をベクトル化します。
            vectorized_f = np.vectorize(f)

            # 関数を適用して 2 次元配列を得ます。
            z = vectorized_f(X, Y)

            # 立方体複体を生成します。
            cubical_complex = gudhi.CubicalComplex(top_dimensional_cells=z)

            # パーシステント図を取得します。
            diagram = cubical_complex.persistence(min_persistence=-4)

            # パーシステント図を表示します。
            gudhi.plot_persistence_diagram(diagram)
            plt.show()

            # パーシステント図をリスト表示します。
            print(diagram)
        
    

　一般に、ある有限単体複体 \(X\) が、3 次元空間 \(\mathbb{R}^3\) の中に区分的に線形に埋め込まれているとします。 つまり、\(X\) は「ポリゴン」です。 ゲームや映画で使われる CG だと思ってください。 このとき、3 次元単位ベクトル v を 1 つ固定するごとに、v 方向の "高さ関数" が定義できます。 よって、3 次元単位ベクトル v を 1 つ固定するごとに、v 方向のパーシステンス図が得られます。 実は、全ての方向のパーシステンス図の情報を保存しておけば、そこから元のポリゴンを完全に復元できることが知られています ([8, 9] を参照)。 このように、パーシステンス図は、形状に関する非常に多くの情報を持っているのです。

距離の話

超局所層理論の話

ライブラリの話