Optimal line such that maximum points are between an upper and lower boundary

Question

I have some 2D data and would like to find a line $y = mx + b$ such that a maximum number of points from the data is captured within the area between $y = mx + b + margin$ and $y = mx + b - margin$. Please see the image for an example (the yellow lines are boundaries, where margin = 20).

If it helps, the values on the y-axis are necessarily integers.

Currently, I have used linear regression (from scikit-learn) to find a line. However, this is not guaranteed to put the most points between the boundaries. I have also looked a bit into support vector regression (SVR), which seems to try to put every point within a certain epsilon of the line, but I'm not sure if this is the right application or if SVR is overkill.

Any Python code would be appreciated, but more generally I would just like a conceptual solution. Does anyone have any insights that would help solve this problem? Thank you!

Edit: the number of points $n$ is in the range of 5000-10000, and theoretically I only need to compute this once.

Answer 1

You can approximate the answer in $O(N log N + kN)$ time where $k$ is the number of sample lines you wish to test. To do so:

1. Choose a uniform set of $k$ points $\theta$ spanning the range [0,π). These are your sample lines.
2. Do a loop starting with the line for $\theta=0$. ($O(N)$)
   1. Get a unit vector $u$ in the direction of $\theta$. ($O(N)$)
   2. Get a vector $\hat u$ perpendicular to this vector. ($O(N)$)
   3. Project all the points onto $\hat u$. ($O(N)$)
   4. Sort the points along $\hat u$. ($O(N)$)
   5. Slide a window the size of your margin along the sorted points and note when the maximal number of points are within this window along with the midpoint of the window when this occurs. ($O(N)$)
3. Return information about the line which had the maximum point count.

The algorithm takes $O(N log N)$ time for initially sorting the points. However, after this we are projecting the points onto a slowly rotating line. The slow rotation means that each reprojection results in the points already being very close to their previously sorted positions:

Most sorting algorithms are able to deal with this situation in nearly-linear time. Note also that as the number of samples increases, the degree to which the points are already sorted increases as well.

The full algorithm looks like this:

#include <algorithm>
#include <cassert>
#include <cmath>
#include <functional>
#include <iostream>
#include <random>
#include <vector>



class FracPoint {
 public:
  float frac = std::numeric_limits<float>::quiet_NaN();
  float x;
  float y;
  FracPoint(float x0, float y0) : x(x0), y(y0) {}
  bool operator<(const FracPoint &o) const {
    return frac<o.frac;
  }
  bool operator>(const FracPoint &o) const {
    return frac>o.frac;
  }
};

struct BestLine {
  int    count;
  double slope;
  double intercept;
};

typedef std::vector<FracPoint> Points;
typedef std::pair<int,double> CountIntercept;



///Projects points to their closest orthogonal point on a line at an angle theta
///to the x-axis. Sets each point's `frac` property to the number of unit
///vectors their projection is from the origin. This can be used to sort points
///along the length of the line.
void GetFractionalDistPointsProjectedToVector(
  Points &pts,
  const double theta
){
  const double tvx = std::cos(theta);
  const double tvy = std::sin(theta);

  //Get a line perpendicular to the one we just generated. This is the line we
  //sort the points along to find the maximal margin.
  const auto vx = -tvy;
  const auto vy = tvx;

  //Do a dot product (v dot s)/(s dot s) and multiply by s to get projection
  for(auto &pt: pts)
    pt.frac = (pt.x*vx+pt.y*vy)/(vx*vx+vy*vy);
}



CountIntercept FindBestIntercept(Points &pts, const double margin){
  Points::iterator upper=pts.begin();
  Points::iterator lower=pts.begin();
  int best_count = 0;
  double best_pos = 0;
  for(;upper!=pts.end();++upper){
    while(upper->frac-lower->frac>margin)
      ++lower;
    const auto count = upper-lower+1;
    if(count>best_count){
      best_count = count;
      best_pos = (upper->frac-lower->frac)/2;
    }
  }
  return CountIntercept(best_count,best_pos);
}



