拟合问题

拟合问题

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

1.1 基于特征值的解法

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

\[\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 .
\]

$\quad$可以看到 $\boldsymbol{A}^{\top} \boldsymbol{A}$ 为一个实对称矩阵，而根据线性代数的矩阵，实对称矩阵总是可以利用特征值分解进行对角化的:

\[\boldsymbol{A}^{\top} \boldsymbol{A}=\boldsymbol{V} \boldsymbol{\Lambda} \boldsymbol{V}^{-1}
\]

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

\[\boldsymbol{x}=\alpha_{1} \boldsymbol{v}_{1}+\ldots+\alpha_{n} \boldsymbol{v}_{n}
\]

那么不难看出:

\[\boldsymbol{V}^{-1} \boldsymbol{x}=\boldsymbol{V}^{\top} \boldsymbol{x}=\left[\alpha_{1}， \ldots， \alpha_{n}\right]^{\top}
\]

这里做一个简单的解释：其中 $\boldsymbol{V}^{-1}$可以表示为

\[\boldsymbol{V}^{-1}=
\begin{bmatrix}
v_1
\\ v_2
\\ \cdot
\\ \cdot
\\ v_n
\end{bmatrix}
\]

则 $\boldsymbol{V}^{-1}x$ 为

\[\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}
\]

于是目标函数变为:

\[\|\boldsymbol{A x}\|_{2}^{2}=\sum_{k=1}^{n} \lambda_{k} \alpha_{k}^{2}
\]

而 $\|\boldsymbol{x}\|=1$ 意味着 $\alpha_{1}^{2}+\ldots+\alpha_{k}^{2}=1$ ，而特征值部分的 $\lambda_{k}$ 又是降序排列的，所以就取 $\alpha_{1}= 0， \ldots， \alpha_{n-1}=0， \alpha_{n}=1$ ，即：

\[\boldsymbol{x}^{*}=\boldsymbol{v}_{n}
\]

这里我们推断出了 $a$ 的具体形式，也就是得到了 $V^{\top}x$ 的具体形式为前面所有维度均为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}
\]

则 $x$ 可取 $x=v_n$：

\[\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}
\]

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

1.2 基于奇异值分解的解法

$\quad$奇异值解法上述问题也可以换一种角度来看，使用奇异值分解 (Singular Value Decomposition， SVD) 来处理。由于任意矩阵都可以进行奇异值分解，于是我们对 $\boldsymbol{A}$ 进行 SVD，可得:

\[\boldsymbol{A}=\boldsymbol{U} \boldsymbol{\Sigma} \boldsymbol{V}^{\top}
\]

$\quad$其中 $\boldsymbol{U}， \boldsymbol{V}$ 为正交矩阵， $\boldsymbol{\Sigma}$ 为对角阵，称为奇异值矩阵，其对角线元系为 $\boldsymbol{A}$ 的奇异值，不妨它们是由大到小排列的。我们可以把 SVD 的结果代回到线性最小二乘中，由于 $\boldsymbol{U}$ 为正交矩阵，在计算二范数时会被消掉。我们会发现它们和特征值法实际上是一致的：

\[\boldsymbol{x}^{\top} \boldsymbol{A}^{\top} \boldsymbol{A x}=\boldsymbol{x}^{\top} \boldsymbol{V} \boldsymbol{\Sigma}^{2} \boldsymbol{V}^{\top} \boldsymbol{x} .
\]

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

2.0 平面拟合

$\quad$我们先看平面的拟合问题。给定一组由 $n$ 个点组成的点云 $\boldsymbol{X}=\left\{\boldsymbol{x}_{1}， \ldots， \boldsymbol{x}_{n}\right\}$ ，其中每个点的坐标取三维欧氏坐标 $\boldsymbol{x}_{k} \in \mathbb{R}^{3}$ 。然后我们寻找一组平面参数 $\boldsymbol{n}， d$ ，使得:

\[\forall k \in[1, n], \boldsymbol{n}^{\top} \boldsymbol{x}_{k}+d=0
\]

$\quad$其中 $\boldsymbol{n} \in \mathbb{R}^{3}$ 为法向量， $d \in \mathbb{R}$ 为截距。显然，上述问题有四维末知量，而每个点提供了三个方程。当我们有多个点时，由于噪声影响，上述方程大概率是无解的 (超定的)。因此，我们往往会求最小二乘解 (Linear Least Square )。

