将平面拟合到三维中的一组点：scipy.optimize.minimize vs scipy.linalg.lstsq

Question

给定3D中的一组点，一般问题是以以下形式找到平面方程的a, b, c系数：

z = a*x + b*y + c

从而得到的平面是该组点的最佳拟合面。

1. 在this SO answer中，函数scipy.optimize.minimize用于解决此问题。

   它依赖于对系数的初始猜测，并最小化一个误差函数，该函数求出每个点到平面表面的距离之和。

2. 在this code（基于this other SO answer）中，scipy.linalg.lstsq函数用于解决相同的问题（当限制为一阶多项式时）。

   它求解方程z = A*C中的C，其中A是点集的x,y坐标的串联，z是点集的z坐标，C是a,b,c系数。

   与上面方法中的代码不同，这个方法似乎不需要对平面系数进行初始猜测。

由于minimize函数需要一个初始猜测，这意味着它可能收敛到最优解，也可能不收敛到最优解（取决于猜测的好坏）。第二种方法是否有类似的警告，或者它是否会返回始终精确的解？

Answer 1

我想用另一种方法来完成答案，以便找到适合R^3中一组点的最佳平面。 实际上，lstsq方法工作得很好，除非在特定情况下（例如）所有点的x坐标为0（或相同）。在这种情况下，lstsq中使用的矩阵的列不是线性独立的。例如：

A = [[ 0   y_0    1]
     [ 0   y_1    1]
     ...
     [ 0   y_k    1] 
     ...
     [ 0   y_N    1]]

为了避免这个问题，可以直接在点集的中心坐标上使用svd。实际上，svd在lstsq中使用，但不在同一个矩阵中使用。

这是一个python示例，给出了coords数组中各点的坐标：

# barycenter of the points
# compute centered coordinates
G = coords.sum(axis=0) / coords.shape[0]

# run SVD
u, s, vh = np.linalg.svd(coords - G)

# unitary normal vector
u_norm = vh[2, :]

使用这种方法，vh矩阵是一个3x3矩阵，其行中包含正交向量。前两个向量在平面上构成正交基，第三个向量是垂直于平面的单位向量。

如果你真的需要a, b, c参数，你可以从法向量得到它们，因为法向量的坐标是(a, b, c)，假设平面的方程是ax + by + cz + d = 0。

Answer 2

最小平方（scipy.linalg.lstsq）保证收敛。实际上，有一个封闭形式的解析解（由(A^T A)^-1 A^Tb给出（其中^T是矩阵转置，^-1是矩阵反转）

然而，标准优化问题一般是不可解的，我们不能保证找到一个最小值。然而，对于给定的方程，找到一些a, b, c，这样z = a*x + b*y + c，我们有一个线性优化问题（约束和目标在我们试图优化的变量中是线性的）。线性优化问题通常是可解的，因此scipy.optimize.minimize应该收敛到最优值。

注意：即使我们做了z = a*x + b*y + d*x^2 + e*y^2 + f，这在我们的约束中也是线性的——我们不必把自己限制在(x,y)的线性空间中，因为我们已经有了这些点(x, y, x^2, y^2)。在我们的算法中，这些看起来就像矩阵A中的点。所以我们可以用最小二乘法得到一个高阶多项式！

暂且不提：实际上，所有不使用精确解析解的解算器通常都停在实际解的某个可接受范围内，因此很少有情况下我们得到了精确的解，但它往往非常接近，以至于我们在实践中接受它为精确解。此外，即使是最小二乘解算器也很少使用解析解，而是使用更快速的方法，如牛顿法。

如果你要改变优化问题，这就不是真的。有一些问题我们通常可以找到一个最优值（其中最大的一类称为凸优化问题——尽管有许多非凸问题我们可以在一定条件下找到一个最优值）。

如果你对学习更多感兴趣，看看博伊德和范登伯格的Convex Optimization。第一章不需要太多的数学背景，概述了一般优化问题，以及它如何与可解优化问题（如线性规划和凸规划）相关。