拟合问题(直线平面拟合)

拟合问题

1.0 线性最小二乘的几种解法

1.1 基于特征值的解法

\quad 特征值解法 代数意义上的线性最小二乘是指给定矩阵 $\mathbf{A}\in {\mathbb{R}}^{m\times n}$，计算 ${\mathbf{x}}^{\ast }\in {\mathbb{R}}^n$，使得:

$${\mathbf{x}}^{\ast }=\arg \underset{\mathbf{x}}{\min }\Vert \mathbf{A}\mathbf{x}{\Vert }_2^2=\arg \underset{\mathbf{x}}{\min }{\mathbf{x}}^{\top }{\mathbf{A}}^{\top }\mathbf{A}\mathbf{x},\text{ s.t., }\Vert \mathbf{x}\Vert =1.$$
\quad 可以看到 ${\mathbf{A}}^{\top }\mathbf{A}$ 为一个实对称矩阵， 而根据线性代数的矩阵， 实对称矩阵总是可以利用特征 值分解进行对角化的:

$${\mathbf{A}}^{\top }\mathbf{A}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{-1}$$
\quad 其中 $\mathbf{\Lambda }$ 为对角特征值矩阵， 不妨认为它们按照从大到小的顺序排列， 记作 ${\lambda }_1,\dots ,{\lambda }_n$ 。$\mathbf{V}$ 为正交 矩阵， 其列向量为每一维特征值对应的特征向量， 记作 ${\mathbf{v}}_1，\dots {\mathbf{v}}_n$ ， 它们构成了一组单位正交基。而 任意 $x$ 总是可以被这组单位正交基线性表出的:

$$\mathbf{x}={\alpha }_1{\mathbf{v}}_1+\dots +{\alpha }_n{\mathbf{v}}_n$$
那么不难看出:

$${\mathbf{V}}^{-1}\mathbf{x}={\mathbf{V}}^{\top }\mathbf{x}={\left[{\alpha }_1，\dots ，{\alpha }_n\right]}^{\top }$$

这里做一个简单的解释：其中 ${\mathbf{V}}^{-1}$可以表示为
$${\mathbf{V}}^{-1}=\left[\begin{array}{c}{v}_1\\ v_2\\ \cdot \\ \cdot \\ v_n\end{array}\right]$$
则 ${\mathbf{V}}^{-1}x$ 为
$${\mathbf{V}}^{-1}x=\left[\begin{array}{c}{v}_{1}x\\ v_{2}x\\ \cdot \\ \cdot \\ v_{n}x\end{array}\right]=\left[\begin{array}{c}{v}_1({\alpha }_1{\mathbf{v}}_1+\dots +{\alpha }_n{\mathbf{v}}_n)\\ v_2({\alpha }_1{\mathbf{v}}_1+\dots +{\alpha }_n{\mathbf{v}}_n)\\ \cdot \\ \cdot \\ v_{n}x\end{array}\right]=\left[\begin{array}{c}{a}_1\\ a_2\\ \cdot \\ \cdot \\ a_n\end{array}\right]$$

于是目标函数变为:
$$\Vert \mathbf{Ax}{\Vert }_2^2=\sum _{k=1}^n{\lambda }_k{\alpha }_k^2$$
而 $\Vert \mathbf{x}\Vert =1$ 意味着 ${\alpha }_1^2+\dots +{\alpha }_k^2=1$ ， 而特征值部分的 ${\lambda }_k$ 又是降序排列的， 所以就取 ${\alpha }_1=0，\dots ，{\alpha }_{n-1}=0，{\alpha }_n=1$ ， 即：
$${\mathbf{x}}^{\ast }={\mathbf{v}}_n$$

这里我们推断出了 $a$ 的具体形式，也就是得到了 $V^{\top }x$ 的具体形式为前面所有维度均为0，只有最后一个维度是$1$，在已有的限制条件下如何能够得到当前的结果。我们有:
$${\mathbf{V}}^{\top }x=\left[\begin{array}{c}{v}_1\\ v_2\\ \cdot \\ \cdot \\ v_n\end{array}\right]x=\left[\begin{array}{c}0\\ 0\\ \cdot \\ \cdot \\ 1\end{array}\right]$$
则 $x$ 可取 $x=v_n$：
$${\mathbf{V}}^{\top }x=\left[\begin{array}{c}{v}_1\\ v_2\\ \cdot \\ \cdot \\ v_n\end{array}\right]v_n=\left[\begin{array}{c}{v}_1{v}_n\\ v_2{v}_n\\ \cdot \\ \cdot \\ v_n{v}_n\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ \cdot \\ \cdot \\ 1\end{array}\right]$$

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

1.2 基于奇异值分解的解法

