ホモロジー群を計算しよう

最近、トポロジーを用いた位相幾何学がいろいろな分野で応用されていると思います(センサーネットワーク、ロボティクス、物理学、DNAや高分子など)。

https://mathsoc.jp/publication/tushin/1604/1604tamaki.pdf

トポロジーというと例えばドーナツの穴の数$g$が集合の性質を特徴づけるといったことが有名だと思いますが、こういったものを決めるのがベッチ数という数です。これはホモロジー群の次元で決まるため、ホモロジー群を計算することは位相幾何学において重要でしょう。


まずは基本概念についておさらいします。

# 単体

複数の点を頂点としてできる図形を単体といいます。単体には向きがあります。また、単体に含まれる点の数を次元と言います。$r$次元の単体を$r$単体と呼ぶこともあります。

例えば$(p_0p_1)$は点$p_0$から点$p_1$へ向かう線分です。

![線分.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/190a0ae7-ed59-f67b-b95f-75740e4fd472.png)

$-$をつけると向きが逆になり、
$(p_1p_0) = -(p_0p_1)$
となります。(これが負符号の定義になります。)
ただし、1点（＝0単体）$p_i$については$-p_i = p_i$と定義します。
２次元以上の場合も含めると、一般的には

```math
\begin{align}
P & = \begin{pmatrix}
0& 1 & \cdots & n \\
i_0 & i_1 & \cdots & i_n
\end{pmatrix} \\
(p_0p_1\cdots p_n) &= \mathrm{sgn}(P)(p_{i_0}p_{i_1}\cdots p_{i_n})
\end{align}
```

と定義されます。

$(p_0p_1p_2)$は点$p_0$、$p_1$、$p_2$を頂点とする三角形でこの順に３点を回る方向に向きがあります。

![三角.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/1abafe19-696b-eb86-a956-fbcca142f7c4.png)

数学的に表すなら

```math
(p_0 p_1\cdots p_n) = \left\{ \sum_i a_i \overrightarrow{p_0p_i} | 0\leq a_i \leq 1, \sum_i a_i \leq 1\right\}
```

に向きがついているものと言えます。

さらに、"何もない状態"というものも考えてこれを$0$と表します。

次にある単体$(p_0p_1\cdots p_n)$からいくつかの点を取り除いたもの$(p_{i_0}p_{i_1}\cdots p_{i_l})$からなる集合$K$を単体的複体といいます。

例えば

```math
K_2 = \{0,p_0,p_1,p_2,(p_0p_1),(p_1p_2),(p_2p_0),(p_0p_1p_2)\}
```

などです。(普通0は要素に含めませんが、ここでは分かりやすくするために要素としてカウントしています。)

様々な図形をポリゴン化して(数学的には同相写像による変換という)、その単体的複体を考えることで様々な図形のホモロジーを計算できます。例えば、先の三角形の単体的複体$K_2$は円盤を同相変換したものになります。

# 鎖群
単体的複体$K$の要素に対して、整数を係数として同じ次元の単体を線形結合したものを鎖といいます。

例えば$K_2$に対して

```math
2(p_0p_1) + (p_1p_2) -3(p_2p_0)
```

などです。
$r$次元の鎖全体の集合を鎖群といい、

```math
C_r(K) = \left\{ \sum_i c_i \sigma_i^{(r)} \mid \{\sigma_i^{(r)}\} \subset K, c_i \in \mathbb{Z} \right\}
```

と表します。鎖群はベクトルのように計算することができます(加群という)。ただし、整数が係数なので係数に分数をかけることができません(つまり割り算ができません)。

## 境界作用素
 
境界作用素

```math
\partial_r:C_r(K) \rightarrow C_{r-1}(K)
```

を$\sigma_r=(p_0p_1\cdots p_r) \in C_r(K)$に対して

```math
\partial_r \sigma_r = \sum_{i=0}^{r}(-1)^i(p_0p_1\cdots \hat{p}_i \cdots p_r) \qquad(\hat{p}_iは取り除いた頂点)
```

で定義します($\partial_r$の$r$は面倒くさいときは省略します)。これは図形の縁を矢印の向きに足し合わせたものになっています。

境界作用素には

```math
\partial_r \circ\partial_{r+1}=0
```

という性質があります。理由は縁は図形の周りをぐるっと囲む(あるいは表面を覆う)ものになっていて、さらにその縁は始点と終点が一致して打ち消すからです。

例えば$K_2$だと