BestLine FindBestLine(Points &pts, const int samples, const double margin){
  int    best_count = 0;
  double best_theta;
  for(int i=0;i<samples;i++){
    const double theta = i*(M_PI/samples);
    //Project points onto the line defined by theta
    GetFractionalDistPointsProjectedToVector(pts, theta);
    //Since points are already mostly sorted from the last projection, use
    //insertion sort's O(N) best case scenario to make them totally sorted.
    // ska_sort(pts.begin(),pts.end(),[](const FracPoint &pt){return pt.frac;});
    std::sort(pts.begin(),pts.end());
    const auto count_intercept = FindBestIntercept(pts, margin);
    if(count_intercept.first>best_count){
      best_count = count_intercept.first;
      best_theta = theta;
    }
  }

  //Project points onto the line defined by the best theta
  GetFractionalDistPointsProjectedToVector(pts, best_theta);
  //Since points are being projected onto a "randomly" chosen line, we use an
  //O(N log N) sort
  std::sort(pts.begin(),pts.end());
  const auto count_intercept = FindBestIntercept(pts, margin);  

  //`count_intercept.second` isn't actually the intercept, it's the distance
  //along a line perpendicular to the one defined theta at which the best
  //intercept is located. Let's find the real intercept now.
  double intercept;

  {
    //Get vector along theta
    const double tvx = std::cos(best_theta);
    const double tvy = std::sin(best_theta);

    //Get perpendicular vector
    auto vx = -tvy;
    auto vy = tvx;

    //Multiply by `counter_intercept.second` to get absolute position of best
    //point.
    vx *= count_intercept.second;
    vy *= count_intercept.second;

    //Get slope of line
    double slope = std::tan(best_theta);
    //Use point-slope equation of line y-y1=m*(x-x1) and reformulate to y=m*x+b
    //to recover the actual intercept
    intercept = -slope*vx+vy;
  }

  return {best_count, std::tan(best_theta), intercept};
}



int main(int argc, char **argv){
  if(argc!=3){
    std::cout<<"Syntax: "<<argv[0]<<" <N> <Samples>"<<std::endl;
    return -1;
  }

  const int N       = std::stoi(std::string(argv[1]));
  const int samples = std::stoi(std::string(argv[2]));

  //Poor way of initializing PRNG. Fine for this application.
  std::mt19937 gen(1234567); 
  std::uniform_real_distribution<> dis(-100, 100);
  auto dice = std::bind(dis,gen);

  Points pts;
  for(int i=0;i<N;i++)
    pts.emplace_back(dice(),dice());

  const auto ret = FindBestLine(pts, samples, 25);

  std::cout<<"Count          = "<<ret.count<<std::endl;
  std::cout<<"Best slope     = "<<ret.slope<<std::endl;
  std::cout<<"Best intercept = "<<ret.intercept<<std::endl;

  return 0;
}

Now, we time the algorithm for 100 samples at various numbers of points, which gives the following:

1000,       0.014
10000,      0.083
50000,      0.445
100000,     0.943
500000,     5.615
1000000,   12.092
5000000,   68.748
10000000, 144.155

which, indeed, is almost linear (slope of 1.028 on a log-log graph with R2=0.9962).

OpenMP can be used as a simple way to parallelize the operation:

  #pragma omp parallel for firstprivate(pts)
  for(int i=0;i<samples;i++){
    const double theta = i*(M_PI/samples);
    //Project points onto the line defined by theta
    GetFractionalDistPointsProjectedToVector(pts, theta);
    //Since points are already mostly sorted from the last projection, use
    //insertion sort's O(N) best case scenario to make them totally sorted.
    // ska_sort(pts.begin(),pts.end(),[](const FracPoint &pt){return pt.frac;});
    std::sort(pts.begin(),pts.end());
    const auto count_intercept = FindBestIntercept(pts, margin);
    #pragma omp critical
    if(count_intercept.first>best_count){
      best_count = count_intercept.first;
      best_theta = theta;
    }
  }

Answer 2

I am not sure what your size for $n$ is or how frequently you need to compute the desired line, but the most simple solution that should be very close to optimal is simply taking the point set of size $n$, constructing $\binom{n}{2}$ lines and then testing each of these lines in $O(n)$ time each to see how good each line is. With this naive algorithm, we have an overall $O(n^3)$ runtime while only really using $O(n)$ space.

If you construct a data structure for half-space range searching with the goal of counting how many points intersect a given half-space, you could reduce the runtime of the above algorithm to $O(n^2 \log n)$ using $O(n^2)$ space. This is because an optimal version of this type of data structure can return the number of intersections for a half-space in 2D in $O(\log n)$ time and you would only need to do two of these operations (one half-space corresponding to each upper and lower margin line) per line. So for $O(n^2)$ lines with $O(\log n)$ amount of work per line, you have $O(n^2 \log n)$ runtime.

My intuition tells me there might be some randomized algorithm that could approximately solve this problem really efficiently, but I would have to think more about it.