四元数ニューラルネットワークとGHR微積分

これは「FOLIO Advent Calendar 2023」6日目の記事です。

---

ニューラルネットワークで取り扱う数値を実数とは異なる数に拡張することは、機械学習や計算科学の発展における魅力的な課題の一つです。実数を用いた数値表現は多くのタスクにおいて十分な結果をもたらしてきましたが、新たな数値体系を導入することで、今までとは異なる問題が解決できるようになったり実数では見られなかった新たな現象が起こる可能性に期待することができるでしょう。例えば数値が取れる値を±1に制限したBinalized Neural Networksはハードウェアとの相性が良くメモリ効率の良い実装が可能であったり、拡大実数$\bar{\mathbb R}={\mathbb R}\cup\{-\infty,\infty\}$を用いた5層のReLUネットワークには任意の深さのReLUネットワークを埋め込むことができたりします。実数と異なる数への拡張は多種多様な方法を考えることができますが、本稿では複素数と四元数への拡張に焦点を当ててその背後にある数学的な概念やアルゴリズムを探っていきたいと思います。

ニューラルネットワーク

ニューラルネットワーク（以下NN）の中心をなす重要なコンセプトは、データを使った順方向の計算である推論と、その後に得られた結果をもとに重みを調整してネットワークを訓練する学習です。推論では、NNはデータを受け取り、層を通じて情報を処理して出力を生成します。学習では、推論の結果と実際の正解を比較し、逆方向の計算によって重みを修正してパラメータを最適化します。

例えば以下のような入力層・中間層・出力層を持つ3層のNNを考えてみましょう。

推論においては以下のように与えられた入力層の値 $x_1, x_2$ を線形変換し非線形な活性化関数に適用することを繰り返すことで最終的な出力を計算します。

$$\begin{matrix} h_i &=& \sigma\left(\displaystyle\sum_{j=1}^2 w^{(1)}_{i,j}x_j+b^{(1)}_i\right), && i = 1,2,3 \\ y_i &=& \sigma\left(\displaystyle\sum_{j=1}^3 w^{(2)}_{i,j}h_j+b^{(2)}_i\right), && i = 1,2 \\ \end{matrix}$$

また学習においては誤差関数 $L$ が与えることでNNの推論した出力と正解である教師データの間の違い（誤差）を測れるようにし、その誤差をできるだけ小さくするようにパラメータを更新するということを行います。パラメータの更新には誤差関数を更新するパラメータで偏微分した偏微分係数が必要になり、これは以下のように計算することができます。

$$\begin{matrix} \frac{\partial L}{\partial w^{(2)}_{i, j}} &=& \frac{\partial L}{\partial y_i}\frac{\partial y_i}{\partial w^{(2)}_{i, j}} \\ \frac{\partial L}{\partial b^{(2)}_i} &=& \frac{\partial L}{\partial y_i}\frac{\partial y_i}{\partial b^{(2)}_i} \\ \frac{\partial L}{\partial w^{(1)}_{i, j}} &=& \sum_{k=1}^2\frac{\partial L}{\partial y_k}\frac{\partial y_k}{\partial h_i}\frac{h_i}{\partial w^{(1)}_{i, j}} \\ \frac{\partial L}{\partial b^{(1)}_i} &=& \sum_{k=1}^2\frac{\partial L}{\partial y_k}\frac{\partial y_k}{\partial h_i}\frac{h_i}{\partial b^{(1)}_j} \\ \end{matrix}$$

以上の計算を丁寧に見ていくと

• 足し算
• 掛け算
• 活性化関数
• 微分
• 合成関数の微分法

のような要素から構成されていることが分かります。そして誤差関数$L$の形によっては積の微分法が必要になることもあります。

逆に言えば上記のような計算を備えているような数であればNNを考えることができるといえるでしょう。実数の他に足し算・掛け算・関数の微分ができるような数といえば複素数が身近な候補として挙がるのではないでしょうか。

複素解析

複素数を使ったNNが実現可能かどうか考えてみましょう。複素数は足し算と掛け算を行うことができます。関数の微分についてはどうでしょうか？

そもそも実数を引数に取り実数を返す関数、すなわち実関数 $f$ が点 $x$ において微分可能であるとは以下の極限が存在することでした。

$$\lim_{h \to 0}\frac{f(x+h) - f(x)}{h}$$

この時、実数には正負があるため $h$ を正の方向から近づける場合と負の方向から近づける場合のいずれの方法でも計算した値が一致することが、極限が存在する条件として要求されていることに注意してください。

同様に複素数を引数に取り複素数を返す関数、すなわち複素関数 $f$ が点 $z$ において微分可能であるとは通常、以下の極限が存在することとして定義されます。

$$\lim_{h \to 0}\frac{f(z+h) - f(z)}{h}$$

実数の場合と同様に $h$ を0に近づける方法（複素数なので複素平面を考えれば"方向"）によらず上記の計算が同じ値になることが要請されていることに注意してください。複素関数$f$がある開集合$D$の全ての点で微分可能である時、$f$は$D$において正則であるといい、$f$は$D$上の正則関数であると言います。

この正則という性質は非常に強力で、1回でも微分可能であれば無限回微分可能であったり、正則関数は局所的に収束ベキ級数で表すことができたり、解析接続によって一意に定義域を拡張することができたり、様々な性質を持っています。中には次のような定理もあります。

【リウヴィルの定理】 複素数全体において正則であるような複素関数 $f$ が有界であれば、すなわち、ある正の実数の定数 $M$ が存在して任意の複素数 $z$ に対して $|f(z)|\leq M$ が成り立つ時、複素関数 $f$ は定数関数となる。

つまり複素数全体で定義される定数関数ではない正則関数は必ず無限遠点に発散していくような部分を持たなければいけません。

ところで我々は複素数を使ってNNを実装することを目指していました。NNで使われる活性化関数にはシグモイド関数のように0から1の間をとる有界な関数が使われることがあります。しかしリウヴィルの定理によれば複素数全体で定義される有界な関数であると同時に微分可能であるような関数は定数関数を除いて存在しません。つまりシグモイド関数のような都合の良い正則関数は存在しないのです。

いやいや、でも最近はReLUのように非有界な活性化関数も使われるから問題ないのでは？と思われるかもしれません。しかし正則関数にはもう一つNNに用いる上で困った性質があります。それは実数値のみを取る正則関数もまた定数関数に限られるということです（この事実は後で証明します）。NNの学習時には誤差関数を使用しますが、誤差関数は出力結果の良し悪しを測るために大小が決まる序数的な性質が求められるので、出力される値は順序の定まっている実数である方が都合が良いでしょう。しかし複素数の場合にこのような誤差関数を用いると、学習を行うために誤差関数を微分可能だと仮定した途端に定数関数に限られてしまうのです。

以上に見てきたような理由から今回は正則関数を扱う分野である複素解析は使いませんが、複素解析自体は数学の中でも非常に古い歴史を持ち、多くの成果を生み出した重要な分野の一つです。興味がある人は以下のような文献から辿ってみてください。

• 複素解析 - Wikipedia
• 複素関数の基礎のキソ
• L.V.アールフォルス『複素解析』

ウィルティンガー微分

複素関数に複素数の意味での微分可能性を定義する（すなわち正則関数に限ってしまう）と非常に強力な性質を持つが故にNNに用いるには都合が悪くなってしまう事がわかりました。そこで、もう少し緩い意味での微分を考えるために複素数$z$を二次元の実ベクトル$x+iy$だと思って微分することを考えてみましょう。

二次元のベクトルといえば$(x, y)$のように成分を括弧とカンマで区切って書かれるもの（いわゆる数ベクトル）を想像するかもしれませんが、複素数$x+iy$も立派な二次元ベクトルです。もちろん複素数がベクトル空間の定義を満たしていることを確認すれば良いのですが、実際に定数倍や足し算といった線形演算において同様の振る舞いをすることを見ても直感的に理解できると思います。（以下 $a,b,c,d\in{\mathbb R}$）

演算	数ベクトル	複素数
定数倍	$a(x, y) = (ax, ay)$	$a(x+iy)=ax+iay$
足し算	$(a, b) + (c, d) = (a+c, b+d)$	$(a+ib) + (c+id) = (a+c)+i(b+d)$

数ベクトルと複素数の大きな違いはベクトル同士の掛け算が定義されているかどうかです。複素数のように掛け算が定まっているベクトル空間は多元環と呼ばれています。

以下では複素関数を実ベクトルを取って実ベクトルを返す関数のように扱います。これにより実変数の偏微分を自然に考えることができますが、上記のように掛け算が定まっていることで単なるベクトルの関数よりもう少し便利に扱えることを見ていきます。

勾配降下法

さて我々の目的は勾配降下法によって誤差関数を最小化するパラメータを求めることにありました。まず一変数の場合に勾配降下法がなぜうまくいくのかを思い出しましょう。実一変数関数$f(x)$において$x$を微小に$\Delta x$だけ変化させた時、テイラー展開から以下のような式が成り立ちます。

$$f(x+\Delta x) = f(x) + \frac{df}{dx}\Delta x + o(\Delta x)$$

ここで $o(\Delta x)$ はランダウの記号で、$o(\Delta x)$ の項が $\Delta x$ よりも十分早く小さくなることを表しています。

この式を少し変形すると

$$f(x+\Delta x) - f(x) = \frac{df}{dx}\Delta x + o(\Delta x)$$

となります。もし $\Delta x$ がある小さな正の定数 $\eta$ を使って $\Delta x = -\eta\frac{df}{dx}$ と表されていたとすると上記式は

$$f(x+\Delta x) - f(x) = -\eta\left|\frac{df}{dx}\right|^2 + o(\Delta x)$$

となり $\eta$ が十分に小さければ右辺は負であることから $x$ を $x+\Delta x$ に変化させることによって $f(x+\Delta x)$ を $f(x)$ より小さくすることができます。

同じことを$f(x+iy)$についても考えてみましょう。$f$の値を小さくするということが意味を持つように$f$は実数値を取る複素関数であるとします。微小変化させた$f((x+\Delta x)+i(y+\Delta y))$と$f(x+iy)$との差分を考えると

$$f((x+\Delta x)+i(y+\Delta y)) = f(x+iy) + \frac{\partial f}{\partial x}\Delta x + \frac{\partial f}{\partial y}\Delta y + o(\Delta x, i\Delta y)$$

となり、$f(x+iy)$を左辺に移項すると

$$f((x+\Delta x)+i(y+\Delta y)) - f(x+iy) = \frac{\partial f}{\partial x}\Delta x + \frac{\partial f}{\partial y}\Delta y + o(\Delta x, i\Delta y)$$

となります。

もし $\Delta x, \Delta y$ がある小さな正の定数 $\eta$ を使って

$$\begin{matrix} \Delta x &=& -\eta\frac{\partial f}{\partial x} \\ \Delta y &=& -\eta\frac{\partial f}{\partial y} \\ \end{matrix}$$

と表されていたとすると先程の式は

$$f(x+\Delta x, y+\Delta y) - f(x, y) = -\eta\left(\left|\frac{\partial f}{\partial x}\right|^2 + \left|\frac{\partial f}{\partial y}\right|^2\right) + o(\Delta x, \Delta y)$$

となり $\eta$ が十分に小さければ右辺は負の実数であることから、左辺を小さくする事ができます。

ここで $\Delta z$ を $\Delta x + i\Delta y$ と置くと、

$$\begin{matrix} \Delta z &=& \Delta x + i\Delta y \\ &=& -\eta\left(\frac{\partial f}{\partial x} + i\frac{\partial f}{\partial y}\right) \\ \end{matrix}$$

となります。この微分演算子をまとめたものには意味がありそうなので

$$\frac{\partial f}{\partial \overline{z}} = \frac{1}{2}\left(\frac{\partial f}{\partial x} + i\frac{\partial f}{\partial y}\right)$$

と定義することにしましょう。

ここで係数 $\frac{1}{2}$ が掛けられているのと微分する変数が共役複素数 $\overline{z}$ になっているのは、$f(z) = c\overline{z}$ の微分を考えると

$$\begin{matrix} \frac{\partial f}{\partial \overline{z}} &=& \frac{1}{2}\left(\frac{\partial f}{\partial x}+i\frac{\partial f}{\partial y}\right) \\ &=& \frac{1}{2}\left(\frac{\partial c\overline{z}}{\partial x}+i\frac{\partial c\overline{z}}{\partial y}\right) \\ &=& \frac{1}{2}\left(\frac{\partial c(x-iy)}{\partial x}+i\frac{\partial c(x-iy)}{\partial y}\right) \\ &=& \frac{1}{2}\left(c +i (-ic)\right) \\ &=& \frac{1}{2}\left(c + c\right) \\ &=& c \end{matrix}$$

となって係数 $c$ が出てくるように直感的な振る舞いになるよう定めています。

もう一つの演算子

$\frac{\partial}{\partial \overline{z}}$ の性質をいくつか見ていきましょう。

まずこれは通常の微分の線形結合で表されているので以下のような線形性を満たします。

$$\frac{\partial (af)}{\partial \overline{z}} = a\frac{\partial f}{\partial \overline{z}},\ \ a \in {\mathbb C}$$

$$\frac{\partial (f+g)}{\partial \overline{z}} = \frac{\partial f}{\partial \overline{z}} + \frac{\partial g}{\partial \overline{z}}$$

積の微分法についても以下のように成り立つことが分かります。

