SLAM之本质矩阵分解得相对变换的R和t

1 前言

对极几何当中的基础矩阵\(F\)与本质矩阵\(E\)是SLAM当中2D-2D恢复出相对变换的重要变量,其相关推导可参看文献[1]。

在相关矩阵计算时,一般会在两幅图像中,根据特征找到对应匹配对后估计出基础矩阵\(F\)或本质矩阵\(E\),如果是直接估计出基础矩阵\(F\),且知道两幅图中的相机参数分别为\(K_1\)、\(K_2\),则可以直接得到本质矩阵\(E\)

(1-1)$$E=K_2^TFK_1$$

对于本质矩阵与相对变换的\(R\)和\(t\)有如下关系(注意这里指的相对变换是从第一幅图变换到第二幅图的位姿相对变换)

(1-2)$$E=t^\wedge R$$

其中本质矩阵\(E\)还具有如下的性质

(1-3)$$rank(E)=2$$

(1-4)$$2EE^TE=tr(EE^T)E$$

 

2 常规分解[2]

平常从\(E\)得到两幅图像相对变换的\(R\)和\(t\),主要通过奇异值分解得到,即对\(E\)进行SVD分解可得

(2-1)$$E=USV^T$$

这里先放置一个常量\(W\)

(2-2)$$W=\begin{bmatrix}
 0&  -1&0 \\ 
 1& 0 &0 \\ 
 0&0  & 1
\end{bmatrix}$$

那么最终可以得到如下的\(R\)和\(t\)

(2-3)$$\begin{align*}
t_1^\wedge &= UWSU^T \\ 
t_2 &= -t1\\ 
R_1 &=UWV^T \\ 
R_2 &= UW^TV^T
\end{align*}$$

上面总共存在4组可能的解\((t_1,R_1)、(t_1,R_2)、 (t_2,R_1)、(t_2,R_2)\),在这4组解中,仅有一组满足所有投影点在两个相机中都为正深度,因此代入匹配点计算出坐标,即可选出正确的一组。

这里需要说明一点,当计算的本质矩阵\(E\)不是严格满足其性质时,其中\(E\)有一个性质是其特征值前两个为相等的非0值,最后一个为0,即公式(2-1)当中的\(S\)应该为如下形式

(2-4)$$S=\begin{bmatrix}
\sigma  & 0 & 0\\ 
 0&  \sigma &0 \\ 
0 & 0 & 0
\end{bmatrix}$$

这时,可以假设一个严格满足本质矩阵的\(E'\)去逼近给的\(E\)矩阵,这里假设\(E'\)为如下形式

(2-5)$$E'=US'V^T$$

那么可以通过最小化\( \left \| E – E' \right \|_F\)来得到最佳的本质矩阵,即最小化如下方程

(2-6)$$\begin{align*}
min  \left \| E – E' \right \|_F &= min  \left \| U^T(E – E')V \right \|_F \\
 &=min  \left \| U^TEV – U^TE'V \right \|_F \\
 &=min  \left \| S – S' \right \|_F \\
 &=min [(\sigma_1 – \sigma)^2 + (\sigma_2 – \sigma)^2 + \sigma_3^2]
\end{align*}
\Rightarrow
\sigma = (\sigma_1+\sigma_2)/2$$

