在Matlab中设置线性规划代码获取4D或更高维度的顶点

Question

假设我有一个4D的点云space，可以设置任意密集。这些点将位于多面体的表面上，此时多面体是未知物体。（只是为了提供一些无用的可视化，这里是我给出的点云的相当任意的投影 - 即绘制 100000x4 数组的前 3 列）：

我的目标是获得这个多面体的顶点，我为获得这些顶点所做的无数尝试之一是计算点云的凸包。这里的问题是生成的顶点数量与指定的容差不同，并且没有先验的方法来检查使用哪一个。我还注意到这里使用 qhull 算法并不是最优的，因为所有给定的点实际上都位于船体上。现在在我的中，我接受的答案建议设置面部损失函数并使用梯度下降。这似乎是一个不错的方法，但我当时无法实现它。我认为这将是一个线性规划问题，不幸的是我对此一无所知。我浏览了一些在线资源，似乎在该设置中经常遇到这样的极值点（凸优化是关于什么的？），这就是为什么我希望在 SO 上提出后续问题是有意义的.

我也希望我把建议的做法贴在这里不要违反礼仪:

Let xi in R4, i = 1, ..., m, be the PCA-reduced data points.

Let F = (a, b) be a face, where a is in R4 with a • a = 1 and b is in R.

We define the face loss function L as follows, where λ+, λ- > 0 are parameters you choose. λ+ should be a very small positive number. λ- should be a very large positive number.

L(F) = sumi(λ+ • max(0, a • xi + b) - λ- • min(0, a • xi + b))

We want to find minimal faces F with respect to the loss function L. The small set of all minimal faces will describe your polytope. You can solve for minimal faces by randomly initializing F and then performing gradient descent using the partial derivatives ∂L / ∂aj, j = 1, 2, 3, 4, and ∂L / ∂b. At each step of gradient descent, constrain a • a to be 1 by normalizing.

∂L / ∂aj = sumi(λ+ • xj • [a • xi + b > 0] - λ- • xj • [a • xi + b < 0]) for j = 1, 2, 3, 4

∂L / ∂b = sumi(λ+ • [a • xi + b > 0] - λ- • [a • xi + b < 0])

Note Iverson brackets: [P] = 1 if P is true and [P] = 0 if P is false.

有人可以提供一些指导，我应该如何将其放入代码中？我想我可以为此使用 Matlab 的 linprog 吗？例如，我不太确定随机启动一张脸是什么意思，如果此时 a 和 b 是未知数。

此外，这里是一个数据集的 link。

Answer 1

以下使用 RANSAC 方法将 2D 平面拟合到 nD 点云上。之后可以通过这些平面的交集找到顶点。

算法思路：

RANSAC 的总体思路非常简单。

Select 个称为 假设内点的随机数据点
从假设的内点生成模型
找到其他与模型很好拟合的点，称之为共识集
重复步骤 1-3，直到找到具有所需幅度的共识集的模型。
基于最大共识集改进模型。

在我们的示例中，我们可以通过以下方式使用它：

Select 三个假设的内点。
这三个点所跨越的平面就是我们的模型。
计算整个数据集的点到平面的正交距离。 Select 那些在预定义公差范围内的点。
重复步骤 1-3 和 select 具有最大共识集的平面。
（我们现在可以找到适合平面的最小二乘法，但在这种情况下似乎没有必要。）
保存找到的平面并从我们需要搜索的点中删除其共识集。
重复上述步骤，直到找到所有平面。
将平面相交以获得顶点。

实现这个所需的代码：

我们需要一些代码来实现这个，这超出了这个答案的范围：

计算我们的数据点与模型平面的距离：distancePointsAffineSpace。
相交平面得到edges/vertices：affineSpaceIntersection.
从最终集中删除重复点：uniquetol。（从版本 R2015a 开始，此功能是内置的。）

实施：

有了这个，RANSAC 部分就非常简单了：

%% // Match planes to dataset X:
%  // Choose 3 Points randomly. Generate plane. Find points within tol.
pointsWithinTolOf = @(Points,tol,Space) ...
                      distancePointsAffineSpace(Points, Space)<tol;
availablePoints = 1:size(X,1);
[foundPlane, pointsOfPlane] = deal(cell(0));
for i = 1:maxNumPlanes
    disp(['Plane #',num2str(i)]);
    bestNumInliers = 0;
    for j = 1:numIterations
        randomPointsIdxs = availablePoints(randperm(numel(availablePoints),3));
        currentPlane = X(randomPointsIdxs,:);
        inliers = find(pointsWithinTolOf(X, PointPlaneDistTol, currentPlane));
        numInliers = numel(inliers);
        if numInliers > bestNumInliers
            bestCurrentPlane = currentPlane;
            bestInliers = inliers;
            bestNumInliers = numInliers;
        end
    end
    foundPlane{i} = bestCurrentPlane;
    pointsOfPlane{i} = bestInliers;
    availablePoints = setdiff(availablePoints, bestInliers);
    if isempty(availablePoints)
        break;
    end
end
numPlanes = numel(foundPlane);

剩下的代码比较长，我就把基本思路写下来，link写到代码中。

找到模型飞机后，我们：

通过与所有模型平面相交并检查相交线上是否有数据点来找到每个平面的边。
通过相交其边来计算每个平面的顶点。

您可以下载此 here 的完整代码。

可视化：

您似乎正在处理由三个挤压三角形和三个挤压四边形构成的 3D 莫比乌斯带的边界。这是投影到 3D（投影到您的 2D 屏幕上）的此条带的 4D 旋转。

在Matlab中设置线性规划代码获取4D或更高维度的顶点

Set up linear programming code in Matlab to obtain vertices in 4D or higher dimensions

matlab

computational-geometry

linear-programming

算法思路：

实现这个所需的代码：

实施：

可视化：