$$\begin{matrix} \frac{\partial (fg)}{\partial \overline{z}} &=& \frac{1}{2}\left(\frac{\partial (fg)}{\partial x}+i\frac{\partial (fg)}{\partial y}\right) \\ &=& \frac{1}{2}\left(\left(\frac{\partial f}{\partial x}g+f\frac{\partial g}{\partial x}\right)+i\left(\frac{\partial f}{\partial y}g+f\frac{\partial g}{\partial y}\right)\right) \\ &=&\frac{1}{2}\left(\frac{\partial f}{\partial x}+i\frac{\partial f}{\partial y}\right)g + f\frac{1}{2}\left(\frac{\partial g}{\partial x}+i\frac{\partial g}{\partial y}\right) \\ &=& \frac{\partial f}{\partial z}g + f\frac{\partial g}{\partial z} \end{matrix}$$

最後に合成関数の微分法についてはどうでしょうか？

$$\frac{\partial (f\circ g)}{\partial \overline{z}} = \frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial x}+i\frac{\partial (f\circ g)}{\partial y}\right)$$

ここで

$$\frac{\partial (f\circ g)}{\partial x} = \frac{\partial (f\circ g)}{\partial g}\frac{\partial g}{\partial x}$$

と評価したくなるところですが$\frac{\partial}{\partial g}$は実変数による微分ではなく先程定義した複素数による微分演算子なのでこのような関係式が成り立つかどうかまだ分かっていません。

議論を進めるために$g(z)$を実部と虚部に分けて$g(z) = u(z) + iv(z)$と表すと、実変数については合成関数の微分法が使えるので

$$\begin{matrix} \frac{\partial (f\circ g)}{\partial \overline{z}} &=& \frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial x}+i\frac{\partial (f\circ g)}{\partial y}\right) \\ &=& \frac{1}{2}\left(\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial x}+\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial x}\right)+i\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial y}+\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial y}\right)\right) \\ &=& \frac{1}{2}\left(\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial x}+i\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial y}\right)+\left(\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial x}+i\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial y}\right)\right) \\ &=&\frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial u}\left(\frac{\partial u}{\partial x}+i\frac{\partial u}{\partial y}\right)+\frac{\partial (f\circ g)}{\partial v}\left(\frac{\partial v}{\partial x}+i\frac{\partial v}{\partial y}\right)\right) \\ &=&\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial \overline{z}}+\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial \overline{z}} \end{matrix}$$

と評価することができます。ここで定義より $\frac{\partial (f\circ g)}{\partial \overline{g}}$ と $\frac{\partial \overline{g}}{\partial z}$ の積は

$$\begin{matrix} \frac{\partial (f\circ g)}{\partial \overline{g}}\frac{\partial \overline{g}}{\partial z} &=& \frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial u}+i\frac{\partial (f\circ g)}{\partial v}\right)\left(\frac{\partial u}{\partial z}-i\frac{\partial v}{\partial z}\right) \\ &=& \frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial z}+\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial z} - i\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial v}{\partial z}-\frac{\partial (f\circ g)}{\partial v}\frac{\partial u}{\partial z}\right)\right) \end{matrix}$$

となり、また

$$\begin{matrix} \overline{\frac{\partial \overline{(f\circ g)}}{\partial \overline{g}}}\frac{\partial g}{\partial z} &=& \frac{1}{2}\overline{\left(\frac{\partial \overline{(f\circ g)}}{\partial u}+i\frac{\partial \overline{(f\circ g)}}{\partial v}\right)}\frac{\partial g}{\partial z} \\ &=& \frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial u}-i\frac{\partial (f\circ g)}{\partial v}\right)\left(\frac{\partial u}{\partial z}+i\frac{\partial v}{\partial z}\right) \\ &=& \frac{1}{2}\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial z}+\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial z} + i\left(\frac{\partial (f\circ g)}{\partial u}\frac{\partial v}{\partial z}-\frac{\partial (f\circ g)}{\partial v}\frac{\partial u}{\partial z}\right)\right) \end{matrix}$$

となります。ここで複素数$z$による微分演算子を新たに

$$\frac{\partial f}{\partial z} = \frac{1}{2}\left(\frac{\partial f}{\partial x}-i\frac{\partial f}{\partial y}\right)$$

と定義すると、上記2段目の式から

$$\overline{\frac{\partial \overline{(f\circ g)}}{\partial \overline{g}}}\frac{\partial g}{\partial z} = \frac{\partial (f\circ g)}{\partial g}\frac{\partial g}{\partial z}$$

と書けることが分かります。以上より

$$\frac{\partial (f\circ g)}{\partial g}\frac{\partial g}{\partial z} + \frac{\partial (f\circ g)}{\partial \overline{g}}\frac{\partial \overline{g}}{\partial z} = \frac{\partial (f\circ g)}{\partial u}\frac{\partial u}{\partial z}+\frac{\partial (f\circ g)}{\partial v}\frac{\partial v}{\partial z}$$

となり $\frac{\partial}{\partial z}$ による合成関数の微分法は

$$\frac{\partial (f\circ g)}{\partial z} = \frac{\partial (f\circ g)}{\partial g}\frac{\partial g}{\partial z} + \frac{\partial (f\circ g)}{\partial \overline{g}}\frac{\partial \overline{g}}{\partial z}$$

と書けることが分かりました。

これまでの議論に出てきた演算子

$$\frac{\partial }{\partial z}, \frac{\partial }{\partial \overline{z}}$$

は ウィルティンガー微分（またはCR微積分） と呼ばれています。

ウィルティンガー微分を用いるメリットの一つは、複素関数の実微分を考える際に$x, y$という実数に言及せず、複素数とその共役複素数の言葉のみで語れるというところにあります。この恩恵はあとで複素NNを実装する際により鮮明になります。

形式的な導出

上記の議論では勾配降下法の観点から $\frac{\partial}{\partial\overline{z}}$ の必要性を、合成関数の微分法の観点から $\frac{\partial}{\partial z}$ の必要性を与えて定義しましたが、以下のような形式的な議論からウィルティンガー微分を導くことも可能です。

引数を二次元ベクトル$x+iy$とみなしたときの複素関数 $f$ の全微分は

$$df = \frac{\partial f}{\partial x}dx + \frac{\partial f}{\partial y}dy$$

と書けます。

また共役複素数 $\overline{z}$ を用いれば $x, y$ は

$$\begin{matrix} x &=& \frac{z + \overline{z}}{2} \\ y &=& \frac{z - \overline{z}}{2i} \\ \end{matrix}$$

と書けることから、全微分

$$\begin{matrix} dx &=& \frac{dz + d\overline{z}}{2} \\ dy &=& \frac{dz - d\overline{z}}{2i} \\ \end{matrix}$$

を考えて $f$ の全微分に代入して整理すると

$$df = \frac{1}{2}\left(\frac{\partial f}{\partial x}-i\frac{\partial f}{\partial y}\right)dz+\frac{1}{2}\left(\frac{\partial f}{\partial x}+i\frac{\partial f}{\partial y}\right)d\overline{z}$$

となり $z$ と $\overline{z}$ を変数とした時の $f$ の全微分が

$$df = \frac{\partial f}{\partial z}dz + \frac{\partial f}{\partial \overline{z}}d\overline{z}$$

と書けることを考えると

$$\begin{matrix} \frac{\partial f}{\partial z} &=& \frac{1}{2}\left(\frac{\partial f}{\partial x}-i\frac{\partial f}{\partial y}\right) \\ \frac{\partial f}{\partial \overline{z}} &=& \frac{1}{2}\left(\frac{\partial f}{\partial x}+i\frac{\partial f}{\partial y}\right) \\ \end{matrix}$$

と対応していることが分かります。

コーシー・リーマンの方程式

複素関数 $f$ が正則関数であることは、 $f$ が $f(z) = u(z) + iv(z)$ と表される時に

$$\begin{matrix} \frac{\partial u}{\partial x} &=& \frac{\partial v}{\partial y} \\ \frac{\partial u}{\partial y} &=& -\frac{\partial v}{\partial x} \\ \end{matrix}$$

という関係式を満たすことの必要十分条件であることが知られています。この関係式はコーシー・リーマンの方程式として知られています。

ところで $\frac{\partial f}{\partial \overline{z}}$ は定義より

$$\begin{matrix} \frac{\partial f}{\partial \overline{z}} &=& \frac{1}{2}\left(\frac{\partial f}{\partial x} + i\frac{\partial f}{\partial y}\right) \\ &=& \frac{1}{2}\left(\left(\frac{\partial u}{\partial x} + i\frac{\partial v}{\partial x}\right) + i\left(\frac{\partial u}{\partial y}+i\frac{\partial v}{\partial y}\right)\right) \\ &=& \frac{1}{2}\left(\left(\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}\right) + i\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)\right) \\ \end{matrix}$$

となることから $\frac{\partial f}{\partial \overline{z}} = 0$ という式とコーシー・リーマンの方程式は同値であることが分かります。

さらに $f$ が実数値のみを取る正則関数だと仮定しましょう。仮定より $v(z)$ は恒等的に値が0になりますからコーシー・リーマンの方程式より

$$\begin{matrix} \frac{\partial u}{\partial x} &=& \frac{\partial v}{\partial y} &=& 0 \\ \frac{\partial u}{\partial y} &=& -\frac{\partial v}{\partial x} &=& 0 \\ \end{matrix}$$

となり $f$ は定数関数となることが分かります。

複素ニューラルネットワーク

ウィルティンガー自動微分

それでは上述のウィルティンガー微分を使って複素NNを実装してみましょう。とはいえ複素NNの各パラメータによる微分を手で導出することは大変ですから自動微分のテクニックを用いることにします。名付けてウィルティンガー自動微分です^[1]。

通常の実数におけるフォワードモード自動微分は二重数という代数的な構造を使って実装することができます（参考）。しかしウィルティンガー微分において同様のフォワードモード自動微分を考える場合、複素数と共役複素数による2つの微分が存在するため二重数だとデータを格納する場所が一つ足りません。そこで二重数を拡張した以下のような数を考えます。名付けて三重複素数です。

三重複素数${\mathbb T}$は複素数に以下の性質を満たす数 $\epsilon, \epsilon'$ を添加して作ります。

$$\begin{matrix} \epsilon, \epsilon' \notin {\mathbb C} \\ \epsilon^2 = \epsilon'^2 = \epsilon\epsilon' = \epsilon'\epsilon = 0 \end{matrix}$$

つまり三重複素数は

