将平面拟合到 3D 中的一组点:scipy.optimize.minimize vs scipy.linalg.lstsq

Fit plane to a set of points in 3D: scipy.optimize.minimize vs scipy.linalg.lstsq

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

z = a*x + b*y + c

使得生成的平面是该组点的最佳拟合

  1. this SO answer, the function scipy.optimize.minimize就是用来解决这个问题的。

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

  2. this code (based on this other SO answer) the scipy.linalg.lstsq 中,函数用于解决相同的问题(当限制为一阶多项式时)。

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

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

由于 minimize 函数需要初始猜测,这意味着它可能会或可能不会收敛到最优解(取决于猜测的好坏)。第二种方法是否有类似的警告,或者它会 return 始终是精确的解决方案吗?

最小二乘法 (scipy.linalg.lstsq) 保证收敛。事实上,有一个封闭形式的解析解(由(A^T A)^-1 A^Tb给出(其中^T是矩阵转置,^-1是矩阵求逆)

但是,标准优化问题通常无法解决 - 我们不能保证找到最小值。然而,对于给定的方程,找到一些 a, b, c 使得 z = a*x + b*y + c,我们有一个线性优化问题(约束和 objective 在我们试图优化的变量中是线性的).线性优化问题一般是可解的,所以scipy.optimize.minimize应该会收敛到最优值。

注意:即使我们这样做 z = a*x + b*y + d*x^2 + e*y^2 + f,这在我们的约束中也是线性的——我们不必将自己限制在 (x,y) 的线性 space,因为我们将已经有了这些要点 (x, y, x^2, y^2)。对于我们的算法,这些看起来就像矩阵 A 中的点。所以我们实际上可以使用最小二乘法得到更高阶的多项式!

一个简短的旁白:实际上,所有不使用精确解析解的求解器通常都停留在实际答案的某个可接受范围内,因此很少见我们得到 exact 解的情况,但它往往非常接近以至于我们在实践中接受它是精确的。此外,即使是最小二乘求解器也很少使用解析解,而是求助于更快的方法,例如牛顿法。

如果你要改变优化问题,这就不是真的了。我们通常可以找到某些 class 类问题的最优值(其中最大的 class 称为凸优化问题——尽管有许多非凸问题我们可以找到特定条件下的最优值)。

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

我想用另一种方法来完成答案,以便找到适合 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,但不在同一个矩阵中。

这是一个 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.