```math
\begin{align}
\partial(p_0p_1p_2) &= (p_0p_1) + (p_1p_2) + (p_2p_0) \\
\partial(p_0p_1) &= p_0 - p_1 = p_0 + p_1 \\
\partial p_0 &= 0
\end{align} 
```

となります。ここで$(p_2p_0)=-(p_0p_2)$を使いました。

また、鎖はベクトルのように扱えるということを思い出してもらうと、境界作用素は線形作用素に対応していて、鎖を基底として行列で表すことができます。
例えば$C_1(K_2)$の基底を

```math
\{ (p_0p_1), (p_1p_2), (p_2p_0) \}
```

$C_0(K_2)$の基底を

```math
\{ p_0, p_1, p_2 \}
```

として

```math
\partial_1 = \begin{pmatrix}
-1 & 0 & 1 \\
1 & -1 & 0 \\
0 & 1 & -1
\end{pmatrix}
```

と書けます。

# ホモロジー

次の２つの鎖群の部分加群を考えます。

```math
\begin{align}
Z_r(K) &= \{ c \in C_r(K) | \partial c = 0 \} \\
B_r(K) &= \{ c \in C_r(K) | \exists c' \in C_{r+1}(K) \quad s.t. \quad c = \partial c' \}
\end{align}
```

$\partial\partial = 0$の性質があるので、$B_r(K) \subset Z_r(K)$となりますが、
$K$の種類によっては$\tilde{\sigma} \in Z_r(K) $で、$\tilde{\sigma} \notin B_r(K) $となるものが存在します。この様な鎖のなす集合がホモロジー群です。

ただ、$\forall b \in B_r(K)$を足しても、$\tilde{\sigma}+b \in Z_r(K)$なので、このような単体は無数にあることになってしまいます。つまり、$b$はちょうど微積分でいうところの積分定数のようなものになっています。

https://qiita.com/hiro949/items/9021c202dcc10792f214

これでは自由度が残るのでこれを無視します。

数学的に言うと剰余類

```math
[z] = \{ z' \in Z_r(K) | z-z' \in B_r(K) \}
```

を考えて、この剰余類の線形結合全体のなす加群をホモロジー群$H_r(K) = Z_r(K)/B_r(K)$と言います。

![鎖複体.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/3533ceb2-1c4f-ba45-991a-d635f7ff8c34.png)

$K$は図形とそのパーツだけでできているので、基本的には$Z_r(K)$の元は大抵どこかの縁になっていて$B_r(K)$に含まれるはずです。では、どういうときにホモロジー群が生じるかというと、例えば周期的境界条件があるときです。ざっくり言うとホモロジー群の元は図形の周期的境界としてふさわしい面や線のことです。

# ホモロジー群の計算方法

ホモロジー群の基底は行列表示した境界演算子$\partial_r$を特異値分解することで得られます。整数係数の特異値分解はスミス標準化と呼ばれていて、特異値は必ず自然数になります。ここで掃き出し法をやる際に、整数係数ゆえに分数を掛けることができない、つまり割り算ができないことに注意してください。ちなみに実数係数の鎖群に対するホモロジー群は特異ホモロジー群と言い、このときは普通のガウスの掃き出し法で計算できます。

```math
\begin{align}
\partial &= U_{r-1} T_{r-1}(V_r)^{-1} \\
T_{r-1} &= \begin{pmatrix}
t^{r-1}_1 & 0 \cdots & 0 & 0 & \cdots & 0\\
0 & \ddots & 0 & \vdots & \ddots & \vdots \\
\vdots & 0 &t^{r-1}_k & 0 & \cdots & 0 \\
0 & \cdots & 0 & 0 & \cdots & 0 \\
\vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
0 & \cdots & 0 & 0 & \cdots & 0
\end{pmatrix} \\
U_{r-1} &= \{ u^{r-1}_1 u^{r-1}_2 \cdots u^{r-1}_k \cdots \} \\
V_{r} &= \{ v^{r}_1 v^{r}_2 \cdots v^{r}_k \cdots \} \\
\end{align}
```

このとき、$u^{r-1}\_1 u^{r-1}\_2 \cdots u^{r-1}\_k$が$B_{r-1}(K)$の基底ベクトル、

$v^{r}\_{k+1} v^{r}\_{k+2} \cdots v^{r}\_{d_r}$が$Z_r(K)$の基底ベクトルとなっています。ここで$d_r = dimC_r(K)$です。

これは、$T_{r-1}$によって前の$r$個の要素は生き残り、後ろの要素は$T_{r-1}$のゼロベクトルの行によって消えてしまうからです。