即当\(S'\)满足如下形式时,有最佳逼近

(2-6)$$S'=\begin{bmatrix}
 (\sigma _1 + \sigma _2)/2& 0 & 0\\ 
0 &  (\sigma _1 + \sigma _2)/2 & 0\\ 
0 & 0 & 0
\end{bmatrix}$$

其中\(E\)的特征值由大到小排列为\(\sigma _1\)、\(\sigma _2\) 、\(\sigma _3\)

因此最终使用(2-7)算式当中的\(S'\)替换(2-3)当中的\(S\)即可。

 

3 快速分解[3]

上面是通过奇异值分解得到\(R\)和\(t\),这里直接求取\(R\)和\(t\)。

先导出一个等式

(3-1)$$\begin{align*}
E^Tt&=(t^\wedge R)^Tt\\
 &= R^T(t^\wedge)^Tt\\ 
 &= R^T(-t^\wedge)t\\ 
 &= -R^T(t^\wedge t)\\ 
 &= 0
\end{align*}$$

从上面公式(3-1)中,可以发现,通过\(E\)可以直接解一个方程\(E^Tt=0\)得到一组非0特殊解\(t\)(真实\(t\)为这个解的\(k\)倍,这个\(k\)目前还不知道,后面可以通过计算反算出真实的\(t\))

从(2-1)可得到\(E\)的模

(3-2)$$\left \| E \right \|_F=\left \| S \right \|_F=\sqrt{2} \sigma$$

而从(2-3)可以得到

(3-3)$$\sqrt{2} \left \| t \right \|_2= \left \| t^\wedge \right \|_F=\left \| WS \right \|_F=\sqrt{2} \sigma$$

结合(3-2)与(3-3)有

(3-4)$$\left \| E \right \|_F=\sqrt{2}  \left \| t \right \|_2$$

从(3-4)可以发现,真实的\(E\)的模为\(t\)的\(\sqrt{2}\)倍,在这个情况下,同时缩放不影响\(R\)的旋转性质,因此不妨取\(E\)的模为\(\sqrt{2}\),\(t\)的模为1,即这时将\(E\)与\(t\)进行缩放,得到

(3-5)$$E'=\sqrt{2}E/ \left \| E \right \|_F$$

(3-6)$$t'=t/ \left \| t \right \|_2$$

则对于公式(1-2)展开有

(3-7)$$\begin{bmatrix}
e_1^T\\ 
e_2^T\\ 
e_3^T
\end{bmatrix}=E'=t'^\wedge R=\begin{bmatrix}
t_1^T\\ 
t_2^T\\ 
t_3^T
\end{bmatrix}R \Leftrightarrow 
\left\{\begin{matrix}
Re_1=t_1\\ 
Re_2=t_2\\ 
Re_3=t_3
\end{matrix}\right. $$

上面公式中的\(e_i\)、\(t_i\)分别表示\(E'\)、\(t'^\wedge \)的第\(i\)行数据。

在(3-7)中需要注意,由于\(E'\)矩阵不可逆(其本身性质决定,它的秩为2)因此并不能直接求得\(R\),但是因为\(E'\)秩为2,因此在(3-7)中的\(e_i\)中必有两个线性无关向量,这里假设选取的两个线性无关的向量为\(e_1\)、\(e_2\),其中\(e_3\)必定可以使用\(e_1\)、\(e_2\)的线性表示。

对于公式(3-7)中的变量进行如下操作

(3-8)$$\begin{align*}
t_1 \times t_2 &= (Re_1)\times (Re_2) \\ 
 &= (Re_1) ^\wedge (Re_2)\\ 
 &= Re_1^\wedge R^T(Re_2)\\
 &=Re_1^\wedge R^TRe_2\\
&=Re_1^\wedge e_2\\
&=R(e_1 \times e_2)
\end{align*}$$

根据叉乘的定理可知\(e_1 \times e_2\)必定与\(e_1\)、\(e_2\)线性无关,因此可以得到一个如下的可逆矩阵\(A\)

(3-9)$$A=\begin{bmatrix}
e_1 & e_2 & e_1 \times e_2
\end{bmatrix}$$

以及一个矩阵\(B\)

(3-10)$$B=\begin{bmatrix}
t_1 & t_2 & t_1 \times t_2
\end{bmatrix}$$

则由(3-7)到(3-10)可以得到如下方程

(3-8)$$ RA=B$$

最终可以得到旋转矩阵

(3-9)$$ R = B A^{-1}$$

由于求得的\(t\)有一个缩放,如果想恢复成原本\(E\)的大小,则可以根据公式(1-2)与(3-9)得到的\(R\)得到\(t\)

(3-10)$$ t^\wedge = ER^T$$

需要说明,由于当\(E=-E\)时也可以通过上面方式得到一个\(R\),因此最终的求解方式是:

①当\(E=E\)时,求得一组\(R_1 , t_1\)

②当\(E=-E\)时,求得一组\(R_2 , t_2\)

这时\(t_1=t_2\)或者\(t_1=-t_2\),这里直接记\(t1=\left | t_1 \right |\),\(t2=-\left | t_1 \right |\)

则最终可能产生4组可能的解\((t_1,R_1)、(t_1,R_2)、 (t_2,R_1)、(t_2,R_2)\),和第2点一样,可以使用正深度法选择出正确的那组解。

可以发现,这个方法分解\(R\)、\(t\)只有解方程与求逆,并且操作矩阵仅为3阶,都有比较简单的常规表达式,因此比起奇异值分解得到变换的\(R\)和\(t\),具有更高效的特点。

 

4 三角测量

当通过以上方法的得到两张图对应的变换\(R\)和\(t\)后,需要恢复出对应点的在空间中的位置时(上面从多组解中筛选出正确的那最解就需要用到),在单目SLAM中,就需要通过三角测量来估计其深度。

配对的两个图像上的特征点,假设归一化后的坐标为\(x_1\)与\(x_2\),这时根据得到的变换,则它们满足

(4-1)$$s_2x_2=s_1Rx_1+t$$

公式中的的\(s_1\)、\(s_2\)分别为两个特征点对应的深度。由于求解变量只有2个,则直接解一个线性方程组即可求得最佳的\(s_1\)、\(s_2\)

 

参考:

[1] 杨述强. 基于单目视觉的场景三维重建与飞行器位姿求解关键技术研究[D]. 国防科学技术大学, 2014.

[2] 高翔 , 张涛. 视觉SLAM十四讲从理论到实践[M]. 北京:电子工业出版社, 2017.3

[3]李国栋, 田国会, 王洪君,等. 弱标定立体图像对的欧氏极线校正框架[J]. 光学精密工程, 2014, 22(7):1955-1961.

SLAM之李群李代数以及其它数学工具

背景 
李群/李代数,因为其特殊的性质,在无人驾驶领域,应用得比较广泛。当在构建数学模型中应用到对应的李群/李代数时,这时当知道一些李群/李代数的性质,以及对应的映射关系时,对一些数学推导可以起到事半功倍的效果。这篇文章主要是对文献[1][2]中用提到的一些公式,进行记录。

定义
1 SO(3)
一个李群SO(3)中的元素\(R\),总能在李代数中找到包含3个元素的向量\(\phi\)与其对应,这里\(R\)是一个\(3 \times 3\)的旋转矩阵.
对于本文的旋转矩阵\(R\),有如下性质 
(1.1)$$
|R|=1\\
R^T = R^{-1}
$$

约定如下的反对称符号 
(1.2)$$
\phi ^\wedge =\begin{bmatrix} 
0 & -\phi_2  & \phi_1\\  
\phi_2 &  0 & -\phi_0\\
-\phi_1 & \phi_0 & 0
\end{bmatrix}_{3 \times 3}
\\
\phi =\begin{bmatrix} 
0 & -\phi_2  & \phi_1\\  
\phi_2 &  0 & -\phi_0\\
-\phi_1 & \phi_0 & 0
\end{bmatrix}^\vee = \begin{bmatrix} 
\phi_0 \\
\phi _1\\
\phi _2
\end{bmatrix}
$$

对于\(R\)与\(\phi\)通过指数映射有如下关系 
(1.3)$$\begin{aligned}
R &= \exp(\phi^\wedge)
\\
\phi^\wedge &= \log(R)
\\
\phi &= \log(R)^\vee\\
(\phi^\wedge)^T &=-\phi^\wedge=(-\phi)^\wedge
\end{aligned}$$

因为一个旋转向量\(\mathbf{a}\)与旋转角\(\theta\)可以表示一个旋转,因此对于如下表达式 
(1.4)$$
\phi = \theta \mathbf{a}
\\
||\mathbf{a}||_2 = 1 \\
\theta = |\phi|
$$
则可以建立联系 
(1.5)$$
R = \cos \theta I + (1-\cos\theta)\mathbf{a}\mathbf{a}^T+ \sin\theta \mathbf{a}^\wedge
$$
其中\(R\)的左右雅克比矩阵\(J_l,J_r\)可表示为 
(1.6)$$\begin{cases}
J_l &=I+\dfrac{1-\cos\theta}{\theta^2}\phi^\wedge+\dfrac{\theta-\sin\theta}{\theta^3}(\phi^\wedge)^2
\\
\\
J_r &=I-\dfrac{1-\cos\theta}{\theta^2}\phi^\wedge+\dfrac{\theta-\sin\theta}{\theta^3}(\phi^\wedge)^2
\\
\\
J_l^{-1}&=I-\dfrac{1}{2}\phi^\wedge+(\dfrac{1}{\theta^2}-\dfrac{1+\cos\theta}{2\theta\sin\theta})(\phi^\wedge)^2
\\
\\
J_r^{-1}&=I+\dfrac{1}{2}\phi^\wedge+(\dfrac{1}{\theta^2}-\dfrac{1+\cos\theta}{2\theta\sin\theta})(\phi^\wedge)^2
\end{cases}
$$

当轴角趋近0时,可简化为如下公式

(1.7)$$\begin{cases}
J_l|_{\theta \rightarrow 0} &\approx I + \dfrac{1}{2}\phi^\wedge+\dfrac{1}{6}(\phi^\wedge)^2
\\
\\
J_r|_{\theta \rightarrow 0} &\approx I – \dfrac{1}{2}\phi^\wedge+\dfrac{1}{6}(\phi^\wedge)^2
\\
\\
J_l^{-1}|_{\theta \rightarrow 0} &\approx I-\dfrac{1}{2}\phi^\wedge+\dfrac{1}{12}(\phi^\wedge)^2
\\
\\
J_r^{-1}|_{\theta \rightarrow 0} &\approx I+\dfrac{1}{2}\phi^\wedge+\dfrac{1}{12}(\phi^\wedge)^2
\end{cases}$$

对于(1.7)忽略二阶小量,可得到如下结果

(1.8)$$\begin{cases}
J_l|_{\theta \rightarrow 0} &\approx I + \dfrac{1}{2}\phi^\wedge
\\
\\
J_r|_{\theta \rightarrow 0} &\approx I – \dfrac{1}{2}\phi^\wedge
\\
\\
J_l^{-1}|_{\theta \rightarrow 0} &\approx I-\dfrac{1}{2}\phi^\wedge
\\
\\
J_r^{-1}|_{\theta \rightarrow 0} &\approx I+\dfrac{1}{2}\phi^\wedge
\end{cases}$$

根据(1.6)定义,\(R\)的左雅克比矩阵与\(R^T\)的右雅克比矩阵相同,\(R\)的右雅克比矩阵与\(R^T\)的左雅克比矩阵相同。

对于一个比较小的量\(\delta \phi\),对旋转矩阵进行左右扰动,有如下结论 
(1.9)$$
\log(\exp(\phi^\wedge)\exp(\delta \phi^\wedge))^\vee \approx J_r^{-1} \delta \phi +\phi
\\
\log(\exp(\delta \phi^\wedge)\exp(\phi^\wedge))^\vee \approx J_l^{-1} \delta \phi +\phi
$$

注意,仔细观察(1.6)与(1.9)式,可以发现\(\theta\)加减\(2n\) \(\pi\),旋转矩阵\(R\)的值不变。

同样,有如下近似 
(1.10)$$
\exp(\delta \phi^\wedge) \exp(\phi^\wedge) \approx \exp((\phi + J_l^{-1}\delta \phi)^\wedge)
\\
\exp(\phi^\wedge)\exp(\delta \phi^\wedge) \approx \exp((\phi + J_r^{-1}\delta \phi)^\wedge)
\\
\exp((\phi +\delta \phi)^\wedge) \approx \exp((J_l\delta \phi)^\wedge)\exp(\phi^\wedge) \approx \exp(\phi^\wedge)\exp((J_r\delta \phi)^\wedge)$$

下面列举一些常见的性质 
(1.11)$$
J_l(\phi) = J_r(-\phi)
\\
J_l = RJ_r
$$

(1.12)$$
R J_r = J_rR
\\
RJ_l=J_lR
$$

(1.13)$$\begin{aligned}
J_l&=RJ_r=J_rR=J_r^T
\\
J_r&=R^TJ_l=J_lR^T=J_l^T\end{aligned}
$$

(1.14)$$
R(J_r\mathbf{x})^\wedge=(J_l\mathbf{x})^\wedge R
$$

(1.15)$$R = I + \phi^\wedge J_l$$

(1.16)$$I = R^T + \phi^\wedge J_r$$

(1.17)$$\theta = \arccos(\frac{\mathbf{Tr}(R)-1}{2})$$

(1.18)$$\mathbf{a}^\wedge = \frac{\theta}{\sin\theta}\frac{R-R^T}{2}$$

(1.19)$$R\phi = \phi$$

(1.20)$$R\phi^\wedge = \phi^\wedge R$$
对于任意3个元素的向量\(\mathbf{x}, \mathbf{y}\),以及任意\(3 \times 3\)矩阵\(W\)有 
(1.21)$$
\mathbf{x}^\wedge \mathbf{x} = 0
$$

(1.22)$$
\mathbf{x}^T\mathbf{x}^\wedge = \mathbf{0}^T
$$

(1.23)$$
\mathbf{x}^T\mathbf{y}^\wedge \mathbf{x} = 0
$$

(1.24)$$
\mathbf{x}^\wedge \mathbf{y}^\wedge =\mathbf{y}\mathbf{x}^T- (\mathbf{x}^T\mathbf{y})I
$$

(1.25)$$
\mathbf{x}^\wedge \mathbf{y} = – \mathbf{y}^\wedge \mathbf{x}
$$

(1.26)$$
(\mathbf{x}^\wedge)^T = -\mathbf{x}^\wedge
$$

(1.27)$$
(\mathbf{x} + \mathbf{y})^\wedge = \mathbf{x}^\wedge + \mathbf{y}^\wedge
$$

(1.28)$$(\mathbf{x}^{\wedge}\mathbf{y})^\wedge= \mathbf{x}^{\wedge} \mathbf{y}^\wedge – \mathbf{y}^{\wedge} \mathbf{x}^\wedge = \mathbf{y}\mathbf{x}^T – \mathbf{x}\mathbf{y}^T$$

(1.29)$$
\mathbf{x} \times \mathbf{y} = \mathbf{x}^\wedge \mathbf{y}
$$

(1.30)$$
(\mathbf{x}^\wedge \mathbf{y}) \cdot \mathbf{z} =(\mathbf{y}^\wedge \mathbf{z}) \cdot \mathbf{x}=(\mathbf{z}^\wedge \mathbf{x}) \cdot \mathbf{y}
$$
上面公式中的\(\times\)表示叉乘.

(1.31)$$
(W\mathbf{x})^\wedge = \mathbf{x}^\wedge(\mathbf{Tr}(W)I – W) – W^T\mathbf{x}^\wedge
$$

(1.32)$$
\mathbf{x}^\wedge W \mathbf{y}^\wedge =\mathbf{Tr}(W^T\mathbf{yx}^T)I-W^T\mathbf{yx}^T +(\mathbf{Tr}(\mathbf{yx}^T) – \mathbf{yx}^T)(W^T-\mathbf{Tr}(W)I)$$


(1.33)$$
(R\mathbf{x})^\wedge = R\mathbf{x}^\wedge R^T =(2\cos \theta + 1)\mathbf{x}^\wedge -\mathbf{x}^\wedge R – R^T\mathbf{x}^\wedge
$$

(1.34)$$
\exp((R\mathbf{x})^\wedge) = R \exp(\mathbf{x}^\wedge)R^T
$$

(1.35)$$
R\exp((R^T\mathbf{x})^\wedge) = \exp(\mathbf{x}^\wedge)R
$$

(1.36)$$
\frac{\partial \mathbf{a}^\wedge \mathbf{b}}{\partial \mathbf{a}}= -\mathbf{b}^\wedge
$$

(1.37)$$
\log(\exp(\phi_1^\wedge)\exp(\phi_2^\wedge))^\vee=\phi_1+\phi_2+\dfrac{1}{2}\phi_1^\wedge\phi_2+\dfrac{1}{12}(\phi_1^\wedge\phi_1^\wedge\phi_2+\phi_2^\wedge\phi_2^\wedge\phi_1)+……
$$

(1.38)$$R(\mathbf{x}\times\mathbf{y})=(R\mathbf{x})\times(R\mathbf{y})$$

(1.39)$$\mathbf{x}\cdot\mathbf{y}=(R\mathbf{x})\cdot(R\mathbf{y})$$

(1.40)$$(R\mathbf{x})\cdot\mathbf{y}=\mathbf{x}\cdot(R^T\mathbf{y})$$

对于车辆局部坐标系下的角速度\(\mathbf{w}_l\),以及全局坐标系下的角速度\(\mathbf{w}_g\),与车辆位姿本身的旋转矩阵\(R\)(局部到全局坐标系表示)有如下关系

(1.41)$$
\begin{cases} R &=\exp(\phi^\wedge)
\\
\dot{R}&=R\mathbf{w}_l^\wedge=\mathbf{w}_g^\wedge R
\\
\mathbf{w}_l&=J_r\dot{\phi}
\\
\mathbf{w}_g&=J_l\dot{\phi}=R\mathbf{w}_l
\\
\dot{\phi} &=J_l^{-1}\mathbf{w}_g=J_r^{-1}\mathbf{w}_l
\\
\\
\dot{J}_l&=\mathbf{w}_g^\wedge J_l + \dfrac{\partial \mathbf{w}_g}{\partial \phi}=\mathbf{w}_g^\wedge J_l-\mathbf{w}_g^\wedge J_l+R\dfrac{\partial \mathbf{w}_l}{\partial \phi}=R\dfrac{\partial \mathbf{w}_l}{\partial \phi}
\\
&=\dot{R}J_r+R\dot{J}_r
\\
&=\mathbf{w}_g^\wedge R J_r+R\dot{J}_r
\\
&=\mathbf{w}_g^\wedge J_l+R\dot{J}_r
\\
\\
\dot{J}_r &=R^T\dfrac{\partial \mathbf{w}_g}{\partial \phi} \end{cases}
$$

2 SE(3)
一个李群SE(3)中的元素\(T\),总能在李代数中找到包含6个元素的向量\(\xi\)与其对应, 其中\(\xi\)表示如下 
(2.1)$$
\xi =
\begin{bmatrix}
\rho\\ 
\phi
\end{bmatrix}
$$
其中\(\rho\)与\(\phi\)均为大小为3的向量,这里约定如下的符号 
(2.2)$$
\xi^\wedge =
\begin{bmatrix}
\phi^\wedge & \rho \\
\mathbf{0}^T & 0
\end{bmatrix}_{4 \times 4}
\\
\xi =
\begin{bmatrix}
\phi^\wedge & \rho \\
\mathbf{0}^T & 0
\end{bmatrix}^\vee
\\
\xi ^ \curlywedge=\begin{bmatrix} 
\phi ^\wedge & \rho^\wedge\\
0 & \phi^\wedge
\end{bmatrix}_{6 \times 6}
\\
\xi=\begin{bmatrix} 
\phi ^\wedge & \rho^\wedge\\
0 & \phi^\wedge
\end{bmatrix}^ \curlyvee
\\
\mathbf{p}^{\odot}= \begin{bmatrix} 
\mathbf{e}\\
\eta
\end{bmatrix}^{\odot}_{4 \times 1}
= \begin{bmatrix} 
\eta I & -\mathbf{e}^\wedge\\
\mathbf{0}^T & \mathbf{0}^T
\end{bmatrix}_{4 \times 6}
\\
\mathbf{p}^{\circledcirc}= \begin{bmatrix} 
\mathbf{e}\\
\eta
\end{bmatrix}^{\circledcirc}_{4 \times 1}
= \begin{bmatrix} 
\mathbf{0} & \mathbf{e}^\wedge\\
-\mathbf{e}^\wedge & \mathbf{0}
\end{bmatrix}_{6 \times 4}
$$
上面公式中的\(\mathbf{e}\)为大小为3的向量,\(\eta\)为一个实数。

下面建立向量与矩阵的转换关系   
(2.3)$$
\begin{aligned} 
T &= \exp(\xi^\wedge)\\
&=
\begin{bmatrix} 
\exp(\phi^\wedge) & J_l\rho\\ 
\mathbf{0}^T & 1
\end{bmatrix}\\
&=
\begin{bmatrix} 
R & \mathbf{t}\\ 
\mathbf{0}^T & 1
\end{bmatrix}_{4 \times 4}\\
&=I +\xi^\wedge + \frac{1-\cos \theta}{\theta^2}(\xi^\wedge)^2+\frac{\theta- \sin\theta}{\theta^3}(\xi^\wedge)^3
\end{aligned}
$$

(2.4)$$
\xi^\wedge = \log(T)$$

(2.5)$$
T^{-1}=\begin{bmatrix} 
R^T & -R^T\mathbf{t}\\ 
\mathbf{0}^T & 1
\end{bmatrix}
$$

同样的\(T\)也有对应的左右雅克比矩阵\(\mathcal{J}_l,\mathcal{J}_r\) 
(2.6)$$
\mathcal{J}_l = \begin{bmatrix} 
J_l & Q_l\\
0 &  J_l
\end{bmatrix}
\\
\mathcal{J}_r = \begin{bmatrix} 
J_r & Q_r\\
0 &  J_r
\end{bmatrix}
\\
\mathcal{J}_l^{-1} = \begin{bmatrix} 
J_l^{-1} & -J_l^{-1}Q_lJ_l^{-1}\\
0 &  J_l^{-1}
\end{bmatrix}
\\
\mathcal{J}_r^{-1}= \begin{bmatrix} 
J_r^{-1} & -J_r^{-1}Q_rJ_r^{-1}\\
0 &  J_r^{-1}
\end{bmatrix}
$$
其中\(Q_l,Q_r\)表达式如下 
(2.7)$$\begin{cases}
Q_l&=\dfrac{1}{2}\rho^\wedge
\\
\\
&+\dfrac{\theta-\sin\theta}{\theta^3}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge+\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&+\dfrac{\theta^2+2\cos\theta-2}{2\theta^4}(\phi^\wedge\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge\phi^\wedge-3\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&+\dfrac{2\theta-3\sin\theta+\theta\cos\theta}{2\theta^5}(\phi^\wedge\rho^\wedge\phi^\wedge\phi^\wedge+\phi^\wedge\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
Q_r&=-\dfrac{1}{2}\rho^\wedge
\\
\\
&+\dfrac{\theta-\sin\theta}{\theta^3}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge-\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&-\dfrac{\theta^2+2\cos\theta-2}{2\theta^4}(\phi^\wedge\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge\phi^\wedge-3\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&+\dfrac{2\theta-3\sin\theta+\theta\cos\theta}{2\theta^5}(\phi^\wedge\rho^\wedge\phi^\wedge\phi^\wedge+\phi^\wedge\phi^\wedge\rho^\wedge\phi^\wedge)
\end{cases}$$

对于(2.7)式,当轴角趋近0时,有

(2.8)$$\begin{cases}Q_l|_{\theta \rightarrow 0}&\approx \dfrac{1}{2}\rho^\wedge
\\
\\
&+\dfrac{1}{6}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge+\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&+\dfrac{1}{24}(\phi^\wedge\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge\phi^\wedge-3\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&+\dfrac{1}{120}(\phi^\wedge\rho^\wedge\phi^\wedge\phi^\wedge+\phi^\wedge\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
Q_r|_{\theta \rightarrow 0}&\approx -\dfrac{1}{2}\rho^\wedge
\\
\\
&+\dfrac{1}{6}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge-\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&-\dfrac{1}{24}(\phi^\wedge\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge\phi^\wedge-3\phi^\wedge\rho^\wedge\phi^\wedge)
\\
\\
&+\dfrac{1}{120}(\phi^\wedge\rho^\wedge\phi^\wedge\phi^\wedge+\phi^\wedge\phi^\wedge\rho^\wedge\phi^\wedge)
\end{cases}$$

对于(2.8)忽略二阶小量有

(2.9)$$\begin{cases}
Q_l|_{\theta \rightarrow 0}&\approx 
\dfrac{1}{6}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge)+\dfrac{1}{2}\rho^\wedge
\\
\\
Q_r|_{\theta \rightarrow 0}&\approx \dfrac{1}{6}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge)-\dfrac{1}{2}\rho^\wedge
\end{cases}$$

(2.10)$$
Q_r(\xi) = Q_l(-\xi) = RQ_l(\xi) + (J_l\rho)^\wedge RJ_l
$$

对于\(T\),其伴随矩阵定义为 
(2.11)$$\begin{aligned}
\mathcal{T}(T) &=\exp(\xi^ \curlywedge)\\
&=I +\xi^\curlywedge \mathcal{J}_l\\
&=\begin{bmatrix}
R & \mathbf{t}^\wedge R\\
0 & R
\end{bmatrix}_{6 \times 6}\\
&=I
\\&+\frac{3\sin\theta – \theta\cos \theta}{2\theta}\xi^\curlywedge
\\&+
\frac{4-\theta\sin\theta- 4\cos \theta}{2\theta^2}(\xi^\curlywedge)^2
\\&+\frac{\sin\theta – \theta\cos\theta}{2\theta^3}(\xi^\curlywedge)^3
\\&+\frac{2-\theta\sin\theta -2\cos\theta}{2\theta^4}(\xi^\curlywedge)^4
\end{aligned}$$

对于一个比较小的量\(\delta \xi\),对其进行左右扰动有如下结论 
(2.12)$$
\log(\exp(\xi^\wedge)\exp(\delta \xi^\wedge))^\vee \approx \mathcal{J}_r^{-1} \delta \xi +\xi
\\
\log(\exp(\delta \xi^\wedge)\exp(\xi^\wedge))^\vee \approx \mathcal{J}_l^{-1} \delta \xi +\xi
$$
同样,有如下近似 
(2.13)$$
\exp(\delta \xi^\wedge) \exp(\xi^\wedge) \approx \exp((\xi + \mathcal{J}_l^{-1}\delta \xi)^\wedge)
\\
\exp(\xi^\wedge)\exp(\delta \xi^\wedge) \approx \exp((\xi + \mathcal{J}_r^{-1}\delta \xi)^\wedge)
\\
\exp((\xi +\delta \xi)^\wedge) \approx \exp((\mathcal{J}_l\delta \xi)^\wedge)\exp(\xi^\wedge) \approx \exp(\xi^\wedge)\exp((\mathcal{J}_r\delta \xi)^\wedge)$$

下面列举一些常见性质 
(2.14)$$
(-\xi)^\wedge = -\xi^\wedge
$$

(2.15)$$
(-\xi)^ \curlywedge = -\xi^ \curlywedge
$$

(2.16)$$
T\xi^\wedge = \xi^\wedge T
$$

(2.17)$$
\mathcal{T}\xi^\curlywedge= \xi^\curlywedge \mathcal{T}
$$

(2.18)$$
\mathcal{T}\xi = \xi
$$

(2.19)$$
(\xi^\wedge)^4 + \theta^2(\xi^\wedge)^2 = \mathbf{0}
$$

(2.20)$$
(\xi^{ \curlywedge})^5+ 2\theta^2(\xi^{ \curlywedge})^3+\theta^4\xi^{ \curlywedge}=\mathbf{0}
$$
(2.21)$$
T^{-1}= \exp((-\xi)^\wedge) = \exp(-\xi^\wedge) =\begin{bmatrix} 
R^T & -R^T\mathbf{t}\\ 
\mathbf{0}^T &1
\end{bmatrix}
$$

(2.22)$$
\mathcal{T}^{-1}(T) =\mathcal{T}(T^{-1}) =\begin{bmatrix} 
R^T & -R^T\mathbf{t}^\wedge\\
0 & R^T
\end{bmatrix}
$$

(2.23)$$
\mathcal{T}(T_1) \mathcal{T}(T_2) = \mathcal{T}(T_1T_2)$$

(2.24)$$
\mathcal{J}_l(\xi) =\mathcal{T}(T) \mathcal{J}_r(\xi)  =\mathcal{J}_r(-\xi)$$

(2.25)$$\begin{aligned}
\mathcal{J}_l&=I\\
&+\frac{4-\theta \sin \theta – 4 \cos \theta}{2\theta^2}\xi^ \curlywedge \\
&+ \frac{4\theta- 5\sin\theta + \theta\cos \theta}{2\theta^3}(\xi^ \curlywedge)^2\\
&+ \frac{2-\theta \sin\theta- 2 \cos\theta}{2\theta^4}(\xi^ \curlywedge)^3\\
&+ \frac{2\theta -3\sin\theta + \theta \cos \theta}{2\theta^5}(\xi ^ \curlywedge)^4
\end{aligned}$$

对于SE(3)中对应的李代数向量\(\mathbf{h},\mathbf{g}\),以及任意大小为4的向量\(\mathbf{p}\),存在如下性质 

(2.26)$$(\mathbf{h}+\mathbf{g})^\wedge = \mathbf{h}^\wedge + \mathbf{g}^\wedge$$

(2.27)$$
(\mathbf{h}+\mathbf{g})^\curlywedge = \mathbf{h}^\curlywedge + \mathbf{g}^\curlywedge
$$

(2.28)$$
\mathbf{h}^\curlywedge \mathbf{g}= – \mathbf{g}^\curlywedge \mathbf{h}
$$

(2.29)$$(\mathbf{h}^\curlywedge \mathbf{g})^\wedge= \mathbf{h}^\wedge \mathbf{g}^\wedge –  \mathbf{g}^\wedge \mathbf{h}^\wedge $$

(2.30)$$
\mathbf{h}^\curlywedge \mathbf{h}= 0
$$

(2.31)$$
\mathbf{g}^\wedge \mathbf{p}=\mathbf{p}^{\odot}\mathbf{g}
$$

(2.32)$$
\mathbf{p}^T \mathbf{g}^\wedge = \mathbf{g}^T\mathbf{p}^{\circledcirc}
$$

(2.33)$$(\mathcal{T}(T)\mathbf{g})^\wedge = T\mathbf{g}^\wedge T^{-1}$$

(2.34)$$
(\mathcal{T}(T)\mathbf{g})^\curlywedge=\mathcal{T}(T)\mathbf{g}^\curlywedge \mathcal{T}(T)^{-1}
$$

(2.35)$$
\exp((\mathcal{T}(T)\mathbf{g})^\wedge ) = T \exp(\mathbf{g}^\wedge) T^{-1}
$$

(2.36)$$
\exp(\mathbf{g}^\wedge ) T= T \exp((\mathcal{T}(T)^{-1}\mathbf{g})^\wedge)
$$

(2.37)$$
\exp((\mathcal{T}(T)\mathbf{g})^\curlywedge) = \mathcal{T}(T) \exp(\mathbf{g}^\curlywedge) \mathcal{T}(T)^{-1}
$$

由于左右雅克比逆在一些算法中, 可能会被使用到且其计算复杂, 当\(\phi\)均接近0时, 这时左右雅克比逆的表达式, 根据上面的公式,直接将二阶小量忽略掉可得到如下公式

(2.38)$$\begin{cases}
\mathcal{J}_l^{-1} &\approx I- \dfrac{1}{2}\begin{bmatrix}
\phi^\wedge & \rho^\wedge-\dfrac{1}{6}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge)
\\
\mathbf{0} & \phi^\wedge
\end{bmatrix}
\\
\\
\mathcal{J}_r^{-1} &\approx I+ \dfrac{1}{2}\begin{bmatrix}
\phi^\wedge & \rho^\wedge+\dfrac{1}{6}(\phi^\wedge\rho^\wedge+\rho^\wedge\phi^\wedge)
\\
\mathbf{0} & \phi^\wedge
\end{bmatrix}
\end{cases}$$

特别地,当\(\xi\)接近0时,同样忽略二阶小量,则上面公式可变为如下形式

(2.39)$$\begin{cases}
\mathcal{J}_l^{-1} &\approx I- \dfrac{1}{2}\begin{bmatrix}
\phi^\wedge & \rho^\wedge
\\
\mathbf{0} & \phi^\wedge
\end{bmatrix}
\\
\\
\mathcal{J}_r^{-1} &\approx I+ \dfrac{1}{2}\begin{bmatrix}
\phi^\wedge & \rho^\wedge
\\
\mathbf{0} & \phi^\wedge
\end{bmatrix}
\end{cases}$$

由于\(\xi\)由\(\phi\)与\(\rho\)组成,而在实际微扰中,\(\phi\)的值更加容易接近0,而\(\rho\)的值通常与空间位置值量级一致(一个远离0的值),因此(2.39)形式的二阶小量忽略往往不成立,这里更建议使用(2.38)公式来进行近似(即只考虑旋转部分的微扰来获取二阶近似公式)。

根据上面的公式,当\(\phi\)的值接近0时,忽略一阶以上小量,可以大致得到如下结果

(2.40)$$
2\mathcal{J}_r^{-1} -\mathcal{T}(T)\approx I
\\
2\mathcal{J}_r+\mathcal{T}(T)\approx 3I
$$

3 Sim(3)
一个李群Sim(3)中的元素\(S\),总能在李代数中找到包含7个元素的向量\(\zeta\)与其对应, 其中\(\zeta\)表示如下 
(3.1)$$
\zeta  =\begin{bmatrix} 
\rho\\
\phi \\
\sigma
\end{bmatrix}
$$
约定如下符号 
(3.2)$$
\zeta ^\wedge =\begin{bmatrix} 
\sigma I + \phi ^\wedge & \rho\\
\mathbf{0}^T & 0
\end{bmatrix}
$$

对应的指数映射如下 
(3.3)$$
\begin{aligned} 
S&=\exp(\zeta ^\wedge)\\
&=\begin{bmatrix} 
e^\sigma \exp(\phi^\wedge) & J_s\rho\\
\mathbf{0}^T & 1
\end{bmatrix}\\
&=
\begin{bmatrix} 
sR & J_s\rho\\
\mathbf{0}^T & 1
\end{bmatrix}\\
&=
\begin{bmatrix} 
sR & \mathbf{t}\\
\mathbf{0}^T & 1
\end{bmatrix}
\end{aligned}
$$
上面公式中\(J_s\)为 
(3.4)$$
J_s=\frac{e^\sigma – 1}{\sigma}I+
\frac{\sigma e^\sigma \sin\theta + (1-e^\sigma \cos\theta)\theta}{\sigma ^2 + \theta^2}a^\wedge +
(\frac{e^\sigma – 1}{\sigma} – \frac{(e^\sigma \cos\theta – 1)\sigma + \theta e^\sigma \sin\theta}{\sigma ^2 + \theta^2})a^\wedge a^\wedge
$$

(3.5)$$
S^{-1}=\begin{bmatrix} 
\frac{R^T}{s} & – \frac{R^T \mathbf{t}}{s}\\
\mathbf{0}^T & 1
\end{bmatrix}
$$

4 矩阵指数
对于一个方阵\(A\),其矩阵指数为 
(4.1)$$
\exp(A) = I + A + \frac{1}{2!}A^2 + \frac{1}{3!}A^3+… = \sum_{n=0}^{\infty}\frac{1}{n!}A^n
$$
当\(A\)是一个趋近0矩阵的小量时,根据公式(4.1),忽略高阶小量有 
(4.2)$$
\exp(A) \approx I + A
$$
对于方阵\(A\),还有如下性质 
(4.3)$$
\exp(A)^{-1}=\exp(-A)$$

(4.4)$$
\exp(A^T)=\exp(A)^T
$$

(4.5)$$
\det(\exp(A))=\exp(\mathbf{Tr}(A))$$

(4.6)$$
\exp(BAB^{-1})=B\exp(A)B^{-1}
$$

如果方阵\(A\)与\(B\)可交换,即\(AB=BA\),那么它们指数函数满足 
(4.7)$$
\exp(A+B)=\exp(A)\exp(B)$$

5 应用
5.1 相机位姿
可以想象,一个相机,当相对原点先旋转\(R\),然后再进行\(\mathbf{t}\)的偏移。因为原点已知,因此如果知道\(R\)和\(\mathbf{t}\),那么我们就可以确定出相机的位置与朝向,也因此,可以使用如下的\(T\)来表示相机的位姿: 
(5.1)$$
T=\begin{bmatrix} 
R &  \mathbf{t}\\  
\mathbf{0}^T &  1 
\end{bmatrix}
$$

仔细观察(5.1)的形式,可以发现和上面SE(3)的形式类似,因此一般会使用\(\xi\)来描述相机的位姿。

注意: 
(1)当使用\(\xi\)来描述相机位姿时 , \(\xi\)当中的\(\phi\)可以用来描述位姿的旋转,但是 \(\rho\)却不是很好理解。毕竟按这种表示方式,偏移还和旋转相关。 
(2)其实,我们也可以直接使用\(\phi\)与\(\mathbf{t}\)来描述相机位姿,

同样的,如果相对原点先旋转\(R\),然后再进行\(s\)的缩放,最后再进行\(t\)的偏移,则根据这几个参数,同样可以确定出相机的位置与朝向,因此也可以使用如下的\(T\)来表示相机的位姿: 
(5.2)$$
T=\begin{bmatrix} 
sR &  \mathbf{t}\\  
\mathbf{0}^T &  1 
\end{bmatrix}
$$
仔细观察(5.2)的形式,可以发现和上面Sim(3)的形式类似,因此一般会使用\(\zeta\)来描述相机的位姿。

对比(5.1)与(5.2)可以发现,(5.1)是(5.2)的特殊形式,对于没有特殊说明,当使用到(5.1)的表示形式时,这里都默认缩放系数\(s=1\)。

5.2 相机光心
空间中一个3D点\(\mathbf{p}\)变换到相机\(T\)坐标系下到位姿为\(\mathbf{p}'\),这里我们约定:变换的方式为先旋转,再缩放,最后再偏移。则3D点的变换可表示为 
(5.3)$$
\mathbf{p}'=sR\mathbf{p}+\mathbf{t}
$$

一般地,我们约定相机在自己相机坐标系下的坐标为原点,如果相机在空间中地坐标为\(\mathbf{c}\)(光心),那么根据公式(5.3)它们有如下关系 
(5.4)$$
\mathbf{0}=sR\mathbf{c}+\mathbf{t}
$$
通过公式(5.4)可以得到光心的表达式 
(5.5)$$
\mathbf{c}=-\frac{1}{s}R^T\mathbf{t}
$$
特别地,当时缩放系数为1时,则光心可表示为 
(5.6)$$
\mathbf{c}=-R^T\mathbf{t}
$$

5.3 点的变换
公式(5.3)可以写成如下形式 
(5.7)$$
\begin{bmatrix} 
\mathbf{p}'\\  

\end{bmatrix}
=
\begin{bmatrix} 
sR & \mathbf{t}\\  
\mathbf{0}^T & 1 
\end{bmatrix}
\begin{bmatrix} 
\mathbf{p}\\  

\end{bmatrix}
= T
\begin{bmatrix} 
\mathbf{p}\\  

\end{bmatrix}
$$
对于一个点\(\mathbf{p}\),先经过\(T_0\)的变换,再经过\(T_1\)的变换,最终得到\(\mathbf{p}'\),则可表达为 
(5.8)$$
\begin{bmatrix} 
\mathbf{p}'\\  

\end{bmatrix}
= T_1T_0
\begin{bmatrix} 
\mathbf{p}\\  

\end{bmatrix}
= T'
\begin{bmatrix} 
\mathbf{p}\\  

\end{bmatrix}
$$
因此中间的变换可表示为 
(5.9)$$
T' = T_1T_0
$$

5.4 位姿变换
我们知道,对于空间中一堆3D点云,我们按照(5.7)让点云在刚体变换\(T\)下进行坐标的变换,可以发现,点云中各个点的相对位置不会改变。
同样地,对于空间中一些相机和一些点云组成的整体,现在让它们进行一个刚体变换\(T\),我们希望变换之后,相机之间,以及相机与点云之间,它们的相对位置信息不会改变。
首先,相机与点的相对位置关系不变可以体现在点在相机坐标系下的坐标不变,也因此,如果相机原始位姿为\(T_0\),经过\(T\)变换后的位姿\(T_1\),对于一个空间原始的空间点\(\mathbf{p}\)有  
(5.10)$$
T_0
\begin{bmatrix} 
\mathbf{p}\\  
1
\end{bmatrix}
=
\begin{bmatrix} 
\mathbf{p}'\\  
1
\end{bmatrix}
=T_1(T\begin{bmatrix} 
\mathbf{p}\\  
1
\end{bmatrix})
=T_1T\begin{bmatrix} 
\mathbf{p}\\  
1
\end{bmatrix}
$$

从上面的公式可以得到 
(5.11)$$
T_1 = T_0T^{-1}
$$
按照(5.11)式定义的相机变换,可以满足相机与空间3D点之间的相对位置关系,而对于相机间的相对位置关系,要使其既满足(5.11)的变换关系,又能使得相对位置关系不因(5.11)的变换而改变,则对于位姿为\(T_i\)与\(T_j\)的相机,它们的相对位置关系\(T_{ij}\),可以定义为如下形式 
(5.12)$$
T_{ij}=T_iT_j^{-1}
$$
公式(5.12)就是一些PGO优化中边的测量形式。
需要注意的是,文献[1]中提到的另一种测量形式如下  
(5.13)$$
T_{ij}=T_i^{-1}T_j
$$
对于(5.13)这种测量形式不建议使用,因为这种测量在优化中,如果有一个顶点不变的情况下,理论上没有问题,但是,往往我们更加注重的是满足位姿与点可共同变换的测量,因此建议尽量使用(5.12)的形式,因为这个公式具有更强的物理意义!

6 矩阵向量运算
下面记录一些矩阵向量的公式[3],所写公式中点矩阵变量/向量均满足维度要求下成立 
(6.1)$$
ABC=A(BC)$$

(6.2)$$
(AB)^T=B^TA^T
$$

(6.3)$$
(A^T)^{-1}=(A^{-1})^T
$$

(6.4)$$
(AB)^{-1}=B^{-1}A^{-1}
$$

(6.5)$$
\mathbf{Tr}(ABC) = \mathbf{Tr}(BCA) = \mathbf{Tr}(CAB)$$

(6.6)$$
\frac{\partial A^{-1}}{\partial \mathbf{x}}=A^{-1}\frac{\partial A }{\partial \mathbf{x}}A^{-1}
$$

(6.7)$$
\frac{\partial AB}{\partial \mathbf{x}}=\frac{\partial A}{\partial \mathbf{x}}B+A\frac{\partial B}{\partial \mathbf{x}}
$$

(6.8)$$\frac{\partial (aB)}{\partial \mathbf{x}}=B \otimes \frac{\partial a}{\partial \mathbf{x}} + a\frac{\partial B}{\partial \mathbf{x}}$$

上面公式中的\(a\)为一个数,\(\otimes\)为 Kronecker product运算符。

(6.9)$$\frac{\partial \mathbf{Tr}(A)}{\partial \mathbf{x}}=\mathbf{Tr}(\frac{\partial A}{\partial \mathbf{x}})$$

(6.10)$$\frac{\partial \mathbf{x}^TA \mathbf{x}}{\partial \mathbf{x}}=(A+A^T)\mathbf{x}$$

(6.11)$$\det(I+\mathbf{a}\mathbf{b}^T)=1+\mathbf{a}^T\mathbf{b}$$

(6.12)$$
\det(A^{-1})=\frac{1}{\det(A)}
$$

(6.13)$$
||A||_F=\sqrt{\mathbf{Tr}(AA^T)}
$$

(6.14)$$
||\mathbf{x}+\mathbf{y}||_2^2 = ||\mathbf{x}||_2^2 + ||\mathbf{y}||_2^2 + 2\mathbf{x}^T\mathbf{y}
$$

(6.15)$$
(I+\mathbf{xy}^T)^{-1}=I-\frac{\mathbf{xy}^T}{1+\mathbf{y}^T\mathbf{x}}
$$
Sherman–Morrison 公式 
(6.16)$$
(A+\mathbf{xy}^T)^{-1}=A^{-1}-\frac{A^{-1}\mathbf{xy}^TA^{-1}}{1+\mathbf{y}^TA^{-1}\mathbf{x}}
$$

Sherman–Morrison–Woodbury公式 
(6.17)$$
\begin{aligned} 
(A+UCV)^{-1}&=A^{-1}-A^{-1}U(C^{-1}+VA^{-1}U)^{-1}VA^{-1}
\\
AB(D+CAB)^{-1}&=(A^{-1}+BD^{-1}C)^{-1}BD^{-1}
\\
(D+CAB)^{-1}CA&=D^{-1}C(A^{-1}+BD^{-1}C)^{-1}
\\
(A+BC)^{-1}&=A^{-1}-A^{-1}B(I+CA^{-1}B)^{-1}CA^{-1}
\end{aligned}
$$

(6.18)$$
\begin{aligned} 
(A+B)^{-1} &= A^{-1} – A^{-1}B(B+BA^{-1}B)^{-1}BA^{-1}
\\
&=A^{-1}-A^{-1}(I+BA^{-1})^{-1}BA^{-1}
\\
&=A^{-1}-(A+AB^{-1}A)^{-1}
\end{aligned}
$$

(6.19)$$
\begin{bmatrix} 
A & 0\\  
C &  D 
\end{bmatrix}
^{-1}
=
\begin{bmatrix} 
A^{-1} & 0\\  
-D^{-1}CA^{-1} &  D^{-1} 
\end{bmatrix}
$$

(6.20)$$
\begin{bmatrix} 
A & B\\  
C & D
\end{bmatrix}^{-1}
=
\begin{bmatrix} 
-C^{-1}D(B-AC^{-1}D)^{-1} & C^{-1}+C^{-1}D(B-AC^{-1}D)^{-1}AC^{-1}\\  
(B-AC^{-1}D)^{-1} & -(B-AC^{-1}D)^{-1}AC^{-1}
\end{bmatrix}
$$

(6.21)$$ \begin{bmatrix} 
A & B\\
C & D
\end{bmatrix}^{-1}
=
\begin{bmatrix} 
-(C-DB^{-1}A)^{-1}DB^{-1} & (C-DB^{-1}A)^{-1}
\\
B^{-1}+B^{-1}A(C-DB^{-1}A)^{-1}DB^{-1} & -B^{-1}A(C-DB^{-1}A)^{-1}
\end{bmatrix}
$$

(6.22)$$
\begin{bmatrix} 
A & B\\  
C & D
\end{bmatrix}^{-1}
=
\begin{bmatrix} 
A^{-1}+A^{-1}B(D-CA^{-1}B)^{-1}CA^{-1} & -A^{-1}B(D-CA^{-1}B)^{-1}
\\
-(D-CA^{-1}B)^{-1}CA^{-1} & (D-CA^{-1}B)^{-1}
\end{bmatrix}
$$

(6.23)$$
\begin{bmatrix} 
A & B\\  
C & D
\end{bmatrix}^{-1}
=
\begin{bmatrix} 
(A-BD^{-1}C)^{-1} &-(A-BD^{-1}C)^{-1}BD^{-1}
\\
-D^{-1}C(A-BD^{-1}C)^{-1} & D^{-1}+D^{-1}C(A-BD^{-1}C)^{-1}BD^{-1}
\end{bmatrix}
$$

如果\(A\)为对称矩阵,则必然存在一个旋转矩阵\(R\)使得 
(6.24)$$
R^TAR=\mathbf{diag}(d_0,d_1,d_2)$$

对于一个复数方阵\(A\),有如下性质 
(6.25)$$
\det(\exp(A))=\exp(\mathbf{Tr}(A))
\\
\exp(A)^{-1}=\exp(-A)
\\
\exp(A^T)=\exp(A)^T
\\
\exp(BAB^{-1})=B\exp(A)B^{-1}
$$

如果方阵\(A\)与\(B\)可交换,即\(AB=BA\),那么它们指数函数满足 
(6.26)$$
\exp(A+B)=\exp(A)\exp(B)$$

对于同维度的向量\(\mathbf{x}\),\(\mathbf{y}\),以及一个向量\(\mathbf{z}\),它们之间满足如下关系

(6.27)$$\begin{aligned}
\mathbf{x}^T \mathbf{yz}
&= \mathbf{y}^T \mathbf{xz}\\
&= \mathbf{z} \otimes \mathbf{x}^T \mathbf{y}\\
&= \mathbf{z} \otimes \mathbf{y}^T \mathbf{x}
\end{aligned}$$

 

(6.28)$$\begin{bmatrix} 
R_2 & t_2
\\
0 & 1
\end{bmatrix}^{-1} \begin{bmatrix} 
R_1 & t_1
\\
0 & 1
\end{bmatrix} = \begin{bmatrix} 
R_2^TR_1 & R_2^T(t_1-t_2)
\\
0 & 1
\end{bmatrix}
$$

(6.29) 在右手坐标系中,大拇指方向指向z轴正向,在x-y平面绕z轴逆时针旋转\(\theta\)角度得到的旋转矩阵为如下的\(R_z\)(旋转前的状态,通过\(R_z\)变换为旋转后的状态,即旋转矩阵乘以原坐标)

$$R_z = \begin{bmatrix} 
\cos \theta & -\sin \theta & 0
\\
\sin \theta & \cos\theta & 0
\\
0 & 0& 1
\end{bmatrix}
$$

(6.30) 在右手坐标系中,大拇指方向指向z轴正向,在x-y平面绕z轴顺时针旋转\(\theta\)角度得到的旋转矩阵为如下 \(R_z\)(旋转前的状态,通过\(R_z\)变换为旋转后的状态)

$$R_z = \begin{bmatrix} 
\cos \theta & \sin \theta & 0
\\
-\sin \theta & \cos\theta & 0
\\
0 & 0& 1
\end{bmatrix}
$$

(6.31) 在右手坐标系中,大拇指方向指向x轴正向,在y-z平面绕x轴逆时针旋转\(\theta\)角度得到的旋转矩阵为如下\(R_x\)(旋转前的状态,通过\(R_x\)变换为旋转后的状态)

$$R_x = \begin{bmatrix} 
1 & 0 & 0
\\
0 & \cos\theta & -\sin\theta
\\
0 & \sin\theta& \cos\theta
\end{bmatrix}
$$

(6.32) 在右手坐标系中,大拇指指向y轴正向,在z-x平面绕y轴逆时针旋转\(\theta\)角度得到的旋转矩阵为如下\(R_y\)(旋转前的状态,通过\(R_y\)变换为旋转后的状态)

$$R_y = \begin{bmatrix} 
\cos \theta & 0 & \sin\theta
\\
0 & 1 & 0
\\
-\sin\theta & 0& \cos\theta
\end{bmatrix}
$$

7 叉乘运算

设向量\(\mathbf{a}=[a_0,a_1,a_2]^T、\mathbf{b}=[b_0,b_1,b_2]^T、\mathbf{c}、\mathbf{d}\)为大小为3的向量。这里定义叉乘  

(7.1)

$$\mathbf{a} \times \mathbf{b} =\begin{bmatrix}
a_1b_2-a_2b_1
\\
a_2b_0-a_0b_2
\\
a_0b_1-a_1b_0
\end{bmatrix}$$

上面公式里以右手规则来确定叉乘后的向量方向:右手食指指向\(\mathbf{a}\)的方向,右手中指指向\(\mathbf{b}\)的方向,右手大拇指指向\(\mathbf{a}\times\mathbf{b}\)的方向。

对于叉乘有如下性质

(7.2)$$\mathbf{a} \times \mathbf{0} = \mathbf{0}$$

(7.3)$$\mathbf{a} \times \mathbf{a} = \mathbf{0}$$

(7.4)$$\mathbf{a} \times \mathbf{b} = -\mathbf{b} \times \mathbf{a}$$

(7.5)$$(k\mathbf{a}) \times \mathbf{b} = \mathbf{a} \times (k\mathbf{b})=k(\mathbf{a}\times\mathbf{b})$$

(7.6)$$\mathbf{a} \times (\mathbf{b} + \mathbf{c})= \mathbf{a} \times \mathbf{b}+\mathbf{a} \times \mathbf{c}$$

(7.7)$$(\mathbf{b} + \mathbf{c})\times \mathbf{a}=\mathbf{b} \times \mathbf{a} + \mathbf{c} \times \mathbf{a}$$

(7.8)$$\mathbf{a} \cdot (\mathbf{b} \times \mathbf{c}) = (\mathbf{a} \times \mathbf{b})\cdot \mathbf{c}=(\mathbf{c}\times\mathbf{a})\cdot \mathbf{b}$$

(7.9)$$\mathbf{a} \times (\mathbf{b} \times \mathbf{c}) = (\mathbf{c}\times \mathbf{a})\times\mathbf{b}+(\mathbf{a}\times \mathbf{b})\times \mathbf{c}=(\mathbf{a} \cdot\mathbf{c})\mathbf{b} – (\mathbf{a} \cdot\mathbf{b})\mathbf{c}$$

(7.10)$$(\mathbf{a} \times \mathbf{b})\times \mathbf{c} = \mathbf{a}\times (\mathbf{b}\times \mathbf{c}) – \mathbf{b}\times(\mathbf{a}\times\mathbf{c})=(\mathbf{a}\cdot\mathbf{c})\mathbf{b}-(\mathbf{b}\cdot\mathbf{c})\mathbf{a}$$

(7.11)$$||\mathbf{a}\times\mathbf{b}||^2+(\mathbf{a}\cdot \mathbf{b})^2=||\mathbf{a}||^2||\mathbf{b}||^2$$

(7.12)$$\dfrac{\partial (\mathbf{a} \times \mathbf{b})}{\partial t}=\dfrac{\partial \mathbf{a}}{\partial t}\times \mathbf{b}+\mathbf{a} \times \dfrac{\partial \mathbf{b}}{\partial t}$$

(7.13)$$|\mathbf{a} \times \mathbf{b}|=|\mathbf{a}||\mathbf{b}|\sin \theta_{a,b}$$

公式里\(\theta_{a,b}\)为向量\(\mathbf{a}\)和\(\mathbf{b}\)的夹角。

8 常见矩阵结果

(8.1)$$\begin{bmatrix} a_{00}&a_{01} \\ a_{10}&a_{11} \end{bmatrix}^{-1}=\dfrac{1}{a_{00}a_{11} – a_{01}a_{10}} \begin{bmatrix} a_{11}&-a_{01} \\ -a_{10}&a_{00} \end{bmatrix}$$

(8.2)$$\begin{bmatrix} a_{00}&a_{01}&a_{02} \\ a_{10}&a_{11}&a_{12} \\ a_{20}&a_{21}&a_{22} \end{bmatrix}^{-1}=\dfrac{1}{a_{00}a_{11}a_{22}-a_{00}a_{12}a_{21}-a_{01}a_{10}a_{22}+a_{01}a_{12}a_{20}+a_{02}a_{10}a_{21}-a_{02}a_{11}a_{20}} \begin{bmatrix} a_{11}a_{22}-a_{12}a_{21}&a_{02}a_{21}-a_{01}a_{22}&a_{01}a_{12}-a_{02}a_{11} \\ a_{12}a_{20}-a_{10}a_{22}&a_{00}a_{22}-a_{02}a_{20}&a_{02}a_{10}-a_{00}a_{12} \\ a_{10}a_{21}-a_{11}a_{20}&a_{01}a_{20}-a_{00}a_{21}&a_{00}a_{11}-a_{01}a_{10} \end{bmatrix}$$

9 一些其它等式

(9.1)$$\begin{cases} \sum_{i=1}^ni&=\dfrac{1}{2}n(n+1) \\ \\ \sum_{i=1}^ni^2&=\dfrac{1}{6}n(n+1)(2n+1) \\ \\ \sum_{i=1}^ni^3&=\dfrac{1}{4}n^2(n+1)^2 \\ \\ \sum_{i=1}^ni^4&=\dfrac{1}{30}n(n+1)(2n+1)(3n^2+3n-1) \\ \\ \sum_{i=1}^ni^5&=\dfrac{1}{12}n^2(n+1)^2(2n^2+2n-1) \\ \\ \sum_{i=1}^ni^6&=\dfrac{1}{42}n(n+1)(2n+1)(3n^4+6n^3-3n+1) \\ \\ \sum_{i=1}^ni^7&=\dfrac{1}{24}n^2(n+1)^2(3n^4+6n^3-n^2-4n+2) \\ \\ \sum_{i=1}^ni^8&=\dfrac{1}{90}n(10n^8+45n^7+60n^6-42n^4+20n^2-3) \end{cases}$$

参考
[1] 高翔,张涛. 视觉SLAM十四讲:从理论到实践[M]. 电子工业出版社,2017.
[2] 高翔,谢晓佳. 机器人学中的状态估计[M]. 西安交通大学出版社,2018
[3] Petersen K B, Pedersen M S. The matrix cookbook[M]. Technical University of Denmark, 2012

 

3D点云配准之刚体变换

前言

3D点云配准,在无人驾驶领域有一定分量。自己转行之后,开始接触这个领域。点云配准的一个非常难的点就是如何快速查找到比较可靠的匹配对。本文这里不展开讨论,在这里,仅叙述在找到匹配对后,如何反算出对应的刚体变换参数(注意,这里仅考虑一个刚体变换)。

 

问题

假如有一个3D点\(p_1\),在经过如下的刚体变换后得到新的位置点\(p_2\)

(1)$$p_2=sRp_1+t  $$

公式(1)中的\(s\)为缩放系数,\(R\)为旋转矩阵,\(t\)为偏移。

现在我们需要求解的问题是,已知有\(n\)组\((p_1,p_2)_i\),它们每一组的匹配对应的权重为\(w_i\),现在需要反求解公式(1)当中的变换参数。

 

求解

对于求解一个优化问题的未知参数,可以按常规方法构建一个最小二乘问题

(2)$$\begin{align*}
min E &=\sum_{i=1}^{n}\left \| w_i(p_{2i} – (sRp_{1i}+t)) \right \|_2^2 \\ 
 &=\sum_{i=1}^{n}w_i^2\left \| p_{2i} – sRp_{1i}-t \right \|_2^2 \\ 
&=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1}) +  \bar{p_2} -sR \bar{p_1}- t \right \|_2^2\\
&=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1})\right \|_2^2 +
\\ &\quad \sum_{i=1}^{n}w_i^2\left \| \bar{p_2} -sR \bar{p_1}- t \right \|_2^2 +
\\ &\quad 2\sum_{i=1}^{n}w_i^2 (p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1}))^T(\bar{p_2} -sR \bar{p_1}- t)
\end{align*}$$

上面公式里的\(\bar{p_1}\)、\(\bar{p_2}\)为引入的消去变量。从公式(2)可以发现最小化目标函数由3部分组成,由于前两项均为非负数,而第3项无法保证值为非负,因此,可以考虑将第3项变成一个与\(s\)、\(R\)、\(t\)无关的常量,进而达到先消去第3部分的表达式目的,这里令它等于0,即

(3)$$\begin{align*}
0 &=2\sum_{i=1}^{n}w_i^2 (p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1}))^T(\bar{p_2} -sR \bar{p_1}- t)\\
  &=2[\sum_{i=1}^{n}w_i^2 (p_{2i} –  \bar{p_2})- sR\sum_{i=1}^{n}w_i^2 (p_{1i}- \bar{p_1})]^T(\bar{p_2} -sR \bar{p_1}- t)\\
