Careercup

2014-05-02 07:49

题目链接

原题：

Given a set of n points (coordinate in 2d plane) within a rectangular space, find out a line (ax+by=c), from which the sum of the perpendicular distances of all the points will be minimum. This can has a general usecase like, in a village there are few house, you have to lay a road, such that sum of all the approach roads from each house will be minimum.

题目：给定二维平面上的n个点，请找到一条直线，使得这些点到这条直线的距离之和最小。

解法：点到直线的距离是|a * x + b * y + c| / sqrt(a * a + b * b)。那么这个距离之和就是一堆绝对值的和。你可能很快会想到线性回归的那幅图片：一堆点均匀分散在一条直线的两侧。很遗憾，这个不是最优解，因为最小二乘满足的条件是方差最小，也就是平方和最小(学过线性代数应该知道p-范数的概念，平方和再开方就是2-范数。而范数，就是线性空间里对长度的度量算子。范数最小的时候，也就是我们用最小二乘期待的结果)。平方和最小和绝对值和最小并非同一概念。那么绝对值的和什么时候最小呢？我找到了一篇不太知名的论文，源头链接在此(里面有人提到了那篇论文的链接，我不直接贴链接是因为那人直接把论文内容贴了出来，有侵权嫌疑)。证明实在太绕了，所以我没工夫仔细推理完，只是偷懒引用了其中的一个结论：满足距离之和最短的直线，一定会穿过n个点中的两点。有这个结论，就可以用两层循环遍历所有的直线组合，然后再用一层循环计算所有点的距离。选取距离之和最小的那条直线作为结果。这样写出来的算法，时间复杂度是O(n^3)的。看来不靠谱，但好歹是个解法吧。如果你用最小二乘法来做，时间复杂度应该是O(n^2)的，我相信得到的就算不是最优解，也应该很接近。所以有时你连算法都拿不准的时候，还不如用次优解来代替最优解，因为其他方面的优势可以作为权衡因素。

代码：

// http://www.careercup.com/question?id=4907555595747328
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>
using namespace std;

struct Line {
    double a;
    double b;
    double c;
    Line(double _a = 0, double _b = 0, double _c = 0): a(_a), b(_b), c(_c) {};
};

struct Point {
    double x;
    double y;
    Point(double _x = 0, double _y = 0): x(_x), y(_y) {};
};

double calcDist(const Point &p, const Line &line)
{
    return abs(line.a * p.x + line.b * p.y + line.c) / sqrt(line.a * line.a + line.b * line.b);
}

void calcLine(const Point &p1, const Point &p2, Line &line)
{
    line.a = p2.y - p1.y;
    line.b = p1.x - p2.x;
    line.c = -((line.a * p1.x + line.b * p1.y) + (line.a * p2.x + line.b * p2.y)) / 2.0;
}

int main()
{
    vector<Point> p;
    int n;
    Line line, min_line;
    int i, j, k;
    double dist, min_dist;
    
    while (cin >> n && n > 0) {
        p.resize(n);
        for (i = 0; i < n; ++i) {
            cin >> p[i].x >> p[i].y;
        }
        
        do {
            if (n == 1) {
                line = Line(0, 1, -p[0].y);
                break;
            } else if (n == 2) {
                min_dist = 0;
                calcLine(p[0], p[1], min_line);
                break;                
            }
            
            min_dist = -1;
            for (i = 0; i < n; ++i) {
                for (j = i + 1; j < n; ++j) {
                    dist = 0;
                    calcLine(p[i], p[j], line);
                    for (k = 0; k < n; ++k) {
                        if (k == i && k == j) {
                            continue;
                        }
                        dist += calcDist(p[k], line);
                    }
                    if (min_dist < 0 || dist < min_dist) {
                        min_dist = dist;
                        min_line = line;
                    }
                }
            }
        } while (0);
        
        if (min_line.a != 0.0) {
            cout << min_line.a << 'x';
        }
        if (min_line.b != 0.0) {
            cout << setiosflags(ios::showpos) << min_line.b << 'y';
        }
        if (min_line.c != 0.0) {
            cout << min_line.c;
        }
        cout << resetiosflags(ios::showpos) << "=0" << endl;
        cout << min_dist << endl;
    }
    
    return 0;
}