四元數Quaternion的基本運算

本文主要介紹四元數Quaternion的一些基本運算法則。四元數的概念，更像是複數的一個推廣，在圖形學和工程學中有大量的應用，在蛋白質結構預測軟體AlphaFold和MEGA-Protein中都大量的使用了四元數的計算。而大部分的四元數的教材中寫的計算法則，經常把各類乘法混在一起使用，閱讀起來非常難受，是以隻好自己總結一下四元數的相關運算。并且跟我們所熟悉的複數運算有一定的對比，更加容易去了解四元數的概念。

技術背景

在前面一篇​​文章​​中我們介紹了歐拉角死鎖問題的一些産生背景，還有基于四元數的求解方案。四元數這個概念雖然重要，但是很少會在通識教育課程中涉及到，更多的是一些圖形學或者是工程學當中才會進行講解。本文主要是面向四元數，相比上一篇文章更加詳細的介紹和總結一下四元數的一些運算法則，還有基于四元數的插值法。

基本運算

說到四元數，很多人可能會覺得有點陌生，但是如果說複數，很多人就都有學習過。我們一般用\(z=x+iy\)這樣的形式去定義一個複數（Complex Number），其中\(x\)是實部，而\(y\)是虛部，\(i\)是虛數機關，并且有\(i^2=-1\)這樣的特性。并且對于一個虛數而言，如果取自然指數（運算規則為：\(e^{i\theta}=cos\theta+i\ sin\theta\)），還能夠得到一個很美的數學公式：\(e^{i\pi}=-1\)，這就是非常著名的歐拉公式。

而四元數Quaternion這個概念的提出，更像是對複數的一個擴充，我們通常把四元數寫成這樣的形式：

\[q=s+ix+jy+kz\]

其中\(s,x,y,z\)都是實數，并滿足這樣的一些運算規則：

\[i^2=j^2=k^2=ijk=-1\\i\times j=k,j\times k=i,k\times i=j\\

j\times i=-k,k\times j=-i,i\times k=-j

\]

以上都是四元數的一些基本定義，接下來我們逐一看一下四元數的一些基本運算。

四元數加法

兩個四元數的加法就是将“實部虛部”對應位置做元素求和：

\[q_1+q_2=(s_1+ix_1+jy_1+kz_1)+(s_2+ix_2+jy_2+kz_2)=(s_1+s_2)+i(x_1+x_2)+j(y_1+y_2)+k(z_1+z_2)\]

可以簡單證明，四元數的加法滿足交換律、結合律和配置設定律，這裡不過多展開介紹。

四元數縮放

在系數縮放這一點上，四元數與複數是一緻的：

\[\lambda q=\lambda s+i\lambda x+j\lambda y+k\lambda z\]

逐一對四元數中的各項元素進行縮放即可。

四元數乘法

四元數的乘法是所有元素之前都要運算一遍：

\[\begin{align*}q_1q_2&=(s_1+ix_1+jy_1+kz_1)*(s_2+ix_2+jy_2+kz_2)\\

&=(s_1s_2-x_1x_2-y_1y_2-z_1z_2)\\

&\ +i(s_1x_2+s_2x_1+y_1z_2-y_2z_1)\\

&\ +j(s_1y_2+s_2y_1+x_2z_1-x_1z_2)\\

&\ +k(s_1z_2+s_2z_1+x_1y_2-x_2y_1)

\end{align*}

\]

這個運算過程是這樣的，我們令\(q=q_1q_2=s+ix+jy+kz\)，那麼這個乘法運算最終組成\(s\)這個元素的，分别是\(s_1*s_2, (ix_1)*(ix_2), (iy_1)*(iy_2), (iz_1)*(iz_2)\)這些項，而\(i^2=j^2=k^2=-1\)，是以最終得到\(s=(s_1s_2-x_1x_2-y_1y_2-z_1z_2)\)。但是直到這裡為止，我們所涉及到的元素乘法隻是在“實部”和相同的“虛數機關”之間的運算，如果一旦涉及到不同的“虛數機關”之間的乘法運算，那就要自動轉化成向量叉乘，比如\(ij=ji=k\)。是以，我們要計算\(q\)中的\(z\)項的時候，隻需要計算\(s_1z_2, s_2z_1, x_1y_2, x_2y_1\)這些項即可，同時注意符号的變換，那麼得到的最終的結果就是如上所示。