&=2[(\sum_{i=1}^{n}w_i^2p_{2i}-\bar{p_2}\sum_{i=1}^{n}w_i^2)-sR(\sum_{i=1}^{n}w_i^2p_{1i}-\bar{p_1}\sum_{i=1}^{n}w_i^2)]^T(\bar{p_2} -sR \bar{p_1}- t)
\end{align*}$$

要满足(3)式,很明显,对于任意的\(s\)、\(R\)、\(t\)均使得公式(3)满足时(相当于公式(3)中的值与\(s\)、\(R\)、\(t\)无关),\(\bar{p_1}\)、\(\bar{p_2}\)值必然满足如下条件

(4)$$\sum_{i=1}^{n}w_i^2p_{2i}-\bar{p_2}\sum_{i=1}^{n}w_i^2=0,\sum_{i=1}^{n}w_i^2p_{1i}-\bar{p_1}\sum_{i=1}^{n}w_i^2=0$$

(5)$$\bar{p_1}=\sum_{i=1}^{n}w_i^2p_{1i} / \sum_{i=1}^{n}w_i^2,\bar{p_2}=\sum_{i=1}^{n}w_i^2p_{2i} / \sum_{i=1}^{n}w_i^2$$

当我们构造的\(\bar{p_1}\)、\(\bar{p_2}\)满足公式(5)时,这时(2)式的最小二乘的目标函数就变成了如下形式

