拟合问题

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​=argxmin​x⊤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 其中 $\mathbf{\Lambda }$ 为对角特征值矩阵， 不妨认为它们按照从大到小的顺序排列， 记作 λ 1 , … , λ n lambda_{1}, ldots, lambda_{n} λ1​,…,λn​ 。$\mathbf{V}$ 为正交 矩阵， 其列向量为每一维特征值对应的特征向量， 记作 v 1 ， … v n boldsymbol{v}_{1}， ldots boldsymbol{v}_{n} v1​，…vn​ ， 它们构成了一组单位正交基。而 任意 $x$ 总是可以被这组单位正交基线性表出的:

x = α 1 v 1 + … + α n v n boldsymbol{x}=alpha_{1} boldsymbol{v}_{1}+ldots+alpha_{n} boldsymbol{v}_{n} x=α1​v1​+…+αn​vn​
那么不难看出:

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= ​v1​v2​⋅⋅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= ​v1​xv2​x⋅⋅vn​x​ ​= ​v1​(α1​v1​+…+αn​vn​)v2​(α1​v1​+…+αn​vn​)⋅⋅vn​x​ ​= ​a1​a2​⋅⋅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​
而 $\Vert \mathbf{x}\Vert =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$ 的具体形式，也就是得到了 V ⊤ x V^{top}x V⊤x 的具体形式为前面所有维度均为0，只有最后一个维度是$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= ​v1​v2​⋅⋅vn​​ ​x= ​00⋅⋅1​ ​
则 $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= ​v1​v2​⋅⋅vn​​ ​vn​= ​v1​vn​v2​vn​⋅⋅vn​vn​​ ​= ​00⋅⋅1​ ​

quad 至此我们看到， 线性最小二乘的最优解即为最小特征值向量。而由于该问题想要解的是 $Ax=0$ 问题， 所以也可以称为零空间解。注意这里的特征值分解是对 A ⊤ A boldsymbol{A}^{top} boldsymbol{A} A⊤A 做的， 而不是直接对 $\mathbf{A}$ 来 做 ( $\mathbf{A}$ 也不能保证一定能对角化， 而 A ⊤ A boldsymbol{A}^{top} boldsymbol{A} A⊤A 是实对称矩阵， 可以对角化)。

1.2 基于奇异值分解的解法

