Fitting two lines to a set of 2D points

Question

Given a set of 2D points (black dots in the picture) I need to choose two lines to somehow represent those points. I probably need to minimize the sum of squares of [distance of x to the closer of two lines]^2, although if any other metric makes this easier, this is fine too.

The obvious but ineffective approach is to try min squares over all 2^n partitions. A practical approach is probably iterative improvement, maybe starting with a random partition.

Is there any research on algorithms to handle this problem?

Hemichordate

You could do a linear regression over all points (and thus find a line dividing the upper from the lower points). Then for both lower and upper set of points, do another linear regression, respectively. – Unthoughtof

Yup, I'm just not sure how well this would work in case the lines intersect. – Asexual

So it would also be possible, if the 2 lines were orthogonal to each other? – Unthoughtof

Also, clearly there cannot always be a unique best solution (e.g. if the points are totally random (white noise). Then it should still yield 2 lines? – Unthoughtof

White noise probably has a unique best solution. Cases like 1 point only or multiple points that are strictly linear are obviously not unique, as the second line doesn't matter at all. – Asexual

Answer 1

Two related concepts from the literature come to mind.

First, my intuition is that there should be a way to interpret this problem as estimating the parameters for a mixture model. I haven't worked out the details, since parameter estimation typically uses expectation--maximization, and I can just describe how that would work: initialize a partition into two parts, then alternately run a regression on each part and reassign points based on their distance to the new regression lines.

Second, assuming that the input is relatively clean, you should be able to get a good initial partition using RANSAC. For some small k, take two disjoint random samples of k points and fit lines through them, then assign every other point. For a (100 − x)% chance of success you'll want to repeat this about ln(100/x) × 2^2k−1 times and take the best one.

Biforate

I once wrote answers like that and got down voted because all the content is actually in the links. It might be a good idea to explain briefly what the terms of the links mean, so the text becomes useful without following the links. – Unthoughtof

also, you might want to add factor analysis (en.wikipedia.org/wiki/Factor_analysis) to your list. – Unthoughtof

@Unthoughtof I described the concrete instantiations of these ideas, viz. "initialize a partition into two parts, then alternately run a regression on each part and reassign points based on their distance to the new regression lines" and "For some small k, take two disjoint random samples of k points and fit lines through them, then assign every other point". I don't see how factor analysis applies here. – Biforate

Just because the question displays an image, it does not necessarily mean it is an image processing problem. It could also be a plot of statistical variables (with X, Y being the variables). If you look at my new comments below the question, you see, this is what I try to find out. – Unthoughtof

Thanks! Expectation-maximization led me to a paper on this topic: cse.sc.edu/~songwang/CourseProj/proj2005/report/guo.pdf – Asexual

Answer 2

I think this can be formulated as an explicit optimization problem:

 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}

We do the assignment of points to a line and the regression (minimization of the sum of the squared residuals) at the same time.

This is a non-convex quadratically constrained mixed-integer quadratic programming problem. Solvers that can handle this include Gurobi, Baron, Couenne, Antigone.

We can reformulate this a bit:

 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)

This makes it a straight convex MIQP model. This will allow more solvers to be used (e.g. Cplex) and it is much easier to solve. Some other formulations are here. Some of the models mentioned do not need the bounds I had to use in the above big-M formulation. It is noted these models deliver proven optimal solutions (for the non-convex models this would require a global solver; the convex models are easier and don't need this). Furthermore, instead of a least squares objective, we can also form an L1 objective. In the latter case we end up with a linear MIP model.

A small test confirms this works:

This problem has 50 points, and needed 1.25 seconds using Cplex's MIQP solver on a slow laptop. This may be a case of a statistical problem where MIP/MIQP methods have something to offer.

Adenoma

Answer 3

In OPL CPLEX if I start with the example curve fitting from Model Buidling

Let me add a few points to the .dat first

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];

then .mod with MIP and absolute value for distance

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);
}

gives

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]

With quadratic some reformulation

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);
}

gives

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]

Lanai