给定一组 3D 点,一般问题是找到平面方程的a, b, c
系数,形式如下:
z = a*x + b*y + c
这样生成的平面就与该点集最合适。
在这个 SO 答案中,函数 scipy.optimize.最小化 用于解决此问题。
它依赖于对系数的初始猜测,并最小化误差函数,该函数将每个点到平面的距离相加。
在此代码中(基于其他 SO 答案),scipy.linalg.lstsq 函数用于解决相同的问题(当限制为一阶多项式时)。
它求解方程
z = A*C
中的C
,其中A
是点集的x,y
坐标的串联,z
是集合的z
坐标,C
是a,b,c
系数。与上述方法中的代码不同,此方法似乎不需要对平面系数进行初始猜测。
由于minimize
函数需要初始猜测,这意味着它可能会也可能不会收敛到最佳解决方案(取决于猜测的好坏)。第二种方法是否有类似的警告,或者它会返回始终精确的解决方案吗?
二乘法(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
中的点。因此,我们实际上可以使用最小二乘法获得高阶多项式!
简短的旁白:实际上,所有不使用精确解析解的求解器通常都停留在实际答案的某个可接受的范围内,因此我们很少得到精确的解,但它往往非常接近,以至于我们在实践中接受它是精确的。此外,即使是最小二乘求解器也很少使用解析解,而是求助于牛顿方法等更快的方法。
如果要更改优化问题,则不会发生这种情况。对于某些类别的问题,我们通常可以找到最优值(其中最大的一类称为凸优化问题 - 尽管有许多非凸问题,我们可以在某些条件下找到最优值)。
如果您有兴趣了解更多信息,请查看Boyd和Vandenberghe的凸优化。第一章不需要太多的数学背景,它概述了一般优化问题以及它与线性和凸规划等可求解优化问题的关系。
用另一种方法完成答案,以便找到适合 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
的。