将平面拟合到 3D 中的一组点：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坐标，@ 987654334@ 是a,b,c 系数。

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

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

Answer 1

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

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

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

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

# 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 中的点。所以我们实际上可以使用最小二乘法得到一个更高阶的多项式！

简短说明：实际上，所有不使用精确解析解的求解器通常都停止在实际答案的某个可接受范围内，因此我们得到的情况很少exact 解决方案，但它往往如此接近，以至于我们在实践中接受它作为精确解决方案。此外，即使是最小二乘求解器也很少使用解析解，而是求助于牛顿法等更快的方法。

如果您要更改优化问题，这将是不正确的。我们通常可以找到某些类别的问题的最佳值（其中最大的一类称为凸优化问题 - 尽管有许多非凸问题我们可以在某些条件下找到最佳值）。

如果您有兴趣了解更多信息，请查看 Boyd 和 Vandenberghe 的 Convex Optimization。第一章不需要太多的数学背景，它概述了一般的优化问题以及它与线性和凸规划等可解决的优化问题的关系。