$${\mathbb T} = \{a + b\epsilon + c\epsilon' | a, b, c \in {\mathbb C}\}$$

のように表される数です。

複素関数 $f$ に対して三重複素数

$$f + \frac{\partial f}{\partial z}\epsilon + \frac{\partial f}{\partial \overline{z}}\epsilon'$$

を対応させると、この対応はフォワードモードのウィルティンガー自動微分を与えます。

実際にHaskellで実装して確かめてみましょう。

import Data.Complex (Complex(..))
import Linear.Conjugate (conjugate)

-- | 実数
type R = Double

-- | 複素数
type C = Complex R

-- | 三重複素数
data T = T C C C deriving Show

complex, gradZ, gradZC :: T -> C
complex (T a b c) = a
gradZ   (T a b c) = b
gradZC  (T a b c) = c

複素数としてHaskellが標準で提供しているData.Complexに Double を適用したものを使います。三重複素数はそのように定義した複素数 C の三つ組として定義し、それぞれの要素を取り出す関数を定義しておきます。

このように定義した三重複素数 T にウィルティンガー微分と整合性のある演算を定義していきましょう。

instance Num T where
  (T z z' zc') + (T w w' wc') = T (z + w) (z' + w') (zc' + wc')
  (T z z' zc') * (T w w' wc') = T (z * w) (z' * w + z * w') (zc' * w + z * wc')
  negate (T z z' zc') = T (negate z) (negate z') (negate zc')
  abs (T z z' zc') = T (abs z) w' wc' where
    w'  = conjugate z / (2 * abs z) * z'  + z / (2 * abs z) * conjugate zc'
    wc' = conjugate z / (2 * abs z) * zc' + z / (2 * abs z) * conjugate z'
  signum (T z _ _) = T (signum z) 0 0
  fromInteger n = T (fromInteger n) 0 0

instance Fractional T where
  recip (T z z' zc') = T (recip z) (-1 * recip (z * z) * z') (-1 * recip (z * z) * zc')
  fromRational r = T (fromRational r) 0 0

線形性や積の微分法を使って定義していることが分かります。絶対値 abs の微分については $|z| = \sqrt{z\overline{z}}$ と表せることから合成関数の微分法を使って

$$\begin{matrix} \frac{\partial \sqrt{f\overline{f}}}{\partial z} &=& \frac{\overline{f}}{2|f|}\frac{\partial f}{\partial z} + \frac{f}{2|f|}\frac{\partial \overline{f}}{\partial z} \\ &=& \frac{\overline{f}}{2|f|}\frac{\partial f}{\partial z} + \frac{f}{2|f|}\overline{\frac{\partial f}{\partial \overline{z}}} \\ \end{matrix}$$

のように計算することができます。2行目への変形は $\frac{\partial \overline{f}}{\partial z} = \overline{\frac{\partial f}{\partial \overline{z}}}$ となることを用いていますが、これは定義より

$$\begin{matrix} \frac{\partial \overline{f}}{\partial z} &=& \frac{1}{2}\left(\frac{\partial\overline{f}}{\partial x}-i\frac{\partial\overline{f}}{\partial y}\right) \\ &=& \frac{1}{2}\left(\overline{\frac{\partial f}{\partial x}}-i\overline{\frac{\partial f}{\partial y}}\right) \\ &=& \overline{\frac{1}{2}\left(\frac{\partial f}{\partial x}+i\frac{\partial f}{\partial y}\right)} \\ &=& \overline{\frac{\partial f}{\partial \overline{z}}} \\ \end{matrix}$$

となることから分かります。

recip についても同様に合成関数の微分法を使って考えることができますが $\frac{1}{z}$ が正則関数であることから $\frac{\partial}{\partial\overline{z}}$ の項が0になり消えるので通常の微分と同様に計算することができています。

instance Floating T where
  pi                 = T pi 0 0
  exp   (T z z' zc') = T (exp z)   (exp z * z')              (exp z * zc')
  log   (T z z' zc') = T (log z)   (z' / z)                  (zc' / z)
  sin   (T z z' zc') = T (sin z)   (cos z * z')              (cos z * zc')
  cos   (T z z' zc') = T (cos z)   (-sin z * z')             (-sin z * zc')
  asin  (T z z' zc') = T (asin z)  (z' / sqrt (1 - z ** 2))  (zc' / sqrt (1 - z ** 2))
  acos  (T z z' zc') = T (acos z)  (-z' / sqrt (1 - z ** 2)) (-zc' / sqrt (1 - z ** 2))
  atan  (T z z' zc') = T (atan z)  (z' / (1 + z ** 2))       (zc' / (1 + z ** 2))
  sinh  (T z z' zc') = T (sinh z)  (cosh z * z')             (cosh z * zc')
  cosh  (T z z' zc') = T (cosh z)  (sinh z * z')             (sinh z * zc')
  asinh (T z z' zc') = T (asinh z) (z' / sqrt (1 + z ** 2))  (zc' / sqrt (1 + z ** 2))
  acosh (T z z' zc') = T (acosh z) (z' / sqrt (z ** 2 - 1))  (zc' / sqrt (z ** 2 - 1))
  atanh (T z z' zc') = T (atanh z) (z' / (1 - z ** 2))       (zc' / (1 - z ** 2))

Floating で定義される関数も全て正則関数なので通常の微分と同様に計算することができます。

それでは簡単な例で実装した三重複素数を使って実際にウィルティンガー微分ができているか見てみましょう。

> z = T (3:+4) 1 0

> z^2 + 2 * z + 1
T (0.0 :+ 32.0) (8.0 :+ 8.0) (0.0 :+ 0.0)

> abs z
T (5.0 :+ 0.0) (0.3 :+ (-0.4)) (0.3 :+ 0.4)

まず $z = (3+4i) + \epsilon$ と定義します。実数のフォワードモード自動微分の時と同様、微分したい変数部分に1を、それ以外の変数部分には0を代入します。

最初の $f(z) = z^2 + 2z + 1$ を評価すると、

$$\begin{matrix} f(3+4i) &=& 32i \\ \frac{\partial f}{\partial z}(3+4i) &=& 8 + 8i \\ \end{matrix}$$

となっています。$\frac{\partial f}{\partial\overline{z}}$ に対応する項は0ですが、この値を求めるには$z$として $(3+4i) + \epsilon'$ を代入して計算する必要があるので、共役複素数による微分係数はこの計算だけでは分かりません。

次に $f(z) = |z|$ を評価すると、

$$\begin{matrix} f(3+4i) &=& 5 \\ \frac{\partial f}{\partial z}(3+4i) &=& 0.3 - 0.4i \\ \end{matrix}$$

となることが分かります。元の数から$\overline{\frac{\partial f}{\partial z}}$ を引くと期待通り絶対値が小さくなるようになっていることが分かりますね。

複素NN

それでは定義した三重複素数を使って学習可能な複素NNを作っていきましょう。今回は線形代数のライブラリとして linear を使います。

NNを実装するために線形層と活性化関数を定義していきましょう。

import GHC.TypeNats (KnownNat, Nat)

import Linear.Matrix ((!*))
import Linear.V (V(..), fromVector)

-- | 行列
type M (m :: Nat) (n :: Nat) a = V m (V n a)

-- | 線形層
linear :: Num a => Matrix m n a -> Vector m a -> Vector n a -> Vector m a
linear w b z = w !* z ^+^ b

-- | 活性化関数の定義
class Activatable a where
  sigmoid :: a -> a

instance Activatable R where
  sigmoid x = 1 / (1 + exp (-x))

instance Activatable C where
  sigmoid z = fmap sigmoid z

instance Activatable T where
  sigmoid (T z@(x :+ y) z' zc') = T (sigmoid z) w' wc' where
    sigmoid' x = sigmoid x * (1 - sigmoid x)
    w'  = (sigmoid' x :+ (-sigmoid' y)) * z'  + (sigmoid' x :+ sigmoid' y) * conjugate zc'
    wc' = (sigmoid' x :+ (-sigmoid' y)) * zc' + (sigmoid' x :+ sigmoid' y) * conjugate z'

活性化関数は実数、複素数、三重複素数それぞれに適用できるようにするために型クラスで定義しています。

活性関数は実数に対しては以下のような関数になっており

$${\rm \sigma_{\mathbb R}}(x) = \frac{1}{1+\exp(-x)}$$

この活性化関数はグラフを書くと以下のようなS字形になっているのでシグモイド関数と呼ばれています。。

このように0と1がハッキリと区別できるので実NNでは分類問題を解くためにこの関数が用いられることも多いです（定量的には $x=0$ を決定境界として区別します）。

シグモイド関数の導関数は以下のように計算することができます。

$$\begin{matrix} \frac{\partial{\rm \sigma_{\mathbb R}}(x)}{\partial x} &=& \frac{\exp(-x)}{(1+\exp(-x))^2} \\ &=& \frac{\exp(-x)}{1+\exp(-x)}\frac{1}{1+\exp(-x)} \\ &=& (1-{\rm \sigma_{\mathbb R}}(x)){\rm \sigma_{\mathbb R}}(x) \end{matrix}$$

複素数に対しては実部と虚部ぞれぞれに上記のシグモイド関数を適用する形で定義されています。

$${\rm \sigma_{\mathbb C}}(x+iy) = {\rm \sigma_{\mathbb R}}(x) + i\ {\rm \sigma_{\mathbb R}}(y)$$

この関数は複素数全体で定義されますが有界なので正則関数ではありません。なのでウィルティンガー微分を考える必要が出てくるわけです。この関数のウィルティンガー自動微分は定義から実部と虚部に言及する形で実装する必要があり少し煩雑になっています。逆に言えば今までの標準的な関数が複素数を実部と虚部に分解することなく複素数を使う形のまま計算できていたのはウィルティンガー微分の恩恵と言えるでしょう。

XOR問題

それでは複素NNを使って具体的な問題を解いてみましょう。通常のNNではXOR問題として

入力1	入力2	出力
0	0	0
1	0	1
0	1	1
1	1	0

という対応を学習することは2層のニューラルネットワークでは解くことができず3層ニューラルネットワークを用いる必要があることが知られています。

実は複素NNはこのXOR問題を2層で解くことができるのです。

NITTA, Tohru. Orthogonality of decision boundaries in complex-valued neural networks. Neural computation, 2004, 16.1: 73-97.

上記論文ではXOR問題を以下のようにエンコードしています。

入力（$z=x+iy$）	出力（$Z=X+iY$）
$-1-i$	$1+0i$
$-1+i$	$0+0i$
$1-i$	$1+i$
$1+i$	$0+i$

この時、入力の $x, y$ を $0, 1$ に、出力の $1+0i, 0+i$ を $0$ に, $0+0i, 1+i$ を $1$ に対応付けることでXOR問題を解くことができるというわけです。

実際に二層複素NNを使ってこのような対応を学習することが可能か実験してみましょう。

import Data.List (cycle, foldl')

import qualified Data.Vector as V
import System.Random (randomRIO, initStdGen)

-- | リストからベクトルに変換する
toV :: KnownNat n => [a] -> V n a
toV = fromJust . fromVector . V.fromList

-- | リストから行列に変換する
toM :: (KnownNat m, KnownNat n) => [[a]] -> M m n a
toM = toV . map toV

main :: IO ()
main = do
  -- パラメータを乱数で初期化
  w <- (uncurry (:+)) <$> randomRIO ((-1, -1), (1, 1)) :: IO C
  b <- (uncurry (:+)) <$> randomRIO ((-1, -1), (1, 1)) :: IO C
  let -- 教師データ
      trainingData = take 1000 $ cycle [
          ((-1):+(-1), 1:+0),
          ((-1):+  1 , 0:+0),
          (  1 :+(-1), 1:+1),
          (  1 :+  1 , 0:+1)
          ]
      -- 学習モデル
      model w b x = fmap sigmoid $ linear w b x
      -- 学習係数
      learningRate = 0.1 :+ 0
      -- パラメータ更新の1ステップ
      updateParam (w, b) (x, y) =
          let -- パラメータを整形
              mw dw = toM @1 @1 [[T w dw 0]]
              vb db = toV @1     [T b db 0]
              vx    = toV @1     [T x  0 0]
              ty    = T y 0 0
              -- 誤差関数の微分係数を計算
              df (dw, db) = gradZ $ (abs (ty - V.head (model (mw dw) (vb db) vx))) ^ 2
              dfdw = df (1, 0)
              dfdb = df (0, 1)
           -- パラメータの更新
           in (w - learningRate * conjugate dfdw, b - learningRate * conjugate dfdb)
      -- 学習データによるパラメータ更新
      (w', b') = foldl' updateParam (w, b) trainingData
      mw = toM @1 @1 [[w']]
      vb = toV @1 [b']
  -- 学習結果を表示
  putStrLn $ concat ["(-1) :+ (-1) | ", show $ model mw vb (toV @1 [(-1):+(-1)])]
  putStrLn $ concat ["(-1) :+   1  | ", show $ model mw vb (toV @1 [(-1):+  1 ])]
  putStrLn $ concat ["  1  :+ (-1) | ", show $ model mw vb (toV @1 [  1 :+(-1)])]
  putStrLn $ concat ["  1  :+   1  | ", show $ model mw vb (toV @1 [  1 :+  1 ])]

学習モデルは

$$f(z) = \sigma_{\mathbb C}\left(wz+b\right)$$

という二層NNで、誤差関数は

$$|y - f(x)|^2$$

という二乗誤差を測る関数になっています。

各パラメータ $w, b$ による微分係数はそれぞれの三重複素数の$\epsilon$に1を対応させて1回ずつ計算することにより計算しています。微分係数のベクトルの計算が一度の計算で求まらないのはフォワードモード自動微分の欠点ですね^[2]。

実際に実行してみると

$ cabal run
(-1) :+ (-1) | Vector [0.9718682830396489 :+ 5.092719254396843e-2]
(-1) :+   1  | Vector [7.407505253473312e-2 :+ 1.96242268684455e-2]
  1  :+ (-1) | Vector [0.9279916101837028 :+ 0.9586299723769597]
  1  :+   1  | Vector [2.8978429197853116e-2 :+ 0.8963086315651159]

となり教師データの通りに対応を学習できていることが分かります（種明かしをすると十分大きな正の実数$a$を使って$w=ai, b=0$とすれば求める解となります）。

ところで出力の $1+0i, 0+i$ を $0$, $0+0i, 1+i$ を $1$ に対応付けるというのはどういうことでしょうか？これが許されるなら出力が二次元で二層の実NNでも同じことができるのではないでしょうか？

実は複素NNは実NNの特殊な場合にあたるので確かに二次元二層実NNでも同様の学習は可能です。しかし複素NNなら一つの数（複素数）で分類できるという違いが重要になります。出力が二次元の二層実NNの場合、出力を一つの実数にするには2つの実数を足し合わせるなどの計算を新たに行う必要がありますが、そうすると集約の計算が三層目に対応して結局三層NNになってしまいます。実数の場合はシグモイド関数の形から0, 1を対応させて分類を行っていたと思いますが、複素数の場合には $\sigma_{\mathbb C}(z)$ の絶対値を図示すると

このように複素平面の領域を4つに分割するような形になっており、この領域それぞれを $0,1$ に対応させることでXOR問題を解いているという仕組みだったのです（領域のどこにあるかを判定するための決定境界は $x=0, y=0$ という直行する2つの直線になります）。

先に触れた通り、複素NNは実NNの特殊な場合にあたります。学習モデルの計算

$$f(z) = \sigma_{\mathbb C}\left(wz+b\right)$$

を実成分に分解して書くと

$$\begin{matrix} f(x+iy) &=& \sigma_{\mathbb C}\left((w_x+iw_y)(x+iy)+(b_x+ib_y)\right) \\ &=& \sigma_{\mathbb C}\left((w_xx-w_yy+b_x)+i(w_yx+w_xy+b_y)\right) \\ &=& \sigma_{\mathbb R}\left(w_xx-w_yy+b_x\right)+i\sigma_{\mathbb R}\left(w_yx+w_xy+b_y\right) \\ \end{matrix}$$

と計算できます。ここで

$$\begin{matrix} z &\mapsto& \left(\begin{matrix}x\\ y\end{matrix}\right) \\ w &\mapsto& \left(\begin{matrix}w_x&-w_y\\ w_y&w_x\end{matrix}\right) \\ b &\mapsto& \left(\begin{matrix}b_x\\ b_y\end{matrix}\right) \\ \end{matrix}$$

という対応を考えると

$$\begin{matrix} f\left(\left(\begin{matrix}x\\ y\end{matrix}\right)\right) &=& \sigma\left(\left(\begin{matrix}w_x&-w_y\\ w_y&w_x\end{matrix}\right)\left(\begin{matrix}x\\ y\end{matrix}\right)+\left(\begin{matrix}b_x\\ b_y\end{matrix}\right)\right) \\ &=& \left(\begin{matrix}\sigma_{\mathbb R}\left(w_xx-w_yy+b_x\right)\\ \sigma_{\mathbb R}\left(w_yx+w_xy+b_y\right)\end{matrix}\right) \\ \end{matrix}$$

となり上記の複素NNが二次元二層実NNに埋め込まれていることが分かります。異なる点は重み $w$ が、二次元二層実NNであれば四成分を自由に設定することができますが、複素NNを埋め込んだ場合にはこの自由度が二成分に制限されてしまう所にあります。ですが、この制約があることで重み$w$は入力の拡大と回転を行っていると自然に解釈することができるようになり、このような操作と親和性のあるデータであれば帰納バイアスとして有利に働く可能性もあるのです。同様の解釈においてバイアス$b$は平行移動に対応しています。

このように複素NNは一つの数による表現力が実NNよりも高いことや、拡大と回転と平行移動という幾何学的な直感を伴うデータの変換により実現されていることなどから、実NNとは異なる表現力を持ち、実NNが苦手とする問題をより簡単に解けるようになることが期待されているのです。

野生の複素数データ

複素NNを実装して実際に動くことはわかりましたが、そもそもどのようなデータが複素数で表されているのでしょうか？

まず信号処理の文脈においてはフーリエ変換やヒルベルト変換を用いて生の信号を複素数に変換して用いることが多々あります。複素数に変換することで振幅と位相の値を分離することができるので直感的にも分かりやすい形式で扱うことができます。

このように位相の情報が得られる場合にデータを複素数で表現する例は多く、例えば衛星から地上の様子を測る画像データでは反射強度の他に位相の情報も記録されており、この強度と位相を複素数にエンコードすることで、1ピクセルが複素数に対応している複素画像（SLC画像）になっています（SAR干渉画像の作成手順と見方 | 国土地理院）。また医療で用いられるMRI（Magnetic Resonance Imaging（磁気共鳴画像））においてもデータは位相情報が加味された複素数で表現されています（第1回：MRイメージングの基礎）。

四元数

さて実数、複素数ときたら次に来るのは四元数ですよね？実は四元数が発見されたのは1843年10月16日なので今年は四元数発見からちょうど180周年なのです。せっかくなので四元数を使ったNNについても考えてみましょう。

四元数は

$$i^2 = j^2 = k^2 = ijk = -1$$

という関係式を満たす元を実数に加えて作られる数です（この関係式をハミルトンがブルーム橋の石に刻みつけてから今年でちょうど180年というわけです）。

集合として書くと

$${\mathbb H} = \{a+bi+cj+dk|a,b,c \in {\mathbb R}\}$$

となります。

複素数を実部と虚部に分けたのと同様に、四元数の実数部分に対応する $a$ はスカラー部、それ以外の $bi+cj+dk$ はベクトル部と呼ばれます。いずれも今日のスカラーとベクトルの起源となる言葉です^[3]。

2つの四元数の足し算は以下のように成分ごとの足し算として定義されます。

$$(a+bi+cj+dk) + (x+yi+zj+wk) = (a+x) + (b+y)i + (c+z)j + (d+w)k$$

一般の四元数の掛け算を見る前に $1,i,j,k$ の積がどのような結果になるか見てみましょう。

$\times$	$1$	$i$	$j$	$k$
$1$	$1$	$i$	$j$	$k$
$i$	$i$	$-1$	$k$	$-j$
$j$	$j$	$-k$	$-1$	$i$
$k$	$k$	$j$	$-i$	$-1$

この表は上述の四元数の関係式から計算することができます。

よく見てみると $ij\neq ji$ となっていることが分かると思います。このように四元数において掛け算は可換ではないのです。

上記の計算と分配法則を使えば四元数の掛け算は以下のように計算することができます。

$$\begin{matrix} (a+bi+cj+dk) \times (x+yi+zj+wk) &=& (ax - by - cz - dw) \\ &+& (ay + bx + cw - dz) i \\ &+& (az - bw + cx + dy) j \\ &+& (aw + bz - cy + dx) k \\ \end{matrix}$$

四元数 $q=a+bi+cj+dk$ に対してその 共役 $\overline{q}$ を

$$\overline{q} = a-bi-cj-dk$$

と定義します。

$q$と$\overline{q}$ の積は

$$\begin{array}{ccc} q\overline{q} &=& (aa + bb + cc + dd) \\ &+& (-ab + ba - cd + dc) i \\ &+& (-ac + bd + ca - db) j \\ &+& (-ad - bc + cb + da) k \\ &=& a^2 + b^2 + c^2 + d^2 \\ \end{array}$$

と計算することができ、これを使えば $q$ の逆元 $q^{-1}$ は

$$q^{-1} = \frac{\overline{q}}{q\overline{q}}$$

と書けることが分かります。

四元数$q$の四元数$\mu$による回転として以下のような操作も知られています。

$$q^\mu = \mu q \mu^{-1}$$

（四元数$q$の肩に$\mu$を書いていますがこれは累乗ではなく上記式によって定義される操作です）
この操作は実際に3次元空間で考る回転操作に対応していることが知られています（四元数と三次元空間における回転 | 高校数学の美しい物語）。3次元空間の回転を表すには座標軸周りの回転の組み合わせで表すオイラー角のような方法もありますが、こういった方法にはジンバルロックと呼ばれる回転の自由度が意図せず失われる現象が起こることが知られています。四元数を使った3次元空間の回転ではジンバルロックが起こらないため3DCGを始め様々な分野で活用されています。さらに二重数と同様の考え方で四元数を2つ組み合わせることで、回転だけでなく 6DoF（位置と姿勢（回転））を同時に表現できる双対四元数 という便利な数も知られています。

回転において特に $\mu$ が虚数単位 $i, j, k$ の何れかであれば（$\eta\in\{i,j,k\}$）

$$q^{\eta} = \eta q \eta^{-1} = \eta q \overline{\eta} = \eta q (-\eta) = -\eta q\eta$$

が成り立ちます。これを用いると

$$\begin{array}{ccl} q-q^i-q^j-q^k &=& q-(-iqi)-(-jqj)-(-kqk) \\ &=& (a+bi+cj+dk)-(a+bi-cj-dk)-(a-bi+cj-dk)-(a-bi-cj+dk) \\ &=& -2a +2bi +2cj +2dk &=& -2\overline{q} \end{array}$$

つまり

$$\overline{q} = \frac{1}{2}\left(-q+q^i+q^j+q^k\right)$$

となり四元数では共役を取る操作が虚数単位による回転を使って計算できることが分かります。

四元数解析

四元数を使ったNNを考えるために、四元数関数の四元数の意味における微分を考えてみましょう。上述のように四元数は積が非可換ですから以下のように2つの微分の定義を考えることができます。

$$\lim_{h\to 0}h^{-1}\left(f(z+h) - f(z)\right)$$

$$\lim_{h\to 0}\left(f(z+h) - f(z)\right)h^{-1}$$

実数、複素数の場合と同様、上記の極限は$h$を0に近づける方向によらず同じ値になる場合 前者は左微分、後者は右微分 と呼ばれています。

このような四元数の意味において微分可能な関数にはどのようなものがあるでしょうか？四元数関数 $f(x+yi+zj+wk)$ を

$$\begin{array}{ccl} f(x+yi+zj+wk) &=& a(x+yi+zj+wk) + b(x+yi+zj+wk)i + c(x+yi+zj+wk)j + d(x+yi+zj+wk)k \\ &=& a(x, y, z, w) + b(x, y, z, w)i + c(x, y, z, w)j + d(x, y, z, w)k \\ \end{array}$$

のように四元数を引数に取る実数値関数$a,b,c,d$を使って分解して考えてみましょう（今後の議論のために$a,b,c,d$は十分に滑らかであると仮定しておきます）。

$f$が左微分可能な場合、まず$x$方向の微分を考えると

$$\begin{array}{ccl} && \displaystyle\lim_{h\to 0}\left(\frac{a(x+h, y, z, w)-a(x,y,z,w)}{h}+\frac{b(x+h, y, z, w)-b(x,y,z,w)}{h}i+\frac{c(x+h, y, z, w)-c(x,y,z,w)}{h}j+\frac{d(x+h, y, z, w)-d(x,y,z,w)}{h}k\right) \\ &=& \frac{\partial a}{\partial x} + \frac{\partial b}{\partial x}i + \frac{\partial c}{\partial x}j + \frac{\partial d}{\partial x}k \\ \end{array}$$

となります。ここで$h$は四元数ではなく実数として考えていることに注意してください。
同様に$y$方向の微分を考えると

$$\begin{array}{ccl} && \displaystyle\lim_{h\to 0}\left(\frac{a(x, y+h, z, w)-a(x,y,z,w)}{ih}+\frac{b(x, y+h, z, w)-b(x,y,z,w)}{ih}i+\frac{c(x, y+h, z, w)-c(x,y,z,w)}{ih}j+\frac{d(x, y+h, z, w)-d(x,y,z,w)}{ih}k\right) \\ &=& \frac{\partial a}{\partial y}i^{-1} + \frac{\partial b}{\partial y}i^{-1}i + \frac{\partial c}{\partial y}i^{-1}j + \frac{\partial d}{\partial y}i^{-1}k \\ &=& -\frac{\partial a}{\partial y}i + \frac{\partial b}{\partial y} - \frac{\partial c}{\partial y}k + \frac{\partial d}{\partial y}j \\ \end{array}$$

となります。
同様に$z$方向の微分を考えると

$$\begin{array}{ccl} && \displaystyle\lim_{h\to 0}\left(\frac{a(x, y, z+h, w)-a(x,y,z,w)}{jh}+\frac{b(x, y, z+h, w)-b(x,y,z,w)}{jh}i+\frac{c(x, y, z+h, w)-c(x,y,z,w)}{jh}j+\frac{d(x, y, z+h, w)-d(x,y,z,w)}{jh}k\right) \\ &=& \frac{\partial a}{\partial z}j^{-1} + \frac{\partial b}{\partial z}j^{-1}i + \frac{\partial c}{\partial z}j^{-1}j + \frac{\partial d}{\partial z}j^{-1}k \\ &=& -\frac{\partial a}{\partial z}j + \frac{\partial b}{\partial z}k + \frac{\partial c}{\partial z} - \frac{\partial d}{\partial z}i \\ \end{array}$$

となります。
同様に$w$方向の微分を考えると

$$\begin{array}{ccl} && \displaystyle\lim_{h\to 0}\left(\frac{a(x, y, z, w+h)-a(x,y,z,w)}{kh}+\frac{b(x, y, z, w+h)-b(x,y,z,w)}{kh}i+\frac{c(x, y, z, w+h)-c(x,y,z,w)}{kh}j+\frac{d(x, y, z, w+h)-d(x,y,z,w+h)}{kh}k\right) \\ &=& \frac{\partial a}{\partial w}k^{-1} + \frac{\partial b}{\partial w}k^{-1}i + \frac{\partial c}{\partial w}k^{-1}j + \frac{\partial d}{\partial w}k^{-1}k \\ &=& -\frac{\partial a}{\partial w}k - \frac{\partial b}{\partial w}j + \frac{\partial c}{\partial w}i + \frac{\partial d}{\partial w} \\ \end{array}$$

となります。

四元数の意味で微分可能というのは上記4つの式が全て一致するということです。すなわち

$$\begin{matrix} \frac{\partial a}{\partial x} + \frac{\partial b}{\partial x}i + \frac{\partial c}{\partial x}j + \frac{\partial d}{\partial x}k &=& \frac{\partial b}{\partial y} - \frac{\partial a}{\partial y}i + \frac{\partial d}{\partial y}j - \frac{\partial c}{\partial y}k \\ &=& \frac{\partial c}{\partial z} - \frac{\partial d}{\partial z}i - \frac{\partial a}{\partial z}j + \frac{\partial b}{\partial z}k \\ &=& \frac{\partial d}{\partial w} + \frac{\partial c}{\partial w}i - \frac{\partial b}{\partial w}j - \frac{\partial a}{\partial w}k \\ \end{matrix}$$

が成り立ち、これより

$$\begin{matrix} \frac{\partial a}{\partial x} &=& \frac{\partial b}{\partial y} &=& \frac{\partial c}{\partial z} &=& \frac{\partial d}{\partial w} \\ \frac{\partial b}{\partial x} &=& -\frac{\partial a}{\partial y} &=& -\frac{\partial d}{\partial z} &=& \frac{\partial c}{\partial w} \\ \frac{\partial c}{\partial x} &=& \frac{\partial d}{\partial y} &=& -\frac{\partial a}{\partial z} &=& -\frac{\partial b}{\partial w} \\ \frac{\partial d}{\partial x} &=& -\frac{\partial c}{\partial y} &=& \frac{\partial b}{\partial z} &=& -\frac{\partial a}{\partial w} \\ \end{matrix}$$

が成り立つ必要があることが分かります。

ここで1番目の式を更に$y$で偏微分することにより

$$\frac{\partial^2 a}{\partial x\partial y} = \frac{\partial^2 c}{\partial y\partial z}$$

また4番目の式を更に$z$で偏微分することにより

$$\frac{\partial^2 c}{\partial y\partial z} = \frac{\partial^2 a}{\partial z\partial w}$$

よって

$$\frac{\partial^2 a}{\partial x\partial y} = \frac{\partial^2 a}{\partial z\partial w}$$

となります。同様に1番目の式を更に$y$で偏微分することから

$$\frac{\partial^2 a}{\partial x\partial y} = \frac{\partial^2 d}{\partial y\partial w}$$

また3番目の式を更に$w$で偏微分することから

$$\frac{\partial^2 d}{\partial y\partial w} = -\frac{\partial^2 a}{\partial z\partial w}$$

が分かり、よって

$$\frac{\partial^2 a}{\partial x\partial y} = -\frac{\partial^2 a}{\partial z\partial w}$$

となります。

以上より

$$\frac{\partial^2 a}{\partial x\partial y} = \frac{\partial^2 a}{\partial z\partial w} = -\frac{\partial^2 a}{\partial z\partial w} = 0$$

となることが分かります。

同様の操作を繰り返すことで異なる変数による2階の偏微分は全て0になることが分かります。

$$\begin{matrix} \frac{\partial^2 a}{\partial x\partial y} &=& \frac{\partial^2 a}{\partial x\partial z} &=& \frac{\partial^2 a}{\partial x\partial w} &=& \frac{\partial^2 a}{\partial y\partial z} &=& \frac{\partial^2 a}{\partial y\partial w} &=& \frac{\partial^2 a}{\partial z\partial w} &=& 0 \\ \frac{\partial^2 b}{\partial x\partial y} &=& \frac{\partial^2 b}{\partial x\partial z} &=& \frac{\partial^2 b}{\partial x\partial w} &=& \frac{\partial^2 b}{\partial y\partial z} &=& \frac{\partial^2 b}{\partial y\partial w} &=& \frac{\partial^2 b}{\partial z\partial w} &=& 0 \\ \frac{\partial^2 c}{\partial x\partial y} &=& \frac{\partial^2 c}{\partial x\partial z} &=& \frac{\partial^2 c}{\partial x\partial w} &=& \frac{\partial^2 c}{\partial y\partial z} &=& \frac{\partial^2 c}{\partial y\partial w} &=& \frac{\partial^2 c}{\partial z\partial w} &=& 0 \\ \frac{\partial^2 d}{\partial x\partial y} &=& \frac{\partial^2 d}{\partial x\partial z} &=& \frac{\partial^2 d}{\partial x\partial w} &=& \frac{\partial^2 d}{\partial y\partial z} &=& \frac{\partial^2 d}{\partial y\partial w} &=& \frac{\partial^2 d}{\partial z\partial w} &=& 0 \\ \end{matrix}$$

よって $a, b, c, d$ はそれぞれ $x, y, z, w$ のみに依存する関数と定数の和に綺麗に分解できることが分かります。

$$\begin{matrix} a(x, y, z, w) &=& a_1(x) + a_2(y) + a_3(z) + a_4(w) + C_a \\ b(x, y, z, w) &=& b_1(x) + b_2(y) + b_3(z) + b_4(w) + C_b \\ c(x, y, z, w) &=& c_1(x) + c_2(y) + c_3(z) + c_4(w) + C_c \\ d(x, y, z, w) &=& d_1(x) + d_2(y) + d_3(z) + d_4(w) + C_d \\ \end{matrix}$$

これをもう一度、一階偏微分の等式に代入してみましょう。

$$\begin{matrix} \frac{\partial a_1}{\partial x}(x) &=& \frac{\partial b_2}{\partial y}(y) &=& \frac{\partial c_3}{\partial z}(z) &=& \frac{\partial d_4}{\partial w}(w) \\ \frac{\partial b_1}{\partial x}(x) &=& -\frac{\partial a_2}{\partial y}(y) &=& -\frac{\partial d_3}{\partial z}(z) &=& \frac{\partial c_4}{\partial w}(w) \\ \frac{\partial c_1}{\partial x}(x) &=& \frac{\partial d_2}{\partial y}(y) &=& -\frac{\partial a_3}{\partial z}(z) &=& -\frac{\partial b_4}{\partial w}(w) \\ \frac{\partial d_1}{\partial x}(x) &=& -\frac{\partial c_2}{\partial y}(y) &=& \frac{\partial b_3}{\partial z}(z) &=& -\frac{\partial a_4}{\partial w}(w) \\ \end{matrix}$$

これは変数が違うにも関わらず等号が成立しているので、それぞれが変数に依存しない異なる定数 $K_1, K_2, K_3, K_4$ であることを示しています。

よって

$$\begin{matrix} a(x, y, z, w) &=& K_1x - K_2y - K_3z - K_4w + C_a \\ b(x, y, z, w) &=& K_2x + K_1y + K_4z - K_3w + C_b \\ c(x, y, z, w) &=& K_3x - K_4y + K_1z + K_2w + C_c \\ d(x, y, z, w) &=& K_4x + K_3y - K_2z + K_1w + C_d \\ \end{matrix}$$

となり、これは $f$ が

$$f(x+yi+zj+wk) = (x+yi+zj+wk)(K_1+K_2i+K_3j+K_4k)+(C_a+C_bi+C_cj+C_dk)$$

と表されることを示しています。

つまり四元数関数は左微分可能だと仮定すると定数もしくは線形関数に限定されてしまうのです。右微分可能性についても同様の議論が成立します。これは四元数NNを考えるには明らかに強すぎる制限でしょう。

ここでの議論は HADDAD, Z. Two remarks on the quaternions. Pi Mu Epsilon Journal, 1981, 7.4: 221-231. を参考にしています。

実ベクトルとしての四元数による微分

四元数関数を四元数の意味で微分しようとすると単純な関数しか扱えない事がわかりました。そこで複素数の場合と同様に多元環として四元数$q$を四次元の実ベクトル$x+yi+zj+wk$だと思って微分することを考えてみましょう。

四元数を引数に取る実数値関数 $f$について四元数$q$とそれを$\Delta q = \Delta x+\Delta yi+\Delta zj+\Delta wk$ だけ微小変化させた $q+\Delta q$ における評価値を考えてテイラー展開すると

$$f(q+\Delta q) = f(q) + \frac{\partial f}{\partial x}\Delta x + \frac{\partial f}{\partial y}\Delta y + \frac{\partial f}{\partial z}\Delta z + \frac{\partial f}{\partial w}\Delta w + o(\Delta x, \Delta y, \Delta z, \Delta w)$$

となり小さな正の定数 $\eta$ を使って $\Delta x = -\eta\frac{\partial f}{\partial x}, \Delta y = -\eta\frac{\partial f}{\partial y}, \Delta z = -\eta\frac{\partial f}{\partial z}, \Delta w = -\eta\frac{\partial f}{\partial w}$ とする、すなわち

$$\begin{matrix} \Delta q &=& -\eta\frac{\partial f}{\partial x} - \eta\frac{\partial f}{\partial y}i - \eta\frac{\partial f}{\partial z}j - \eta\frac{\partial f}{\partial w}k \\ &=& -\eta\left(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}i + \frac{\partial f}{\partial z}j + \frac{\partial f}{\partial w}k\right) \\ \end{matrix}$$

とおくと

$$f(q+\Delta q) - f(q) = -\eta\left(\left|\frac{\partial f}{\partial x}\right|^2 + \left|\frac{\partial f}{\partial y}\right|^2 + \left|\frac{\partial f}{\partial z}\right|^2 + \left|\frac{\partial f}{\partial w}\right|^2\right) + o(\Delta x, \Delta y, \Delta z, \Delta w)$$

となり $\eta$が十分に小さければ右辺は負の実数になることから、左辺を小さくできる ことが分かります。

ここで $\Delta q$ に現れる微分演算子の部分を使って、四元数$q$による微分演算子 $\frac{\partial}{\partial \overline{q}}$ を

$$\frac{\partial f}{\partial \overline{q}} = \frac{1}{4}\left(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}i + \frac{\partial f}{\partial z}j + \frac{\partial f}{\partial w}k\right)$$

と定義します（ここで四元数は積が非可換なので$i,j,k$の位置を変えた"右微分"も同様に考えられることに注意）。

係数 $\frac{1}{4}$ と微分する変数が共役四元数になっている訳は $f(q) = c\overline{q}$ という線形関数の微分を考えると

$$\begin{matrix} \frac{\partial f}{\partial \overline{q}} &=& \frac{1}{4}\left(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}i + \frac{\partial f}{\partial z}j + \frac{\partial f}{\partial w}k\right) \\ &=& \frac{1}{4}\left(\frac{\partial c\overline{q}}{\partial x} + \frac{\partial c\overline{q}}{\partial y}i + \frac{\partial c\overline{q}}{\partial z}j + \frac{\partial c\overline{q}}{\partial w}k\right) \\ &=& \frac{1}{4}\left(\frac{\partial c(x-yi-zj-wk)}{\partial x} + \frac{\partial c(x-yi-zj-wk)}{\partial y}i + \frac{\partial c(x-yi-zj-wk)}{\partial z}j + \frac{\partial c(x-yi-zj-wk)}{\partial w}k\right) \\ &=& \frac{1}{4}\left(c + (c(-i))i - (c(-j))j - (c(-k))k\right) \\ &=& \frac{1}{4}\left(c + c + c + c\right) \\ &=& c \\ \end{matrix}$$

と自然に係数が出てくるように設定されています。

このように定義した微分演算子 $\frac{\partial}{\partial \overline{q}}$ はどのような性質を満たすでしょうか。

まず線形性については四元数$a$と四元数関数$f, g$について

$$\frac{\partial (af)}{\partial \overline{q}} = a\frac{\partial f}{\partial \overline{q}}$$

$$\frac{\partial (f+g)}{\partial \overline{q}} = \frac{\partial f}{\partial \overline{q}} + \frac{\partial g}{\partial \overline{q}}$$

が成り立ちます（ただし左微分を考えているので$\frac{\partial (fa)}{\partial\overline{q}} = \frac{\partial f}{\partial\overline{q}}a$ が一般には成り立たないことに注意）。

それでは積の微分法についてはどうでしょうか？実は四元数の積が非可換であることから今までのような単純な積の微分法は成り立ちません。このことを簡単な例で見てみましょう。$q\overline{q}$ の微分を考えると

$$\begin{matrix} \frac{\partial q\overline{q}}{\partial \overline{q}} &=& \frac{1}{4}\left(\frac{\partial (q\overline{q})}{\partial x} + \frac{\partial (q\overline{q})}{\partial y}i + \frac{\partial (q\overline{q})}{\partial z}j + \frac{\partial (q\overline{q})}{\partial w}k\right) \\ &=& \frac{1}{4}\left(\frac{\partial \left(x^2+y^2+z^2+z^2\right)}{\partial x} + \frac{\partial \left(x^2+y^2+z^2+z^2\right)}{\partial y}i + \frac{\partial \left(x^2+y^2+z^2+z^2\right)}{\partial z}j + \frac{\partial \left(x^2+y^2+z^2+z^2\right)}{\partial w}k\right) \\ &=& \frac{1}{4}\left(2x + 2yi + 2zj + 2wk\right) \\ &=& \frac{1}{2}q \\ \end{matrix}$$

となりますが、

$$\begin{matrix} \frac{\partial q}{\partial \overline{q}} &=& \frac{1}{4}\left(\frac{\partial q}{\partial x} + \frac{\partial q}{\partial y}i + \frac{\partial q}{\partial z}j + \frac{\partial q}{\partial w}k\right) \\ &=& \frac{1}{4}\left(\frac{\partial (x+yi+zj+wk)}{\partial x} + \frac{\partial (x+yi+zj+wk)}{\partial y}i + \frac{\partial (x+yi+zj+wk)}{\partial z}j + \frac{\partial (x+yi+zj+wk)}{\partial w}k\right) \\ &=& \frac{1}{4}(1+ii+jj+kk) \\ &=& -\frac{1}{2} \\ \end{matrix}$$

$$\begin{matrix} \frac{\partial \overline{q}}{\partial \overline{q}} &=& 1 \\ \end{matrix}$$

より

$$\frac{\partial q}{\partial \overline{q}}\overline{q} + q\frac{\partial \overline{q}}{\partial \overline{q}} = -\frac{1}{2}\overline{q} + q \neq \frac{1}{2}q = \frac{\partial q\overline{q}}{\partial \overline{q}}$$

となり単純な積の微分法が成立しないことが分かります。

実は同様に単純な方法では合成関数の微分法も成立しないようになっています。これでは誤差逆伝播法が使えないので四元数NNを学習させることが出来ません。

GHR微積分

実は前述の四元数による微分 $\frac{\partial}{\partial\overline{q}}$ を少し拡張することで、積の微分法と合成関数の微分法が考えられる微分法が知られています。それがGHR微積分です。

XU, Dongpo, et al. Enabling quaternion derivatives: The generalized HR calculus. Royal Society open science, 2015, 2.8: 150255.

上記論文において（左）GHC微積分は以下のような微分演算子によって定義されています。

$$\frac{\partial f}{\partial q^\mu} = \frac{1}{4}\left(\frac{\partial f}{\partial x} - \frac{\partial f}{\partial y}i^\mu - \frac{\partial f}{\partial z}j^\mu - \frac{\partial f}{\partial w}k^\mu\right)$$

$$\frac{\partial f}{\partial \overline{q^\mu}} = \frac{1}{4}\left(\frac{\partial f}{\partial x} + \frac{\partial f}{\partial y}i^\mu + \frac{\partial f}{\partial z}j^\mu + \frac{\partial f}{\partial w}k^\mu\right)$$

ここで $q^\mu$ は回転であり $\mu$ は0ではない四元数です。

GHR微積分では以下のような積の微分法と合成関数の微分法が成り立ちます。

積の微分法

$$\frac{\partial (fg)}{\partial q^\mu} = f\frac{\partial g}{\partial q^\mu} + \frac{\partial f}{\partial q^{g\mu}}g$$

$$\frac{\partial (fg)}{\partial \overline{q^\mu}} = f\frac{\partial g}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial \overline{q^{g\mu}}}g$$

合成関数の微分法

$$\frac{\partial (f\circ g)}{\partial q^\mu} = \frac{\partial f}{\partial g^\nu}\frac{\partial g^\nu}{\partial q^\mu} + \frac{\partial f}{\partial g^{\nu i}}\frac{\partial g^{\nu i}}{\partial q^\mu} + \frac{\partial f}{\partial g^{\nu j}}\frac{\partial g^{\nu j}}{\partial q^\mu} + \frac{\partial f}{\partial g^{\nu k}}\frac{\partial g^{\nu k}}{\partial q^\mu}$$

$$\frac{\partial (f\circ g)}{\partial \overline{q^\mu}} = \frac{\partial f}{\partial g^\nu}\frac{\partial g^\nu}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial g^{\nu i}}\frac{\partial g^{\nu i}}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial g^{\nu j}}\frac{\partial g^{\nu j}}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial g^{\nu k}}\frac{\partial g^{\nu k}}{\partial \overline{q^\mu}}$$

