四元数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 + i x + j y + k z

$q=s+ix+jy+kz$

其中 $s,x,y,z$ 都是实数，并满足这样的一些运算规则：

i^{2} = j^{2} = k^{2} = i j k = - 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

$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} + i x_{1} + j y_{1} + k z_{1}) + (s_{2} + i x_{2} + j y_{2} + k z_{2}) = (s_{1} + s_{2}) + i (x_{1} + x_{2}) + j (y_{1} + y_{2}) + k (z_{1} + z_{2})

$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)$

可以简单证明，四元数的加法满足交换律、结合律和分配律，这里不过多展开介绍。

四元数缩放

在系数缩放这一点上，四元数与复数是一致的：

λ q = λ s + i λ x + j λ y + k λ z

$\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中都大量的使用了四元数的计算。而大部分的四元数的教材中写的计算法则，经常把各类乘法混在一起使用，阅读起来非常难受，因此只好自己总结一下四元数的相关运算。并且跟我们所熟悉的复数运算有一定的对比，更加容易去理解四元数的概念。