Раздел «Образование».FIVTLecturesTerm4Lecture05:
<< Список лекций ФИВТ, 4-й семестр, 2009 г.

Лекция 5. Задача кластеризации точек в метрическом пространстве

О задаче кластеризации

Алгоритмы кластеризации часто следуют одной из двух идей:

Качество кластеризации — понятие, определяемое в каждой задаче по-своему.

Для случая точках в метрическом пространстве пространстве, задача кластеризации могла бы выглядеть следующим образом:

Задача: Покрытие точек шарами

Даны точки в метрическом пространстве. Найдите минимальное количество шаров радиуса R, которые покрывают все данные точки.

В этом случае качество кластера (шара) можно определить как количество точек, которые он покрывает.

В некоторых задачах мера "плохости кластера" есть отношение его размера на количества точек в кластере, где размер может быть определен как объём шара, внутрь которого помещаются все точки кластера, либо как диаметр множества точек кластера (расстояние между двумя самыми удаленными точками кластера), либо как среднеквадратическое отклонение точек от центра кластера (если в пространстве точек определено понятие центра множества точек), либо как-нибудь ещё.

Кластеризация может быть шагом различных эвристических алгоритмов. Например, при приближённом решении Задачи Коммивояжера имеет смысл разбить точки на группы близких точек — кластеров, чтобы решать Задачу Коммивояжера отдельно в каждом кластере, а затем на множестве центров кластеров.

Задача "Лучшая кластеризация"

Задача 046: Лучшая кластеризация

Даны точки в трёхмерном евклидовом пространстве. Разбейте множество точек на множество непересекающихся подмножеств — кластеров — размера больше 1 так, чтобы величина D была минимальна.

D = сумма D_i, i = 1,2,...,K;

D_i = дисперсия точек в i-м кластере (средний квадрат расстояния точек кластера от центра кластера).

Количество кластеров K можно выбрать любым.

Рассмотрим алгоритм основанный на идее слияния.

  1. Отберём некоторую хорошую часть возможных рёбер. Для каждой точки найдём не более 20 соседей. Не будем брать ребра, которые слишком длинные (например, которые в 2 раза длиннее, чем средняя длина ребра для концов рёбер или которые в 2 раза длинее чем самое короткое найденное на текущий момент ребро исходящее из первой вершины ребра). Отбор некоторой части рёбер позволит уменьшить время работы алгоритма.
  2. Отсортируем отобранные рёбра.
  3. Будем перебирать рёбра в порядке возрастания длины и выполнять слияние кластеров, которым принадлежать концы ребра.
    • Если это один и тот же кластер, то переходим к следующему ребру
    • Если ухудшает или не слишком сильно улучшает целевую функцию D, то переходим к следующему ребру
  4. Избавимся от одноточечных кластеров. Например, присоединим их к ближайшему кластеру, а лучше к такому кластеру, присоединение к которому не сильно ухудшает целевую функцию D.

#include <cstdio>
#include <vector>
#include <map>
#include <algorithm>
using namespace std;

#define tr(c,it) tr2(c,it, c.begin())
#define tr2(c, it, it2)  for( typeof(c.begin()) it = it2; it != c.end(); it++ )
#define all(c) c.begin(),c.end()
#define pb push_back
#define mp make_pair

typedef vector<Point>::iterator Pi;
typedef pair<double, pair<Pi, Pi> > Edge;
inline double sqr(double a) {return a*a;}

class Cluster;
class Point;

class Point {
public:
  double x,y,z;
  double min_dist;
  int neib_count;
  Cluster *cluster;

  Point(double _x = 0, double _y = 0, double _z = 0): 
    x(_x), y(_y), z(_z), min_dist(1E+255), neib_count(0) {};

  Point operator *(const double &r) {
    return Point(x * r, y * r, z * r);
  }

  Point operator +(const Point &r) {
    return Point(x + r.x, y + r.y, z + r.z);
  }
  
  &Point operator +=(const Point &r) {
    x += r.x; y += r.y; z += r.z;
    return *this;
  }
};

inline double dist2(const Point &p, const Point &q) {
  return sqr(p.x - q.x) + sqr(p.y - q.y) + sqr(p.z - q.z);
}

class Cluster {

public:
  Point center;
  double sigma2;
  int size;
  double min_dist;

