四元數與旋轉--計算、推導與實現

上一篇說到四元數的旋轉表示，即從軸角法表示變化成的單位四元數的形式，即：

\[p=(\cos\frac{\theta}{2}, \boldsymbol{n}\sin\frac{\theta}{2}) \]

但是並沒有說它是怎么來的。區區一個結果怎么能行呢，所以這里我就來細細說一說四元數的運算法則。

為了方便表示，這里規定兩個四元數作為我們接下來演示用的小白鼠：

\[q_1 = (s, \boldsymbol{u}) = a_1+b_1i+c_1j+d_1k \]

\[q_2 = (t, \boldsymbol{v}) = a_2+b_2i+c_2j+d_2k \]

四元數基本運算

加減

加減法很簡單，沒什么好說的，只要把對應位置的元素加減就行：

\[q_1\pm q_2=(a_1 \pm a_2)+(b_1 \pm b_2)i+(c_1 \pm c_2)j+(d_1 \pm d_2)k \]

寫成矢量、標量部分分開的表示形式也可以：

\[q_1 \pm q_2=(s \pm t, \boldsymbol{u} \pm \boldsymbol{v}) \]

乘法

之前給出過有關 \(i, j, k\) 的運算律：

1. \(i^2=j^2=k^2=-1\);
2. \(ij=k, jk=i, ki=j\);
3. \(ij=-ji, jk=-kj, ki=-ik\).

根據這些，我們把兩個四元數當做多項式來計算，就可以得出結果。由於結果太長了，我們換個寫法：假定 \(q=q_1q_2=[w, x, y, z]\), 那么

\[w = a_1a_2-b_1b_2-c_1c_2-d_1d_2 \]

\[x = a_1b_2+a_2b_1+c_1d_2-c_2d_1 \]

\[y = a_1c_2+a_2c_1+d_1b_2-d_2b_1 \]

\[z = a_1d_2+a_2d_1+b_1c_2-b_2c_1 \]

好家伙，這是真的長。有沒有更簡潔的表達方式呢？答案是有的。把四元數寫成矢量、標量部分之后，我們會發現，計算結果的標量部分，也就是上面的 \(w\)，可以寫成：

\[w=st-\boldsymbol{u}\cdot\boldsymbol{v} \]

而矢量部分的計算，也就是上面的 \(x, y, z\)，每一個式子都有四項，這四項可以按照是否含有標量部分因子被分為兩部分：含有標量因子的組合起來，我們看到它是 \(s\boldsymbol{v}+t\boldsymbol{u}\). 而不含有標量因子的部分，看上去像是行列式的形式，這里直接說結論吧，這部分可以寫成 \(\boldsymbol{u\times v}\). 於是，計算結果可以表示成

\[q_1q_2=(st-\boldsymbol{u}\cdot\boldsymbol{v}, s\boldsymbol{v}+t\boldsymbol{u}+\boldsymbol{u}\times\boldsymbol{v}) \]

這時，我們請出老朋友二元復數，我們可以看到很多相似點，四元數的乘積比普通的復數多了一項：\(\boldsymbol{u}\times\boldsymbol{v}\). 這一項決定了四元數的乘法沒有交換律，甚至不具有反交換律。有趣的是，如果矢量部分是同向或反向的，這一項就是 0，而對於一般的復數，虛部只有一個單位 \(i\)，它們只能共線。

模

按照直覺，四元數 \(q\) 的模 \(||q||\)，應該就是把四個成員取平方加起來再開方，也就是常說的二范數。實際上也是如此。這部分就這樣，公式就不寫了。是不是很無聊？這可能是本文最無聊的一個小節了 2333。

共軛

四元數有三個虛部，怎么求共軛呢？答案很簡單，就是把三個虛部都反過來。舉個例子，\(q=(s, \boldsymbol{u})\)，那么 \(q\) 的共軛為

\[q^\star=(s, -\boldsymbol{u}) \]

這個定義符合在二元復數里我們對共軛的一般認知，也是就是，滿足共軛的運算律。比如說：

\[qq^\star=q^\star q=||q||^2 \]

\[q+q^\star=2s, q-q^\star=2\boldsymbol{u} \]

具體證明按照上面的乘法運算律算一算就知道了，不算復雜。

旋轉的四元數表示究竟是怎么來的

旋轉的分解

為了搞清楚這個問題，我們首先要將旋轉操作相關的量用四元數表示出來。其實只有三個相關的向量：旋轉前向量 \(\boldsymbol{u}\), 旋轉后向量 \(\boldsymbol{v}\), 和旋轉軸 \(\boldsymbol{n}\) 和一個角度，即旋轉角 \(\theta\). 旋轉前后的向量可以用純四元數表示，即轉前：\(u=(0, \boldsymbol{u})\)，轉后：\(v=(0, \boldsymbol{v})\). 注意，這里的轉軸的表示向量 \(\boldsymbol{n}\) 是一個單位向量。

要解決這種這種繞固定軸旋轉的問題，一個自然的思路是將被旋轉的向量分解成平行於軸向的 \(\boldsymbol{u}_{\parallel}\) 和垂直於軸向的 \(\boldsymbol{u}_\perp\).

顯然，旋轉對於平行分量是無作用的，只作用於垂直分量。暫且不考慮與轉軸平行的分量，我們將目光放在垂直分量上。把 \(\boldsymbol{u}_\perp\) 旋轉 \(\theta\) 角度，得到的結果就是 \(\boldsymbol{v}_\perp\)，所以后者相對前者的分量，也就是在 \(\boldsymbol{u}_\perp\) 和 \(\boldsymbol{n}\times\boldsymbol{u}_\perp\) 這兩個垂直方向上的分量，我們是知道的，我們有：

