拟合问题(直线平面拟合)
拟合问题
1.0 线性最小二乘的几种解法
1.1 基于特征值的解法
quad
特征值解法 代数意义上的线性最小二乘是指给定矩阵
A
∈
R
m
×
n
boldsymbol{A} in mathbb{R}^{m times n}
A∈Rm×n，计算
x
∗
∈
R
n
boldsymbol{x}^{*} in mathbb{R}^{n}
x∗∈Rn，使得:
x
∗
=
arg
min
x
∥
A
x
∥
2
2
=
arg
min
x
x
⊤
A
⊤
A
x
,
s.t.,
∥
x
∥
=
1.
boldsymbol{x}^{*}=arg min _{boldsymbol{x}}|boldsymbol{A} boldsymbol{x}|_{2}^{2}=arg min _{boldsymbol{x}} boldsymbol{x}^{top} boldsymbol{A}^{top} boldsymbol{A} boldsymbol{x}, quad text { s.t., }|boldsymbol{x}|=1 .
x∗=argxmin∥Ax∥22=argxminx⊤A⊤Ax, s.t., ∥x∥=1.
quad
可以看到
A
⊤
A
boldsymbol{A}^{top} boldsymbol{A}
A⊤A 为一个实对称矩阵， 而根据线性代数的矩阵， 实对称矩阵总是可以利用特征 值分解进行对角化的:
A
⊤
A
=
V
Λ
V
−
1
boldsymbol{A}^{top} boldsymbol{A}=boldsymbol{V} boldsymbol{Lambda} boldsymbol{V}^{-1}
A⊤A=VΛV−1
quad
其中
Λ
boldsymbol{Lambda}
Λ 为对角特征值矩阵， 不妨认为它们按照从大到小的顺序排列， 记作
λ
1
,
…
,
λ
n
lambda_{1}, ldots, lambda_{n}
λ1,…,λn 。
V
boldsymbol{V}
V 为正交 矩阵， 其列向量为每一维特征值对应的特征向量， 记作
v
1
，
…
v
n
boldsymbol{v}_{1}， ldots boldsymbol{v}_{n}
v1，…vn ， 它们构成了一组单位正交基。而 任意
x
x
x 总是可以被这组单位正交基线性表出的:
x
=
α
1
v
1
+
…
+
α
n
v
n
boldsymbol{x}=alpha_{1} boldsymbol{v}_{1}+ldots+alpha_{n} boldsymbol{v}_{n}
x=α1v1+…+αnvn
那么不难看出:
V
−
1
x
=
V
⊤
x
=
[
α
1
，
…
，
α
n
]
⊤
boldsymbol{V}^{-1} boldsymbol{x}=boldsymbol{V}^{top} boldsymbol{x}=left[alpha_{1}， ldots， alpha_{n}right]^{top}
V−1x=V⊤x=[α1，…，αn]⊤
这里做一个简单的解释：其中
V
−
1
boldsymbol{V}^{-1}
V−1可以表示为
V
−
1
=
[
v
1
v
2
⋅
⋅
v
n
]
boldsymbol{V}^{-1}= begin{bmatrix} v_1 \ v_2 \ cdot \ cdot \ v_n end{bmatrix}
V−1=
v1v2⋅⋅vn
则V
−
1
x
boldsymbol{V}^{-1}x
V−1x 为
V
−
1
x
=
[
v
1
x
v
2
x
⋅
⋅
v
n
x
]
=
[
v
1
(
α
1
v
1
+
…
+
α
n
v
n
)
v
2
(
α
1
v
1
+
…
+
α
n
v
n
)
⋅
⋅
v
n
x
]
=
[
a
1
a
2
⋅
⋅
a
n
]
boldsymbol{V}^{-1}x= begin{bmatrix} v_1 x \ v_2x \ cdot \ cdot \ v_nx end{bmatrix}= begin{bmatrix} v_1 (alpha_{1} boldsymbol{v}_{1}+ldots+alpha_{n} boldsymbol{v}_{n}) \ v_2(alpha_{1} boldsymbol{v}_{1}+ldots+alpha_{n} boldsymbol{v}_{n}) \ cdot \ cdot \ v_nx end{bmatrix} = begin{bmatrix} a_1 \ a_2 \ cdot \ cdot \ a_n end{bmatrix}
V−1x=
v1xv2x⋅⋅vnx
=
v1(α1v1+…+αnvn)v2(α1v1+…+αnvn)⋅⋅vnx
=
a1a2⋅⋅an
于是目标函数变为:
∥
A
x
∥
2
2
=
∑
k
=
1
n
λ
k
α
k
2
|boldsymbol{A x}|_{2}^{2}=sum_{k=1}^{n} lambda_{k} alpha_{k}^{2}
∥Ax∥22=k=1∑nλkαk2
而
∥
x
∥
=
1
|boldsymbol{x}|=1
∥x∥=1 意味着
α
1
2
+
…
+
α
k
2
=
1
alpha_{1}^{2}+ldots+alpha_{k}^{2}=1
α12+…+αk2=1 ， 而特征值部分的
λ
k
lambda_{k}
λk 又是降序排列的， 所以就取
α
1
=
0
，
…
，
α
n
−
1
=
0
，
α
n
=
1
alpha_{1}= 0， ldots， alpha_{n-1}=0， alpha_{n}=1
α1=0，…，αn−1=0，αn=1 ， 即：
x
∗
=
v
n
boldsymbol{x}^{*}=boldsymbol{v}_{n}
x∗=vn
这里我们推断出了
a
a
a 的具体形式，也就是得到了
V
⊤
x
V^{top}x
V⊤x 的具体形式为前面所有维度均为0，只有最后一个维度是
1
1
1，在已有的限制条件下如何能够得到当前的结果。我们有:
V
⊤
x
=
[
v
1
v
2
⋅
⋅
v
n
]
x
=
[
0
0
⋅
⋅
1
]
boldsymbol{V}^{top}x= begin{bmatrix} v_1 \ v_2 \ cdot \ cdot \ v_n end{bmatrix}x= begin{bmatrix} 0 \ 0 \ cdot \ cdot \1 end{bmatrix}
V⊤x=
v1v2⋅⋅vn
x=
00⋅⋅1
则x
x
x 可取
x
=
v
n
x=v_n
x=vn：
V
⊤
x
=
[
v
1
v
2
⋅
⋅
v
n
]
v
n
=
[
v
1
v
n
v
2
v
n
⋅
⋅
v
n
v
n
]
=
[
0
0
⋅
⋅
1
]
boldsymbol{V}^{top}x= begin{bmatrix} v_1 \ v_2 \ cdot \ cdot \ v_n end{bmatrix}v_n= begin{bmatrix} v_1v_n \ v_2v_n \ cdot \ cdot \ v_nv_n end{bmatrix}= begin{bmatrix} 0 \ 0 \ cdot \ cdot \1 end{bmatrix}
V⊤x=
v1v2⋅⋅vn
vn=
v1vnv2vn⋅⋅vnvn
=
00⋅⋅1
quad
至此我们看到， 线性最小二乘的最优解即为最小特征值向量。而由于该问题想要解的是
A
x
=
0
A x= 0
Ax=0 问题， 所以也可以称为零空间解。注意这里的特征值分解是对
A
⊤
A
boldsymbol{A}^{top} boldsymbol{A}
A⊤A 做的， 而不是直接对
A
boldsymbol{A}
A 来 做 (
A
boldsymbol{A}
A 也不能保证一定能对角化， 而
A
⊤
A
boldsymbol{A}^{top} boldsymbol{A}
A⊤A 是实对称矩阵， 可以对角化)。
1.2 基于奇异值分解的解法
quad
奇异值解法 上述问题也可以换一种角度来看， 使用奇异值分解 (Singular Value Decomposition， SVD) 来处理。由于任意矩阵都可以进行奇异值分解， 于是我们对
A
boldsymbol{A}
A 进行 SVD， 可得:
A
=
U
Σ
V
⊤
boldsymbol{A}=boldsymbol{U} boldsymbol{Sigma} boldsymbol{V}^{top}
A=UΣV⊤
quad
其中
U
，
V
boldsymbol{U}， boldsymbol{V}
U，V 为正交矩阵，
Σ
boldsymbol{Sigma}
Σ 为对角阵， 称为奇异值矩阵， 其对角线元系为
A
boldsymbol{A}
A 的奇异值， 不妨它们 是由大到小排列的。我们可以把 SVD 的结果代回到线性最小二乘中， 由于
U
boldsymbol{U}
U 为正交矩阵， 在计算 二范数时会被消掉。我们会发现它们和特征值法实际上是一致的：
x
⊤
A
⊤
A
x
=
x
⊤
V
Σ
2
V
⊤
x
.
boldsymbol{x}^{top} boldsymbol{A}^{top} boldsymbol{A x}=boldsymbol{x}^{top} boldsymbol{V} boldsymbol{Sigma}^{2} boldsymbol{V}^{top} boldsymbol{x} .
x⊤A⊤Ax=x⊤VΣ2V⊤x.
quad
于是， 类似于特征值的解法， 我们取
x
boldsymbol{x}
x 为
V
boldsymbol{V}
V 的最后一列即可。总而言之， 如果我们对一组点进行平面拟合， 只需要把所有点的坐标排成 矩阵
A
boldsymbol{A}
A ， 然后求
A
boldsymbol{A}
A 的最小奇异值对应的右奇异值向量， 或者求
A
⊤
A
boldsymbol{A}^{top} boldsymbol{A}
A⊤A 的最小特征值对应的特征向量即可。
2.0 平面拟合
quad
我们先看平面的拟合问题。给定一组由
n
n
n 个点组成的点云
X
=
{
x
1
，
…
，
x
n
}
boldsymbol{X}=left{boldsymbol{x}_{1}， ldots， boldsymbol{x}_{n}right}
X={x1，…，xn} ， 其中每个点的 坐标取三维欧氏坐标
x
k
∈
R
3
boldsymbol{x}_{k} in mathbb{R}^{3}
xk∈R3 。然后我们寻找一组平面参数
n
，
d
boldsymbol{n}， d
n，d ， 使得:
∀
k
∈
[
1
,
n
]
,
n
⊤
x
k
+
d
=
0
forall k in[1, n], boldsymbol{n}^{top} boldsymbol{x}_{k}+d=0
∀k∈[1,n],n⊤xk+d=0
quad
其中
n
∈
R
3
boldsymbol{n} in mathbb{R}^{3}
n∈R3 为法向量，
d
∈
R
d in mathbb{R}
d∈R 为截距。显然，上述问题有四维末知量， 而每个点提供了三个方程。当我们有多个点时，由于噪声影 响，上述方程大概率是无解的 (超定的)。因此， 我们往往会求最小二乘解 (Linear Least Square )。
使其误差最小化:
min
n
,
d
∑
k
=
1
n
∥
n
⊤
x
k
+
d
∥
2
2
.
min _{boldsymbol{n}, d} sum_{k=1}^{n}left|boldsymbol{n}^{top} boldsymbol{x}_{k}+dright|_{2}^{2} .
n,dmink=1∑n
n⊤xk+d
22.
quad
如果取齐次坐标， 还可以再化简一下该问题。一个三维空间点的齐次坐标是四维的， 但实际处理 当中只需在末尾加上 1 即可， 我们记作:
x
~
=
[
x
⊤
,
1
]
⊤
∈
R
4
.
tilde{boldsymbol{x}}=left[boldsymbol{x}^{top}, 1right]^{top} in mathbb{R}^{4} .
x~=[x⊤,1]⊤∈R4.
quad
于是
n
~
=
[
n
⊤
，
d
]
⊤
∈
R
4
tilde{boldsymbol{n}}=left[boldsymbol{n}^{top}， dright]^{top} in mathbb{R}^{4}
n~=[n⊤，d]⊤∈R4 也是一个齐次向量， 上述方程可以写为 :
min
n
~
∑
k
=
1
n
∥
x
~
k
⊤
n
~
∥
2
2
.
min _{tilde{n}} sum_{k=1}^{n}left|tilde{boldsymbol{x}}_{k}^{top} tilde{boldsymbol{n}}right|_{2}^{2} .
n~mink=1∑n
x~k⊤n~
22.
quad
下标 2 表示取二范数， 即欧氏空间的常规范数， 上标二表示取欧氏范数的平方和。上述问题是求 和形式的线性最小二乘， 我们还可以把它写成矩阵形式。将点云的所有点写在一个矩阵中， 记作:
X
~
=
[
x
~
1
,
…
,
x
~
n
]
tilde{boldsymbol{X}}=left[tilde{boldsymbol{x}}_{1}, ldots, tilde{boldsymbol{x}}_{n}right]
X~=[x~1,…,x~n]
quad
那么该问题中的求和号也可以省略掉:
min
n
~
∥
X
~
⊤
n
~
∥
2
2
min _{tilde{boldsymbol{n}}}left|tilde{boldsymbol{X}}^{top} tilde{boldsymbol{n}}right|_{2}^{2}
n~min
X~⊤n~
22
quad
这个问题即线性代数中的解方程问题: 给定任意一个矩阵
A
boldsymbol{A}
A (注意不是方阵)， 我们希望找一 个非零向量
x
boldsymbol{x}
x ， 使得
A
x
boldsymbol{A x}
Ax 能够最小化。当然， 如果
x
boldsymbol{x}
x 取为零， 该乘积自然是零， 但是我们不想找这 种平凡的解， 所以要给
x
boldsymbol{x}
x 加上约束
x
≠
0
boldsymbol{x} neq 0
x=0 。同时， 如果
x
boldsymbol{x}
x 乘上非零常数
k
k
k， 那么
A
x
boldsymbol{A} boldsymbol{x}
Ax 也会被放大
k
k
k 倍， 平方之后就是
k
2
k^{2}
k2 倍。所以我们不想讨论
x
boldsymbol{x}
x 的长度， 只想知道它的方向， 于是又设定
∥
x
∥
=
1
|boldsymbol{x}|=1
∥x∥=1 。
quad
对于
A
boldsymbol{A}
A， 我们就不施加什么约束了。在点云平面提取问题中， 上面的
A
boldsymbol{A}
A 为一个
R
n
×
4
mathbb{R}^{n times 4}
Rn×4 的矩阵， 而
x
boldsymbol{x}
x 则为
R
4
mathbb{R}^{4}
R4 中的单位向量。我们可以问， 取什么样的
x
boldsymbol{x}
x 时，
A
x
boldsymbol{A x}
Ax 能达到最大值或者最小值。
quad
这里的
A
boldsymbol{A}
A 其实就是前面方程中的
X
~
tilde{boldsymbol{X}}
X~，这里的
x
boldsymbol{x}
x 其实就是前面的
n
~
tilde{boldsymbol{n}}
n~。
3.0 直线拟合
quad
直线拟合。我们仍然设点集为
X
X
X 由
n
n
n 个三维点组成。不过,我们可以用若干种不同的方式来描述一根直线,例如把直线视为两个平面的交线,或者使用直线上一点再加上直线方向矢量来描述直线。后者更直观一些。
quad
设直线上的点
x
x
x 满足方程:
x
=
d
t
+
p
x=dt+p
x=dt+p
quad
其中
d
,
p
∈
R
3
,
t
∈
R
boldsymbol{d}, boldsymbol{p} in mathbb{R}^{3}, t in mathbb{R}
d,p∈R3,t∈R。
d
boldsymbol{d}
d 为直线的方向, 满足
∥
d
∥
=
1
|boldsymbol{d}|=1
∥d∥=1 ，
p
boldsymbol{p}
p 为直线
l
l
l 上某个点，
t
t
t 为直线参数。
我们想求的是
d
boldsymbol{d}
d 和
p
p
p , 共 6 个末知数。显然, 当给定的点集较大时, 这依然是一个超定方程, 需要 构造最小二乘问题进行求解。
对于任意一个不在
l
l
l 上的点
x
k
x_{k}
xk , 我们可以利用勾股定理, 计算它离直线垂直距离的平方:
f
k
2
=
∥
x
k
−
p
∥
2
−
∥
(
x
k
−
p
)
⊤
d
∥
2
,
f_{k}^{2}=left|boldsymbol{x}_{k}-boldsymbol{p}right|^{2}-left|left(boldsymbol{x}_{k}-boldsymbol{p}right)^{top} boldsymbol{d}right|^{2},
fk2=∥xk−p∥2−
(xk−p)⊤d
2,
这里做个解释，这里假设直线外有一点
x
k
x_k
xk，做该点到直线的垂线，同时假设直线上有一点
p
p
p，且该点不是刚刚直线外那一点的垂点
O
O
O，那么此时由直线上的点，垂点，直线外的点
x
k
x_k
xk 三个点构成了一个直角三角形，求
∣
∣
x
k
O
∣
∣
2
||x_kO||_2
∣∣xkO∣∣2，那么根据勾股定理知道斜边
∥
x
k
−
p
∥
2
left|boldsymbol{x}_{k}-boldsymbol{p}right|^{2}
∥xk−p∥2，也可算出斜边在直线上的投影
∥
(
x
k
−
p
)
⊤
d
∥
2
left|left(boldsymbol{x}_{k}-boldsymbol{p}right)^{top} boldsymbol{d}right|^{2}
(xk−p)⊤d
2，即可得到点到直线的距离。
然后构造最小二乘问题, 求解
d
boldsymbol{d}
d 和
p
p
p :
(
d
,
p
)
∗
=
arg
min
d
,
p
∑
k
=
1
n
f
k
2
,
s.t.
∥
d
∥
=
1.
(boldsymbol{d}, boldsymbol{p})^{*}=arg min _{boldsymbol{d}, boldsymbol{p}} sum_{k=1}^{n} f_{k}^{2}, quad text { s.t. }|boldsymbol{d}|=1 .
(d,p)∗=argd,pmink=1∑nfk2, s.t. ∥d∥=1.
由于每个点的误差项已经取了平方形式, 在此只需求和即可。
接下来我们分离
d
boldsymbol{d}
d 部分和
p
boldsymbol{p}
p 部分。先考虑
∂
f
k
2
∂
p
frac{partial f_{k}^{2}}{partial boldsymbol{p}}
∂p∂fk2 , 得到:
∂
f
k
2
∂
p
=
−
2
(
x
k
−
p
)
+
2
(
x
k
−
p
)
⊤
d
⏟
标量,
=
d
⊤
(
x
k
−
p
)
d
,
=
(
−
2
)
(
I
−
d
d
⊤
)
(
x
k
−
p
)
.
begin{array}{l} frac{partial f_{k}^{2}}{partial boldsymbol{p}}=-2left(boldsymbol{x}_{k}-boldsymbol{p}right)+2 underbrace{left(boldsymbol{x}_{k}-boldsymbol{p}right)^{top} boldsymbol{d}}_{text {标量, }=boldsymbol{d}^{top}left(boldsymbol{x}_{k}-boldsymbol{p}right)} boldsymbol{d}, \ =(-2)left(boldsymbol{I}-boldsymbol{d} boldsymbol{d}^{top}right)left(boldsymbol{x}_{k}-boldsymbol{p}right) text {. } \ end{array}
∂p∂fk2=−2(xk−p)+2标量, =d⊤(xk−p)
(xk−p)⊤dd,=(−2)(I−dd⊤)(xk−p).
于是整体的目标函数关于
p
p
p 导数为:
∂
∑
k
=
1
n
f
k
2
∂
p
=
∑
k
=
1
n
(
−
2
)
(
I
−
d
⊤
)
(
x
k
−
p
)
,
=
(
−
2
)
(
I
−
d
d
⊤
)
∑
k
=
1
n
(
x
k
−
p
)
.
begin{aligned} frac{partial sum_{k=1}^{n} f_{k}^{2}}{partial boldsymbol{p}} & =sum_{k=1}^{n}(-2)left(boldsymbol{I}-boldsymbol{d}^{top}right)left(boldsymbol{x}_{k}-boldsymbol{p}right), \ & =(-2)left(boldsymbol{I}-boldsymbol{d} boldsymbol{d}^{top}right) sum_{k=1}^{n}left(boldsymbol{x}_{k}-boldsymbol{p}right) . end{aligned}
∂p∂∑k=1nfk2=k=1∑n(−2)(I−d⊤)(xk−p),=(−2)(I−dd⊤)k=1∑n(xk−p).
为了求最小二乘的极值, 令它等于零, 则得到:
p
=
1
n
∑
k
=
1
n
x
k
,
boldsymbol{p}=frac{1}{n} sum_{k=1}^{n} boldsymbol{x}_{k},
p=n1k=1∑nxk,
说明
p
p
p 应该取点云的中心。于是, 我们可以先确定
p
p
p , 然后再考虑
d
∘
d_{circ}
d∘ 此时
p
p
p 已经被求解出来, 不 妨记
y
k
=
x
k
−
p
boldsymbol{y}_{k}=boldsymbol{x}_{k}-boldsymbol{p}
yk=xk−p , 视
y
k
boldsymbol{y}_{k}
yk 为已知量, 对误差项进行简化:
f
k
2
=
y
k
⊤
y
k
−
d
⊤
y
k
y
k
⊤
d
f_{k}^{2}=boldsymbol{y}_{k}^{top} boldsymbol{y}_{k}-boldsymbol{d}^{top} boldsymbol{y}_{k} boldsymbol{y}_{k}^{top} boldsymbol{d}
fk2=yk⊤yk−d⊤ykyk⊤d
不难看出第一个误差项并不含
d
boldsymbol{d}
d , 如何取
d
boldsymbol{d}
d 并不影响它, 可以舍去。而第二项求最小化相当 于去掉负号后, 求最大化:
d
∗
=
arg
max
d
∑
k
=
1
n
d
⊤
y
k
y
k
⊤
d
=
∑
k
=
1
n
∥
y
k
⊤
d
∥
2
2
boldsymbol{d}^{*}=arg max _{boldsymbol{d}} sum_{k=1}^{n} boldsymbol{d}^{top} boldsymbol{y}_{k} boldsymbol{y}_{k}^{top} boldsymbol{d}=sum_{k=1}^{n}left|boldsymbol{y}_{k}^{top} boldsymbol{d}right|_{2}^{2}
d∗=argdmaxk=1∑nd⊤ykyk⊤d=k=1∑n
yk⊤d
22
如果记:
A
=
[
y
1
⊤
⋯
y
n
⊤
]
,
boldsymbol{A}=left[begin{array}{c} boldsymbol{y}_{1}^{top} \ cdots \ boldsymbol{y}_{n}^{top} end{array}right],
A=
y1⊤⋯yn⊤
,
quad
那么该问题变为:
d
∗
=
arg
max
d
∥
A
d
∥
2
2
.
boldsymbol{d}^{*}=arg max _{boldsymbol{d}}|boldsymbol{A} boldsymbol{d}|_{2}^{2} .
d∗=argdmax∥Ad∥22.
quad
这个问题与平面拟合很相似, 无非是把最小化变成了最大化。对于平面拟合, 我们求这个问题 的最小化; 对于直线拟合, 则求其最大值。按照前文的讨论, 取
d
boldsymbol{d}
d 为最小特征值或者奇异值向量, 就得到该问题的最小化解；反之, 求取最大特征值向量时，就得到最大化解。于是该问题的解应该 取
d
boldsymbol{d}
d 为
A
boldsymbol{A}
A 的最大右奇异向量, 或者
A
⊤
A
boldsymbol{A}^{top} boldsymbol{A}
A⊤A 的最大特征值对应的特征向量。