\quad 奇异值解法 上述问题也可以换一种角度来看， 使用奇异值分解 (Singular Value Decomposition， SVD) 来处理。由于任意矩阵都可以进行奇异值分解， 于是我们对 $\mathbf{A}$ 进行 SVD， 可得:
$$\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^{\top }$$
\quad 其中 $\mathbf{U}，\mathbf{V}$ 为正交矩阵， $\mathbf{\Sigma }$ 为对角阵， 称为奇异值矩阵， 其对角线元系为 $\mathbf{A}$ 的奇异值， 不妨它们 是由大到小排列的。我们可以把 SVD 的结果代回到线性最小二乘中， 由于 $\mathbf{U}$ 为正交矩阵， 在计算 二范数时会被消掉。我们会发现它们和特征值法实际上是一致的：

$${\mathbf{x}}^{\top }{\mathbf{A}}^{\top }\mathbf{Ax}={\mathbf{x}}^{\top }\mathbf{V}{\mathbf{\Sigma }}^2{\mathbf{V}}^{\top }\mathbf{x}.$$
\quad 于是， 类似于特征值的解法， 我们取 $\mathbf{x}$ 为 $\mathbf{V}$ 的最后一列即可。总而言之， 如果我们对一组点进行平面拟合， 只需要把所有点的坐标排成 矩阵 $\mathbf{A}$ ， 然后求 $\mathbf{A}$ 的最小奇异值对应的右奇异值向量， 或者求 ${\mathbf{A}}^{\top }\mathbf{A}$ 的最小特征值对应的特征向量即可。

2.0 平面拟合

\quad 我们先看平面的拟合问题。给定一组由 $n$ 个点组成的点云 X = { x 1 ， … ， x n } \boldsymbol{X}=\left\{\boldsymbol{x}_{1}， \ldots， \boldsymbol{x}_{n}\right\} X={x1​，…，xn​} ， 其中每个点的 坐标取三维欧氏坐标 ${\mathbf{x}}_k\in {\mathbb{R}}^3$ 。然后我们寻找一组平面参数 $\mathbf{n}，d$ ， 使得:
$$\forall k\in [1,n],{\mathbf{n}}^{\top }{\mathbf{x}}_k+d=0$$
\quad 其中 $\mathbf{n}\in {\mathbb{R}}^3$ 为法向量， $d\in \mathbb{R}$ 为截距。显然，上述问题有四维末知量， 而每个点提供了三个方程。当我们有多个点时，由于噪声影 响，上述方程大概率是无解的 (超定的)。因此， 我们往往会求最小二乘解 (Linear Least Square )。

使其误差最小化:

$$\underset{\mathbf{n},d}{\min }\sum _{k=1}^n{\Vert {\mathbf{n}}^{\top }{\mathbf{x}}_k+d\Vert }_2^2.$$
\quad 如果取齐次坐标， 还可以再化简一下该问题。一个三维空间点的齐次坐标是四维的， 但实际处理 当中只需在末尾加上 1 即可， 我们记作:

$$\tilde{\mathbf{x}}={\left[{\mathbf{x}}^{\top },1\right]}^{\top }\in {\mathbb{R}}^4.$$
\quad 于是 $\tilde{\mathbf{n}}={\left[{\mathbf{n}}^{\top }，d\right]}^{\top }\in {\mathbb{R}}^4$ 也是一个齐次向量， 上述方程可以写为 :

$$\underset{\tilde{n}}{\min }\sum _{k=1}^n{\Vert {\tilde{\mathbf{x}}}_k^{\top }\tilde{\mathbf{n}}\Vert }_2^2.$$
\quad 下标 2 表示取二范数， 即欧氏空间的常规范数， 上标二表示取欧氏范数的平方和。上述问题是求 和形式的线性最小二乘， 我们还可以把它写成矩阵形式。将点云的所有点写在一个矩阵中， 记作:

$$\tilde{\mathbf{X}}=\left[{\tilde{\mathbf{x}}}_1,\dots ,{\tilde{\mathbf{x}}}_n\right]$$
\quad 那么该问题中的求和号也可以省略掉:

$$\underset{\tilde{\mathbf{n}}}{\min }{\Vert {\tilde{\mathbf{X}}}^{\top }\tilde{\mathbf{n}}\Vert }_2^2$$
\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$ 倍。所以我们不想讨论 $\mathbf{x}$ 的长度， 只想知道它的方向， 于是又设定 $\Vert \mathbf{x}\Vert =1$ 。
\quad 对于$\mathbf{A}$， 我们就不施加什么约束了。在点云平面提取问题中， 上面的 $\mathbf{A}$ 为一个 ${\mathbb{R}}^{n\times 4}$ 的矩阵， 而 $\mathbf{x}$ 则为 ${\mathbb{R}}^4$ 中的单位向量。我们可以问， 取什么样的 $\mathbf{x}$ 时， $\mathbf{Ax}$ 能达到最大值或者最小值。

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

3.0 直线拟合

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

我们想求的是 $\mathbf{d}$ 和 $p$ , 共 6 个末知数。显然, 当给定的点集较大时, 这依然是一个超定方程, 需要 构造最小二乘问题进行求解。
对于任意一个不在 $l$ 上的点 $x_k$ , 我们可以利用勾股定理, 计算它离直线垂直距离的平方:
$$f_k^2={\Vert {\mathbf{x}}_k-\mathbf{p}\Vert }^2-{\Vert {\left({\mathbf{x}}_k-\mathbf{p}\right)}^{\top }\mathbf{d}\Vert }^2,$$