quad 奇异值解法 上述问题也可以换一种角度来看， 使用奇异值分解 (Singular Value Decomposition， SVD) 来处理。由于任意矩阵都可以进行奇异值分解， 于是我们对 $\mathbf{A}$ 进行 SVD， 可得:
A = U Σ V ⊤ boldsymbol{A}=boldsymbol{U} boldsymbol{Sigma} boldsymbol{V}^{top} A=UΣV⊤
quad 其中 $\mathbf{U}，\mathbf{V}$ 为正交矩阵， $\mathbf{\Sigma }$ 为对角阵， 称为奇异值矩阵， 其对角线元系为 $\mathbf{A}$ 的奇异值， 不妨它们 是由大到小排列的。我们可以把 SVD 的结果代回到线性最小二乘中， 由于 $\mathbf{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 于是， 类似于特征值的解法， 我们取 $\mathbf{x}$ 为 $\mathbf{V}$ 的最后一列即可。总而言之， 如果我们对一组点进行平面拟合， 只需要把所有点的坐标排成 矩阵 $\mathbf{A}$ ， 然后求 $\mathbf{A}$ 的最小奇异值对应的右奇异值向量， 或者求 A ⊤ A boldsymbol{A}^{top} boldsymbol{A} A⊤A 的最小特征值对应的特征向量即可。

2.0 平面拟合

quad 我们先看平面的拟合问题。给定一组由 $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 。然后我们寻找一组平面参数 $\mathbf{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\in \mathbb{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,dmin​k=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~min​k=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 这个问题即线性代数中的解方程问题: 给定任意一个矩阵 $\mathbf{A}$ (注意不是方阵)， 我们希望找一 个非零向量 $\mathbf{x}$ ， 使得 $\mathbf{Ax}$ 能够最小化。当然， 如果 $\mathbf{x}$ 取为零， 该乘积自然是零， 但是我们不想找这 种平凡的解， 所以要给 $\mathbf{x}$ 加上约束 $\mathbf{x}\ne 0$ 。同时， 如果 $\mathbf{x}$ 乘上非零常数 $k$， 那么 $\mathbf{A}\mathbf{x}$ 也会被放大 $k$ 倍， 平方之后就是 k 2 k^{2} k2 倍。所以我们不想讨论 $\mathbf{x}$ 的长度， 只想知道它的方向， 于是又设定 $\Vert \mathbf{x}\Vert =1$ 。
quad 对于$\mathbf{A}$， 我们就不施加什么约束了。在点云平面提取问题中， 上面的 $\mathbf{A}$ 为一个 R n × 4 mathbb{R}^{n times 4} Rn×4 的矩阵， 而 $\mathbf{x}$ 则为 R 4 mathbb{R}^{4} R4 中的单位向量。我们可以问， 取什么样的 $\mathbf{x}$ 时， $\mathbf{Ax}$ 能达到最大值或者最小值。

quad 这里的 $\mathbf{A}$ 其实就是前面方程中的 X ~ tilde{boldsymbol{X}} X~，这里的 $\mathbf{x}$ 其实就是前面的 n ~ tilde{boldsymbol{n}} n~。

3.0 直线拟合

quad 直线拟合。我们仍然设点集为 $X$ 由 $n$ 个三维点组成。不过,我们可以用若干种不同的方式来描述一根直线,例如把直线视为两个平面的交线,或者使用直线上一点再加上直线方向矢量来描述直线。后者更直观一些。
quad 设直线上的点 $x$ 满足方程:
$$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。$\mathbf{d}$ 为直线的方向, 满足 $\Vert \mathbf{d}\Vert =1$ ， $\mathbf{p}$ 为直线$l$ 上某个点， $t$ 为直线参数。

我们想求的是 $\mathbf{d}$ 和 $p$ , 共 6 个末知数。显然, 当给定的点集较大时, 这依然是一个超定方程, 需要 构造最小二乘问题进行求解。
对于任意一个不在 $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$，且该点不是刚刚直线外那一点的垂点 $O$，那么此时由直线上的点，垂点，直线外的点 x k x_k xk​ 三个点构成了一个直角三角形，求 ∣ ∣ x k O ∣ ∣ 2 ||x_kO||_2 ∣∣xk​O∣∣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，即可得到点到直线的距离。

然后构造最小二乘问题, 求解 $\mathbf{d}$ 和 $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,pmin​k=1∑n​fk2​, s.t. ∥d∥=1.
由于每个点的误差项已经取了平方形式, 在此只需求和即可。
接下来我们分离 $\mathbf{d}$ 部分和 $\mathbf{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)⊤d​​d,=(−2)(I−dd⊤)(xk​−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=1n​fk2​​​=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=n1​k=1∑n​xk​,
说明 $p$ 应该取点云的中心。于是, 我们可以先确定 $p$ , 然后再考虑 d ∘ d_{circ} d∘​ 此时 $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⊤yk​yk⊤​d
不难看出第一个误差项并不含 $\mathbf{d}$ , 如何取 $\mathbf{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∗=argdmax​k=1∑n​d⊤yk​yk⊤​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 这个问题与平面拟合很相似, 无非是把最小化变成了最大化。对于平面拟合, 我们求这个问题 的最小化; 对于直线拟合, 则求其最大值。按照前文的讨论, 取 $\mathbf{d}$ 为最小特征值或者奇异值向量, 就得到该问题的最小化解；反之, 求取最大特征值向量时，就得到最大化解。于是该问题的解应该 取 $\mathbf{d}$ 为 $\mathbf{A}$ 的最大右奇异向量, 或者 A ⊤ A boldsymbol{A}^{top} boldsymbol{A} A⊤A 的最大特征值对应的特征向量。