最小包含球

説明

http://citeseer.ist.psu.edu/welzl91smallest.html

d 次元の点たちをすべて含むような球を move-to-front heuristics で求める．

最小包含球問題はナイーブには次の全探索アルゴリズムで解くことができる:

def min_disk(P, R):
  if ps is empty:
    D := ball enclosing R
  else:
    choose p <- P
    D := min_disk(P-{p}, R)
    if p is not in D:
      D := min_disk(P-{p}, R+{p})
  return D

どのように p を取るかによって計算量は大きく変わる．最も悪いケース，すなわちすべての点の組について包含球を求めることになる場合の計算量は O( n^{d+1} ) であるが，ランダムに選んだとすると再帰的探索が動く確率がその点が包含球の境界である確率なので O(d/n) となり，全体で O(2^d n) と評価できる．

ランダムに取ってもこれだけうまくいくのだから，良い選び方をすればより効率が良くなるのではないか，と考えるのは自然である．そのような点の選び方の基準として move-to-front heuristics がある．これは包含球の境界になりそうなものから取っていく heuristics である．具体的には p is not in D が成立したとき，p は境界候補なので P の先頭に持っていく( move-to-front )．こうすることによって多くのケースで効率よく解くことができるとわかっている．

細かいところ: 新しい点が付け加えられたときに最小包含球がどのように変化するかは，差分ベクトルを gram schmidt で直交化して得られる新しい方向ベクトルの分だけ中心ベクトルをずらせばよい．この方法を取ると縮退するケースを考慮しなくてよいので実装が簡単になる．

計算量

平均 O(2^d n)．

使い方

コンパイル前に意味があるもの．

void add(const point& p);

新しい点を min_ball に含めるべき点列に加える．

void compile();

加えられている点から最小包含球を求める．

コンパイル後に参照して意味があるもの:

point center;

最小包含球の中心．

number radius;

最小包含球の半径の二乗．

void support(OUT out);

最小包含球を張る点を OUT にコピーする．

表面にあるすべての点ではないので注意．そちらは半径と中心から再計算すること．

TODO

• 長すぎるのでコンパクトに．

ソースコード

struct min_ball {
  point center;
  number radius2;
  void add(const point& p) {
    ps.push_back(p);
  }
  void compile() {
    m = 0;
    center.resize(dim, 0);
    radius2 = -1;
    make_ball(ps.end());
  }
  template <class OUT>
  void support(OUT out) {
    copy(ps.begin(), supp_end, out);
  }
  min_ball() {
    for (int i = 0; i <= dim; ++i) {
      c[i].resize(dim, 0);
      v[i].resize(dim, 0);
    }
  }
private:
  list<point> ps;
  list<point>::iterator supp_end;
  int m;
  point v[dim+1], c[dim+1];
  number z[dim+1], r[dim+1];
  void pop() { --m; }
  void push(const point& p) {
    if (m == 0) {
      c[0] = p; r[0] = 0;
    } else {
      number e = dist2(p, c[m-1]) - r[m-1];
      point delta = p - c[0];
      v[m] = p - c[0];
      for (int i = 1; i < m; ++i)
        v[m] -= v[i] * dot(v[i], delta) / z[i];
      z[m] = dot(v[m], v[m]);
      c[m] = c[m-1] + e*v[m]/z[m]/2;
      r[m] = r[m-1] + e*e/z[m]/4;
    }
    center  = c[m];
    radius2 = r[m]; ++m;
  }
  void make_ball(list<point>::iterator i) {
    supp_end = ps.begin();
    if (m == dim + 1) return;
    for (list<point>::iterator k = ps.begin(); k != i; ) {
      list<point>::iterator j = k++;
      if (dist2(*j, center) > radius2) {
        push(*j); make_ball(j); pop();
        move_to_front(j);
      }
    }
  }
  void move_to_front(list<point>::iterator j) {
    if (supp_end == j) ++supp_end;
    ps.splice (ps.begin(), ps, j);
  }
};

Verified

• PKU 2069 (Super Star), 31 MS, 280 KB