使用Eigen和SVD找到一组点的最佳拟合平面

Question

在一个3D空间中给定N个点集，我正在尝试使用SVD和Eigen寻找最佳拟合平面。

我的算法如下：

1. 将数据点围绕（0,0,0）居中。 2. 形成点坐标的3xN矩阵。 3. 计算矩阵的SVD。 4. 将对应于最小奇异值的最小奇异向量设置为平面的法线。 5. 将平面到原点的距离设置为normal∙centroid。

我无法想出如何使用 Eigen's SVD Module来找到对应于点坐标矩阵的最小奇异值的最小奇异向量。

到目前为止，我有以下代码（算法的步骤1、2和5）：

Eigen::Matrix<float, 3, 1> mean = points.rowwise().mean();
const Eigen::Matrix3Xf points_centered = points.colwise() - mean;

int setting = Eigen::ComputeThinU | Eigen::ComputeThinV;
Eigen::JacobiSVD<Eigen::Matrix3Xf> svd = points_centered.jacobiSvd(setting);

Eigen::Vector3d normal = **???**

double d = normal.dot(mean);

- dustymax

---

你目前尝试了哪些代码？ - Avi Ginsburg

已将代码添加到问题中。 - dustymax

请查看此演示文稿的第97-99页。 - ggael

@ggael，我想使用SVD，我读到在某些情况下应该更精确。 - dustymax

Answer 1

表示 U = svd.matrixU()，向量 U.col(0) 和 U.col(1) 定义了你的平面的基础，而 U.col(2) 则是你的平面的法线。

U.col(0) 还定义了具有最大标准差的方向。

- billx

---

请问能否再详细解释一下为什么要使用标志 ComputeFullU 而不是 ComputeThinU？根据 https://eigen.tuxfamily.org/dox/classEigen_1_1JacobiSVD.html，`ComputeThinU` 将计算前三个奇异向量，因为我们有一个 3xN 的矩阵。而第三个正是我们需要的。所以我有点困惑。 - vicky Lin

你能否再详细解释一下为什么使用ComputeFullU标志而不是ComputeThinU标志？根据https://eigen.tuxfamily.org/dox/classEigen_1_1JacobiSVD.html，`ComputeThinU`将计算前三个奇异向量，因为我们有一个3xN的矩阵。而第三个奇异向量正是我们所需要的。所以我在这里有些困惑。 - undefined

这是一个相当老的答案。现在我重新阅读了一下，我确实不明白为什么需要 ComputeFullU。 - billx

Answer 2

你的问题基本上是如何使用Eigen JacobiSVD模块进行最小二乘拟合。这里有一个更有帮助的示例link。最小二乘拟合的基本思想是，首先将所有N-1个点与N个点之间的向量差计算出来，然后尝试将所有这样的N-1个向量近似为两个基向量的线性组合，这两个基向量定义了2D平面。

- Zhuoran He

---

这个链接是关于“正常”的最小二乘拟合，而问题是关于总体最小二乘法的。 - ggael

我意识到我的算法并不是最优的。把质心放在最佳拟合平面上，而不是N个点之一，更有意义。我刚刚使用Lagrange乘数法从数学上推导出，问题中给出的算法已经最小化了N个点到平面距离的平方和。这是一个与坐标无关的结论。在推导之后，算法中就没有更多的最小二乘要做了。实现该算法的正确代码在billx的答案中给出。 - Zhuoran He

我的示例执行普通的最小二乘拟合，其中通过使用最小二乘公式x=(A'*A)^(-1)*A'*b解决了一个超定线性方程组Ax=b。该实现更为巧妙。首先对 A 进行奇异值分解，即 A=U*S*V'，然后解决方案为 x=V*S^(-1)*U'*b。这避免了求解 A'*A 的逆矩阵，因为它具有与 A 本身相比的平方条件数。 - Zhuoran He

这个问题比我之前想象的更有趣。一个看起来像最小二乘的问题最终并不需要进行任何最小二乘计算。让我来提高一下难度！ - Zhuoran He