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

6
给定一组2D点(图片中的黑点),我需要选择两条线来代表这些点。 我可能需要最小化[到两条线中距离较近的x的距离]^ 2之和的平方和,尽管如果任何其他指标能使此问题更容易,则也可以接受。
明显但无效的方法是尝试在所有2 ^ n个分区上进行min二次项。 一个实用的方法可能是迭代改进,也许从随机分区开始。
是否有研究处理此问题的算法?

two lines


1
你可以对所有点进行线性回归(从而找到一条将上部点与下部点分开的直线)。然后分别对下部和上部点集进行另一次线性回归。 - BitTickler
2
是的,我只是不确定如果这些线相交时会有多好的效果。 - Rafał Dowgird
那么如果这两条线互相垂直,也是可能的吗? - BitTickler
此外,显然并不总是存在唯一的最佳解决方案(例如,如果点是完全随机的(白噪声)。那么它仍应该产生2条线?) - BitTickler
白噪声可能有一个独特的最佳解决方案。像仅有1个点或多个严格线性的点这样的情况显然是不唯一的,因为第二条线完全没有影响。 - Rafał Dowgird
3个回答

6

我认为这可以被表述为一个明确的优化问题:

 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),并且更容易解决。一些其他的公式在这里。其中一些提到的模型不需要我上面使用的边界。需要注意的是,这些模型提供了被证明的最优解(对于非凸模型,这需要全局求解器进行求解;而凸模型则更容易,不需要这样做)。此外,我们还可以形成一个L1目标函数,而不是最小二乘目标函数。在后一种情况下,我们得到了一个线性MIP模型。

一个小测试证实了这个方法的可行性:

enter image description here

这个问题有50个点,在一台缓慢的笔记本电脑上,使用Cplex的MIQP求解器需要1.25秒钟。这可能是一个统计问题,MIP/MIQP方法可能会有所帮助。


3

我能想到两个相关概念。

首先,我的直觉是应该有一种方法将这个问题解释为估计混合模型的参数。我还没有解决细节问题,因为参数估计通常使用期望最大化算法,我只能描述一下它是如何工作的:初始化一个包含两个部分的分区,然后交替对每个部分进行回归,并根据点到新回归线的距离重新分配点。

其次,假设输入相对清洁,您应该能够使用RANSAC获得良好的初始分区。对于一些小k,取两个不相交的k个点的随机样本,并通过它们适应线条,然后指定每个其他点。为了获得(100-x)%的成功率,您需要重复约ln(100/x)× 22k-1次,并选择最好的一个。


我曾经写过类似的答案,但因为所有内容都在链接中而被踩了。简要解释一下链接术语可能是个好主意,这样即使不跟随链接,文本也能变得有用。 - BitTickler
另外,您可能想将因子分析(https://en.wikipedia.org/wiki/Factor_analysis)添加到您的列表中。 - BitTickler
@BitTickler,我描述了这些想法的具体实例,即“将分区初始化为两个部分,然后交替在每个部分上运行回归并根据它们到新回归线的距离重新分配点”和“对于一些小的k,取两个不相交的随机样本k个点,并通过它们拟合线,然后分配其他点”。我不明白因子分析在这里如何适用。 - David Eisenstat
仅仅因为问题展示了一张图片,并不意味着它一定是一个图像处理问题。它也可能是统计变量的绘图(其中X,Y是变量)。如果您查看我在问题下面的新评论,您会发现这就是我想要找出的内容。 - BitTickler
1
谢谢!期望最大化算法让我了解到这个主题的论文:https://cse.sc.edu/~songwang/CourseProj/proj2005/report/guo.pdf - Rafał Dowgird

1
如果我从模型构建中的示例Model Building开始,使用 OPL CPLEX,那么请让我先对 .dat 文件进行一些补充curve fitting
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];

然后使用MIP进行.mod计算,并使用绝对值计算距离。
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]

enter image description here


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接