$$\frac{\partial (f\circ g)}{\partial q^\mu} = \frac{\partial f}{\partial \overline{g^\nu}}\frac{\partial \overline{g^\nu}}{\partial q^\mu} + \frac{\partial f}{\partial \overline{g^{\nu i}}}\frac{\partial \overline{g^{\nu i}}}{\partial q^\mu} + \frac{\partial f}{\partial \overline{g^{\nu j}}}\frac{\partial \overline{g^{\nu j}}}{\partial q^\mu} + \frac{\partial f}{\partial \overline{g^{\nu k}}}\frac{\partial \overline{g^{\nu k}}}{\partial q^\mu}$$

$$\frac{\partial (f\circ g)}{\partial \overline{q^\mu}} = \frac{\partial f}{\partial \overline{g^\nu}}\frac{\partial \overline{g^\nu}}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial \overline{g^{\nu i}}}\frac{\partial \overline{g^{\nu i}}}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial \overline{g^{\nu j}}}\frac{\partial \overline{g^{\nu j}}}{\partial \overline{q^\mu}} + \frac{\partial f}{\partial \overline{g^{\nu k}}}\frac{\partial \overline{g^{\nu k}}}{\partial \overline{q^\mu}}$$

ここで $\nu$ は0ではない任意の四元数です。