よって$v^{r}\_{k+1} v^{r}\_{k+2} \cdots v^{r}\_{d_r}$のうち、$u^{r}\_1 u^{r}\_2 \cdots u^{r}\_k$の線形結合で表せないものの剰余類がホモロジー群の基底になります。

## 注意点
境界演算子を作用させると重み$t^r\_i$がかかります。よって$t^r\_i>1$の場合は
$v^{r}\_i = \sum_i c\_{ij} u^{r}\_j$のように表せたとしても、実際に$\partial$を作用させて出てくるのは

```math
\partial \left( \sum_j c_{ij} u^{r}_j \right) = t^{r-1}_i v^{r-1}_i
```

です。よって$w \not\equiv 0 \mod  t^{r-1}_i$に対してのみ $wv^{r-1}_i$はホモロジー群の元になります。つまり、ホモロジー群をベクトル空間のように直和分解したとき、$[v^{r-1}_i]$を基底とする１次元空間は$\mathbb{Z}\_{t^{r-1}\_i}$に同型になります。この特異値$t$をねじれ率と言います。

# 一般の位相空間のホモロジー

同相写像で単体と同相な集合は同相写像で写すことで同様に扱えます。そのような同相写像が存在しない場合も含めてより広く位相空間を扱うには特異ホモロジー群というものを使う必要があります。

https://qiita.com/hiro949/items/89bb7bcab7a24c2c09ab

# 実装

では以上の計算過程をpython３コードで計算してみましょう。

単体は頂点番号のリスト形式で与えられます。単体的複体はさらにそのリストKで与えられます。
これを入力として、リストK中の$r$次元鎖を$K$の要素の順番で基底としてホモロジーを計算します。

例えば先の三角形の例：


![三角.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/1abafe19-696b-eb86-a956-fbcca142f7c4.png)

つまり、

```math
K = \{0,p_0,p_1,p_2,(p_0p_1),(p_1p_2),(p_2p_0),(p_0p_1p_2)\}
```
ならば、

```input.py
K = [[0],[1],[2],[0,1],[1,2],[2,0],[0,1,2]]
```

と表します。(各点$p$のインデックスをそのままリストにすればよいです。)

今回は11種類の図形に対応する単体的複体を用意しました。
全部やると多すぎるので、特に有名な$K_9$(トーラス)と$K_{10}$(クラインの瓶)をデモしてみたいと思います。

```main.py
import numpy as np
from HomologyGroup import calc_HomologyGroupList

# simplicical-complex inputs for tests
K0 = [[0]]
K1 = [[0],[1]]
K2 = [[0],[1],[0,1]]
# circumference
K3 = [[0],[1],[2],[0,1],[1,2],[2,0]]
# disk
K4 = [[0],[1],[2],[0,1],[1,2],[2,0],[0,1,2]]
# spherical
K5 = [[0],[1],[2],[3],\
      [0,1],[0,2],[0,3],[1,2],[1,3],[2,3],\
      [0,1,2],[0,1,3],[0,2,3],[1,2,3]]
# sphere
K6 = [[0],[1],[2],[3],\
      [0,1],[0,2],[0,3],[1,2],[1,3],[2,3],\
      [0,1,2],[0,1,3],[0,2,3],[1,2,3],[0,1,2,3]]
# Mobius band
K7 = [[0],[1],[2],[3],[4],[5],\
     [0,1],[0,2],[0,5],[1,2],[1,3],[1,4],[1,5],[2,3],[2,4],[3,4],[3,5],[4,5],\
    [0,1,2],[2,1,4],[2,4,3],[3,4,5],[3,5,1],[5,0,1]]
# projective space
K8 = [[0],[1],[2],[3],[4],[5],\
     [0,1],[0,2],[0,3],[0,4],[0,5],[1,2],[1,3],[1,4],[1,5],[2,3],[2,4],[2,5],[3,4],[3,5],[4,5],\
    [0,1,2],[0,2,3],[0,3,5],[0,4,1],[0,5,4],[1,3,5],[1,4,3],[1,5,2],[2,4,3],[2,5,4]]
# torus
K9 = [[0],[1],[2],[3],[4],[5],[6],[7],[8],\
    [0,1],[0,2],[0,3],[0,4],[0,6],[0,7],[1,2],[1,4],[1,5],[1,7],[1,8],\
    [2,3],[2,5],[2,6],[2,8],[3,4],[3,5],[3,7],[3,8],\
    [4,5],[4,6],[4,8],[5,6],[5,7],[6,7],[6,8],[7,8],\
    [0,1,4],[1,4,5],[1,2,5],[2,5,6],[0,6,2],[0,6,4],\
    [3,4,5],[3,7,5],[5,6,7],[6,7,8],[4,8,6],[3,4,8],\
    [0,3,7],[0,1,7],[1,7,8],[1,2,8],[2,8,3],[0,3,2]]
# Klein bottle
K10 = [[0],[1],[2],[3],[4],[5],[6],[7],[8],\
    [0,1],[0,2],[0,3],[0,4],[0,6],[0,7],[1,2],[1,5],[1,6],[1,7],[1,8],\
    [2,3],[2,4],[2,5],[2,8],[3,4],[3,5],[3,7],[3,8],\
    [4,5],[4,6],[4,8],[5,6],[5,7],[6,7],[6,8],[7,8],\
    [0,4,2],[2,4,5],[1,2,5],[1,5,6],[0,1,6],[0,6,4],\
    [3,5,4],[3,7,5],[5,7,6],[6,7,8],[4,6,8],[3,4,8],\
    [0,7,3],[0,1,7],[1,8,7],[1,2,8],[2,3,8],[0,3,2]]

H = calc_HomologyGroupList(K9)
betti = np.array([ len(hi) for hi in H ])
Euler = np.sum( np.array([ pow(-1,i) for i in range(len(betti)) ])*betti )

for i, hi in enumerate(H):
    print("H{0} = {1}\ntor={2}".format(i,hi[0].T,hi[1]))

print(Euler)
```

