定义问题

给定一组3D点 {pi} <script type="math/tex" id="MathJax-Element-1">\{p_i\}</script>，我们想找到这组3D点满足的平面参数，即平面的法向量 n <script type="math/tex" id="MathJax-Element-2">n</script>和中心q<script type="math/tex" id="MathJax-Element-3">q</script> :

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

注意： n <script type="math/tex" id="MathJax-Element-6">n</script>表示平面的法向量，并且为单位向量： nTn=1<script type="math/tex" id="MathJax-Element-7">n^Tn = 1</script>

求解问题

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

<script type="math/tex; mode=display" id="MathJax-Element-8"> \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} </script>
<script type="math/tex; mode=display" id="MathJax-Element-9"> where A(q) \triangleq [\cdots, p_i-q, \cdots] </script>

第一步先固定平面法向量求 q <script type="math/tex" id="MathJax-Element-10">q</script>，则（3）对q<script type="math/tex" id="MathJax-Element-11">q</script>求偏导并且令导数为零：

<script type="math/tex; mode=display" id="MathJax-Element-12"> \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} </script>

求解（4）式可以得到最优的平面中心 q∗ <script type="math/tex" id="MathJax-Element-13">q^*</script> ：

<script type="math/tex; mode=display" id="MathJax-Element-14"> \begin{equation} q^* = \frac{1}{|\{p_i\}|} \sum_i p_i. \end{equation} \tag{5} </script>
这个结果和我们直观理解一致，并且这个结果与平面法向量的取值无关。

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

第二步求解平面的法向量 n <script type="math/tex" id="MathJax-Element-15">n</script>， 我们可以将（5）式带入（3）式中得如下结果：

<script type="math/tex; mode=display" id="MathJax-Element-16"> \begin{equation} \text{cost}(n; q^* ) \triangleq n^\intercal A( q^* ) A( q^* )^\intercal n = n^\intercal B( q^* ) n \end{equation} \tag{6} </script>

其中 B(q∗)≜A(q∗)A(q∗)⊺ <script type="math/tex" id="MathJax-Element-17">B(q^* ) \triangleq A( q^* ) A( q^* )^\intercal</script>，由于 A(q) <script type="math/tex" id="MathJax-Element-18">A(q)</script>是一个 3×|{pi}| <script type="math/tex" id="MathJax-Element-19"> 3 \times |\{p_i\}| </script>维的矩阵，所以 B(q) <script type="math/tex" id="MathJax-Element-20">B(q)</script>是一个 3×3 <script type="math/tex" id="MathJax-Element-21"> 3 \times 3 </script> 维的正定矩阵，现在问题可以表示为如下形式：

<script type="math/tex; mode=display" id="MathJax-Element-22"> \begin{align} n^* = \arg & \min_n & & n^\intercal B( q^* ) n \\ & \text{s.t.} & & n^\intercal n=1 \nonumber \end{align} \tag{7} </script>

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

<script type="math/tex; mode=display" id="MathJax-Element-23"> \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} </script>

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

代码

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

PCA

下面简单说说主成分分析（PCA），主要是是为了备忘。
记 x1,...,xp <script type="math/tex" id="MathJax-Element-25">x_1, ..., x_p</script> 为p个原始特征，设新的特征 ξi,i=1,...,p <script type="math/tex" id="MathJax-Element-26">\xi_i, i = 1, ..., p</script>是原始特征的线性组合

<script type="math/tex; mode=display" id="MathJax-Element-27"> \xi_i = \sum_{j = 1}^p = a_{ij}x_j = \mathbf{a_i^Tx} \tag{9} </script>

为了统一 ξi <script type="math/tex" id="MathJax-Element-28">\xi_i</script> 的尺度，不妨要求线性组合系数的模为1，即：

<script type="math/tex; mode=display" id="MathJax-Element-29"> \mathbf{a_i^T a_i = 1} \tag{10} </script>

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

<script type="math/tex; mode=display" id="MathJax-Element-30"> \mathbf{\xi = A^T x} \tag{11} </script>

其中， ξ <script type="math/tex" id="MathJax-Element-31">\xi</script>是由新特征 ξi <script type="math/tex" id="MathJax-Element-32">\xi_i</script>组成的向量， A <script type="math/tex" id="MathJax-Element-33">A</script> 是特征变换矩阵。要求解的是最优的正交变换矩阵A<script type="math/tex" id="MathJax-Element-34">A</script>，它使得新特征 ξi <script type="math/tex" id="MathJax-Element-35">\xi_i</script> 的方差达到极致。正交变换保证了新特征之间不相关，而新特征的方差越大，则样本在该维特征上的差异就越大，因而这一特征就越重要。
考虑第一个新特征 ξi <script type="math/tex" id="MathJax-Element-36">\xi_i</script>：