証明は論文を参照していただければと思いますが、一例として積の微分法が先程の例で成り立っていることを見てみましょう。GHR微積分の積の微分法によると $q\overline{q}$ の微分は

$$\begin{array}{ccc} \frac{\partial q\overline{q}}{\partial \overline{q}} &=& q\frac{\partial \overline{q}}{\partial \overline{q}} + \frac{\partial q}{\partial \overline{q^{\overline{q}}}}\overline{q} \\ &=& q + \frac{1}{4}\left(\frac{\partial q}{\partial x} + \frac{\partial q}{\partial y}i^{\overline{q}} + \frac{\partial q}{\partial z}j^{\overline{q}} + \frac{\partial q}{\partial w}k^{\overline{q}}\right)\overline{q} \\ &=& q + \frac{1}{4}\left(1 + i\overline{q}i\overline{q}^{-1} + j\overline{q}j\overline{q}^{-1} + k\overline{q}k\overline{q}^{-1}\right)\overline{q} \\ &=& q + \frac{1}{4}\left(\overline{q} + i\overline{q}i + j\overline{q}j + k\overline{q}k\right) \\ &=& q + \frac{1}{4}\left(-2 \overline{\overline{q}}\right) \\ &=& q - \frac{1}{2}q \\ &=& \frac{1}{2}q \\ \end{array}$$