(6)$$\begin{align*}
min E &=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1})\right \|_2^2 + \sum_{i=1}^{n}w_i^2\left \| \bar{p_2} -sR \bar{p_1}- t \right \|_2^2\\
&=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1})\right \|_2^2+ \left \| \bar{p_2} -sR \bar{p_1}- t \right \|_2^2 \sum_{i=1}^{n}w_i^2
\end{align*}$$

从公式(6)可以发现,目标最小值由两部分组成,第一部分只和\(R\)、\(s\)有关,而第二项当得到\(R\)、\(s\)后,可以马上确定出使第2项为0的最佳\(t\),因此求解这类问题可以先优化求解第一部分,然后再根据第二部分确定出最佳\(t\)。

注意:这里之所以可以分成两部分求解,因为公式的第二部分总可以找到一个\(t\)使得其表达式值为0对于其它可类似拆分的优化问题,如果不能确保第另一项为0,则不可采用这种方式求解,而需要采用来回迭代求解的方式进行试解

 

根据公式(6)则展开第一部分有

(7)$$\begin{align*}
min E_1 &=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}- sR(p_{1i}- \bar{p_1})\right \|_2^2 \\
&= \sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}\right \|_2^2 + \sum_{i=1}^{n}w_i^2\left \|  sR(p_{1i}- \bar{p_1})\right \|_2^2 – 2\sum_{i=1}^{n}w_i^2(p_{2i} –  \bar{p_2})^T sR(p_{1i}- \bar{p_1}) \\
&= \sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}\right \|_2^2 + \sum_{i=1}^{n}w_i^2   s^2(p_{1i}- \bar{p_1})^TR^TR(p_{1i}- \bar{p_1}) – 2s\sum_{i=1}^{n}w_i^2(p_{2i} –  \bar{p_2})^T R(p_{1i}- \bar{p_1})\\
&=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}\right \|_2^2 + \sum_{i=1}^{n}w_i^2   s^2(p_{1i}- \bar{p_1})^T(p_{1i}- \bar{p_1}) – 2s\sum_{i=1}^{n}tr(w_i^2R(p_{1i}- \bar{p_1})(p_{2i}- \bar{p_2})^T)\\
&=\sum_{i=1}^{n}w_i^2\left \| p_{2i} –  \bar{p_2}\right \|_2^2 + s^2\sum_{i=1}^{n}w_i^2\left \| p_{1i} –  \bar{p_1}\right \|_2^2 – 2s \times tr(R \sum_{i=1}^{n}w_i^2(p_{1i}- \bar{p_1})(p_{2i}- \bar{p_2})^T)
\end{align*}$$

