モチベーション

最近Computer Vision関連の領域で、回転・並進等価性を維持しながら畳み込むことのできるNeural Networkが多く発案されている。

グラフやボクセルのこのようなネットワークについては今まで有識者の方々が図を使って親切に解説してくださっている記事があったので、なんとなくは理解できていた。

今回このシリーズの3D点群バージョン(Tensor field networks: Rotation- and translation-equivariant neural networks for 3D point clouds)を自力で読もうとしたところ、球面調和関数とかいうガッチガチの数学用語が沢山出てきて全く意味が分からなかった。

そこで、Equivariant Neural Networksの理解を深めるためにも、まず球面調和関数について勉強することとした。

球面調和関数の定義(Wikiから引用)

一般的に球面調和関数は以下のいずれかを示すとされている。

$n$次元ラプラス方程式の解となる斉次多項式を単位球面で制限することで得られる関数
1.の$n=3$の場合、もしくは$(r,\theta,\phi)$で書いたラプラス方程式の変数分離解を記述するのに用いることができる関数$Y^n_k(\theta,\phi)$

ここで$n$次元ラプラス方程式の解とは$u=u(x_1,...,x_n)$に対して

$\Delta u = \sum \frac{\partial^2 u}{\partial x_i^2} = 0$

を満たすものであり、斉次多項式（あるいは同次多項式）は非零の項の次数がすべておなじである多項式のことである。(例はWikiを参照)

ここで、より細かい定義を見ていく。

$(x_1,...,x_n) \in \mathbb{R}^n$とし、$\phi:\mathbb{R}^n \rightarrow \mathbb{C}$が二階微分可能であるとする。

$\Delta \phi = \sum_i \frac{\partial^2}{\partial x_i^2} \phi$としたとき、$\Delta$をラプラス作用素という。

$\mathbb{R}^n$上の多項式$p(x_1,...,x_n)$で$\Delta p = 0$を満たす$p$ （つまり$n$次ラプラス方程式の解になるもの）を調和多項式という。

このとき、ラプラス作用素は回転行列$R$について$\Delta p (R(x)) = R (\Delta p (x))$を満たす。（ラプラシアンの回転不変性）
ので、調和多項式の定義は座標の取り方に依存しない。

さらに、$p$が斉次多項式であるとき、$p$を単位球面

$S^{n-1} = \{(x_1,...,x_n) \in \mathbb{R}^n | \sum_i x_i^2 =1 \}$

に制限した写像（つまり、入力を単位球面上の$(x_1,...,x_n)$のみとした関数・写像）をk次の球面調和関数という。

ラプラシアンの回転不変性の証明

$\mathbf{x}=(x_1,...,x_n)^T$としたとき、$\mathbf{y} = R\mathbf{x} + \mathbf{z}$（$R$は回転行列、$\mathbf{z}$は平行移動に対応するベクトル）について、

$\Delta u = \sum_i \frac{\partial^2}{\partial x_i^2} u = \sum_j \frac{\partial^2}{\partial y_j^2} u$

であることを示す。

$u_{x_i} = \sum_j \left\{ \frac{\partial u}{\partial y_j} \frac{\partial y_j}{\partial x_i} \right\}$

ここで$\frac{\partial u}{\partial y_j} = u_{y_j}$と表すことができ、$\frac{\partial y_j}{\partial x_i}$は回転行列$R$の$ji$成分である。

$u_{x_i} = \sum_j u_{y_j} R_{ji}$

同様に二階微分も求める。

$u_{x_i x_i} = \sum_j \left\{ \frac{\partial u_{x_i}}{\partial y_j} \frac{\partial y_j}{\partial x_i} \right\} = \sum_j \left\{ \frac{\partial}{\partial y_j}(\sum_k u_{y_k} R_{ki}) R_{ji} \right\} = \sum_j R_{ji}^2 u_{y_j y_j}$

したがってこれを足し合わせると

$\sum_i u_{x_i x_i} = \sum_i \sum_j R_{ji}^2 u_{y_j y_j} = \sum_j (\sum_i R_{ji}^2) u_{y_j y_j}$

回転行列$R$は直行行列であるため各列のノルムは1である。したがって$\sum_i R_{ji}^2=1$が成立するのでラプラシアンの回転不変性が示された。

球面調和関数の例

以下は3次元の球面調和関数について考える。$\mathbb{R}^3$を極座標$(r,\theta,\phi)$で表すと、以下が球面調和関数になることが知られている。

