В универе дали задание, написать алгоритм кластеризации точек методом k-means. Задание дано по принципу "получил - делай" 

Проблема заключается в том, что я не имею и малейшего представления что это такое, с чем это едят и тем более как это готовят =) Буду очень благодарен, если кто-нибудь поделится информацией или исходниками. 

гугл любит это делать  

1) выбираем произвольно K центров тяжести
2) кластеризуем, то есть проходимся по всем точкам и "прилепляем" их к ближайшему кластеру
3) рассчитываем центры тяжестей полученных кластеров
4) если центр тяжести изменился с последнего пересчета, то переходим на пункт 2 и продолжаем действия

Silent, а можно поподробнее?

K-means для чайников, учебное пособие в картинках:
http://ru.wikipedia.org/wiki/K-means

Спасибо

Вообще к-средних не алгоритм кластеризации. Этo скорее алгоритм оптимизации кластеров. Оn не может дать ни оценки сколько кластеров нужно ни каk выбирать начальные кластеры. А этo задачи на порядок сложнее чем оптимизация кластеров.

А что самое страшное, таk это то, что оn может никогда не сойтись либo может взять очень много итераций. А каждая итерация o(n*k) где n - количество точек, k - количествo центроидов.

neutrino

Оценка оптимального числа кластеров делается тем же способом, что и для других методов кластеризации - взглядом на график поведения остаточной дисперсии вдоль оси числа кластеров. Исключение из этих методов и из схемы определения - пожалуй, только лямбда-метод Загоруйко на основе построения соединяющего все классифицируемые точки минимального связного подграфа-остова (сама трудоемкость алгоритма построения минимального связного подграфа - квадратична по числу точек) и деления наиболее "нагруженных" ребер графа ("нагрузка" определяется специфическим образом, но это, к счастью, ничего не усугубляет).

Трудозатраты к-средних на каждой итерации тоже можно снизить, как это ни неочевидно. Ссылки на 2 статьи - в моей заметке.

neutrino! Все вроде как верно.
Но человек не ищет хороший алгоритм кластеризации. Ему конкретно задали сваять K-means, что бы по этому поводу не думала мировая научная общественность.

А в отношении самого алгоритма...
Очень забавные результаты можно наблюдать на "регулярных" множествах.

VictorTsaregorodtsev, Спасибо. Я почитал немного информации на Вашем сайте. Сурьезно. Я кстати тоже написал один алгоритм кластеризации. Без претензий на новшество (может такое есть). Вся идея сводится к тому, что я ищу пики функции плотности распределения данных. А функция плотности строится следующим образом. В каждой точке данных (в гиперкубе, где количество измерений = количеству признаков классификации) я строю гауссиан. То есть добавляю функцию гауссиана с мат ожиданием в самой точке и с некоторой дисперсией, что влияет также и на некоторую округу этой точки. Теперь остается найти максимумы получившейся функции плотности. Для этого я придумал следующее.
Можно найти такое расстояние между гауссианами, что если сближать их ближе этого расстояния функция суммы этих гауссианов будет иметь своим максимумом середину этого расстояния. Конечно тут все зависит от конкретных параметров гауссианов (дисперсий). Так вот, допустим мы нашли это расстояние - а. Далее мы делим весь наш гиперкуб на маленькие гиперкубы с диагональю а. Легко показать что для каждого такого гиперкуба может быть только 2 ситуации:
- есть один пик 
- пиков нет

Значит если мы будем идти из серидины каждого гиперкуба по градиенту, то обязательно наткнемся на максимум, если он есть.

Алгоритм работает очень быстро.

neutrino