  // ссылка на кластер, в который был "влит" данный;
  // равна this, если кластер реальный, 
  // то есть существует (не объединён с другим кластером)
  Cluster *parent;

  // номер кластера
  int id;

  // конструктор одноточечного кластера
  Cluster(const Point &p) {
    init();
    center = p;
    size = 1;
  }

  Cluster() {
    init();
  }

  void init() {
    sigma2 = 0;
    size = 0;
    id = 0;
    parent = this;
    min_dist = 1E+255;
  };

  // ленивое вычисление реального кластера,
  // внутри которого находится данный
  // (см. алгоритм Крускала вычисления остовного дерева)
  Cluster* root() {
    return parent = ( (parent == this) ?  this : parent->root() );
  }
};

Cluster join(Cluster &c1, Cluster &c2) {
  Cluster c;
  c.size = c1.size + c2.size;
  c.center = (c1.center*c1.size + c2.center*c2.size) * (1.0/ (c.size));
  double d1 = dist2(c.center, c1.center);
  double d2 = dist2(c.center, c2.center);
  c.sigma2 = ((c1.sigma2 + d1) * c1.size + (c2.sigma2 + d2) * c2.size) / c.size;
  return c;
}

void register_joined_cluster(Cluster &c, Cluster &p, Cluster &q) {
  Cluster *c2 = new Cluster(c);
  q.parent = p.parent = c2;
  c2->parent = c2;
}


int main () {
  int n;
  scanf("%d", &n);
  vector<Point> points(n);
  vector<Edge> edges;

  tr(points, p) {
    scanf("%lf%lf%lf", &p->x, &p->y, &p->z);
  };

  tr(points, p) {
    tr2(points, q, p+1) {
      // обходим все пары точек
      double d = dist2(*p, *q);
      if ( p->neib_count > 20 && q->neib_count > 20 ) continue;
      
      // Здесь можно написать различные проверки на хорошесть ребра p--q
      // Например:
      if ( p->min_dist * 1.8 <= d && q->min_dist * 1.8 <= d ) continue;
      
      if (p->min_dist > d) p->min_dist = d;
      if (q->min_dist > d) q->min_dist = d;
      p->neib_count++;
      q->neib_count++;
      // но в массив ребёр добавляем лишь перспективные,
      // прошедшие выполненные выше проверки
      edges.pb( mp(d, mp(p, q) ) );
    }
    p->cluster = new Cluster(*p);
  }

  sort( all(edges) );
  
  // первый проход по рёбрам
  tr( edges, e ) {
    Cluster *p = e->second.first->cluster->root();
    Cluster *q = e->second.second->cluster->root();

    if ( p == q ) continue;
    Cluster c = join(*p, *q);
    // TODO: проверям, хороший ли получился кластер
    // и регистрируем его, если хороший;
    // на первом проходе условие хорошести можно сделать жестким;
    // если оба кластера одноточечные, то нужно выполнять слияние;
    // ... 
  }
    
  // второй проход по рёбрам
  tr( edges, e ) {
    Cluster *p = e->second.first->cluster->root();
    Cluster *q = e->second.second->cluster->root();
    if ( p == q ) continue;
    // здесь условие объединения может быть уже помягче
    
    // здесь также нужно собрать одноточечные кластера --
    // в конце работы алгоритма их быть не должно
  }
  
  // TODO: разобраться с одноточечными кластерами
  
  // Пронумеруем кластера и выведем результат:
  int cnt = 0;

  tr(points, p) {
    if ( p->cluster->root()->id == 0 )  p->cluster->root()->id  = ++cnt;
  }

  printf("%d\n", cnt);

  tr(points, p) {
    printf("%d\n", p->cluster->root()->id);
  }

  return 0;
}

Важные замечания

Метаэвристики

Идеи (метаэвристики), используемые в эвристических алгоритмах

Сами алгоритмы кластеризации могут играть роль шага огрубления (local + scale) данных для решения других задач на графах. Вместо задачи на точкам мы переходим к задаче на кластерах, которых существенно меньше и к которым можно применять более сложные (с большей ассимптотикой по времени) алгоритмы. После получения грубого решения на кластерах можно вернуться к задаче на точках и уточнить грубое решение.

Алгоритм Крускала

В основе лежит алгоритм Крускала.