使其误差最小化:

\[\min _{\boldsymbol{n}, d} \sum_{k=1}^{n}\left\|\boldsymbol{n}^{\top} \boldsymbol{x}_{k}+d\right\|_{2}^{2} .
\]

$\quad$如果取齐次坐标，还可以再化简一下该问题。一个三维空间点的齐次坐标是四维的，但实际处理当中只需在末尾加上 1 即可，我们记作:

\[\tilde{\boldsymbol{x}}=\left[\boldsymbol{x}^{\top}, 1\right]^{\top} \in \mathbb{R}^{4} .
\]

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

\[\min _{\tilde{n}} \sum_{k=1}^{n}\left\|\tilde{\boldsymbol{x}}_{k}^{\top} \tilde{\boldsymbol{n}}\right\|_{2}^{2} .
\]

$\quad$下标 2 表示取二范数，即欧氏空间的常规范数，上标二表示取欧氏范数的平方和。上述问题是求和形式的线性最小二乘，我们还可以把它写成矩阵形式。将点云的所有点写在一个矩阵中，记作:

\[\tilde{\boldsymbol{X}}=\left[\tilde{\boldsymbol{x}}_{1}, \ldots, \tilde{\boldsymbol{x}}_{n}\right]
\]

$\quad$那么该问题中的求和号也可以省略掉:

\[\min _{\tilde{\boldsymbol{n}}}\left\|\tilde{\boldsymbol{X}}^{\top} \tilde{\boldsymbol{n}}\right\|_{2}^{2}
\]

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

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

3.0 直线拟合

$\quad$直线拟合。我们仍然设点集为 $X$ 由 $n$ 个三维点组成。不过,我们可以用若干种不同的方式来描述一根直线,例如把直线视为两个平面的交线,或者使用直线上一点再加上直线方向矢量来描述直线。后者更直观一些。
$\quad$设直线上的点 $x$ 满足方程:

\[x=dt+p
\]

$\quad$其中$\boldsymbol{d}, \boldsymbol{p} \in \mathbb{R}^{3}, t \in \mathbb{R}$。$\boldsymbol{d}$ 为直线的方向, 满足 $\|\boldsymbol{d}\|=1$ ， $\boldsymbol{p}$ 为直线$l$ 上某个点， $t$ 为直线参数。

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

\[f_{k}^{2}=\left\|\boldsymbol{x}_{k}-\boldsymbol{p}\right\|^{2}-\left\|\left(\boldsymbol{x}_{k}-\boldsymbol{p}\right)^{\top} \boldsymbol{d}\right\|^{2},
\]

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

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

\[(\boldsymbol{d}, \boldsymbol{p})^{*}=\arg \min _{\boldsymbol{d}, \boldsymbol{p}} \sum_{k=1}^{n} f_{k}^{2}, \quad \text { s.t. }\|\boldsymbol{d}\|=1 .
\]

由于每个点的误差项已经取了平方形式, 在此只需求和即可。
接下来我们分离 $\boldsymbol{d}$ 部分和 $\boldsymbol{p}$ 部分。先考虑 $\frac{\partial f_{k}^{2}}{\partial \boldsymbol{p}}$ , 得到:

\[\begin{array}{l}
\frac{\partial f_{k}^{2}}{\partial \boldsymbol{p}}=-2\left(\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$ 导数为:

\[\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}
\]

为了求最小二乘的极值, 令它等于零, 则得到:

\[\boldsymbol{p}=\frac{1}{n} \sum_{k=1}^{n} \boldsymbol{x}_{k},
\]

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

\[f_{k}^{2}=\boldsymbol{y}_{k}^{\top} \boldsymbol{y}_{k}-\boldsymbol{d}^{\top} \boldsymbol{y}_{k} \boldsymbol{y}_{k}^{\top} \boldsymbol{d}
\]

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

\[\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}
\]

如果记:

\[\boldsymbol{A}=\left[\begin{array}{c}
\boldsymbol{y}_{1}^{\top} \\
\cdots \\
\boldsymbol{y}_{n}^{\top}
\end{array}\right],
\]

$\quad$那么该问题变为:

\[\boldsymbol{d}^{*}=\arg \max _{\boldsymbol{d}}\|\boldsymbol{A} \boldsymbol{d}\|_{2}^{2} .
\]

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