首页 文章

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

提问于
浏览
6

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

z = a*x + b*y + c

这样得到的平面最适合那组点 .

它依赖于系数的初始猜测,并最小化误差函数,该函数将每个点与平面表面的距离相加 .

它解决方程 z = A*C 中的 C ,其中 A 是该组点的 x,y 坐标的串联, z 是该组的 z 坐标, Ca,b,c 系数 .

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

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

2 回答

  • 8

    保证最小二乘( 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 中的点 . 所以我们实际上可以使用最小二乘法获得更高阶的多项式!

    A brief aside: 实际上,所有不使用精确解析解的求解器一般都停留在实际答案的某个可接受的范围内,因此我们很少得到确切的解,但它往往非常接近我们在实践中完全接受它 . 此外,即使最小二乘解算器很少使用解析解,而是采用像牛顿方法更快的方法 .

    如果您要更改优化问题,则不会这样 . 存在某些类型的问题,我们通常可以找到最佳值(这些问题中最大的类称为凸优化问题 - 尽管存在许多非凸问题,我们可以在某些条件下找到最佳值) .

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

  • 2

    我想用另一种方法完成答案,以找到适合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 .

相关问题