Алгоритмы / Еще раз о поиске наибольшего сгущения в облаке точек

В очередной раз мне попалась задача – найти в облаке точек место их наибольшего сгущения. На этот раз ситуация была такой:есть некоторое количество (можно считать, что не более 16 миллионов) измерений набора параметров. Число параметров в наборе – от 2 до 5.

измерение параметров может быть относительно успешным – тогда их результат будет неподалеку от истинного (параметры и тип распределения неизвестны), либо не успешным – тогда результат будет случайным (опять-таки с неизвестными параметрами распределения). Определить по одиночному измерению, было ли оно успешным, нельзя.

Можно считать, что точка сгущения существует. Если их несколько с близким качеством (которое формально так и не определяется), можно выдать любую.

ответ нужно дать за один проход по исходным данным: перевычислять их или сохранять целиком – дорого.

И, как обычно, хочется, чтобы алгоритм выглядел попроще, а работал побыстрее.

Для простоты можно считать, что измеренные значения параметров – вещественные числа от 0 до 1. Результат хочется получить с разрешением 3 десятичных знака для случая 5 параметров, и 4 или больше – для 4 и менее параметров.
Понятно, что большое число исходных точек делает сомнительными методы, основанные на обработке пар точек. Даже эффективный поиск расстояния до ближайшей точки потребует заметных усилий по реализации.
Построение гистограмм тоже непросто: заранее неизвестен даже желаемый размер ячейки. Возможно, например, что все измерения дадут результаты в шарике диамером 0.1, а успешные будут локализованы в области размером 0.003. И надо будет выдать положение этой области (в N-мерном пространстве!)
Можно было бы построить гистограммы распределения отдельных параметров, понадеявшись, что сгущение N-мерных точек можно будет восстановить по положениям сгущений проекций точек на отдельные координаты, но лучше оставить этот вариант в качестве запасного: терять информацию о взаимосвязях параметров опасно, а при проекциях могут возникнуть паразитные области сгущения.
Вариант, который мне показался наиболее перспективным – использование k-d дерева. Если мы построим бинарное дерево, каждый узел которого соответствует делению области пространства на две части по одной из координат (координаты перебираются по циклу – x1,x2,…,xk,x1,x2,…), и найдем количество точек, попавших в каждую из областей, то можно сделать так.
Пусть точек N. Выберем значение K, меньшее N, например, K=[sqrt(N)] или K=[N2/3]. Найдем область наименьшего объема (т.е. лежащую на наибольшей глубине дерева), содержащую не менее K точек. Если таких областей несколько, возьмем ту, в которой точек больше всего. После этого начнем делить ее пополам (спускаться по дереву дальше), каждый раз выбирая половину, в которой точек больше. Когда дошли до листа (например, до отдельной точки), выдаем его в качестве ответа.
Можно найти примеры, в которых этот алгоритм упустит главную область сгущения, а выберет вместо нее какую-нибудь локальную аномалию, или в качестве ответа выдаст точку, находящуюся далеко от центра сгущения, но можно ожидать, что в большинстве случаев его ответы будут более-менее адекватными. К сожалению, строить k-d дерево дорого, да и памяти оно занимает немало. Я готов выделить 8 байт на точку, но тратить больше было бы нежелательно.
Чтобы сэкономить время и память, я решил построить k-d дерево неявно.
Пусть у нас есть точка (x1,x2,…,xk). Запишем ее координаты в двоичной системе
x1=0.x11x12 x13…
x2=0.x21x22 x23…
x3=0.x31x32 x33…
…
И постоим 64-битное целое число
P= x11 x21… xk1 x12… xk2…
Это число представляет путь по дереву на глубину 64 уровня. Кроме того, оно позволяет восстановить координаты точки с требуемой точностью.
Если мы отсортируем коды P для всех точек по возрастанию, то получим точки в порядке обхода дерева. Для любых двух точек Pa и Pb легко найти самый глубокий узел, общими потомками которого они являются: достаточно найти старший ненулевой бит в числе Pa^Pb. И все необходимые фрагменты задачи на таком представлении решаются в несколько строчек:
Поиск уровня, на котором находится самая маленькая точка, содержащая не менее K точек:
ulong mask=ulong.MaxValue;
K--;
for(int i=K;i<m_np;i++) mask=Math.Min(mask,m_Arr[i]^m_Arr[i-K]);

Здесь m_Arr – массив кодов, m_np – число заполненных элементов в нем, mask – число, старший бит которого определяет искомый уровень.
Поиск области на найденном уровне, содержащей максимальное количество точек (мы знаем, что их не меньше K):
mask=InvMask(mask);
int ms=0,me=0;
for(int i=0;i<m_np-K;i++){
ulong a=m_Arr[i];
int h=i+K;
if(((a^m_Arr[h])&mask)==0){
while(hce) {
ms=h;
ce=InvMask(ce);
while(ms>0 && ((m_Arr[ms-1]^samp)&ce)==0) ms--;
} else {
me=h;
cb=InvMask(cb);
while(me<m_np-1 && ((m_Arr[me+1]^samp)&cb)==0) me++;
}

Здесь спуск идет не обязательно на один уровень.
Таким образом, задача поиска сгущения решается.
Вся программа уложилась в 100 строк, исходник (на C#) можно найти здесь.
Реальная задача имеет несколько более сложную формулировку. Диапазон изменения параметров заранее неизвестен, кроме того, измерений в выборке может быть больше запланированных 16 миллионов, причем «хорошие» измерения могут идти ближе к концу. Но небольшие модификации алгоритма позволяют справиться с этими проблемами. Например, если место в массиве кончилось, а точки продолжают поступать, мы можем отсортировать и проредить массив уже набранных кодов – на положение и качество сгущения это не повлияет.