我学习《视觉SLAM十四讲》中遇到了SVD法求解ICP问题,比较困惑,如何从最小误差中构造SVD分解,如何得到R呐?,在此做出详解。读完这个你就悟了!
1. 问题描述
在ICP算法中,我们有两个点云:
- 源点云 P = { p i } \mathbf{P} = \{\mathbf{p}_i\} P={pi},其中 p i \mathbf{p}_i pi 是第 i i i个点的坐标。
- 目标点云 Q = { q i } \mathbf{Q} = \{\mathbf{q}_i\} Q={qi},其中 q i \mathbf{q}_i qi 是第 i i i个点的坐标。
目标是找到一个旋转矩阵 ( R ) 和一个平移向量 ( T ),使得通过这个变换,源点云 ( \mathbf{P} ) 中的点与目标点云 Q \mathbf{Q} Q中的点对齐,即:
min R , T ∑ i = 1 N ∥ q i − ( R p i + T ) ∥ 2 \min_{R, T} \sum_{i=1}^{N} \| \mathbf{q}_i - (R\mathbf{p}_i + T) \|^2 minR,T∑i=1N∥qi−(Rpi+T)∥2
2. 去除质心
首先,为了简化问题,我们将两个点云的质心对齐。质心去除的步骤如下:
- 计算源点云的质心
p
centroid
\mathbf{p}_\text{centroid}
pcentroid和目标点云的质心
q
centroid
\mathbf{q}_\text{centroid}
qcentroid:
p centroid = 1 N ∑ i = 1 N p i \mathbf{p}_\text{centroid} = \frac{1}{N} \sum_{i=1}^{N} \mathbf{p}_i pcentroid=N1∑i=1Npi
q centroid = 1 N ∑ i = 1 N q i \mathbf{q}_\text{centroid} = \frac{1}{N} \sum_{i=1}^{N} \mathbf{q}_i qcentroid=N1∑i=1Nqi - 将两个点云分别移动到质心为原点的坐标系:
p i ′ = p i − p centroid \mathbf{p}_i' = \mathbf{p}_i - \mathbf{p}_\text{centroid} pi′=pi−pcentroid
q i ′ = q i − q centroid \mathbf{q}_i' = \mathbf{q}_i - \mathbf{q}_\text{centroid} qi′=qi−qcentroid
此时优化问题变为:
min R ∑ i = 1 N ∥ q i ′ − R p i ′ ∥ 2 \min_{R} \sum_{i=1}^{N} \| \mathbf{q}_i' - R\mathbf{p}_i' \|^2 minR∑i=1N∥qi′−Rpi′∥2
3. 误差函数展开
为了最小化误差函数,我们将其展开:
E ( R ) = ∑ i = 1 N ∥ q i ′ − R p i ′ ∥ 2 E(R) = \sum_{i=1}^{N} \| \mathbf{q}_i' - R\mathbf{p}_i' \|^2 E(R)=∑i=1N∥qi′−Rpi′∥2
展开误差函数:
E
(
R
)
=
∑
i
=
1
N
[
(
q
i
′
−
R
p
i
′
)
T
(
q
i
′
−
R
p
i
′
)
]
E(R) = \sum_{i=1}^{N} \left[ (\mathbf{q}_i' - R\mathbf{p}_i')^T(\mathbf{q}_i' - R\mathbf{p}_i') \right]
E(R)=∑i=1N[(qi′−Rpi′)T(qi′−Rpi′)]
进一步展开得到:
E ( R ) = ∑ i = 1 N [ q i ′ T q i ′ − 2 q i ′ T R p i ′ + p i ′ T R T R p i ′ ] E(R) = \sum_{i=1}^{N} \left[ \mathbf{q}_i'^T \mathbf{q}_i' - 2\mathbf{q}_i'^T R \mathbf{p}_i' + \mathbf{p}_i'^T R^T R \mathbf{p}_i' \right] E(R)=∑i=1N[qi′Tqi′−2qi′TRpi′+pi′TRTRpi′]
由于 R T R = I R^T R = I RTR=I(因为 R R R是正交矩阵),所以第三项为 p i ′ T p i ′ v \mathbf{p}_i'^T \mathbf{p}_i'v pi′Tpi′v。因为我们只关心与 ( R ) 相关的项,所以可以忽略常数项,简化为:
E ( R ) = − 2 ∑ i = 1 N q i ′ T R p i ′ E(R) = -2 \sum_{i=1}^{N} \mathbf{q}_i'^T R \mathbf{p}_i' E(R)=−2∑i=1Nqi′TRpi′
4. 构建协方差矩阵 H H H
最小化误差 ( E® ) 等价于最大化:
∑ i = 1 N q i ′ T R p i ′ \sum_{i=1}^{N} \mathbf{q}_i'^T R \mathbf{p}_i' ∑i=1Nqi′TRpi′
我们可以将线性项 ∑ i = 1 N q i ′ T R p i ′ \sum_{i=1}^{N} \mathbf{q}_i'^T R \mathbf{p}_i' ∑i=1Nqi′TRpi′写成矩阵形式:
∑ i = 1 N q i ′ T R p i ′ = trace ( R ∑ i = 1 N p i ′ q i ′ T ) \sum_{i=1}^{N} \mathbf{q}_i'^T R \mathbf{p}_i' = \text{trace} \left( R \sum_{i=1}^{N} \mathbf{p}_i' \mathbf{q}_i'^T \right) ∑i=1Nqi′TRpi′=trace(R∑i=1Npi′qi′T)
这里用到了矩阵迹的性质: trace ( A B ) = trace ( B A ) \text{trace}(\mathbf{A} \mathbf{B}) = \text{trace}(\mathbf{B} \mathbf{A}) trace(AB)=trace(BA),所以可以重新排列为:
trace ( R ∑ i = 1 N p i ′ q i ′ T ) \text{trace} \left( R \sum_{i=1}^{N} \mathbf{p}_i' \mathbf{q}_i'^T \right) trace(R∑i=1Npi′qi′T)
因此,我们可以构建协方差矩阵 H H H:
H = ∑ i = 1 N p i ′ q i ′ T H = \sum_{i=1}^{N} \mathbf{p}_i' \mathbf{q}_i'^T H=∑i=1Npi′qi′T
于是,最大化线性项的过程等价于最大化 trace ( R H ) \text{trace}(R H) trace(RH)。
几何上,协方差矩阵 ( H ) 捕捉了去质心后的源点云和目标点云的点对之间的分布关系。通过最大化 ( \text{trace}(R H) ),我们实际上是在寻找一个旋转矩阵 ( R ),使得源点云在目标点云坐标系下的投影尽可能地接近。
5. 最优化问题中的 SVD 解法
通过 SVD 分解,我们知道:
trace ( R H ) = trace ( R U Σ V T ) \text{trace}(RH) = \text{trace}(RU \Sigma V^T) trace(RH)=trace(RUΣVT)
利用矩阵迹的性质(即 trace ( A B C ) = trace ( C A B ) = trace ( B C A ) \text{trace}(ABC) = \text{trace}(CAB) = \text{trace}(BCA) trace(ABC)=trace(CAB)=trace(BCA)),我们可以将上式改写为:
trace ( R H ) = trace ( R U Σ V T ) = trace ( V T R U Σ ) \text{trace}(RH) = \text{trace}(RU \Sigma V^T) = \text{trace}(V^T RU \Sigma) trace(RH)=trace(RUΣVT)=trace(VTRUΣ)
注意到 trace ( V T R U Σ ) \text{trace}(V^T RU \Sigma) trace(VTRUΣ) 的最大化主要依赖于矩阵 V T R U V^T RU VTRU的结构。因为 V T R U V^T RU VTRU 是一个正交矩阵,而对于给定的对角矩阵 Σ \Sigma Σ,当 V T R U V^T RU VTRU 是一个对角矩阵时,迹达到最大值。这个对角矩阵应该是单位矩阵 I I I(因为这时各项奇异值直接相加),因此我们可以设置:
R
=
V
U
T
R = VU^T
R=VUT
(视觉SLAM十四讲里面给的是
R
=
U
V
T
R = UV^T
R=UVT,我严重怀疑这本书印错了!!!)
6. 结论:用 SVD 解出 R R R
通过上面的推导过程,得出的旋转矩阵 R R R是 V U T VU^T VUT。这就是为什么在得到 H H H之后,通过 SVD 分解 H H H 得到的两个正交矩阵 U U U和 V V V,可以通过计算 R = V U T R = VU^T R=VUT得到旋转矩阵的原因。
5. 几何意义
从几何角度来看,矩阵 H H H捕捉了两个点集之间的协方差关系。SVD 分解将这个关系解构成旋转和尺度变换(通过奇异值)。通过 R = V U T R = VU^T R=VUT的组合,实际上是在寻找最优的旋转矩阵 R R R,使得旋转后的点集与原点集的协方差最大程度匹配。
总结
- 最大化 trace ( R H ) \text{trace}(RH) trace(RH):目标是找到一个旋转矩阵 ( R ),使得目标函数最大化。
- SVD 的应用:通过对 H H H进行 SVD 分解,获得两个正交矩阵 U U U和 V V V,从而构建出最优的旋转矩阵 R = V U T R = VU^T R=VUT。
- 几何意义:SVD 解构了点集的协方差关系,通过旋转矩阵最大化了这些点集之间的对齐程度。
拓展:
为什么
V
T
R
U
V^T RU
VTRU是正交阵?为什么正交阵和对角阵相乘时候,当正交阵是单位阵时候,正交阵和对角阵乘积的迹最大?
在进行矩阵优化和分解时,尤其是在奇异值分解 (SVD) 和旋转矩阵的上下文中,了解 ( V^T R U ) 是正交矩阵的原因非常重要。以下是详细的解释。
证明1
首先我们来证明 ( V^T R U ) 是正交矩阵。
步骤 1: 计算 ( V T R U ) T ( V T R U ) (V^T R U)^T (V^T R U) (VTRU)T(VTRU)
( V T R U ) T ( V T R U ) (V^T R U)^T (V^T R U) (VTRU)T(VTRU)
首先计算 ( V T R U ) T (V^T R U)^T (VTRU)T:
( V T R U ) T = U T R T V (V^T R U)^T = U^T R^T V (VTRU)T=UTRTV
接下来计算乘积:
( V T R U ) T ( V T R U ) = ( U T R T V ) ( V T R U ) (V^T R U)^T (V^T R U) = (U^T R^T V) (V^T R U) (VTRU)T(VTRU)=(UTRTV)(VTRU)
利用矩阵乘法的结合律:
( U T R T V ) ( V T R U ) = U T ( R T ( V V T ) R ) U (U^T R^T V) (V^T R U) = U^T (R^T (V V^T) R) U (UTRTV)(VTRU)=UT(RT(VVT)R)U
因为 V V V 是正交矩阵,所以 V V T = I V V^T = I VVT=I,因此:
U T ( R T I R ) U = U T ( R T R ) U U^T (R^T I R) U = U^T (R^T R) U UT(RTIR)U=UT(RTR)U
由于 R R R 是旋转矩阵,它是正交的,所以 ( R^T R = I ):
U T I U = I U^T I U = I UTIU=I
所以:
( V T R U ) T ( V T R U ) = I (V^T R U)^T (V^T R U) = I (VTRU)T(VTRU)=I
这说明 V T R U V^T R U VTRU 是正交矩阵。
证明2
为了使用公式证明为什么当 A A A) 是正交矩阵, B B B 是对角矩阵时, tr ( A B ) \text{tr}(AB) tr(AB) 的最大值在 A A A 为单位矩阵时达到,我们可以通过迹的性质和正交矩阵的特性来推导。
1. 定义和性质
设 A A A是一个 n × n n \times n n×n 的正交矩阵, B B B 是一个 n × n n \times n n×n 的对角矩阵。我们记 B = diag ( b 1 , b 2 , … , b n ) B = \text{diag}(b_1, b_2, \ldots, b_n) B=diag(b1,b2,…,bn)。
正交矩阵 A A A满足 A T A = I A^T A = I ATA=I,即 ( A ) 的列向量是单位正交的。
2. 迹的定义
矩阵 A B AB AB的迹为:
tr ( A B ) = ∑ i = 1 n ( A B ) i i \text{tr}(AB) = \sum_{i=1}^{n} (AB)_{ii} tr(AB)=∑i=1n(AB)ii
其中 ( A B ) i i (AB)_{ii} (AB)ii表示矩阵 A B AB AB的第 i i i 个对角元素。
3. 矩阵乘积的迹
迹的一个重要性质是:
tr ( A B ) = tr ( B A ) \text{tr}(AB) = \text{tr}(BA) tr(AB)=tr(BA)
因此我们可以等价地考虑 tr ( B A ) \text{tr}(BA) tr(BA),其中 B B B 是对角矩阵:
tr ( B A ) = ∑ i = 1 n ( B A ) i i \text{tr}(BA) = \sum_{i=1}^{n} (BA)_{ii} tr(BA)=∑i=1n(BA)ii
由于 B B B 是对角矩阵,其对角元素为 b 1 , b 2 , … , b n b_1, b_2, \ldots, b_n b1,b2,…,bn,所以:
( B A ) i i = b i ⋅ A i i (BA)_{ii} = b_i \cdot A_{ii} (BA)ii=bi⋅Aii
因此:
tr ( A B ) = tr ( B A ) = ∑ i = 1 n b i ⋅ A i i \text{tr}(AB) = \text{tr}(BA) = \sum_{i=1}^{n} b_i \cdot A_{ii} tr(AB)=tr(BA)=∑i=1nbi⋅Aii
4. 最大化 tr ( A B ) \text{tr}(AB) tr(AB)
要最大化 tr ( A B ) \text{tr}(AB) tr(AB),我们需要最大化上式中的每一项,即 b i ⋅ A i i b_i \cdot A_{ii} bi⋅Aii。
由于 A A A 是正交矩阵,满足 ∑ i = 1 n A i j 2 = 1 \sum_{i=1}^{n} A_{ij}^2 = 1 ∑i=1nAij2=1(即每一列的平方和为 1),因此 A i i A_{ii} Aii 的取值范围为 [ − 1 , 1 ] [-1, 1] [−1,1]。
特别地,当 A A A是单位矩阵时, A i i = 1 A_{ii} = 1 Aii=1 对于所有 i i i 成立,因此此时:
tr ( A B ) = ∑ i = 1 n b i \text{tr}(AB) = \sum_{i=1}^{n} b_i tr(AB)=∑i=1nbi
这是可能得到的最大值,因为如果 A A A不是单位矩阵,那么 A i i A_{ii} Aii会小于 1,从而使得每个 b i ⋅ A i i b_i \cdot A_{ii} bi⋅Aii 的贡献变小,总和 tr ( A B ) \text{tr}(AB) tr(AB) 也会减小。
5. 证明总结
- 迹的计算: tr ( A B ) = ∑ i = 1 n b i ⋅ A i i \text{tr}(AB) = \sum_{i=1}^{n} b_i \cdot A_{ii} tr(AB)=∑i=1nbi⋅Aii。
- 优化目标:我们希望最大化 tr ( A B ) \text{tr}(AB) tr(AB)。
- 单位矩阵的作用:当 A A A是单位矩阵时, A i i = 1 A_{ii} = 1 Aii=1,因此 tr ( A B ) = ∑ i = 1 n b i \text{tr}(AB) = \sum_{i=1}^{n} b_i tr(AB)=∑i=1nbi,达到最大值。
- 非单位矩阵的影响:如果 A A A不是单位矩阵,则 A i i A_{ii} Aii的绝对值小于 1,从而导致 tr ( A B ) \text{tr}(AB) tr(AB)减小。
因此,我们通过公式证明了当 A A A是单位矩阵时, tr ( A B ) \text{tr}(AB) tr(AB) 达到最大值。