のように計算することができ積の微分法を用いない直接の計算と結果が一致することが分かります。

GHR微積分には右微分も存在しますが、実数値を取る四元数関数の場合は左微分と右微分は一致します。本稿ではNNの誤差関数という実数値関数しか考えないため今後の実装では左微分のみを考えます。さらに四元数$q$による積の微分法や合成関数の微分法には共役四元数$\overline{q}$による微分の項が現れないため、今後の実装では$\frac{\partial}{\partial q^\mu}$のみを考えます。

他にも平均値の定理が成り立ったり、関数の虚数単位 $j, k$ の項が0であるとき（つまり複素関数だと考えられるとき）はウィルティンガー微分に一致したりと色々面白い性質が成り立つので興味があれば元論文も是非読んでみてください。

四元数ニューラルネットワーク

それでは自動微分のテクニックを用いてGHR微分を実装してみましょう。名付けてGHR自動微分です。

まずは四元数を定義しておきます。

import Linear.Quaternion (Quaternion(..))
import Linear.V3 (V3(..))

-- | 四元数
type H = Quaternion Double

-- | 1と虚数単位
e, i, j, k :: H
e = Quaternion 1 (V3 0 0 0)
i = Quaternion 0 (V3 1 0 0)
j = Quaternion 0 (V3 0 1 0)
k = Quaternion 0 (V3 0 0 1)

-- | 四元数のスカラー倍
(*|) :: Double -> H -> H
a *| q = fmap (*a) q

五重四元数

$\frac{\partial}{\partial q}$の積の微分法や合成関数の微分法を計算するためには一般の四元数 $\mu$ に関する回転による微分 $\frac{\partial}{\partial q^\mu}$ を計算する必要があります。これは以下のように $\frac{\partial}{\partial q}, \frac{\partial}{\partial q^i}, \frac{\partial}{\partial q^j}, \frac{\partial}{\partial q^k}$ の値から再現することができます。

定義から

$$\begin{matrix} \frac{\partial f}{\partial q} &=& \frac{1}{4}\left(\frac{\partial f}{\partial x}-\frac{\partial f}{\partial y}i-\frac{\partial f}{\partial z}j-\frac{\partial f}{\partial w}k\right) \\ \frac{\partial f}{\partial q^i} &=& \frac{1}{4}\left(\frac{\partial f}{\partial x}-\frac{\partial f}{\partial y}i+\frac{\partial f}{\partial z}j+\frac{\partial f}{\partial w}k\right) \\ \frac{\partial f}{\partial q^j} &=& \frac{1}{4}\left(\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}i-\frac{\partial f}{\partial z}j+\frac{\partial f}{\partial w}k\right) \\ \frac{\partial f}{\partial q^k} &=& \frac{1}{4}\left(\frac{\partial f}{\partial x}+\frac{\partial f}{\partial y}i+\frac{\partial f}{\partial z}j-\frac{\partial f}{\partial w}k\right) \\ \end{matrix}$$

なので

$$\begin{matrix} \frac{\partial f}{\partial x} &=& \frac{\partial f}{\partial q} + \frac{\partial f}{\partial q^i} + \frac{\partial f}{\partial q^j} + \frac{\partial f}{\partial q^k} \\ \frac{\partial f}{\partial y} &=& \left(\frac{\partial f}{\partial q} + \frac{\partial f}{\partial q^i} - \frac{\partial f}{\partial q^j} - \frac{\partial f}{\partial q^k}\right)i \\ \frac{\partial f}{\partial z} &=& \left(\frac{\partial f}{\partial q} - \frac{\partial f}{\partial q^i} + \frac{\partial f}{\partial q^j} - \frac{\partial f}{\partial q^k}\right)j \\ \frac{\partial f}{\partial w} &=& \left(\frac{\partial f}{\partial q} - \frac{\partial f}{\partial q^i} - \frac{\partial f}{\partial q^j} + \frac{\partial f}{\partial q^k}\right)k \\ \end{matrix}$$

と計算できます。GHR微積分の定義から $\mu = \mu_a+\mu_bi+\mu_cj+\mu_dk$ と置くと

$$\begin{array}{ccl} 4\frac{\partial f}{\partial q^\mu} &=& \frac{\partial f}{\partial x} - \frac{\partial f}{\partial y}i^\mu - \frac{\partial f}{\partial z}j^\mu - \frac{\partial f}{\partial w}k^\mu \\ &=& \frac{\partial f}{\partial q} + \frac{\partial f}{\partial q^i} + \frac{\partial f}{\partial q^j} + \frac{\partial f}{\partial q^k} \\ &-& \left(\frac{\partial f}{\partial q} + \frac{\partial f}{\partial q^i} - \frac{\partial f}{\partial q^j} - \frac{\partial f}{\partial q^k}\right)i\mu i\mu^{-1} \\ &-& \left(\frac{\partial f}{\partial q} - \frac{\partial f}{\partial q^i} + \frac{\partial f}{\partial q^j} - \frac{\partial f}{\partial q^k}\right)j\mu j\mu^{-1} \\ &-& \left(\frac{\partial f}{\partial q} - \frac{\partial f}{\partial q^i} - \frac{\partial f}{\partial q^j} + \frac{\partial f}{\partial q^k}\right)k\mu k^{-1} \\ &=& \frac{\partial f}{\partial q}\left(1-i\mu i\mu^{-1}-j\mu j\mu^{-1}-k\mu k\mu^{-1}\right) \\ &-& \frac{\partial f}{\partial q^i}\left(1-i\mu i\mu^{-1}+j\mu j\mu^{-1}+k\mu k\mu^{-1}\right) \\ &-& \frac{\partial f}{\partial q^i}\left(1+i\mu i\mu^{-1}-j\mu j\mu^{-1}+k\mu k\mu^{-1}\right) \\ &-& \frac{\partial f}{\partial q^i}\left(1+i\mu i\mu^{-1}+j\mu j\mu^{-1}-k\mu k\mu^{-1}\right) \\ &=& \frac{\partial f}{\partial q}\left(\mu-i\mu i-j\mu j-k\mu k\right)\mu^{-1} \\ &-& \frac{\partial f}{\partial q^i}\left(\mu-i\mu i+j\mu j+k\mu k\right)\mu^{-1} \\ &-& \frac{\partial f}{\partial q^i}\left(\mu+i\mu i-j\mu j+k\mu k\right)\mu^{-1} \\ &-& \frac{\partial f}{\partial q^i}\left(\mu+i\mu i+j\mu j-k\mu k\right)\mu^{-1} \\ &=& \frac{\partial f}{\partial q}\left(2\overline{\mu}+2\mu\right)\mu^{-1} \\ &-& \frac{\partial f}{\partial q^i}\left(-2\overline{\mu}-2i\mu i\right)\mu^{-1} \\ &-& \frac{\partial f}{\partial q^j}\left(-2\overline{\mu}-2j\mu j\right)\mu^{-1} \\ &-& \frac{\partial f}{\partial q^k}\left(-2\overline{\mu}-2k\mu k\right)\mu^{-1} \\ &=& 2\left(\frac{\partial f}{\partial q}\left(\overline{\mu}+\mu\right) + \frac{\partial f}{\partial q^i}\left(\overline{\mu}+i\mu i\right) + \frac{\partial f}{\partial q^j}\left(\overline{\mu}+j\mu j\right) + \frac{\partial f}{\partial q^k}\left(\overline{\mu}+k\mu k\right)\right)\mu^{-1} \\ &=& 4\left(\frac{\partial f}{\partial q}\mu_a + \frac{\partial f}{\partial q^i}(\mu_a+\mu_bi) + \frac{\partial f}{\partial q^j}(\mu_a+\mu_cj) + \frac{\partial f}{\partial q^k}(\mu_a+\mu_dk)\right)\mu^{-1} \\ \end{array}$$

と計算することができます。
よって $\frac{\partial f}{\partial q^\mu}$ は $\frac{\partial}{\partial q}, \frac{\partial}{\partial q^i}, \frac{\partial}{\partial q^j}, \frac{\partial}{\partial q^k}$ を用いて

$$\frac{\partial f}{\partial q^\mu} = \left(\frac{\partial f}{\partial q}\mu_a + \frac{\partial f}{\partial q^i}(\mu_a+\mu_bi) + \frac{\partial f}{\partial q^j}(\mu_a+\mu_cj) + \frac{\partial f}{\partial q^k}(\mu_a+\mu_dk)\right)\mu^{-1}$$

と表すことができます。

以上より関数の値$f$に加えて$\frac{\partial}{\partial q}, \frac{\partial}{\partial q^i}, \frac{\partial}{\partial q^j}, \frac{\partial}{\partial q^k}$の値を保持しておけば自動微分が実装できそうです。これを実現するために5つの四元数を要素として持つ数を定義しましょう。名付けて五重四元数です。

-- | 五重四元数
data Q = Q H H H H H deriving Show

quaternion, gradQ, gradQi, gradQj, gradQk :: Q -> H
quaternion (Q a b c d e) = a
gradQ      (Q a b c d e) = b
gradQi     (Q a b c d e) = c
gradQj     (Q a b c d e) = d
gradQk     (Q a b c d e) = e

Num, Fractional

五重四元数Qに対してNumのインスタンスを実装するには積の微分法と合成関数の微分法を計算する必要があります。積の微分法はこれまでの議論で実装できそうですが、合成関数の微分法については $\frac{\partial g^\nu}{\partial q^\mu}$ を $\frac{\partial g}{\partial q}, \frac{\partial g}{\partial q^i}, \frac{\partial g}{\partial q^j}, \frac{\partial g}{\partial q^k}$ からどのように計算するかが問題になりそうです。

一般に $\mu$ による回転の微分の定数 $\nu$ による回転は以下のように計算することができます。

$$\begin{array}{ccl} \left(\frac{\partial f}{\partial q^\mu}\right)^\nu &=& \frac{1}{4}\left(\frac{\partial f}{\partial x} - \frac{\partial f}{\partial y}i^\mu - \frac{\partial f}{\partial z}j^\mu - \frac{\partial f}{\partial w}k^\mu\right)^\nu \\ &=& \frac{1}{4}\left(\nu\frac{\partial f}{\partial x}\nu^{-1} - \nu\frac{\partial f}{\partial y}\nu^{-1}\nu\mu i\mu^{-1}\nu^{-1} - \nu\frac{\partial f}{\partial z}\nu^{-1}\nu\mu j\mu^{-1}\nu^{-1} - \nu\frac{\partial f}{\partial w}\nu^{-1}\nu\mu k\mu^{-1}\nu^{-1}\right)^\nu \\ &=& \frac{1}{4}\left(\frac{\partial f^\nu}{\partial x} - \frac{\partial f^\nu}{\partial y}i^{\nu\mu} - \frac{\partial f^\nu}{\partial z}j^{\nu\mu} - \frac{\partial f^\nu}{\partial w}k^{\nu\mu}\right) \\ &=& \left(\frac{\partial f^\nu}{\partial q^{\nu\mu}}\right) \end{array}$$