公式(7)可以看出,这里求最小值又可以分成三部分,而第一部分为常数,因此实际最小化就变成了后面两部分。

这里设置一个如下的中间变量\(H\)

(8)$$ H=\sum_{i=1}^{n}w_i^2(p_{1i}- \bar{p_1})(p_{2i}- \bar{p_2})^T$$

以及变量\(a\)

(9)$$ a = \sum_{i=1}^{n}w_i^2\left \| p_{1i} –  \bar{p_1}\right \|_2^2$$

则(7)式最终等价求如下函数的最小形式

(10)$$\begin{align*}
min E_1'&=  as^2 – 2s \times tr(RH)\\
 &= a(s-\frac{tr(RH)}{a})^2-\frac{tr(RH)^2}{a}\\ 
\end{align*}$$

分析公式(10)可知,要使整体表达式最小,则前面二次方的部分要最小,后面被减部分要最大,从表达式可知,总存在一个\(s\)使得二次方表达式部分取得最小的0值,这时可以得到

(11)$$s=\frac{tr(RH)}{a}$$

而对于(10)式,由于要使得最小,则后面的负数项则必须最大,即(10)要求得最小值,等价于(12)算式求得最大值

(12)$$max E= tr(RH)^2$$

如果尺度\(s\)有大小范围限制,要求最小值,其等价于(13)算式求得最大值(本文这里仅把\(s\)当做一个纯优化求解的可正可负的变量)。