需要注意的是，​

​四元數與複數的最大的一點不同，複數乘法是有交換律的，而四元數沒有​

​。舉個例子說，我們可以計算一下\(q_1,q_2\)的對易：

\[\begin{align*}[q_1, q_2]&=q_1q_2-q_2q_1\\

&=(s_1s_2-x_1x_2-y_1y_2-z_1z_2)+i(s_1x_2+s_2x_1+y_1z_2-y_2z_1)+j(s_1y_2+s_2y_1+x_2z_1-x_1z_2)+k(s_1z_2+s_2z_1+x_1y_2-x_2y_1)-(s_2s_1-x_2x_1-y_2y_1-z_2z_1)-i(s_2x_1+s_1x_2+y_2z_1-y_1z_2)-j(s_2y_1+s_1y_2+x_1z_2-x_2z_1)-k(s_2z_1+s_1z_2+x_2y_1-x_1y_2)\\

&=2i(y_1z_2-y_2z_1)+2j(x_2z_1-x_1z_2)+2k(x_1y_2-x_2y_1)

\end{align*}

\]

那麼也就是說，這兩個四元數\(q_1,q_2\)之間是非對易的，也就是不可交換的。但是，​

​四元數的運算是滿足結合律和配置設定率的​

​。

由于上面的這種四元數乘法展開，寫起來過于繁雜，我們考慮對其進行一定的簡化。如果我們假定2個純虛數（\(s=0\)）：

\[a=ix_1+jy_1+kz_1\\b=ix_2+jy_2+kz_2

\]

其實類似于這種形式的四元數，實際上就是三維空間中的向量，那麼這兩者的點積和叉積有：

\[a\cdot b=x_1x_2+y_1y_2+z_1z_2\\a\times b=(y_1z_2-y_2z_1)i+(z_1x_2-z_2x_1)j+(x_1y_2-x_2y_1)k

\]

需要注意的是，這裡的叉積是向量叉積，跟四元數中的“虛數機關”相比，最大的一點不同就是：在向量叉積中，\(i\times i=0\)，但是在四元數的乘法中，\(i\times i=-1\)（非常重要）。

那麼在有了以上的兩個公式之後，我們就可以對四元數的乘法表示做一個簡化：

\[q_1q_2=s_1s_2-a\cdot b+s_ab+s_ba+a\times b\]

實四元數和純四元數

對于一個實四元數而言，就是取\(x=y=z=0\)：

\[q_r=s\]

對于一個純四元數而言，就是取\(s=0\)：

\[q_i=ix+jy+kz\]

四元數共轭

對四元數的所有“虛部”取負數，即是四元數的共轭：

\[q^*=s-ix-jy-kz\]

機關四元數

四元數的模的定義跟複數是一緻的：

\[|q|=\sqrt{s^2+x^2+y^2+z^2}=\sqrt{qq*}\]

而機關四元數的定義即是模為1的四元數：

\[s^2+x^2+y^2+z^2=1\]

如果給定的一個四元數不是機關四元數，那麼我們可以對其進行規範化：