関数`calc_HomologyGroupList`は単体的複体の２次元リスト（ノードに番号を振ってその番号を羅列したもののリスト）を与えると、各次元のホモロジー群の基底ベクトルとねじれ率のセットを与えます。また、$b_r = dim(H_r(K))$を$r$次元のベッチ数と言い、

```math
\chi = \sum_i (-1)^i b_i
```

をオイラー数と言います。オイラー数と穴の数$g$には

```math
\chi = 2-2g
```

の関係があります。ただし、これは向き付け可能な場合(ねじれのない場合)です。向き付け不可能な場合は別の関係式が成り立ちます。

`calc_HomologyGroupList`を与えるコードは次のようになります。

なお、スミス標準化を与えるコード`smith_normalize.py`の中身は

https://qiita.com/yuji0001/items/64dc97cd4dcebf83d0a8

で紹介されているものを使いました。

```HomologyGroup.py
from math import gcd
import numpy as np
import numpy.linalg as LA
from sympy.combinatorics import Permutation
import smith_normalize as smith

#リスト形式の1単体の符号の判定
def perm_sgn(l1,l2):
    return Permutation([ l2.index(i) for i in l1 ]).signature()

#リスト形式の単体に境界作用素を作用させる。
def simplex_boudary(s,C1):
    chain = np.zeros(len(C1))
    for i in range(len(s)):
        sb = s.copy()
        sb.pop(i)
        isb = [ set(s) for s in C1 ].index(set(sb))
        chain[isb] = pow(-1,i)*perm_sgn(C1[isb],sb)
    return chain

#i次元の境界作用素の行列表示    
def calc_ith_boundary(C1,C2):
    D = np.zeros( (len(C2),len(C1)) ).astype(int)
    for j,s in enumerate(C2):
        D[j] = simplex_boudary(s,C1)
    return D.T

def calc_boundary_operators(SimplicicalComplex):
    size = list(set([ len(s) for s in SimplicicalComplex ]))
    ChainGroupList = [ [ s for s in SimplicicalComplex if len(s) == size[i] ] for i in range(len(size)) ]
    BoundaryList = [np.zeros( (1,len(ChainGroupList[0])) ).astype(int)]
    for i in range(len(size)-1):
        BoundaryList.append( calc_ith_boundary(ChainGroupList[i],ChainGroupList[i+1]) )
    BoundaryList.append(np.array([]))
    return BoundaryList


#Bの基底を掃き出し法で階段行列形式に直す
def colR(M):
    if(M.shape[0]*M.shape[1]==0): return M
    for j in range(M.shape[1]):
        i = list(map(lambda a:np.where(a!=0)[0][0],M.T))[j]
        for j2 in range(M.shape[1]):
            if(j2==j): continue
            c = np.sign(M[i,j2])*(abs(M[i,j2])//M[i,j])
            M[:,j2] -= c*M[:,j]
    return M

def colR_same_tor(M,tortion):
    tor = tortion.copy()
    if(M.shape[0]*M.shape[1]==0): return M
    for j in range(M.shape[1]):
        i = list(map(lambda a:np.where(a!=0)[0][0],M.T))[j]
        for j2 in range(M.shape[1]):
            gt = gcd(tor[j],tor[j2])
            # avoid mixing non-trivial tortion
            if(j2==j or (tor[j2]>tor[j]) or (tor[j]>1 and tor[j2]>1 and gt==1) ): continue
            c = np.sign(M[i,j2])*(abs(M[i,j2])//M[i,j])
            tor[j2] = gt
            M[:,j2] -= c*M[:,j]
    return M,tor

#ホモロジー群の基底の計算
def calc_cohomology(D):
    r = LA.matrix_rank(D)
    invU,S,V = smith.Smith_Normalization(D)
    U = LA.inv(invU)
    Z = V[:,r:]
    tortion = np.diag(S[:r,:r])
    B = U[:,:r]
    Z = colR(Z)
    B, tortion = colR_same_tor(B,tortion)
    return B,tortion,Z

def calc_ith_homology(D1,D2):
    if( (D1==0).all() ): Z1 = np.identity(D1.shape[1]).astype(int)
    else: B0, tortion0, Z1 = calc_cohomology(D1)
    if( len(D2)==0 ):
        B1 = np.array([[]])
        tortion1 = np.array([])
    else: B1,tortion1, Z2 = calc_cohomology(D2)
    Z_nonzero = list(map(lambda z:np.where(z!=0)[0][0],Z1.T))
    B_nonzero = list(map(lambda z:np.where(z!=0)[0][0],B1.T))
    non_trivial = [ [i,1] for i,z in enumerate(Z_nonzero) if not(z in B_nonzero) ]
    non_trivial += [ [np.where(Z_nonzero == ib )[0][0],tortion1[i]] for i,ib in enumerate(B_nonzero) if tortion1[i]>1 ]
    return Z1[:,[ x[0] for x in non_trivial ]],[ x[1] for x in non_trivial ]

def calc_HomologyGroupList(SimplicicalComplex):
    BoundaryList = calc_boundary_operators(SimplicicalComplex)
    HomologyGroupList = []
    for i in range( len(BoundaryList)-1 ):
        HomologyGroupList.append( calc_ith_homology(BoundaryList[i],BoundaryList[i+1]) )
    return HomologyGroupList
```