(13)$$max E= sign(s)tr(RH)$$

接下来,进一步对\(H\)进行奇异值分解,即

(14)$$H=USV^T$$

则求迹时可变成如下等式
(15)$$\begin{align*}
tr(RH) &= tr(RUSV^T) \\ 
&= tr(V^TRUS)\\ 
&= tr(AS)
\end{align*}$$

上面的\(A=V^TRU\),很明显\(A\)为一个正交矩阵,由于\(S\)由奇异值分解时求解实对称矩阵特征值而来,因此对角矩阵\(S\)上的元素均不小于0,这里假设\(S\)对角上的元素为\(\sigma_0 \)、\(\sigma_1 \)、\(\sigma_2\),而\(A\)上的对角元素为\(a_{00}\)、\(a_{11}\)、\(a_{22}\),由于\(A\)的正交性,这几个元素取值范围在[-1,1]之间,则最终(15)式可变成

(16)$$tr(RH)=tr(AS)=a_{00}\sigma_0 +a_{11}\sigma_1 +a_{22}\sigma_2$$

从上面可知,(12)算式要求得最大值,则公式(16)的迹必须取极大或者极小值。由于上面\(\sigma_0 \)、\(\sigma_1 \)、\(\sigma_2\)均不小于0,因此(16)式要取得极值,则必然\(a_{00}\)、\(a_{11}\)、\(a_{22}\)同号时方可取极值。明显地,当\(a_{00}\)、\(a_{11}\)、\(a_{22}\)均取最大值1或者最小值-1时,这时可以使得\(tr(RH) \)的值取得最大,这时正好有一个唯一与之对应的矩阵\(\pm E\),由于取正负单位阵得到的结果使得(12)得到的最大值一致,因此其目前有2组可能解,即