这里做个解释，这里假设直线外有一点 $x_k$，做该点到直线的垂线，同时假设直线上有一点 $p$，且该点不是刚刚直线外那一点的垂点 $O$，那么此时由直线上的点，垂点，直线外的点 $x_k$ 三个点构成了一个直角三角形，求 $||x_{k}O|{|}_2$，那么根据勾股定理知道斜边 ${\Vert {\mathbf{x}}_k-\mathbf{p}\Vert }^2$，也可算出斜边在直线上的投影 ${\Vert {\left({\mathbf{x}}_k-\mathbf{p}\right)}^{\top }\mathbf{d}\Vert }^2$，即可得到点到直线的距离。

然后构造最小二乘问题, 求解 $\mathbf{d}$ 和 $p$ :

$$(\mathbf{d},\mathbf{p}{)}^{\ast }=\arg \underset{\mathbf{d},\mathbf{p}}{\min }\sum _{k=1}^n{f}_k^2,\text{ s.t. }\Vert \mathbf{d}\Vert =1.$$
由于每个点的误差项已经取了平方形式, 在此只需求和即可。
接下来我们分离 $\mathbf{d}$ 部分和 $\mathbf{p}$ 部分。先考虑 $\frac{\partial f_k^2}{\partial \mathbf{p}}$ , 得到:
$$\begin{array}{l}\frac{\partial f_k^2}{\partial \mathbf{p}}=-2\left({\mathbf{x}}_k-\mathbf{p}\right)+2\underset{\text{标量, }={\mathbf{d}}^{\top }\left({\mathbf{x}}_k-\mathbf{p}\right)}{\underbrace{{\left({\mathbf{x}}_k-\mathbf{p}\right)}^{\top }\mathbf{d}}}\mathbf{d},\\ =(-2)\left(\mathbf{I}-\mathbf{d}{\mathbf{d}}^{\top }\right)\left({\mathbf{x}}_k-\mathbf{p}\right)\text{. }\end{array}$$
于是整体的目标函数关于 $p$ 导数为:

$$\begin{array}{rl}\frac{\partial {\sum }_{k=1}^n{f}_k^2}{\partial \mathbf{p}}& =\sum _{k=1}^n(-2)\left(\mathbf{I}-{\mathbf{d}}^{\top }\right)\left({\mathbf{x}}_k-\mathbf{p}\right),\\ & =(-2)\left(\mathbf{I}-\mathbf{d}{\mathbf{d}}^{\top }\right)\sum _{k=1}^n\left({\mathbf{x}}_k-\mathbf{p}\right).\end{array}$$
为了求最小二乘的极值, 令它等于零, 则得到:

$$\mathbf{p}=\frac{1}{n}\sum _{k=1}^n{\mathbf{x}}_k,$$
说明 $p$ 应该取点云的中心。于是, 我们可以先确定 $p$ , 然后再考虑 $d_{\circ }$ 此时 $p$ 已经被求解出来, 不 妨记 ${\mathbf{y}}_k={\mathbf{x}}_k-\mathbf{p}$ , 视 ${\mathbf{y}}_k$ 为已知量, 对误差项进行简化:

$$f_k^2={\mathbf{y}}_k^{\top }{\mathbf{y}}_k-{\mathbf{d}}^{\top }{\mathbf{y}}_k{\mathbf{y}}_k^{\top }\mathbf{d}$$
不难看出第一个误差项并不含 $\mathbf{d}$ , 如何取 $\mathbf{d}$ 并不影响它, 可以舍去。而第二项求最小化相当 于去掉负号后, 求最大化:

$${\mathbf{d}}^{\ast }=\arg \underset{\mathbf{d}}{\max }\sum _{k=1}^n{\mathbf{d}}^{\top }{\mathbf{y}}_k{\mathbf{y}}_k^{\top }\mathbf{d}=\sum _{k=1}^n{\Vert {\mathbf{y}}_k^{\top }\mathbf{d}\Vert }_2^2$$
如果记:

$$\mathbf{A}=\left[\begin{array}{c}{\mathbf{y}}_1^{\top }\\ \cdots \\ {\mathbf{y}}_n^{\top }\end{array}\right],$$
\quad 那么该问题变为:

$${\mathbf{d}}^{\ast }=\arg \underset{\mathbf{d}}{\max }\Vert \mathbf{A}\mathbf{d}{\Vert }_2^2.$$
\quad 这个问题与平面拟合很相似, 无非是把最小化变成了最大化。对于平面拟合, 我们求这个问题 的最小化; 对于直线拟合, 则求其最大值。按照前文的讨论, 取 $\mathbf{d}$ 为最小特征值或者奇异值向量, 就得到该问题的最小化解；反之, 求取最大特征值向量时，就得到最大化解。于是该问题的解应该 取 $\mathbf{d}$ 为 $\mathbf{A}$ 的最大右奇异向量, 或者 ${\mathbf{A}}^{\top }\mathbf{A}$ 的最大特征值对应的特征向量。