### 出力の見方

出力は各次元のホモロジー群の基底となるホモロジー類の代表元（２次元リスト）とそのねじれ率（１次元リスト）の組がタプル型で出力されます。i次元の代表元のリストは入力として与えた単体的複体のリストのうちi次元単体であるものを入力のリストの順番に並べたとき、その係数のリストとなっています。


# 計算例
いくつか例をやってみましょう。

## トーラス

トーラスを同相変換すると、次の右図のようになります。
この図は展開図のようになっていますが、同じ番号の点が複数あって、これらの部分は周期的境界条件を考えています。また、各単体の向きも共有している辺などの境界線で向きが食い違わないようにとる必要があります。(そうならないように十分細かくメッシュをとる必要があって、例えば次図はトーラスにおいてそうならないような最小のメッシュの取り方になっています。)

![トーラス.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/205d4da5-225d-334b-bbaa-f4945a9bbb8a.png)

コードを実行すると出力は以下のようになります：

![トーラスホモロジー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/5f948ed6-85ae-48d2-3937-6066ccb42816.png)

ねじれ率はすべて１なので、係数は整数であり、出力を鎖の形で書くと上記のようになります。
また、図の青、緑線は$H_1$の基底であり、縦横に１本ずつ出ています。

## クラインの壺

ねじれがある例としてクラインの壺も観てみましょう。

![クライン.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/76f47ca7-0c78-291a-f53e-1038a1f3040d.png)

結果は以下の通りです。

![クラインホモロジー.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2445529/7d9ef10c-6507-bb2d-4457-aabb4bf9c5b3.png)

図のピンク、緑の線は１次元ホモロジー群の基底で、
ピンクの基底はねじれ率が２だから$\mathbb{Z}_2$係数となっています。
これがどういうことを意味しているか見てみましょう。
ピンクの線は適当な$B_1(K)$の元を付け足して、オレンジの線に書き直すことができます。実はこの線は周期的境界条件のために水色の線と同じものになっています。したがって、このピンクの線は１本ならホモロジー群の元ですが、２本つなぐと始点と終点が一致してしまい、$B_1$の元になってしまいます。

# 参考
- サンプルコード：https://github.com/hiro949/homology
- 河田敬義 位相幾何学 (現代数学演習叢書 2)
- 中原　幹夫 著・訳 佐久間　一浩 訳　理論物理学のための幾何学とトポロジー１ [原著第2版］