我认为这可以被表述为一个明确的优化问题:
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}
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模型。
一个小测试证实了这个方法的可行性:
这个问题有50个点,在一台缓慢的笔记本电脑上,使用Cplex的MIQP求解器需要1.25秒钟。这可能是一个统计问题,MIP/MIQP方法可能会有所帮助。
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];
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]