Find two lines around which points were randomly generated

Question

Given a list of points that were randomly generated around two lines, find two new lines that match the original lines as closely as possible. Here's the function definition:

getLinesOfBestFit(x: numpy.array, y: numpy.array) -> Tuple[Line, Line]

Line = namedtuple('Line', [('m', float), ('b', float)])

I saw that there is a question about this problem already on StackOverflow, however, I found it difficult to understand. Please help me with understanding. I would like a solution that doesn't rely too heavily on libraries other than numpy so that I can understand the principles being applied.

I can use polyfit to find the line of best fit for a set of points as follows:

m, b = numpy.polyfit(x, y, deg=1)

However, before I can do this approach I need to split the set of input points into two sets of points, one set for each line. Naively, I can achieve that by considering all possible lines, then finding the two lines that minimise the distance of all points from any line. This approach is very inefficient though.

Answer 1

This falls into the general class of clustering problems. After the points are clustered it is straight forward to fit the two lines; however, it is possible to directly formulate it as an unconstrained least squares problem in terms of the distance of each point to the closest of the two candidate lines.

Given data $(x_i, y_i)$ formulate formulate a residual sum of squares as $$ J(b_1, m_1, b_2, m_2) = \sum_i \operatorname{min}( |y_i - (m_1 x_i - b_1)|^2, |y_i - (m_2 x_i - b_2)|^2) $$ then solve for the estimates $\hat b_1, \hat m_1, \hat b_2, \hat m_2 = \arg \min J(b_1, m_1, b_2, m_2)$

A potential issue with this approach is that the residual function is non-differentiable for certain values of the parameters. (If this will be a practical issue is probably data specific). One way around this would be to use a smooth-min function, but in practice using a conditional may work.

The following code implements the naive (non smooth minimum) approach.

import numpy as np
from scipy.optimize import least_squares

# Generate example data for two lines
b1,m1 = 2, 4
b2,m2 = 3,-3
n_samples = 1000
x1 = np.random.rand(n_samples)
y1 = b1 + m1*x1
x2 = np.random.rand(n_samples)
y2 = b2 + m2*x2
data =  np.hstack([[x1,y1], [x2,y2]])

# Define a residual function to minimize
def residual(x, data):
""" 
Input
    x (array): m1, b1, m2, b2
    data (array 2 x N): first row x-values and second row y-values of   data
Returns: 
    min_err (array N): min_err[k] is the y-distance of each point in data[:, k] to
                       either b=x[0] m=x[1]  or b=x[2] m=x[3], whichever is smaller
"""
    l1_err = data[0] * x[1] + x[0] - data[1]
    l2_err = data[0] * x[3] + x[2] - data[1]
    min_err = np.where(np.abs(l1_err) < np.abs(l2_err), l1_err, l2_err)
    return min_err

# Set initial guess of intercepts and slopes for two lines
x0 = 0, 0, 1, 1 # m1, b1, m2, b2
result = least_squares(residual, x0, args=(data,)) #Finds the correct 
print(result.x)
# >>> 3, -3, 2, 4

You could likely use some heuristic to get good initial estimates for x0, e.g. find the approximate "corners" of the points and find the lines going through those points at a cross.

This approach could be generalized to higher dimensions, but the parameterization of the lines would be a bit different (each line would require a point and a direction vector), and one would likely want to compute the orthogonal distance to the each line unless there is some reason to prefer one coordinate to another.

As a final note, given the specific form of the problem, the non-convex optimization procedures noted in the answers to the question you linked, may be more robust than this approach.