使用PCA拟合平面(Plane Fitting using PCA)

定义问题

给定一组3D点$\{p_i\}$，我们想找到这组3D点满足的平面参数，即平面的法向量$n$和中心$q$ :

$$\begin{equation} n^\intercal (p_i-q)=0, \forall i \tag{1} \end{equation}$$
当然在实际数据中（1）式一般是不可能严格满足的，因此可以定义如下函数：
$$\begin{equation} \text{dist}(p_i;n,q) \triangleq n^\intercal (p_i-q) \tag{2} \end{equation}$$
上式表示点到面的符号距离，即就是这个距离可能为正也可能为负

注意：$n$表示平面的法向量，并且为单位向量： $n^Tn = 1$

求解问题

现在这个问题可以转化为一个最小二乘问题：

$$\begin{align} \text{cost}(n,q) & \triangleq \sum_i \text{dist}^2(p_i;n,q) \nonumber \\ & = \sum_i (n^\intercal (p_i-q))^2 \nonumber \\ & = n^\intercal [ \cdots, p_i-q, \cdots ] [ \cdots, p_i^\intercal-q^\intercal, \cdots ]^\intercal n \nonumber \\ & = n^\intercal A(q) A(q)^\intercal n \tag{3} \end{align}$$
$$where A(q) \triangleq [\cdots, p_i-q, \cdots]$$

第一步先固定平面法向量求$q$，则（3）对$q$求偏导并且令导数为零：

$$\begin{equation} \mathbf{0} = \frac{\partial \text{cost}(n,q)}{\partial q} \equiv \sum_i (2nn^\intercal q - 2nn^\intercal p_i). \end{equation} \tag{4}$$

求解（4）式可以得到最优的平面中心$q^*$ ：

$$\begin{equation} q^* = \frac{1}{|\{p_i\}|} \sum_i p_i. \end{equation} \tag{5}$$
这个结果和我们直观理解一致，并且这个结果与平面法向量的取值无关。

或者可以从另外一个角度去理解，直观使用（5）式就可以求出一组3D点的中心，（4）式从另外一个角度验证了（5）式的正确性

第二步求解平面的法向量$n$， 我们可以将（5）式带入（3）式中得如下结果：

$$\begin{equation} \text{cost}(n; q^* ) \triangleq n^\intercal A( q^* ) A( q^* )^\intercal n = n^\intercal B( q^* ) n \end{equation} \tag{6}$$

其中$B(q^* ) \triangleq A( q^* ) A( q^* )^\intercal$，由于$A(q)$是一个$3 \times |\{p_i\}|$维的矩阵，所以$B(q)$是一个$3 \times 3$ 维的正定矩阵，现在问题可以表示为如下形式：

$$\begin{align} n^* = \arg & \min_n & & n^\intercal B( q^* ) n \\ & \text{s.t.} & & n^\intercal n=1 \nonumber \end{align} \tag{7}$$

上式的形式和PCA的目标函数非常相似：

$$\begin{align} \max_{a_1} & & a_1^\intercal \sum a_1 \\ \text{s.t.} & & a_1^\intercal a_1=1 \nonumber \end{align} \tag{8}$$

唯一的区别就在于PCA求的是目标函数的最大值，而我们的问题求的是最小值。所以只需要对$B(q^*)$进行分解，最小特征值对应的特征向量就是平面的法向量。当然我们也可以从直观角度去理解，对于一组3D的点平面的法向量一定是最不重要的那个投影向量（理论情况下所有平面上的点在法向量上的投影为零）。

代码

使用PCA拟合平面是我在看冯晨的《Fast Plane Extraction in Organized Point Clouds Using Agglomerative Hierarchical Clustering》论文时看到的，大家有兴趣的可以看看，代码地址，运行结果如下图所示

PCA

下面简单说说主成分分析（PCA），主要是是为了备忘。
记$x_1, ..., x_p$ 为p个原始特征，设新的特征$\xi_i, i = 1, ..., p$是原始特征的线性组合