$Y^m_k (\theta,\phi) = (-1)^{(m + |m|)/2} \sqrt{\frac{2k+1}{4\pi} \frac{(k-|m|)!}{(k+|m|)!}} P^{|m|}_k (cos \theta) e^{i m \phi}$

ここで$P^m_k$はルジャンドルの倍多項式と呼ばれるもので、ルジャンドルの倍微分方程式の解である。

ルジャンドルの倍多項式：$P^m_k(t) = \frac{1}{2^k}(1-t^2)^{m/2} {\displaystyle \Sigma_{j=0}^{\lfloor (k-m)/2 \rfloor}} \frac{(-1)^j (2k-2j)!}{j! (k-j)! (k-2j-m)!} t^{k-2j-m} (k\in \mathbb{Z} ,m\in \mathbb{Z} , k \geq m \geq 0 )$
ルジャンドルの倍微分方程式：$(1-t^2)y''(t) - 2ty'(t) + \left\{k(k+1) - \frac{m^2}{1-t^2} \right\}y(t) = 0$

この式は前述の3D点群バージョンの論文に沢山出てきていたので、この式が理解できればある程度は読み進められそうだ。

$Y^m_k$を用いることで、$p$における3次ラプラス方程式の解を変数分離形で与えることができる。

$p(r,\theta,\phi) = r^k \left(AY^m_k (\theta,\phi) + BY^{-m}_k (\theta,\phi) \right)$

($Y^m_k$の式の中身をよく見ると$\theta$と$\phi$が既に変数分離されていることがわかる)

証明

$p=p(x,y,z)$が変数分離できると仮定して$p=R(r)\Theta(\theta)\Phi(\phi)$とおくと

$ \Delta p = \frac{\partial^2 p}{\partial r^2} + \frac{\partial^2 p}{\partial \theta^2} + \frac{\partial^2 p}{\partial \phi^2} \
= \frac{1}{r^2} \left\{\Delta_r R(r)\Theta(\theta)\Phi(\phi) + R(r)\Delta_\theta \Theta(\theta)\Phi(\phi) + R(r)\Theta(\theta)\Delta_\phi \Phi(\phi)\right\}$　

ここで$\Delta_r = \frac{\partial}{\partial r}(r^2 \frac{\partial}{\partial r}), \Delta_\theta = \frac{1}{sin \theta} \frac{\partial}{\partial \theta} (sin \theta \frac{\partial}{\partial \theta}),\Delta_\phi = \frac{1}{sin^2 \theta} \frac{\partial^2}{\partial \phi^2}$(cf.極座標ラプラシアン)

したがってラプラス方程式全体を$p/r^2$で割ると$\frac{1}{R}\Delta_r R + \frac{1}{\Theta}\Delta_\theta \Theta + \frac{1}{\Phi}\Delta_\phi \Phi = 0$となる。

$p$が斉次多項式であることを思い出すと、$k$次の斉次多項式には$p(\lambda x_1,...,\lambda x_n)=\lambda^k p(x_1,...,x_n)$となる性質がある(詳しくは斉次多項式のWiki)ので

$p = r^k Y(\theta, \phi)$

という形式で$r$と$(\theta,\phi)$に分解できる。

$R(r)=r^k$とすれば$\frac{1}{R}\Delta_r R = \frac{1}{R(r)} \frac{\partial}{\partial r}(r^2 \frac{\partial}{\partial r} R(r)) = \frac{k(k+1)r^k}{r^k}=k(k+1)$

したがって$\frac{1}{\Theta}\Delta_\theta \Theta + \frac{1}{\Phi}\Delta_\phi \Phi = -k(k+1)$となり、全体に$sin^2 \theta$を掛けると

$\frac{sin \theta}{\Theta}\frac{\partial}{\partial \theta}(sin \theta \frac{\partial \Theta}{\partial \theta})+ \frac{1}{\Phi}\frac{\partial^2 \Phi}{\partial \phi^2} = -k(k+1)sin^2 \theta$

ここで$\frac{1}{\Phi}\frac{\partial^2 \Phi}{\partial \phi^2} = (\phi\text{について一定})$と考えると、$\Phi$は$\Phi'' = -m^2 \Phi$の一般的な解となるので

$\Phi(\phi) = Ae^{i m \phi} + Be^{-i m \phi}$　と表すことができる。($i$は虚数単位、$m$は整数、$A$と$B$は実数)

