将两条线拟合到一组二维点

Question

给定一组二维点（图中的黑点）我需要选择两条线以某种方式表示这些点。我可能需要最小化 [距离 x 到两条线中较近的距离]^2 的平方和，尽管如果有任何其他指标使这更容易，这也很好。

明显但无效的方法是在所有 2^n 分区上尝试最小二乘法。一种实用的方法可能是迭代改进，也许从随机分区开始。

是否有研究处理这个问题的算法？

Answer 1

我想到了文献中的两个相关概念。

首先，我的直觉是应该有一种方法可以将这个问题解释为估计 mixture model. I haven't worked out the details, since parameter estimation typically uses expectation--maximization 的参数，我可以描述它是如何工作的：将一个分区初始化为两部分，然后交替运行对每个部分进行回归，并根据它们与新回归线的距离重新分配点。

其次，假设输入相对干净，您应该能够使用 RANSAC 获得良好的初始分区。对于一些小的 k，取两个不相交的 k 点随机样本，并通过它们拟合线，然后每隔一个点分配一次。为了获得 (100 − x)% 的成功机会，您需要重复 ln(100/x) × 2^2k−1 次并选择最好的一个。

Answer 2

我认为这可以表述为一个显式优化问题：

 min sum(j, r1(j)^2 + r2(j)^2)                 (quadratic)
 subject to
     r1(j) = (y(j) - a0 - a1*x(j)) * δ(j)      (quadratic)
     r2(j) = (y(j) - b0 - b1*x(j)) * (1-δ(j))  (quadratic) 
     δ(j) ∈ {0,1}

我们同时进行点到线的分配和回归（残差平方和的最小化）。

这是一个非凸二次约束混合整数二次规划问题。可以处理此问题的求解器包括 Gurobi、Baron、Couenne、Antigone。

我们可以重新表述一下：

 min sum(j, r(j)^2)                      (convex quadratic)
 subject to
     r(j) = y(j) - a0 - a1*x(j) + s1(j)  (one of these will be relaxed)
     r(j) = y(j) - b0 - b1*x(j) + s2(j)  (all linear)
     -U*δ(j) <= s1(j) <= U*δ(j)
     -U*(1-δ(j)) <= s2(j) <= U*(1-δ(j))
     δ(j) ∈ {0,1}
     s1(j),s2(j) ∈ [-U,U]
     U = 1000                            (reasonable bound, constant)

这使它成为一个直 凸 MIQP 模型。这将允许使用更多求解器（例如 Cplex）并且更容易求解。其他一些公式是 here。提到的一些模型不需要我在上面的 big-M 公式中必须使用的边界。值得注意的是，这些模型提供了经过验证的最佳解决方案（对于非凸模型，这需要全局求解器；凸模型更容易，不需要这个）。此外，除了最小二乘 objective，我们还可以形成 L1 objective。在后一种情况下，我们最终得到一个线性 MIP 模型。

一个小测试证实了这一点：

这道题有 50 分，在慢速笔记本电脑上使用 Cplex 的 MIQP 求解器需要 1.25 秒。这可能是一个统计问题的案例，其中 MIP/MIQP 方法可以提供一些东西。

Answer 3

在 OPL CPLEX 中，如果我从示例开始 curve fitting from Model Buidling

先给.dat补充几点

n=24;
 
x = [1,2,3,4,5,0.0, 0.5, 1.0, 1.5, 1.9, 2.5, 3.0, 3.5, 4.0, 4.5,
     5.0, 5.5, 6.0, 6.6, 7.0, 7.6, 8.5, 9.0, 10.0];
y = [10,20,30,40,50,1.0, 0.9, 0.7, 1.5, 2.0, 2.4, 3.2, 2.0, 2.7, 3.5,
     1.0, 4.0, 3.6, 2.7, 5.7, 4.6, 6.0, 6.8, 7.3];

然后 .mod 使用 MIP 和距离的绝对值

execute
{
  cplex.tilim=10;
}

int n=...;
range points=1..n;
float x[points]=...;
float y[points]=...;

int nbLines=2;

range lines=1..nbLines;

dvar float a[lines];
dvar float b[lines];

// distance between a point and a line
dvar float dist[points][lines];
// minimal distance to the closest line
dvar float dist2[points];

dvar float+ obj;
minimize obj; 

subject to
{
obj==sum(i in points) dist2[i];

forall(i in points,j in lines) dist[i][j]==abs(b[j]*x[i]+a[j]-y[i]);

forall(i in points) dist2[i]==min(l in lines ) dist[i][l];

}

// which line for each point ?
int whichline[p in points]=first({l | l in lines : dist2[p]==dist[p][l]});

execute
{

writeln("b = ",b);
writeln("a = ",a);
writeln("which line = ",whichline);
}

给予

b =  [0.6375 10]
a =  [0.58125 0]
which line =  [2 2 2 2 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

二次方程式

int n=...;
range points=1..n;
float x[points]=...;
float y[points]=...;

int nbLines=2;

range lines=1..nbLines;

dvar float a[lines];
dvar float b[lines];

dvar float distance[points][lines];
dvar boolean which[points]; // 1 means 1, 0 means 2

minimize sum(i in points,j in lines) distance[i][j]^2; 

subject to
{
forall(i in points,j in 0..1) (which[i]==j) => (distance[i][2-j]==b[2-j]*x[i]+a[2-j]-y[i]);

}



execute
{

writeln("b = ",b);
writeln("a = ",a);
writeln("which = ",which);
}

给予

b =  [0.61077 10]
a =  [0.42613 0]
which =  [0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]

将两条线拟合到一组二维点

Fitting two lines to a set of 2D points

algorithm

optimization

curve-fitting

computational-geometry