Работает. Но в подобном подходе (в рамках восстановления многомерной и многомодовой плотности вероятности приемами непараметрической статистики) тоже очень уж много возникает вычислительных затрат, а именно:
1. Оптимизация параметра размытости гауссианы. Если решать задачу более-менее точно, то обычно всё укладывается в несколько проб разных значений и прогнозирование минимума (околооптимума) алгоритмом параболической аппроксимации, пробу в этом спрогнозированном значении и последующем уточнении минимума по трем точкам, образующим выпуклую вверх или вниз (в зависимости от того, как у нас критерий считается) параболу. Но - это же пробы на полной выборке...
2. Расчет градиента. Опять-таки приходится для расчета градиента в точке перепахивать всю выборку (отстрел точек, отстоящих более чем на три расстояния а от текущей точки, т.к. вылетающих за три сигмы и почти не влияющих, сильно расчет не ускорит, т.к. надо будет всё равно считать расстояние между двумя точками (которое входит в показатель экспоненты и при оценивании плотности, т.е. как было, так и остается).
Поэтому насчет быстро... Может быть, но на не сильно больших базах данных. Несколько прогонов выборки, как никак. А для к-средних, при коррекции по репрезентативным подвыборкам, решение может сойтись и за один полный прогон выборки (но для к-средних придется пробовать разное число кластеров для определения нужного их числа, т.е. может в итоге получиться баш на баш по времени вычислений, но к-средние-то запрограммировать проще и быстрее).

neutrino, Еще подумал. Ваш алгоритм может быть быстрее из-за того, что к-средним для каждой точки требуется n попарных вычислений расстояний между векторами и одно добавление вектора (точки) к вектору поправок для ближайшего ядра кластера, а Вашему требуется всего одно вычисление расстояния. Т.е. то, что выиграли на этапе вычисления градиента, можем не растерять за несколько шагов движения (с перерасчетом градиента) по нему.

neutrino, тьфу, чего-то сегодня торможу. 

В предыдущем посте - переменная k вместо n должна быть.

Да и найти-то Вашим способом можно только локальные максимумы плотности распределения вероятности, а как разбивать классифицируемые точки на кластеры? Только однократным исполнением прогона через к-средние (естественно, только прогона точек без последующей коррекции положений ядер). А это еще плюс к общему времени расчетов.
Но в таком гибриде не совсем хорошо сходится используемое к-средними кратчайшее расстояние между точками и соответствие именно восстановленной Вами плотности распределения вероятности - граница между кластерами может пройти там, где не находится локальный минимум плотности. У к-средних границы будут в виде ячеек Вороного, т.е. в виде кусочно-линейных замкнутых или разомкнутых (на краях облака данных) границ, а в Вашем алгоритме границы между кластерами должны быть в виде гладких кривых, но вот как их быстро и для возможности быстрого пользования в дальнейшем восстановить? Как бы не пришлось для каждой приходящей на кластеризацию точки (не входившей ранее в обучающее множество) использовать всё исходное обучающее множество для восстановления плотности (это - общее больное место непараметрического подхода, нельзя в нем отказаться от ранее использованной обучающей выборки и получить более компактную модель или редуцированно-квантованную выборку).

Утверждение об очень быстрой работе алгоритма, основанного на аналитической аппроксимации плотности распределения в многомерном пространстве, не то чтобы вызывает недоверие, а просто не понятно без приведения хотя бы самых общих идей реализации. Попытаюсь объяснить, что именно непонятно.
Не совсем понятно, но пусть. Предположим, что у нас есть способ быстрой оценки вклада «атома», без точного вычисления этого вклада. Ну и что? Ведь рядышком с далекой точкой может находиться очень много других точек. Все они тоже очень далеки и влияние каждой в отдельности очень мало. А сумма?

А ничего страшного. Нас же никто не заставляет аппроксимировать именно нормальными слагаемыми. Будем считать, что распределение «атома» задается обрезанным нормальным распределением. Тогда вроде все хорошо.

Не хорошо. Получаемая сумма не будет дифференцируемой везде функцией. А мы собрались искать чего-то, вычисляя градиенты.

Выход можно искать в построении подходящей функции окна, взяв за начальное приближение нормальную плотность, если уж она нам так нравится. Скажем, можно ее умножить на что-то, похожее на треугольник. Не будем вдаваться в детали, а только перечислим необходимые свойства результата.

Функция окна должна быть
1. Отличной от нуля в ограниченной области определения
2. Всюду дифференцируемой хотя бы один раз
3. Относительно быстро вычисляемой

Удовлетворить всем трем требованиям как-то можно. Но ведь в методе, о котором идет речь, вроде как используется многомерная нормальная плотность без изменений. Хотя это не факт – ведь подробности неизвестны.

В принципе, можно пойти по тому пути, что число слагаемых в аппроксимирующей функции существенно меньше общего числа точек.
Но, во-первых, из намека на алгоритм я не понял, что такой путь выбран.
Во-вторых, такой путь предполагает какой-то способ предварительного поиска областей сгущения точек. Саму такую идею лично я считаю правильной. И реализовать ее можно достаточно эффективно. Но хочу заметить, что поиск таких точек – это и есть, по существу, та самая задача поиска начальных центров кластеризации, необходимость решения которой расценивается, как недостаток метода K-means.


Это сообщение отредактировал(а) Sefko - 14.3.2009, 19:30

Приветствую!

Ого! Сколько понаписали то!

Сначала в общем. У меня вычислительные затраты идут на пожалуй две вещи:
1) построение функции плотности распределения. Если характеристик - c, точек - n а размер гауссиана в одном измерении - m, то сложность построения - O(n*m^c) обычно характеристик несколько, а для моего гиперкуба со стороной 100 пикселей я выбрал размер гауссиана - 15. Этот процесс довольно быстр но одновременно является самой сложной (по времени) частью алгоритма.
2) поиск пиков. Сложность: О((R/m)^c), где R - разрешение гиперкуба в одном измерении. Опять же у меня - 100. Слишком большое разрешение задавать бесполезно, т.к. точности потом все равно нужно добиваться единичным или несколькими прогонами k-means. На самом деле R/m - это я немного слукавил. Так происходит когда дисперсию гауссиана выбирают относительно большую, таким образом то самое максимальное расстояние между пиками должно быть m. Минус этой ситуации в том, что функция плотности в таком случае не будет гладкой. Но как показывает опыт, это не влияет на результат. По крайней мере на моих данных.