$\frac{1}{\Phi}\frac{\partial^2 \Phi}{\partial \phi^2}=-m^2$を代入して全体に$\Theta$をかけると

$\left\{k(k+1)sin^2 \theta -m^2 \right\}\Theta + sin \theta \frac{\partial}{\partial \theta}(sin \theta \frac{\partial \Theta}{\partial \theta}) = 0$

$t=cos\theta$と変数変換し$y(t) = \Theta (arccos \theta)$とすると

$\frac{\partial \Theta}{\partial \theta} = \frac{\partial \Theta}{\partial t}\frac{\partial t}{\partial \theta} = -sin\theta y'(t)$

$\frac{\partial}{\partial \theta}(sin \theta \frac{\partial \Theta}{\partial \theta}) = \frac{\partial}{\partial \theta}\left\{(t^2 -1) y'(t)\right\} = 2t(-sin\theta)y'(t) + (t^2-1)(-sin\theta)y''(t)$

この結果を代入して全体を$sin^2 \theta = 1 - t^2$で割ると

$(1-t^2)y''(t) - 2ty'(t) + \left\{k(k+1) - \frac{m^2}{1-t^2} \right\}y(t) = 0$

これはルジャンドルの倍微分方程式である。

したがって$\Theta (\theta) = P^{|m|}_k(cos\theta)$

これで$p=R(r)\Theta(\theta)\Phi(\phi)$が球面調和関数の条件を満たした。

$(-1)^{(m + |m|)/2} \sqrt{\frac{2k+1}{4\pi} \frac{(k-|m|)!}{(k+|m|)!}}$は規格化定数である。

球面調和関数の性質

内積・直交性

球面調和関数のない席は以下のように定義することができる。

$\langle f|g \rangle_{S^{n-1}} := \int_{S^{n-1}} f(\mathbf{x})g(\mathbf{x})dS$

3次元空間の場合、球面座標$(r,\theta,\phi)$に対し

$\mathrm{d}S = sin\theta \mathrm{d}\theta \mathrm{d}\phi$

が成立する。

$k$次球面調和関数全体のなすベクトル空間を$\mathcal{H}_k$とおくと、$k \neq j$に対し$\mathcal{H}_k$と$\mathcal{H}_j$は前述のない席に関して直行する。

$$
\langle Y^m_k | Y^s_j \rangle_{S^{n-1}} = \begin{cases}
1 & (if \quad m=s,k=j) \\
0 & (otherwise)
\end{cases}
$$

対称性

球面調和関数のパリティ変換($\mathbf{r} \rightarrow -\mathbf{r} (\theta \rightarrow \pi - \theta, \phi \rightarrow \phi+\pi)$)について以下が成立する

$Y^m_k (\pi - \theta, \phi + \pi) = (-1)^m Y^m_k (\theta, \phi)$

加法定理

$\mathbf{x} = (1,\theta,\phi),\mathbf{y}=(1,\theta',\phi')$について内積を計算すると、$\mathbf{x}$と$\mathbf{y}$のなす角$\alpha$を用いて以下のように表される

$\mathbf{x} \cdot \mathbf{y} = sin \theta sin \theta' cos(\phi - \phi') + cos\theta cos\theta' = cos \alpha$

この式を整理していくと下式が得られる(参考)

$P_k (cos \alpha) = P_k (\mathbf{x} \cdot \mathbf{y}) = \frac{4\pi}{2k + 1} \sum^k_{m=-k} Y^m_k (\theta', \phi')*Y^m_k(\theta,\phi)$

これは三角関数の加法定理($cos(\theta - \theta') = cos\theta' cos\theta + sin\theta sin\theta'$)を一般化したものだと考えることができる

また、ウンゼルトの定理より、$\mathbf{x}=\mathbf{y}$の時は$P_k(\mathbf{x}\cdot\mathbf{y})=1$となる。

この記事のまとめ

$k$次球面調和関数$p$とは、斉次多項式かつ調和多項式である$p$を単位球面上に制限したものである。
ラプラシアンの回転不変性より$\Delta p (R(x)) = R(\Delta p(x))$が成立する。
$Y^m_k(\theta,\phi)$がよく球面多項式の例として挙げられる。
内積、対称性、加法定理など、ベクトルに近い性質も持ち合わせている。

随時更新予定です。