四元数Quaternion的基本运算

技术背景

在前面一篇文章中我们介绍了欧拉角死锁问题的一些产生背景，还有基于四元数的求解方案。四元数这个概念虽然重要，但是很少会在通识教育课程中涉及到，更多的是一些图形学或者是工程学当中才会进行讲解。本文主要是面向四元数，相比上一篇文章更加详细的介绍和总结一下四元数的一些运算法则，还有基于四元数的插值法。

基本运算

说到四元数，很多人可能会觉得有点陌生，但是如果说复数，很多人就都有学习过。我们一般用\(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度。这就是旋转的特殊性，在很多情况下，欧拉角的旋转跟绕轴旋转是不能画上等号的。

总结概要

本文主要介绍四元数Quaternion的一些基本运算法则。四元数的概念，更像是复数的一个推广，在图形学和工程学中有大量的应用，在蛋白质结构预测软件AlphaFold和MEGA-Protein中都大量的使用了四元数的计算。而大部分的四元数的教材中写的计算法则，经常把各类乘法混在一起使用，阅读起来非常难受，因此只好自己总结一下四元数的相关运算。并且跟我们所熟悉的复数运算有一定的对比，更加容易去理解四元数的概念。