\[q'=\frac{q}{\sqrt{s^2+x^2+y^2+z^2}}\]

四元數的逆

對于一個機關四元數而言，因為有\(qq^*=1\)，是以機關四元數的逆就是其共轭四元數。如果是對于更加一般的場景，我們可以這樣考慮：

\[q(q^{-1}*|q|^2)=|q|^2\\qq^*=|q|^2\\

q^{-1}=\frac{q^*}{|q|^2}

\]

比較特殊地，對于機關四元數\(q^{-1}=q^*\)。

四元數的二進制表示

類似于複數的二進制形式，通常一個四元數也可以被表示成如下的二進制形式：

\[q=s+v\hat{q}=[s,v\hat{q}]\]

其中\(v=[x,y,z],\hat{q}=[i,j,k]\)。關于此處的乘法描述，其實有一定的不嚴謹性，因為它既不是點積，也不是叉積，也不是外積，而是普通的元素乘。這種元素乘的概念在計算機領域是很常用的，但是在數學上其實并不是很常用。在這種二進制描述下，四元數的乘法形式會略有調整：

\[q_1q_2=[s_1s_2-(v_1\hat{q}_1)\cdot (v_2\hat{q}_2), s_1v_2\hat{q}_2+s_2v_1\hat{q}_1+(v_1\hat{q}_1)\times(v_2\hat{q}_2)]\]

四元數點積

上面的章節中提到過四元數的普通乘法，但其實四元數也像普通的向量一樣可以進行點積運算：

\[q_1\cdot q_2=s_1s_2+v_1\cdot v_2\]

這也是受益于四元數的二進制表示，使得我們在書寫結果的時候可以更加的簡練。

四元數的指數

我們先來回顧一下複數\(z=x+iy\)的指數計算，根據泰勒展開公式\(f(x)=\sum_n\frac{f^{(n)}(x_0)}{n!}(x-x_0)^n\)（比較特殊地，\(e^x=\sum_{k=0}^\infty\frac{x^k}{k!}\)）對\(e^z\)在\(y=0\)處的展開有：

\[\begin{align*}e^{z}&=e^{x+iy}=e^xe^{iy}\\

&=e^x\left(

1+iy-\frac{1}{2!}y^2-\frac{i}{3!}y^3+\frac{1}{4!}y^4+\frac{i}{5!}y^5-\frac{1}{6!}y^6-\frac{i}{7!}y^7+...

\right)

\end{align*}

\]

對比一下常用的三角函數的泰勒展開式（相關證明見參考連結2）：

\[sin\ x=x-\frac{1}{3!}x^3+\frac{1}{5!}x^5-\frac{1}{7!}x^7+...\\cos\ x=1-\frac{1}{2!}x^2+\frac{1}{4!}x^4-\frac{1}{6!}x^6+...

\]

代入可得：

\[e^{z}=e^xe^{iy}=e^x\left(

cos\ y+i\ sin\ y

\right)

\]

也即，對于一個純虛數\(i\theta\)而言，其指數為：\(e^{i\theta}=cos\theta+i\ sin\theta\)。那麼類似的，對于一個二進制表示的四元數\(q=s+v\hat{q}\)有：

\[\begin{align*}e^q&=e^{s+v\hat{q}}\\

&=e^s\left(

1+v\hat{q}-\frac{1}{2!}(v\hat{q})^2-\frac{1}{3!}(v\hat{q})^3+\frac{1}{4!}(v\hat{q})^4+\frac{1}{5!}(v\hat{q})^5-\frac{1}{6!}(v\hat{q})^6-\frac{1}{7!}(v\hat{q})^7+...

\right)

\end{align*}

\]

這裡有一點不同的是，我們計算四元數的幂次的時候需要謹慎，可以先手動計算一下：

\[\begin{align*}(v\hat{q})^2&=0-(v\hat{q})\cdot (v\hat{q})+0+0+(v\hat{q})\times(v\hat{q})=-|v|^2\\

(v\hat{q})^3&=-|v|^2(v\hat{q})\\

(v\hat{q})^4&=|v|^4\\

(v\hat{q})^5&=|v|^4(v\hat{q})\\

(v\hat{q})^6&=-|v|^6\\

(v\hat{q})^7&=-|v|^6(v\hat{q})\\

&...

\end{align*}

\]

代入四元數的指數部分進行計算可得：

\[\begin{align*}e^q&=e^{s+v\hat{q}}\\

&=e^s\left(

1+v\hat{q}-\frac{1}{2!}|v|^2-\frac{1}{3!}|v|^2(v\hat{q})+\frac{1}{4!}|v|^4+\frac{1}{5!}|v|^4(v\hat{q})-\frac{1}{6!}|v|^6-\frac{1}{7!}|v|^6(v\hat{q})+...

\right)\\

&=e^s\left(

cos|v|+\frac{v\hat{q}}{|v|}sin|v|

\right)

\end{align*}

\]

這就是四元數的指數運算。

四元數的指數表示

區分于上一個章節中的四元數的指數運算，這個章節我們是要用一個指數形式去表示任意給定的一個四元數。因為在上一個章節中我們發現，一個四元數的指數形式是另外一個四元數，是以，理論上說我們可以用一個指數形式來表示任意的一個四元數。我們首先還是參考一下複數的指數表示：

\[z=x+iy=\sqrt{x^2+y^2}\left(\frac{x}{\sqrt{x^2+y^2}}+i\frac{y}{\sqrt{x^2+y^2}}\right)=\sqrt{x^2+y^2}e^{i\frac{y}{|y|}\ arccos\left(\frac{x}{\sqrt{x^2+y^2}}\right)}\]

類似地，一個四元數可以表示為：

\[q=s+v\hat{q}=\sqrt{s^2+|v|^2}\left(\frac{s}{\sqrt{s^2+|v|^2}}+\frac{v\hat{q}}{|v|}\frac{|v|}{\sqrt{s^2+|v|^2}}

\right)=\sqrt{s^2+|v|^2}e^{\hat{q}\frac{v}{|v|}\ arccos\left(\frac{s}{\sqrt{s^2+|v|^2}}\right)}

\]

比較有意思的是，如果我們取\(q=s+ix+jy+kz\)中的\(y=0,z=0\)時，我們發現\(v=\pm|v|,\hat{q}=i\)，這樣一來，四元數的指數表示形式就和複數的指數表示形式完全對應上了。

四元數的對數

在上一個章節中，如果我們把一個四元數表示成一個指數的形式，就會很大程度上友善我們去計算一個四元數\(q=s+v\hat{q}\)的對數：

\[\begin{align*}log(q)&=log\left(\sqrt{s^2+|v|^2}e^{\hat{q}\frac{v}{|v|}\ arccos\left(\frac{s}{\sqrt{s^2+|v|^2}}\right)}\right)\\

&=log\left(\sqrt{s^2+|v|^2}\right)+\hat{q}\frac{v}{|v|}\ arccos\left(\frac{s}{\sqrt{s^2+|v|^2}}\right)

\end{align*}

\]

這樣就得到了四元數的對數的二進制表示形式。

四元數的幂次

了解了四元數的指數和對數的計算子產品之後，我們可以計算一個四元數的幂次。正是由于四元數的指數表示形式，使得我們可以将四元數的幂次簡單的轉化成乘法的表示形式：

\[q^t=\left[\sqrt{s^2+|v|^2}e^{\hat{q}\frac{v}{|v|}\ arccos\left(\frac{s}{\sqrt{s^2+|v|^2}}\right)}\right]^t=\left(s^2+|v|^2\right)^{\frac{t}{2}}e^{\hat{q}\frac{vt}{|v|}\ arccos\left(\frac{s}{\sqrt{s^2+|v|^2}}\right)}

\]

那麼這就得到了四元數的幂次表達形式。

歐拉角旋轉四元數

在上一篇​​文章​​中我們提到過，每一個四元數其實都可以對應于三維空間的一個向量旋轉，一個四元數\(q\)作用在一個空間向量\(v\)上就會旋轉得到一個新的空間向量：

\[v'=qvq^*\]

而如果給定了三維空間中的旋轉歐拉角，在四元數中就可以表示為相應的旋轉四元數。比如繞\(X\)軸旋轉\(\beta\)角所對應的四元數為：\(q=cos\frac{\beta}{2}+i\ sin\frac{\beta}{2}\)，繞\(Y\)軸旋轉\(\alpha\)角所對應的四元數為：\(q=cos\frac{\alpha}{2}-j\ sin\frac{\alpha}{2}\)，繞\(Z\)軸旋轉\(\gamma\)角所對應的四元數為：\(q=cos\frac{\gamma}{2}+k\ sin\frac{\gamma}{2}\)。而通常使用的\(ZXY\)順規可表示為：

\[q=\left(cos\frac{\alpha}{2}-j\ sin\frac{\alpha}{2}\right)\left(cos\frac{\beta}{2}+i\ sin\frac{\beta}{2}\right)\left(cos\frac{\gamma}{2}+k\ sin\frac{\gamma}{2}\right)\]

關于更多的旋轉四元數的内容，可以閱讀一下參考連結3中的内容。

向量變換四元數

這個問題的定義是比較清晰的，如果給定空間中的兩個不同的向量，能否直接獲得這兩個向量之間變換的四元數呢？如果用公式來表示就是：已知\(\textbf{v}_1,\textbf{v}_2\)兩個空間向量，求\(q\)使得\(\textbf{v}_2=q\textbf{v}_1q^*\)。關于這個問題的求解，在參考連結3中也是有介紹的，這裡再簡單提一下計算方法：

\[\textbf{u}=\textbf{v}_1\times\textbf{v}_2\\cos\theta=\frac{\textbf{v}_1\cdot\textbf{v}_2}{|\textbf{v}_1||\textbf{v}_2|}\\

q=cos\frac{\theta}{2}+i\ sin\frac{\theta}{2}\textbf{u}\cdot i+j\ sin\frac{\theta}{2}\textbf{u}\cdot j+k\ sin\frac{\theta}{2}\textbf{u}\cdot k

\]

這個算法的本質，其實就是先用向量叉乘找到旋轉軸，然後計算兩個向量之間的夾角，最後再使用四元數的繞旋轉軸旋轉指定角度的公式計算，就可以得到對應的空間向量變換的四元數。

旋轉四元數和變換四元數的不對等場景

這裡為了區分歐拉角旋轉的四元數和繞軸旋轉的四元數，我們将其分别稱呼為旋轉四元數與變換四元數。其實要證明二者是不對等的也很容易，隻要舉出一個例子來證明二者作用的結果不相等即可。其實在參考連結3中已經給出了類似的案例，這裡我們再舉一個簡單例子來進行探讨。首先我們給定一個初始空間向量：

\[\textbf{v}_1=i\]

我們使其先繞\(Y\)軸旋轉90度，再繞\(Z\)軸旋轉90度，那麼利用四元數進行計算，得到的結果是：

\[q_{rot}=(cos\frac{\gamma}{2}+k\ sin\frac{\gamma}{2})(cos\frac{\beta}{2}-j\ sin\frac{\beta}{2})\\\textbf{v}_2=q\textbf{v}_1q^*=k

\]

也就是說，在這個歐拉角對應的旋轉變換下，我們将一個向量從\(X\)軸正方向變換到了\(Z\)軸的正方向上。那麼如果我們考慮使用向量繞軸旋轉的公式來計算的話：

\[\textbf{u}=\textbf{v}_1\times\textbf{v}_2=i\times k=-j\\cos\theta=\frac{\textbf{v}_1\cdot\textbf{v}_2}{|\textbf{v}_1||\textbf{v}_2|}=0\\

q_{trans}=cos\frac{\theta}{2}+i\ sin\frac{\theta}{2}\textbf{u}\cdot i+j\ sin\frac{\theta}{2}\textbf{u}\cdot j+k\ sin\frac{\theta}{2}\textbf{u}\cdot k=

\frac{\sqrt{2}}{2}-j\frac{\sqrt{2}}{2}

\]

顯然，\(q_{rot}\neq q_{trans}\)，隻是在這個案例下有：

\[q_{rot}\textbf{v}_1q_{rot}^*=q_{trans}\textbf{v}_1q_{trans}^*=\textbf{v}_2\]

我們隻需要換一個初始空間向量，比如我們将初始向量設定為\(Z\)軸的負方向，即：\(\textbf{v}_1=-k\)，那麼再看看兩個變換所對應的結果：

\[q_{rot}\textbf{v}_1q_{rot}^*=\frac{1}{4}(1+i-j+k)(-k)(1-i+j-k)=j\\q_{trans}\textbf{v}_1q_{trans}^*=\frac{1}{2}(1-j)(-k)(1+j)=i

\]

顯然這兩個變換之後的結果是不同的。其實道理很簡單，\(i\)在\(q_{rot}\)的作用下，先繞\(Y\)軸轉90度，剛好就到了\(Z\)軸正方向上，此時再繞\(Z\)軸旋轉是不會動的。而\(q_{trans}\)，就隻是繞\(Y\)軸旋轉90度的一個操作，因為從客觀的變換來說，\(q_{rot}\)作用在\(i\)上就等價于隻繞\(Y\)軸旋轉了90度。這就是旋轉的特殊性，在很多情況下，歐拉角的旋轉跟繞軸旋轉是不能畫上等号的。

總結概要