\[\boldsymbol{v}_\perp =\boldsymbol{u}_\perp \cos\theta + \boldsymbol{n}\times\boldsymbol{u}_\perp\sin\theta \]

對比一下我們之前說過的四元數乘積形式，不難發現，如果我們用 \(q=(\cos\theta, \boldsymbol{n}\sin\theta)\) 來左乘 \(u_{\perp}=(0, \boldsymbol{u}_\perp)\)，會得到：

\[qu_{\perp}=(0-\boldsymbol{n}\sin\theta\cdot\boldsymbol{u}_\perp, \boldsymbol{u}_\perp\cos\theta+\boldsymbol{n}\times\boldsymbol{u}_\perp\sin\theta)=(0, \boldsymbol{v}_\perp) \]

代數小 trick

我們已經很接近最終結果了，這個結論和最終結論的差異就在於，被乘的是整個 \(q\) 而不是它的分量。要做的就是找到一種運算，可以不影響平行分量，只作用於垂直分量。看看在本文開始提出的最終結果，如果把它取平方的話，會發現

\[p^2=(\cos\theta, \boldsymbol{n}\sin\theta) \]

這就是我們說的作用於垂直分量的四元數 \(q\).

再進一步，我們把旋轉寫成一個很好看的形式：

\[v=u_{\parallel}+qu_{\perp}=pp^{-1}u_{\parallel}+ppu_{\perp} \]

因為 \(p\) 是單位四元數，所以 \(p^{-1}=p^\star\) （不理解的話可以參考二元復數的性質）注意 \(p\) 的矢量部分，是和轉軸平行的，那么把下面的引理一說，大家就明白我要做什么了：

引理：

假設 \(u=(0, \boldsymbol{u})\) 是一個純四元數，而 \(q=(\alpha, \beta\boldsymbol{n})\)，其中 \(\boldsymbol{n}\) 是一個單位向量，\(\alpha, \beta \in R\)， 那么，根據 \(\boldsymbol{n}\) 和 \(\boldsymbol{u}\) 之間的關系，有以下結論：

1. 若 \(\boldsymbol{n} \parallel \boldsymbol{u}\)， 則 \(qu = uq\);
2. 若 \(\boldsymbol{n} \perp \boldsymbol{u}\)，則 \(qu=uq^\star\).

證明同樣省略，直接計算就得到了。

有沒有發現，旋轉后的兩項剛好分別符合引理中的兩個結論。於是：

\[\array{v&=&pp^\star u_\parallel+ppu_\perp\;\;\\ &=&pu_\parallel p^\star+pu_\perp p^\star\\ &=&p(u_\parallel+u_\perp)p^\star\;\;\\ &=&pup^\star\qquad\qquad\;} \]

這就是我們的結論啦。

四元數運算的 Python 實現

四元數這個東西說年輕也不算年輕，那么是不是有相關的庫呢？我在 GitHub 上面翻了翻，還真有，而且已經有 300+ stars 了。

> 鏈接在這 <

具體來說，這個庫是為 numpy 提供了一個四元數的 dtype，而對於簡單的單個四元數的運算也是完全可以勝任的。注意，如果你對運行速度有要求（或者單純有強迫症，不想每次運行都看到 warning），最好安裝一下 numba，這是一個提高 Python 運行速度的工具。

這個庫的安裝很簡單，因為它已經加入 PyPI 了。注意這個包是基於 numpy 的所以要先安裝 numpy.

pip install numpy-quaternion

有了這個工具，我們就可以輕松地應對四元數的運算了：

import numpy as np
import quaternion

q1 = np.quaternion(1, 2, 3, 4)
q2 = np.quaternion(5, 6, 7, 8)

print("q1 = ", q1)
print("\nq2 = ", q2)
print("\nq1 + q2 = ", q1 + q2) # quaternion(6, 8, 10, 12)
print("\nq2 - q1 = ", q2 - q1) # quaternion(4, 4, 4, 4)
print("\nq2 * q1 = ", q1 * q2) # quaternion(-60, 12, 30, 24)
print("\nThe conjugate of q1 is ", q1.conjugate()) # quaternion(1, -2, -3, -4)
print("\nThe square of q1's norm is ", q1.norm()) # 30.0 這里輸出的實際上並不是模，而是模的平方，算是 bug 吧:-/

arr = np.array([[q1, q2]])
print("\nThere is an example of array of quaternions:\n", arr)
print("\nIts transpose times itself is:\n", arr.T * arr)
print("\nIts elementwise exponential in base e is:\n", np.exp(arr))

注意幾個點：

• 現行版本的內置函數 norm 是有問題的，會返回模的平方，如果需要求模記得取平方根；
• numpy 的一些逐元素執行的運算，如上面例子中的 numpy.exp 在這個庫中也是實現了的。
• 這里我沒有講四元數的指數運算，不過其實很簡單 w 感興趣的可以自行了解或者在評論區留言。
• 四元數乘法沒有交換律，這一點在矩陣運算的時候也有體現。

結尾

例行總結好無聊啊，不想寫，那就寫一下本文的參考吧。

本文的求解思路參考了 Krasjet 的四元數與三維旋轉一文的思路。文章講解十分詳實，推薦閱讀。

就到這里啦，有什么問題或發現什么錯誤可以留言或者私信呀。

祝福每一個在努力學習的人 : )