Оправдала себя еще и следующая эвристика: 
При поиске пика в каждом маленьком гиперкубе после разбиения делается так: определяем градиент в центре (либо подсчитан при построении функции плотности, либо тупо смотрим в у какого соседя значение выше) и скачим по этому направлению до границы гиперкуба. От туда градиентным поиском. Дело в том, что в большинстве гиперкубов на этом поиск и заканчивается, т.к. гиперкубов с пиками не так много.

Я понимаю, что звучит странно, но метод работает хорошо для данных в несколько и несколько десятков тысяч. На миллионах не проверял.
Ключевой момент - выбор дисперсии гауссианов. От этого зависит и результат и скорость работы.

Теперь по порядку.

То есть как раз дисперсию гауссианов? Да к сожалению я не автоматизировал оптимизацию этого параметра. Вот сейчас думаю как бы это сделать.


А аналитически никак нельзя?


А если при построении функции плотности? Ведь частичные производные считаются просто. Я не делал так, но думаю много времени это не добавит.


Как раз на очень больших базах получается куда быстрее k-means. ведь большого числа итеративных циклов зависимых от количества данных у меня нет. Единственное - построение функции плотности. Я экспериментально доказал на нескольких базах, что разрешение гиперкуба (или поля) с одной стороной в 100 пикселей отвечает любому требованию. Результат потом обрабатывается 1-2 итерациями k-means. Да выбор дисперсии гауссианов сильно влияет на результат. 
 

В том то и дело, что может и не может. В наших задачах мне пришлось от него отказаться. Слишком медленно работал а результат был не лучше чем у этого алгоритма. Да, k-means запрограммировать что плюнуть. 


При кластеризации иногда пользуются немного другим методом подсчета расстояния до центроида - по аналогии магнетизма. Все точки влияют на все центроиды, но чем дальше точка тем меньше она влияет (вроде магнитная сила убывает квадратично от расстояния). Как раз это здесь и происходит. Разбиение точек на кластеры может быть реализовано по разному. Можно просто прогоном k-means (так сделал я), можно с учетом каких-то дополнительных ограничений. Я думаю здесь можно использовать несколько иной метод. От каждой точки нужно идти по градиенту. К какому пику попали тот и кластер. Но я этот метод еще особо не обдумывал. Выглядит не быстрым.
 

Да, да здесь я с Вами согласен. У меня проскакивала мысль по этому поводу, но пока подхода у меня нет. Да и времени :(


Ну я вроде как написал уже, если это недостаточная аргументация я не берусь спорить. Вообще я не какой-то там специалист. Но мой алгоритм дал мне хорошие результаты за десятую от времени k-means.
 

В этом случае много точек далеко, не что иное как отдельный кластер, что прекрасно найдется этим алгоритмом. Ведь ему не надо знать изначально сколько кластеров вам надо.


Сумма гауссианов еще какая дифференциируемая функция! Да еще и сколько раз дифферениируемая! Градиенты как раз тут работают замечательно.


Правильно поняли. Я не выбрал такого подхода. Но вот интересно, если точек действительно много, то можно брать, скажем 10% от всех точек. По идеи они должны отражать начальное распределение.


VictorTsaregorodtsev, 
Sefko, 
Большое спасибо за интерес и ваши ответы!

Добавлено @ 23:09
Sefko, Не могли бы вы поподробнее написать про функцию окна. Что-то я все не до конца понял.

Это сообщение отредактировал(а) neutrino - 15.3.2009, 10:46