(17)$$ V^TRU = A = \pm E \Leftrightarrow  R= \pm VU^T$$

根据旋转矩阵的性质,其行列式为1,因此可以从(17)的两组解中选择出满足旋转矩阵性质的解,即

(18)$$ R=\left | UV \right | VU^T$$

结合公式(11),则最终,当得到\(R\)、\(s\)后,可根据(6)式得到\(t\)

(19)$$t = \bar{p_2} -sR \bar{p_1}$$

需要注意,上面的求解过程认为\(s\)为没有限制的参数,当\(s \neq 0\)且为大小限制的参数时,则优化目标函数由(12)式变化成了(13)式。根据上面的方法可以得到最佳\(R\)的结果如下

(20)$$R=sign(s)VU^T$$

对于(20)的结果,如果\(\left | R \right | =1\),则这时求得的旋转矩阵则为最终解。

当(20)式得到的\(\left | R \right |  \neq 1\)时,明显这个结果不是真实的解,这时不能使用(18)算式进行修正,因为这样求解时得到的结果完全与优化目标结果相反。

当限制\(s>0\)时,通过公式(10)可以发现,要求最优值,则要求

(21)$$s=\frac{tr(RH)}{a}>0 \Rightarrow tr(RH)>0$$

然后求解满足如下条件的最大值

(22)$$maxE'=tr(RH)$$

而直接指定缩放系数\(s\)的值,并且这个\(s>0\)时,根据(10)式可知最终也是求解(22)的最大值,因为这时的\(s\)已知,因此这种情况没有(21)的限制。

对于(22)式,最终都会归结到(16)式,由于其元素值本身的性质,在(16)式中,总能找到大于0的结果,因此最终的优化结果均能满足(21)式的要求,因此不管是直接给出一个大于0的缩放系数,还是限制待求解的缩放系数为正,最终优化目标均变成(22)式,根据文献当(16)式中的\(A\)满足如下条件时,可以取得最佳值

(23)$$A=\begin{bmatrix}
1 & 0 & 0\\ 
0 & 1 & 0\\ 
0 & 0 & -1
\end{bmatrix}$$

这时,\(R\)结果如下

(24)$$ R=VAU^T$$

不失一般性,当限制缩放系数为正数时,则最终的旋转矩阵结果如下

(25)$$ R=VZU^T$$

其中\(Z\)取值如下

(26)$$Z=\begin{bmatrix}
1 & 0 & 0\\ 
0 &  1& 0\\ 
 0&0  & \left | UV \right |
\end{bmatrix}$$

步骤

\(s\)为无限制参数

 最终求解带权重的变换矩阵的步骤如下:

第1步:根据公式(5)计算\(\bar{p_1}\)、\(\bar{p_2}\)

第2步:根据公式(8)计算\(H\)

第3步:根据公式(14)对\(H\)进行奇异值分解

第4步:根据公式(18)求得旋转矩阵\(R\)

第5步:根据公式(9)计算\(a\),如果忽略尺度,则跳过本步骤

第6步:根据公式(11)计算尺度\(s\),如果忽略尺度,则跳过本步骤