（2行目から3行目の変形は$\nu$が定数で$x,y,z,w$に依存しないことを前提としていることに注意）

これを使えば、合成関数の微分法は（$\nu$は任意なので$\nu=1$とすると）

$$\begin{matrix} \frac{\partial (f\circ g)}{\partial q^\mu} &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q^\mu} + \frac{\partial f}{\partial g^i}\frac{\partial g^i}{\partial q^\mu} + \frac{\partial f}{\partial g^j}\frac{\partial g^j}{\partial q^\mu} + \frac{\partial f}{\partial g^k}\frac{\partial g^k}{\partial q^\mu} \\ &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q^\mu} + \frac{\partial f}{\partial g^i}\left(\frac{\partial g}{\partial q^{i\mu}}\right)^i + \frac{\partial f}{\partial g^j}\left(\frac{\partial g}{\partial q^{j\mu}}\right)^j + \frac{\partial f}{\partial g^k}\left(\frac{\partial g}{\partial q^{k\mu}}\right)^k \\ &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q^\mu} - \frac{\partial f}{\partial g^i}i\frac{\partial g}{\partial q^{i\mu}}i - \frac{\partial f}{\partial g^j}j\frac{\partial g}{\partial q^{j\mu}}j - \frac{\partial f}{\partial g^k}k\frac{\partial g}{\partial q^{k\mu}}k \\ \end{matrix}$$

と書けるので

$$\begin{matrix} \frac{\partial (f\circ g)}{\partial q} &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q} - \frac{\partial f}{\partial g^i}i\frac{\partial g}{\partial q^i}i - \frac{\partial f}{\partial g^j}j\frac{\partial g}{\partial q^j}j - \frac{\partial f}{\partial g^k}k\frac{\partial g}{\partial q^k}k \\ \frac{\partial (f\circ g)}{\partial q^i} &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q^i} - \frac{\partial f}{\partial g^i}i\frac{\partial g}{\partial q}i - \frac{\partial f}{\partial g^j}j\frac{\partial g}{\partial q^k}j - \frac{\partial f}{\partial g^k}k\frac{\partial g}{\partial q^j}k \\ \frac{\partial (f\circ g)}{\partial q^j} &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q^j} - \frac{\partial f}{\partial g^i}i\frac{\partial g}{\partial q^k}i - \frac{\partial f}{\partial g^j}j\frac{\partial g}{\partial q}j - \frac{\partial f}{\partial g^k}k\frac{\partial g}{\partial q^i}k \\ \frac{\partial (f\circ g)}{\partial q^k} &=& \frac{\partial f}{\partial g}\frac{\partial g}{\partial q^k} - \frac{\partial f}{\partial g^i}i\frac{\partial g}{\partial q^j}i - \frac{\partial f}{\partial g^j}j\frac{\partial g}{\partial q^i}j - \frac{\partial f}{\partial g^k}k\frac{\partial g}{\partial q}k \\ \end{matrix}$$

のように計算することができます。

以上からQのNumのインスタンスは以下のように実装することができます。

-- | ∂/∂q^{,i,j,k} を使った ∂/∂q^μ の計算
ghr :: H -> H -> H -> H -> H -> H
ghr mu@(Quaternion mu_a (V3 mu_b mu_c mu_d)) q' qi' qj' qk' =
  (q' * a + qi' * b + qj' * c + qk' * d) * recip mu
  where
    a = mu_a *| e
    b = mu_a *| e + mu_b *| i
    c = mu_a *| e + mu_c *| j
    d = mu_a *| e + mu_d *| k

-- | GHR微積分における合成関数の微分法
ghrChainRule :: (H, H, H, H) -> (H, H, H, H) -> (H, H, H, H)
ghrChainRule (dfdq, dfdqi, dfdqj, dfdqk) (q', qi', qj', qk') =
  let p'  = dfdq * q'  - dfdqi * i * qi' * i - dfdqj * j * qj' * j - dfdqk * k * qk' * k
      pi' = dfdq * qi' - dfdqi * i * q'  * i - dfdqj * j * qk' * j - dfdqk * k * qj' * k
      pj' = dfdq * qj' - dfdqi * i * qk' * i - dfdqj * j * q'  * j - dfdqk * k * qi' * k
      pk' = dfdq * qk' - dfdqi * i * qj' * i - dfdqj * j * qi' * j - dfdqk * k * q'  * k
   in (p',pi', pj', pk')

instance Num Q where
  (Q q q' qi' qj' qk') + (Q p p' pi' pj' pk') = Q (q + p) (q' + p') (qi' + pi') (qj' + pj') (qk' + pk')
  (Q q q' qi' qj' qk') * (Q p p' pi' pj' pk') = Q (q * p) r' ri' rj' rk' where
    r'  = q * p'  + ghr  p      q' qi' qj' qk' * p
    ri' = q * pi' + ghr (p * i) q' qi' qj' qk' * p
    rj' = q * pj' + ghr (p * j) q' qi' qj' qk' * p
    rk' = q * pk' + ghr (p * k) q' qi' qj' qk' * p
  negate (Q q q' qi' qj' qk') = Q (negate q) (negate q') (negate qi') (negate qj') (negate qk')
  abs (Q q q' qi' qj' qk') = Q (abs q) p' pi' pj' pk' where
    dfdq  = conjugate q / (4 * abs q)
    dfdqi = -i * dfdq * i
    dfdqj = -j * dfdq * j
    dfdqk = -k * dfdq * k
    (p',pi', pj', pk') = ghrChainRule (dfdq, dfdqi, dfdqj, dfdqk) (q', qi', qj', qk')
  signum (Q q q' qi' qj' qk') = Q (signum q) 0 0 0 0
  fromInteger n = Q (fromInteger n) 0 0 0 0

* の実装に積の微分法、absの実装に合成関数の微分法を用いています。absは数式で書くと

$$|q| = \sqrt{x^2+y^2+z^2+w^2}$$

となることより、$\frac{\partial|q|}{\partial q^\mu}$ は

$$\begin{matrix} \frac{\partial|q|}{\partial q^\mu} &=& \frac{1}{4}\left(\frac{\partial|q|}{\partial x} - \frac{\partial|q|}{\partial y}i^\mu - \frac{\partial|q|}{\partial z}j^\mu - \frac{\partial|q|}{\partial w}k^\mu\right) \\ &=& \frac{1}{4}\left(\frac{x}{|q|} - \frac{y}{|q|}i^\mu - \frac{z}{|q|}j^\mu - \frac{w}{|q|}k^\mu\right) \\ &=& \mu\frac{1}{4}\left(\frac{x}{|q|} - \frac{y}{|q|}i - \frac{z}{|q|}j - \frac{w}{|q|}k\right)\mu^{-1} \\ &=& \mu\left(\frac{1}{4}\frac{\overline{q}}{|q|}\right)\mu^{-1} \\ \end{matrix}$$

と計算できるので、あとは合成関数の微分法と組み合わせて実装しています。

同様に Fractional のインスタンスは以下のように実装することができます。

instance Fractional Q where
  recip (Q q q' qi' qj' qk') = Q (recip q) p' pi' pj' pk' where
    dfdq_mu mu = -(recip q) * (view _e (recip q * mu) *| recip mu)
    dfdq  = dfdq_mu e
    dfdqi = dfdq_mu i
    dfdqj = dfdq_mu j
    dfdqk = dfdq_mu k
    (p',pi', pj', pk') = ghrChainRule (dfdq, dfdqi, dfdqj, dfdqk) (q', qi', qj', qk')
  fromRational r = Q (fromRational r) 0 0 0 0

まず一般に $\frac{\partial q}{\partial q^\mu}$ は

$$\begin{matrix} \frac{\partial q}{\partial q^\mu} &=& \frac{1}{4}\left(\frac{\partial q}{\partial x} - \frac{\partial q}{\partial y}i^\mu - \frac{\partial q}{\partial z}j^\mu - \frac{\partial q}{\partial w}k^\mu\right) \\ &=& \frac{1}{4}\left(1 - ii^\mu - jj^\mu - kk^\mu\right) \\ &=& \frac{1}{4}\left(\mu - i\mu i - j\mu j - k\mu k\right)\mu^{-1} \\ &=& \frac{1}{4}\left(2\mu + 2\overline{\mu}\right)\mu^{-1} \\ &=& \Re(\mu)\mu^{-1} \\ \end{matrix}$$

（ここで $\Re(\mu)$ は $\mu$ のスカラー部を表す記号とします）
と計算することができます。

$\frac{\partial q^{-1}}{\partial q}$ を計算するために式

$$1 = qq^{-1}$$

を考えます。両辺を微分すると、

$$\begin{matrix} 0 &=& \frac{\partial qq^{-1}}{\partial q^\mu} \\ &=& q\frac{\partial q^{-1}}{\partial q^\mu} + \frac{\partial q}{\partial q^{q^{-1}\mu}}q^{-1} \\ &=& q\frac{\partial q^{-1}}{\partial q^\mu} + \frac{\partial q}{\partial q^{q^{-1}\mu}}q^{-1}\mu\mu^{-1} \\ &=& q\frac{\partial q^{-1}}{\partial q^\mu} + \Re\left({q^{-1}\mu}\right)\mu^{-1} \\ \end{matrix}$$

となることより $\frac{\partial q^{-1}}{\partial q}$ は

$$\begin{matrix} \frac{\partial q^{-1}}{\partial q^\mu} = -q^{-1}\Re\left({q^{-1}\mu}\right)\mu^{-1} \end{matrix}$$

と計算できることが分かりました。Fractional の recip の実装はこの計算に合成関数の微分法を組み合わせた形になっています。

四元数指数関数

次に Floating の実装について考えてみましょう。Floating を実装するには四元数の三角関数 sin, cos, tan や指数関数／対数関数 exp, log などを実装する必要があります。

四元数における指数関数 $\exp$ は

$$\exp(q) = \sum_{n=0}^\infty \frac{1}{n!}q^n$$

と定義されます。これは $q = a + bi + cj + dk$ とすると $a$ は他の項と可換なので

$$\exp(q) = \exp(a)\exp(bi + cj + dk)$$

と分解できます。ここで ${\rm v} = bi + cj + dk$ と置き

$${\rm v}^2 = (bi + cj + dk)^2 = -b^2-c^2-d^2 = -|{\rm v}|^2$$

に注意すると

$$\begin{matrix} {\rm v}^2 &=& -|{\rm v}|^2 \\ {\rm v}^3 &=& -|{\rm v}|^2{\rm v} \\ {\rm v}^4 &=& |{\rm v}|^4 \\ {\rm v}^5 &=& |{\rm v}|^4{\rm v} \\ {\rm v}^6 &=& -|{\rm v}|^6 \\ {\rm v}^7 &=& -|{\rm v}|^6{\rm v} \\ \vdots \end{matrix}$$

より

$$\begin{matrix} \exp({\rm v}) &=& \sum_{n=0}^\infty \frac{1}{n!}{\rm v}^n \\ &=& 1 + \frac{\rm v}{1!} - \frac{|{\rm v}|^2}{2!} - \frac{|{\rm v}|^2{\rm v}}{3!} + \frac{|{\rm v}|^4}{4!} + \frac{|{\rm v}|^4{\rm v}}{5!} - \frac{|{\rm v}|^6}{6!} - \cdots \\ &=& 1 + \frac{|{\rm v}|\rm v}{1!|{\rm v}|} - \frac{|{\rm v}|^2}{2!} - \frac{|{\rm v}|^3{\rm v}}{3!|{\rm v}|} + \frac{|{\rm v}|^4}{4!} + \frac{|{\rm v}|^5{\rm v}}{5!|{\rm v}|} - \frac{|{\rm v}|^6}{6!} - \cdots \\ &=& \left(1 - \frac{|{\rm v}|^2}{2!} + \frac{|{\rm v}|^4}{4!} - \frac{|{\rm v}|^6}{6!} + \cdots\right) + \frac{\rm v}{|{\rm v}|}\left(\frac{|{\rm v}|}{1!} - \frac{|{\rm v}|^3}{3!} + \frac{|{\rm v}|^5}{5!} - \cdots\right) \\ &=& \cos(|{\rm v}|) + \frac{\rm v}{|{\rm v}|}\sin(|{\rm v}|) \end{matrix}$$

と計算できるので、以上より四元数の指数関数は

$$\exp(q) = \exp(a)\left(\cos(|{\rm v}|) + \frac{\rm v}{|{\rm v}|}\sin(|{\rm v}|)\right)$$

のように計算できることが分かりました。

同様の級数を用いた考え方を使って四元数の三角関数も計算できることが知られています。

さっそくこれらの定義式を使って Floating を実装していきたいところですが、この実装は非常に骨が折れそうです。試しに指数関数のGHR微分を計算してみましたが、実装しやすい簡潔な式にはなりそうにありませんでした^[4]（僕の計算が十分でない可能性もあるので簡潔な式で表せたらぜひ教えてください🙏）。上に挙げたGHR微積分の原論文には指数関数の微分として

$$\frac{\partial \exp(q)}{\partial q^\mu}\mu = \sum_{n=0}^\infty\frac{1}{n!}\sum_{m=1}^nq^{n-m}\Re(q^{m-1}\mu)$$

という表示が与えられていました。近似計算で良ければこの式を有限項で打ち切って計算するといった方法も考えられそうです。

実数でも複素数でも指数関数の微分は指数関数になりました。これは指数関数が微分演算子の固有関数になっていることを意味しています。四元数において指数関数の微分（GHR微積分）が元の指数関数から形を変えてしまうということは、指数関数はもはや微分演算子の固有関数になっていないということです。逆にGHR微積分における微分演算子の固有関数はどのようなものになっているのかも気になるところですね👀

