3
$\begingroup$

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).

linear regression showing boundaries

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.

| cite | improve this question | | | | |
$\endgroup$
  • $\begingroup$ I wonder what this is for, or what motivates this question? $\endgroup$ – Richard Oct 28 '19 at 1:36
  • $\begingroup$ Is the margin given? Or is it also an optimization variable? $\endgroup$ – Wolfgang Bangerth Oct 28 '19 at 3:51
1
$\begingroup$

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:

Rotating least-squares animation

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;
    }
  }
| cite | improve this answer | | | | |
$\endgroup$
  • $\begingroup$ The animation looks very familiar... Did you take it from a CrossValidated.SE question related to PCA? Please put the reference for the animation if you did not create it by yourself... When, I see the animation, I thought, what if you calculate PCA of the point cloud find just first principal axis and then your line that you are looking for is that principal axis and margin would be the standard deviation? $\endgroup$ – Alone Programmer Oct 28 '19 at 14:42
1
$\begingroup$

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.

| cite | improve this answer | | | | |
$\endgroup$
  • $\begingroup$ Thanks! I've added this in the question, but n is 5000-10000 and theoretically I only need to compute this once. I'll try something along these lines. $\endgroup$ – jasli Oct 26 '19 at 20:40
  • $\begingroup$ @jasli yeah, if you only have to compute this quantity a small number of times, the brute force solution should be able to get you the answer within a few minutes I would think if implemented well (given the number of points you have). If you took advantage of multiple cores on your system, it could be reduced further since this is embarrassingly parallel. $\endgroup$ – spektr Oct 27 '19 at 14:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.