第7步:根据公式(19)计算\(t\)

 

\(s\)限制为正数

第1步:根据公式(5)计算\(\bar{p_1}\)、\(\bar{p_2}\)

第2步:根据公式(8)计算\(H\)

第3步:根据公式(14)对\(H\)进行奇异值分解

第4步:根据公式(25)求得旋转矩阵\(R\)

第5步:根据公式(9)计算\(a\),如果忽略尺度,则跳过本步骤

第6步:根据公式(11)计算尺度\(s\),如果忽略尺度,则跳过本步骤

第7步:根据公式(19)计算\(t\)

 

PS:可以发现,权重系数均是以\(w_i^2\)形式出现,如果我们设置的权重系数\(w_i\)本身为实际权重的平方项,则这时可以直接使用\(w_i\)代替\(w_i^2\)

数据拟合愚见

以下为个人愚见,希望对有缘人有帮助,如有不足,欢迎拍砖!

1、数据拟合问题,可简单转换成方程组求解的问题,且也可简单分为线性与非线性问题。

2、对于线性问题,可简单看成一个超定方程的求解,这时采用最小二乘算法求解即可。
比如:需要拟合的方程为$$y=ax_1+bx_2+cx_3+…$$
其中\(y,x_1,x_2,x_3,…\)为采集的数据,需要拟合的变量为\(a,b,c,…\)
如果把上面的例子化成一个方程组可表示如下的超定线性方程组
$$\\y_1=x_{1,1}+bx_{2,1}+cx_{3,1}+…\\y_2=x_{1,2}+bx_{2,2}+cx_{3,2}+…\\y_3=x_{1,3}+bx_{2,3}+cx_{3,3}+…\\…$$
PS:对于一些表面上是非线性问题却可转换成线性拟合的问题,建议先转换成线性拟合得到结果。有些人可能会说,转换后拟合的数据会放大误差,对于这类会放大误差的问题,还是建议先转换成线性拟合,得到一个解,然后将这个解作为初始值带入非线性拟合问题当中进行求解,也不失为一个很好的方法。
例子:需要拟合的方程为$$y=asin(x_1)+bx_{2}^{2}+cx_1x_2$$
其中\(y,x_1,x_2\)为采集的数据,\(a,b,c\)为需要拟合的变量,这时公式可进行简单的变换
$$
\begin{align}
z_1&=sin(x_1)\\
z_2&=x_{2}^{2}\\
z_3&=x_1x_2\\
y&=az_1+bz_2+cz_3 
\end{align}
$$
对上面的公式,可以很快求解出\(a,b,c\)

3、非线性问题的拟合
    可能我们大部分遇到的都是非线性拟合问题,非线性拟合问题又可简单看成一种简单的循环“迭代问题”。这里先简单说一下在求解非线性问题采用算法时注意的地方:(1)智能算法(比如:遗传算法、粒子群算法……)一般都可简单说成是全局优化算法,当然这里所谓的全局优化算法只能说是相对的,毕竟目前还查不到有什么算法能确定出任意的非线性拟合问题的解。(2)剩下的就是一些局部最优算法了(比如:Broyden形式的求解算法、LM算法……),这类算法的特点是在解的附近时能很快收敛到解。也因此,配合好全局与局部算法是一个不错的选择!(3)为了加快找到解的速度,可考虑采用多线程求解方式【如果你的算法允许】。即使你不会多线程,只要你会编程,网上找一下你所熟悉编程语言的多线程操作例子,不出半个小时,你准能搞定!

4、非线性拟合维度降维

可降维的尽量降低拟合变量的维数。比如拟合表达式为$$y=a+bsin(cx)+dcos(cx)$$,需要拟合的变量为\(a,b,c,d\)。乍一看需要拟合的变量为4个,其实如果\(c\)变量已知,上面的表达式就是一个求解线性方程组问题,也就是说,一旦\(c\)确定了,能很快确定出\(a,b,d\)。这时相当于就只对\(c\)变量进行查找最优值即可,也就从4个变量变成了一个变量的拟合。这对求解效率会大大的提高!

5、非线性拟合时可能遇到的问题
(1)拟合公式里含有积分或者微分等项
解:不管是积分或者是微分,它只是一个函数,自己实现即可
比如:需要拟合的公式为
$$y=a\times sin(\frac{b}{x_1})+c\int_{0}^{x_2}cos(ax^2+bx+c)dx$$
其中\(y,x_1,x_2\)为采集的数据,\(a,b,c\)为需要拟合的变量,这时我们可考虑做如下变换
$$f(a,b,c,y)=c\int_{0}^{y}cos(ax^2+bx+c)dx$$
$$y=a\times sin(\frac{b}{x_1}) + f(a,b,c,x_2)$$
当中的\(f(a,b,c,y)\)函数,我们可以自己写一个积分算法进行求解,这个和公式中直接调用正弦函数\(sin(x)\)是一样的,只是程序内部已经帮我们将正弦函数的算法实现。
同样的,如果有微分项,一样可以使用类似的方法进行处理。
PS:即使有一些软件提供积分或者微分等算法,如果那种算法你不熟悉或者怕被做手脚,加之一些算法的局限性,建议自己编程实现这类积分或者微分求解的问题。

(2)拟合的公式为隐函数的拟合
解:对于这个问题,说白了就是转移一种评价方式的问题,这时进行简单的变换即可
比如:需要拟合的方程为
$$y=a\times sin(x)+\frac{b}{cos(y)}+c\sqrt{xy}$$
其中\(x,y\)为采集的数据,\(a,b,c\)为需要拟合的变量,
正常情况下,程序内部的评价方式可能是\(y-{y}'\),其中一个是真实值一个是拟合值,因为这种评价方式比较困难,这时我们可考虑做如下变换
$$f=a\times sin(x)+\frac{b}{cos(y)}+c\sqrt{xy}-y$$
这时的评价方式即变成\(f-0\),这样就完美解决!

6、拟合的评价方式问题(损失函数)
这里的评价方式只针对求解过程中的评价方式,也因此只讨论非线性问题的评价问题。
所谓的评价方式,就是指用什么方式来解答:什么时候拟合结束,怎样的解才算更优?好的评价函数可以有一定降噪作用,这里不细说,仅列举下面的回答!
这个没有统一的答案,仁者见仁智者见智吧!不过,我还是说下我的想法,有些在求解时喜欢采用各种假设检验来检验解的优度,因为我是一个追求程序效率的人,我并不建议在程序中采用这种方式,因为在调用各种假设检验算法时必定会消耗一部分计算时间,因此一个简单高效的评价方法不失为一种好的方式。
常见的,我们会使用如下形式,即
$$S=\sum_{i=1}^{n}(y_i-{y_i}')^2$$

7、拟合初值的选取问题
因为初值对拟合速度甚至结果有一定影响,因此这里就简单说一下确定数学模型后,拟合非线性问题时,初值的选取的问题。
1、如果已知数学模型,有一定物理意义,则建议根据物理意义选取。
2、当无法确定初值时,且你的数学模型有导数(如果求导模型很复杂甚至没有导数,则可进行简单的差分构造),则可以采用如下的办法进行

步骤: 
(1)求出拟合函数的一阶导数【如果有必要可求更高阶导数】 
(2)使用已知数据求出近似点的一阶导数 
(3)代入一阶导数函数以及原函数求得初值近似值 例子: 
已知一组数据\(x,y\)满足如下关系式,求拟合数据\(a,b,c,d\)的初始近似值
$$y=a+b(x-c)^d$$
步骤: 
(1)$${y}'=bd(x-c)^{(d-1)}$$
(2)因为已知\(x,y\)数据,则根据相应算法求得一组\(x,{y}'\)的近似值,这里记\(f={y}'\)
(3)将\(x,{y}'\) 代入(1)式的方程得到如下三个方程进而求解出\(b,c,d\)
$$
\begin{align}
f_1&=bd(x_1-c)^{d-1} \\
f_2&=bd(x_2-c)^{d-1} \\
f_3&=bd(x_3-c)^{d-1} 
\end{align}
$$
(4)取任意一组\(x,y\)然后将\(b,c,d\)一起代入原方程\(y = a + b (x – c) ^d\) 进而可以求得近似值 \(a\) 
(5)至此\( a,b,c,d\)初始近似值确定完毕!

8、常微分方程拟合

之所以在这里单独列一个点,主要是提供一个求解这类问题的一个快速处理方式:
(1)对于微分项,使用对应差分公式代替,然后按正常流程拟合出对应拟合参数。(由于这个时候的微分使用的差分代替,是比较粗糙的求解的,这个拟合值仅能作为初始值参考)
(2)然后采用使用第(1)点得到的拟合参数作为初始值,常微分方程函数主体采用对应常微分方程求解算法(比如龙格库塔法),然后正常进行非线性拟合。
PS:这种求解方式是先简单粗暴找到一个邻近的解,然后再使用比较精确的算法在这个解附近查找最优值.

MathSword 小型数值计算软件

简介

这是一款专门针对数据处理而开发的小型数值计算软件,其本身小巧、绿色、免费,且易于使用。可解决常规数值计算的一些计算问题。

其本身以命令行输入调用函数方式进行数据处理(目前支持大约1760个函数),对于一些特殊功能开发了对应处理窗口,以方便调用。

PS:

1. 您如果有什么好的建议,欢迎给我发邮件 shikang999@126.com

2. 程序运行中出现找不到dll模块的提示,可自行下载 Microsoft Visual c++ 2022(vc2022运行库) x86/x64库进行安装。

下载

MathSword.x64

MathSword.x86

教程

1. 基础操作

2. 线性/非线性规划求解

3. 非线性拟合

4. 非线性方程组求解

5. 数学公式搜索

6. 神经网络拟合

7. 曲线数据提取

8. 启发式求解器

9. 分段多项式拟合

10. 数据降噪与滤波器设计

功能

(1)常规的矩阵运算

(2)线性/非线性 方程组求解

(3)线性/非线性 规划

(4)线性稀疏矩阵求解

(5)常微分方程组求解

(6)一维、二维、多维数据插值

(7)一维、二维、三维定积分

(8)一维、二维、三维作图(曲线、云图)

(9)非线性数据拟合、公式查找、网络建模计算

(10)数据变换(傅里叶、Haar小波…)

(11)数据平滑、微分

(12)多项式计算

(13)大数、复数计算

(14)一些简单的符号计算

(15) 大部分特殊函数计算(贝塞尔函数、椭圆积分函数、雅克比函数、超几何函数、误差函数、Airy函数、Zeta函数等等)

(16)曲线扣图、曲线数据快速处理

(17)一些其它数据处理的相关计算