三维平面的表达方式有很多种，通常采用的形式为$$Ax+By+Cz+D=0$$其中\(\begin{bmatrix}A\\B\\C \end{bmatrix}\)为平面法向量。

已知：采样点集\(\{(x_1,y_1),(x_2,y_2),…,(x_n,y_n)\}\),欲求一平面\(p\),使得所有采样点到\(p\)的距离的平方和\(f\)最小。
$$f=\sum_{i=1}^n (\frac{|Ax_i+By_i+Cz_i+D|}{\sqrt{A^2+B^2+C^2}})^2$$
因为\(A,B,C\)的关联性，我们可以令\(||(A,B,C)^T||=1。\)那么
$$f = \sum(Ax_i+By_i+Cz_i+D)^2$$
这与二维直线\(Ax+By+C=0\)的情况一模一样，具体解法可以参考我的另外一篇帖子。
主要思路是：
首先解D。令\(\frac {\partial f }{\partial D}=\sum2(Ax_i+By_i+Cz_i+D) =0\),可以得到\(D = -\frac{\sum (Ax_i+By_i+Cz_i)}{n} \)。
令\(Y_i = X_i -\bar X\),则 \(f=\sum(Ay_i+By_i+Cz_i)^2\)
$$(f=\begin{bmatrix}A\\B\\C\end{bmatrix}^T\begin{bmatrix}\sum x_2 & \sum xy & \sum xz\\ \sum xy & \sum y^2 &  \sum yz\\ \sum xz & \sum yz & \sum z^2\end{bmatrix}\begin{bmatrix}A\\B\\C\end{bmatrix}$$
因为限定条件\(||(A,B,C)^T||=1。\)以及中间的矩阵为实对称矩阵，上述表达式取最小值时的\(A,B,C\)值对应中间矩阵的最小特征值的最小特征向量。采用Eigen的EigenSolver可以编程求解。
详细说明可以参考这里。