三维平面的表达方式有很多种，通常采用的形式为

$$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可以编程求解。
详细说明可以参考这里。