活性化関数

気を取り直して最後に活性化関数 Activatable のインスタンスを実装していきましょう。四元数における活性化関数も複素数の時と同様に各成分に実数のシグモイド関数を適用する形で実装します。シグモイド関数には指数関数が使われますが、四元数関数の指数関数ではなく、あくまで実数の指数関数なので実装もなんとかできそうです。

instance Activatable H where
  sigmoid q = fmap sigmoid q

instance Activatable Q where
  sigmoid (Q q@(Quaternion x (V3 y z w)) q' qi' qj' qk') = Q (sigmoid q) p' pi' pj' pk' where
    sigmoid' x = sigmoid x * (1 - sigmoid x)
    dfdq  = (sigmoid' x *| e - sigmoid' y *| i - sigmoid' z *| j - sigmoid' w *| k) / 4
    dfdqi = -i * dfdq * i
    dfdqj = -j * dfdq * j
    dfdqk = -k * dfdq * k
    (p',pi', pj', pk') = ghrChainRule (dfdq, dfdqi, dfdqj, dfdqk) (q', qi', qj', qk')

実装した四元数関数の活性化関数は数式で書くと以下のようになります。

$$\sigma_{\mathbb H}(x+yi+zj+wk) = \sigma_{\mathbb R}(x)+\sigma_{\mathbb R}(y)i+\sigma_{\mathbb R}(z)j+\sigma_{\mathbb R}(w)k$$

この関数の微分は以下のように計算することができます。

$$\begin{matrix} \frac{\partial\sigma_{\mathbb H}}{\partial q^\mu} &=& \frac{1}{4}\left(\frac{\partial\sigma_{\mathbb R}(x)}{\partial x} - \frac{\partial\sigma_{\mathbb R}(y)}{\partial y}i^\mu - \frac{\partial\sigma_{\mathbb R}(z)}{\partial z}j^\mu - \frac{\partial\sigma_{\mathbb R}(w)}{\partial w}k^\mu\right) \\ &=& \left(\frac{\partial\sigma_{\mathbb H}}{\partial q}\right)^\mu \end{matrix}$$

各虚数単位の係数となる導関数は実数値関数のため四元数による回転と可換になることを使っています。$\frac{\partial\sigma_{\mathbb H}}{\partial q}$ は

$$\begin{matrix} \frac{\partial\sigma_{\mathbb H}}{\partial q} &=& \frac{1}{4}\left(\frac{\partial\sigma_{\mathbb R}(x)}{\partial x} - \frac{\partial\sigma_{\mathbb R}(y)}{\partial y}i - \frac{\partial\sigma_{\mathbb R}(z)}{\partial z}j - \frac{\partial\sigma_{\mathbb R}(w)}{\partial w}k\right) \\ &=& \frac{1}{4}\left((1-\sigma_{\mathbb R}(x))\sigma_{\mathbb R}(x) - (1-\sigma_{\mathbb R}(y))\sigma_{\mathbb R}(y)i - (1-\sigma_{\mathbb R}(z))\sigma_{\mathbb R}(z)j - (1-\sigma_{\mathbb R}(w))\sigma_{\mathbb R}(w)k\right) \\ \end{matrix}$$

と計算することができます。

四元数NNの学習

以上の用意した関数を使って四元数NNを学習させてみましょう。

今回は単純に以下のような対応を学習させることにします。

$$x+yi+zj+wk \mapsto -y+xi-wj+zk,\ \ x,y,z,w \in\{-1, 1\}$$

四元数NNは3層用意し中間層には活性化関数を適用し出力層はそのまま出力する形にします。

main :: IO ()
main = do
  -- パラメータを乱数で初期化
  let randomQuaterion = (\(x, y, z, w) -> Quaternion x (V3 y z w)) <$> randomRIO ((-1, -1, -1, -1), (1, 1, 1, 1)) :: IO H
  w1 <- randomQuaterion
  w2 <- randomQuaterion
  w3 <- randomQuaterion
  w4 <- randomQuaterion
  b1 <- randomQuaterion
  b2 <- randomQuaterion
  b3 <- randomQuaterion
  let -- 教師データ
      combinations = do
          x <- [-1, 1]
          y <- [-1, 1]
          z <- [-1, 1]
          w <- [-1, 1]
          pure $ Quaternion x (V3 y z w)
      trainingData = take 10000 $ cycle $ map (\q@(Quaternion x (V3 y z w)) -> (q, Quaternion (-y) (V3 x (-w) z))) combinations
      -- 学習モデル
      model w1 w2 b1 b2 x = linear w2 b2 $ fmap sigmoid $ linear w1 b1 x
      -- 学習係数
      learningRate = 0.1 *| e
      -- パラメータ更新の1ステップ
      updateParam (w1, w2, w3, w4, b1, b2, b3) (x, y) =
          let -- パラメータを整形
              mw1 dw1 dw2 = toM @2 @1 [[Q w1 dw1 0 0 0], [Q w2 dw2 0 0 0]]
              vb1 db1 db2 = toV @2     [Q b1 db1 0 0 0,   Q b2 db2 0 0 0]
              mw2 dw3 dw4 = toM @1 @2 [[Q w3 dw3 0 0 0,   Q w4 dw4 0 0 0]]
              vb2 db3     = toV @1     [Q b3 db3 0 0 0]
              vx = toV @1 [Q x 0 0 0 0]
              ty = Q y 0 0 0 0
              -- 誤差関数の微分係数を計算
              loss (dw1, dw2, dw3, dw4, db1, db2, db3) = (abs (ty - view _1 (model (mw1 dw1 dw2) (mw2 dw3 dw4) (vb1 db1 db2) (vb2 db3) vx))) ^ 2
              dfdw1 = gradQ $ loss (1, 0, 0, 0, 0, 0, 0)
              dfdw2 = gradQ $ loss (0, 1, 0, 0, 0, 0, 0)
              dfdw3 = gradQ $ loss (0, 0, 1, 0, 0, 0, 0)
              dfdw4 = gradQ $ loss (0, 0, 0, 1, 0, 0, 0)
              dfdb1 = gradQ $ loss (0, 0, 0, 0, 1, 0, 0)
              dfdb2 = gradQ $ loss (0, 0, 0, 0, 0, 1, 0)
              dfdb3 = gradQ $ loss (0, 0, 0, 0, 0, 0, 1)
           -- パラメータの更新
           in ( w1 - learningRate * conjugate dfdw1
              , w2 - learningRate * conjugate dfdw2
              , w3 - learningRate * conjugate dfdw3
              , w4 - learningRate * conjugate dfdw4
              , b1 - learningRate * conjugate dfdb1
              , b2 - learningRate * conjugate dfdb2
              , b3 - learningRate * conjugate dfdb3
              )
      -- 学習データによるパラメータ更新
      (w1', w2', w3', w4', b1', b2', b3') = foldl' updateParam (w1, w2, w3, w4, b1, b2, b3) trainingData
      mw1 = toM @2 @1 [[w1'], [w2']]
      mw2 = toM @1 @2 [[w3',   w4']]
      vb1 = toV @2 [b1', b2']
      vb2 = toV @1 [b3']
  -- 学習結果を表示
  forM_ combinations $ \q@(Quaternion x (V3 y z w)) ->
    putStrLn $ concat [
      show x, " ",
      show y, " ",
      show z, " ",
      show w, " : ",
      show . view _1 $ model mw1 mw2 vb1 vb2 (toV @1 [q])
      ]

線形層 linear の実装は抽象化のお陰で複素数の時の実装と同じものを使っています。もう少しライブラリとして作り込めばさらに繰り返し部分を排除できるかもしれませんが、今回はこれで最後の実装なので良しとしましょう。

学習時にはパラメータの共役四元数による微分 $\frac{\partial f}{\partial \overline{q}}$ を用いる必要がありますが、四元数関数 $f$ が実数値のみを取る場合に

$$\overline{\frac{\partial f}{\partial q}} = \frac{\partial f}{\partial \overline{q}}$$

が成り立つことを用いて $\frac{\partial f}{\partial q}$ を使って学習を行っています。

実行して途中の誤差関数の値を確認すると、以下のグラフのように学習が進むに連れて誤差が減少していることが確認できました。

実行した結果を見てみましょう。

$ cabal run
-1.0 -1.0 -1.0 -1.0 : Quaternion 1.0031153244170459 (V3 (-1.006197974510911) 0.9869804545007982 (-1.036277345270371))
-1.0 -1.0 -1.0 1.0 : Quaternion 1.0033954624231927 (V3 (-1.0087185397264407) (-1.0038004527756106) (-1.0304861601600945))
-1.0 -1.0 1.0 -1.0 : Quaternion 1.033901941925051 (V3 (-0.9977533502068777) 1.0028839875276538 0.9945971196354323)
-1.0 -1.0 1.0 1.0 : Quaternion 1.0373527412774388 (V3 (-1.001605582508933) (-1.0065179525339232) 1.0073390494493857)
-1.0 1.0 -1.0 -1.0 : Quaternion (-0.9776054834613533) (V3 (-1.0118800578512797) 0.9777776104880347 (-0.9611666284484391))
-1.0 1.0 -1.0 1.0 : Quaternion (-0.9647189035737046) (V3 (-0.9784101835957897) (-1.014642699091501) (-0.9803623423551261))
-1.0 1.0 1.0 -1.0 : Quaternion (-0.9983142129467947) (V3 (-1.0030911750955684) 1.0179038268959517 0.9293817383623729)
-1.0 1.0 1.0 1.0 : Quaternion (-0.9961571570391524) (V3 (-0.9708249842126295) (-1.005721036536831) 0.9639181649188818)
1.0 -1.0 -1.0 -1.0 : Quaternion 0.9443458800939036 (V3 0.9903424798278815 1.007030552675442 (-0.9932629960716303))
1.0 -1.0 -1.0 1.0 : Quaternion 0.9463624390300247 (V3 1.001199162862839 (-1.0011544864972952) (-1.007352964520582))
1.0 -1.0 1.0 -1.0 : Quaternion 1.0351668151999411 (V3 1.0152543250535957 0.9923390405807191 1.0101324026659984)
1.0 -1.0 1.0 1.0 : Quaternion 1.0429114086715854 (V3 1.0114749454601668 (-0.9891407344009309) 1.0116133081673537)
1.0 1.0 -1.0 -1.0 : Quaternion (-1.0034733782284597) (V3 1.0041762611202794 1.0092104655777823 (-0.9007108601232559))
1.0 1.0 -1.0 1.0 : Quaternion (-0.9794129177244678) (V3 1.003624546949018 (-1.0145015842124137) (-0.9375120134583156))
1.0 1.0 1.0 -1.0 : Quaternion (-1.0178119171861661) (V3 1.0163151148456766 1.0180004206865816 0.9535234667336302)
1.0 1.0 1.0 1.0 : Quaternion (-0.999746475988771) (V3 1.000810906053783 (-0.9906703132702619) 0.9784125614126059)

想定通りに対応を学習できていますね👏

本稿で実装した全てのコードはこちらのGistにアップロードしています。

おわりに

今回はフォワードモード自動微分を実装しましたが実用の観点からはリバースモード自動微分の実装も重要でしょう。あと複素数、四元数とくれば次は八元数の世界でどのような微分法を考えることができるのかも気になるところです。八元数は普段は中々お目にかかれませんが超ひも理論の超対称性の記述に応用例があるらしく（「八元数と超ひも理論」）、物理学と機械学習の交流が盛んな今面白いテーマかもしれません。数の一般化は八元数にとどまらず、クリフォード代数に基づくNN（[2302.06594] Geometric Clifford Algebra Networks）や$C^*$-環上のNN（[2206.09513] $C^*$-algebra Net: A New Approach Generalizing Neural Network Parameters to $C^*$-algebra）なども提案されておりまだまだ面白い発展が続きそうです。
本稿を書いたきっかけは2022年12月26日にGHR微積分を用いた誤差逆伝播法の実装について言及された論文（[2212.13082] Quaternion Backpropergation）がarXivにアップロードされていて、偶然目に止まりGHR微積分というものを初めて知りました。この論文ではNNの誤差逆伝播法の計算式が明示的に導出されていますが、本稿では更にGHR微積分を自動微分の枠組みで実装することで誤差逆伝播法に限らない一般的な最適化問題に適用できる形で実現することができました。GHR微積分の積の微分法と合成関数の微分法は四元数の積の非可換性から通常のものとは異なる形になるので最初は自動微分で実装できるのか不安なところもありましたが、なんとか形になったので満足です。昨年暮れに論文を見つけてから週末にコツコツと調べたり実装したり執筆をしたりしていたら（間が空いた期間もありましたが）まるまる1年経ってしまいました😅 次の記事はもう少し早く書けるといいなぁ。

脚注

1. 実装までは載っていませんがウィルティンガー微分で自動微分を実装するアイデアはこちらの論文でも提案されています
   GUO, Chu; POLETTI, Dario. Scheme for automatic differentiation of complex loss functions with applications in quantum physics. Physical Review E, 2021, 103.1: 013309.

2. しかしフォワードモード自動微分には実装が簡単であり計算コストも低いというメリットがあります。フォワードモード自動微分は1回の計算で方向微分を求めることができるので、この性質を利用して乱択した方向ベクトルとその方向微分を使ってパラメータ更新を行うことにより、誤差逆伝播法の逆方向の計算コストを抑えることができると言った学習アルゴリズムも提案されています
   BAYDIN, Atılım Güneş, et al. Gradients without backpropagation. arXiv preprint arXiv:2202.08587, 2022.

3. このあたりの歴史についてはこちらの講義ノートの 「1.1.8 ベクトルとテンソルの概念に関する簡単な歴史」に詳しく書いてあります

4. ぐぬぬ...