$$\xi_i = \sum_{j = 1}^p = a_{ij}x_j = \mathbf{a_i^Tx} \tag{9}$$

为了统一$\xi_i$ 的尺度，不妨要求线性组合系数的模为1，即：

$$\mathbf{a_i^T a_i = 1} \tag{10}$$

将（9）式写成矩阵的形式：

$$\mathbf{\xi = A^T x} \tag{11}$$

其中，$\xi$是由新特征$\xi_i$组成的向量， $A$ 是特征变换矩阵。要求解的是最优的正交变换矩阵$A$，它使得新特征$\xi_i$ 的方差达到极致。正交变换保证了新特征之间不相关，而新特征的方差越大，则样本在该维特征上的差异就越大，因而这一特征就越重要。
考虑第一个新特征$\xi_i$：

$$\xi_1 = \sum_{j = 1}^p = a_{1j}x_j = \mathbf{a_1^Tx} \tag{12}$$

它的方差为

$$\begin{align} var(\xi_1) & = E[\xi_1^2] - E^2[\xi_1]\\ & = E[ \mathbf{ a_1^T xx^T a_1 }] - E[ \mathbf{a_1^Tx} ] E[ \mathbf{a_1^Tx} ] \\ & = \mathbf{a_1^T} E[ \mathbf{xx^T} ] \mathbf{a_1} - \mathbf{a_1^T} E[x] E[x^T] \mathbf{a_1} \\ & = \mathbf{a_1^T} \sum \mathbf{a_1} \tag{13} \end{align}$$

其中， $\sum$ 是 $x$ 的协方差矩阵， 可以用样本来估计， $E[.]$ 是数学期望。要在（10）式约束下最大化$\xi_1$ 的方差，这等价于求下面拉格朗日函数的极值：

$$f(a_1) = \mathbf{ a_1^T \sum a_1 - v(a_1^T a_1 - 1) } \tag{14}$$

其中，$v$ 是拉格朗日乘子， 将（14）式对$\mathbf{a_1}$ 求导并令其为零，得到最优解$\mathbf{a_1}$ 满足

$$\sum \mathbf{a_1} = v \mathbf{a_1} \tag{15}$$

这是协方差矩阵$\sum$ 的特征值， 即$\mathbf{a_1}$ 是矩阵$\sum$ 的特征向量， $v$ 是对应的特征值。 把（15）式带入（13）式中，可得：

$$var(\xi_1) = \mathbf{a_1^T \sum a_1} = v \mathbf{a_1^T a_1} = v \tag{16}$$

因此，最后的$\mathbf{a_1}$ 是矩阵$\sum$ 最大特征值对应的特征向量。 $\xi_1$ 称为第一主成分， 它在原始特征的所有线性组合中方差是最大的。

协方差矩阵$\sum$ 共有 $p$ 个特征值$\lambda_i , i = 1, ..., p$ （包括可能相等的特征值和可能为零的特征值），把他们降序排列。按照上面相同的方式，可以得出有对应这些特征向量构造的$p$ 个主成分$\xi_i, i = i,...,p$ ，全部主成分的方差之和为：

$$\sum_{i = 1}^p var(\xi_i) = \sum_{i = 1}^p \lambda_i \tag{17}$$

变换矩阵$A$ 的各个列向量是由$\sum$ 的正交特征向量组成的， 因此$A^T = A^{-1}$ ， 即$A$ 是正交矩阵。 从$\mathbf{\xi}$ 到$\mathbf{x}$ 的逆变换为：

$$\mathbf{x = A \xi} \tag{18}$$

通常把主成分零均值化，即

$$\mathbf{\xi = A^T(x - \mu)} \\ \mathbf{x = A \xi + \mu} \tag{20}$$

这种平移并不影响主成分的方向。

参考资料

[1] Chen Feng, Fast Plane Extraction in Organized Point Clouds Using Agglomerative Hierarchical Clustering
[2] 主成分分析（PCA）原理详解
[3] 深入理解拉格朗日乘子法（Lagrange Multiplier) 和KKT条件
[4] 张学工, 模式识别（第三版）, 主成分分析, 163-164