<script type="math/tex; mode=display" id="MathJax-Element-37"> \xi_1 = \sum_{j = 1}^p = a_{1j}x_j = \mathbf{a_1^Tx} \tag{12} </script>

它的方差为

<script type="math/tex; mode=display" id="MathJax-Element-38"> \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} </script>

其中， ∑ <script type="math/tex" id="MathJax-Element-39">\sum</script> 是 x <script type="math/tex" id="MathJax-Element-40">x</script> 的协方差矩阵， 可以用样本来估计， E[.]<script type="math/tex" id="MathJax-Element-41">E[.]</script> 是数学期望。要在（10）式约束下最大化 ξ1 <script type="math/tex" id="MathJax-Element-42">\xi_1</script> 的方差，这等价于求下面拉格朗日函数的极值：

<script type="math/tex; mode=display" id="MathJax-Element-43"> f(a_1) = \mathbf{ a_1^T \sum a_1 - v(a_1^T a_1 - 1) } \tag{14} </script>

其中， v <script type="math/tex" id="MathJax-Element-44">v</script> 是拉格朗日乘子， 将（14）式对a1<script type="math/tex" id="MathJax-Element-45"> \mathbf{a_1}</script> 求导并令其为零，得到最优解 a1 <script type="math/tex" id="MathJax-Element-46"> \mathbf{a_1}</script> 满足

<script type="math/tex; mode=display" id="MathJax-Element-47"> \sum \mathbf{a_1} = v \mathbf{a_1} \tag{15} </script>

这是协方差矩阵 ∑ <script type="math/tex" id="MathJax-Element-48">\sum</script> 的特征值， 即 a1 <script type="math/tex" id="MathJax-Element-49"> \mathbf{a_1}</script> 是矩阵 ∑ <script type="math/tex" id="MathJax-Element-50">\sum</script> 的特征向量， v <script type="math/tex" id="MathJax-Element-51">v</script> 是对应的特征值。 把（15）式带入（13）式中，可得：

<script type="math/tex; mode=display" id="MathJax-Element-52"> var(\xi_1) = \mathbf{a_1^T \sum a_1} = v \mathbf{a_1^T a_1} = v \tag{16} </script>

因此，最后的 a1 <script type="math/tex" id="MathJax-Element-53"> \mathbf{a_1}</script> 是矩阵 ∑ <script type="math/tex" id="MathJax-Element-54">\sum</script> 最大特征值对应的特征向量。 ξ1 <script type="math/tex" id="MathJax-Element-55">\xi_1</script> 称为第一主成分， 它在原始特征的所有线性组合中方差是最大的。

协方差矩阵 ∑ <script type="math/tex" id="MathJax-Element-56">\sum</script> 共有 p <script type="math/tex" id="MathJax-Element-57">p</script> 个特征值λi,i=1,...,p<script type="math/tex" id="MathJax-Element-58">\lambda_i , i = 1, ..., p</script> （包括可能相等的特征值和可能为零的特征值），把他们降序排列。按照上面相同的方式，可以得出有对应这些特征向量构造的 p <script type="math/tex" id="MathJax-Element-59">p</script> 个主成分ξi,i=i,...,p<script type="math/tex" id="MathJax-Element-60">\xi_i, i = i,...,p</script> ，全部主成分的方差之和为：

<script type="math/tex; mode=display" id="MathJax-Element-61"> \sum_{i = 1}^p var(\xi_i) = \sum_{i = 1}^p \lambda_i \tag{17} </script>

变换矩阵 A <script type="math/tex" id="MathJax-Element-62">A</script> 的各个列向量是由∑<script type="math/tex" id="MathJax-Element-63">\sum</script> 的正交特征向量组成的， 因此 AT=A−1 <script type="math/tex" id="MathJax-Element-64">A^T = A^{-1}</script> ， 即 A <script type="math/tex" id="MathJax-Element-65">A</script> 是正交矩阵。 从ξ<script type="math/tex" id="MathJax-Element-66">\mathbf{\xi}</script> 到 x <script type="math/tex" id="MathJax-Element-67">\mathbf{x}</script> 的逆变换为：

<script type="math/tex; mode=display" id="MathJax-Element-68"> \mathbf{x = A \xi} \tag{18} </script>

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

<script type="math/tex; mode=display" id="MathJax-Element-69"> \mathbf{\xi = A^T(x - \mu)} \\ \mathbf{x = A \xi + \mu} \tag{20} </script>

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

参考资料

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