Множество точек. Оболочка из ближайших соседей.

[B@R5uk]

Сразу оговорюсь, вопрос ниже НЕ относится к задачам классификации объектов (во всяком случае не на прямую) и НЕ является попыткой построить алгоритм, который быстро находит ближайших соседей в больших множествах (чем меня забрасывает поисковик, при попытке что-то нагуглить). Задача находится на стадии формализации, поэтому приветствуются не только рекомендации к решению, но даже и просто список ключевых слов для поиска статей/книжек.

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

Выбирается одна точка P, назовём её целевой, и требуется построить некоторое подмножество

$$N\subset S\setminus\{P\}$$

соседей этой точки P, которые с одной стороны достаточно блики к ней, а с другой — окружают её со всех сторон. Как формализовать это в формулах я ещё пытаюсь понять.

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

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

Я пытался было думать в сторону жадного алгоритма с отсеканием очередных добавляемых точек, если они попадают в конечно порождённый конус, образованный векторами направлений на точки, уже добавленные во множество N соседей. Несмотря на то, что на плоскости такой подход скорее всего будет работать, с ростом размерности возникнут вычислительные сложности в связи с комбинаторно быстро растущим числом таких выпуклых конусов, которые будут заполнять пространство вокруг целевой точки P.

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

Посоветуйте, пожалуйста, идеи и ключевые слова для поиска статей на тему задачи.

[B@R5uk]

Ну, хорошо. Попробую подойти к задаче с обратной стороны.

Поскольку для задачи Таммеса мне нужен граф соприкосновений (contact graph), то попробую его и строить. По определению, рёбра графа соприкосновений соединяют вершины, расстояния между которыми минимально возможные (и все равные). Поскольку я хочу строить приближение графа соприкосновений для произвольного случая положения точек, то мне требование минимальности нужно дополнить чем-то ещё. Плоский случая наиболее простой, поэтому буду его решать в первую очередь. На сфере (и на плоскости в ограниченной области) рёбра минимальной длины никогда не пересекаются. Это и положу в основу алгоритма построения графа. Получается два правила:

1. На каждом шаге в граф добавляется следующее по минимальности длины ребро, если ничего не мешает;
2. Каждое новое ребро не должно пересекаться с уже имеющимися.

Пока запрограммировал это в лоб. Картинки на выходе выглядят примерно так:

Получающийся граф является своего рода "триангуляцией" множества точек. Правда, в отличие от триангуляции Делоне, которая стремиться строить наиболее "толстые" треугольники, эта стремится к наиболее коротким рёбрам. Оболочка соседей для каждой точки — это точки, имеющие с рассматриваемой общее ребро. Например, для точки 2 на рисунке ими являются точки 5, 4, 3, 12, 17 и 18. Изначальная идея заключалась в том, чтобы для каждой точки строить оболочку соседей, а потом, на основе этой информации, строить весь граф. Предположительно, оболочка будет содержать значительно меньше точек, чем всё множество (хотя, конечно, можно привести вырожденные случаи), поэтому можно ожидать время работы алгоритма порядка ${}_{{}_.}O\lbig(n^2\rbig)$ от числа точек n.

У текущего алгоритма (реализованного в лоб) сложность 4-ой степени по времени и 2-ой по памяти. В нём, для каждого добавленного ребра, я отмечаю все пары точек, видимость которых друг из друга преграждается этим ребром (что блокирует последующее добавление пересекающихся рёбер). Мера расстояния между точками используется только для определения порядка, в котором добавляются рёбра. Для определения видимости точек используются только вектора направленности из точки на точку. По этой причине алгоритм без труда обобщается на пространства произвольной кривизны (просто оперируем с векторами направленности в двумерном линейном касательном пространстве к исходному пространству в точке наблюдения). Четвёртая степень сложности возникает из-за того, что общее число рёбер порядка квадрата (хотя конечное число рёбер в графе может быть заметно меньше), и для каждого добавляемого ребра в худшем случае прямая, содержащая это ребро, разбивает множество точек пополам (с одной стороны от прямой и с другой), а нужно проверить каждую пару по одной с каждой стороны.

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

Может быть существует какой-нибудь другой критерий, эквивалентный не пересекаемости рёбер в двумерном случае, но более легко считаемый в произвольном случае, например, с помощью средств линейной алгебры?

Код программы прилагаю. Файл test_contact_graph.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160335.html
%   B@R5uk, 2025, May

clc
clearvars
format compact
%format long
format short

%   Число точек
num_points = 30;
%   Минимально допустимое расстояние между точками
thresh_dist = 0.05;
%   Порог отсечки на вероятность добавления очередной точки
thresh_area = 0.1;

%   Функция поворота двумерного вектора на 90 градусов против часовой стрелки
func_rot90 = @(v) v ([2 1]) .* [-1 1];
%   Функция дробной части числа
func_frac = @(x) x - floor (x);

%   Задание конфигурации точек
switch 3
    case 1
        coords = test_contact_graph_pgen (num_points, thresh_dist, thresh_area);
        num_points = size (coords, 1);
    case 2
        coords = test_contact_graph_pgen (floor (0.8 * num_points), thresh_dist, thresh_area);
        coords = (0.6 + (sqrt (2) - 0.6) * coords (:, 1)) * [1 1] .* [cos(0.5 * pi * coords (:, 2)) sin(0.5 * pi * coords (:, 2))];
        coords = coords (1 > coords (:, 1) & 1 > coords (:, 2), :);
        coords = [0.2 * test_contact_graph_pgen(num_points - size (coords, 1), 2 * thresh_dist, thresh_area); coords];
        num_points = size (coords, 1);
    case 3
        golden = (sqrt (5) + 1) / 2;
        coords = func_frac ((0.5 : num_points)' * golden);
        %coords = func_frac ((1 : num_points)' * sqrt (2));
        coords = [coords (0.5 : num_points)' / num_points];
end

%   Расстояние между точками
dists = repmat (permute (coords, [3, 1, 2]), num_points, 1);
dists = sqrt (sum ((permute (dists, [2, 1, 3]) - dists) .^ 2, 3));
dists = dists + diag (inf (num_points, 1));

%   Количество пар точек
num_pairs = num_points * (num_points - 1) / 2;
%   Таблица информации о потенциальных рёбрах: индексы вершин, длина, флаг наличия
edges = zeros (num_pairs, 4);
m = 1;
for k = 1 : num_points - 1
    for l = k + 1 : num_points
        edges (m, :) = [k, l, dists(k, l) 0];
        m = m + 1;
    end
end
%   Сортировка рёбер в порядке возрастания длины
edges = sortrows (edges, 3);
%   Таблица связности вершин рёбрами
connectivity = zeros (num_points);
%   Таблица видимости пар вершин (на пути нет ребра)
visibility = ones (num_points);
%   Количество рёбер
num_edges = 0;
%   Количество потенциально невидимых пар
num_vis_count = 0;
%   Перебираем все потенциальные рёбра
for k = 1 : num_pairs
    %   Флаг наличия ребра (сбрасывается, если оно неудачное)
    flag = 1;
    %   Индексы соединяемых вершин
    vrtx_1 = edges (k, 1);
    vrtx_2 = edges (k, 2);
    %   Если между вершинами нет видимости
    if 0 == visibility (vrtx_1, vrtx_2)
        %   Переход к следующему ребру
        continue
    end
    %   Увеличиваем счётчик рёбер
    num_edges = num_edges + 1;
    %   Добавляем индексы рёбер в таблицу связности вершин
    connectivity (vrtx_1, vrtx_2) = num_edges;
    connectivity (vrtx_2, vrtx_1) = num_edges;
    %   Отмечаем наличие ребра в таблице информации о рёбрах
    edges (k, 4) = 1;
    %   Разбиваем множество точек на две половины, расположенные по разные
    %   стороны прямой 0, содержащей ребро. Для этого находим вектор,
    %   ортогональный прямой 0, и считаем скалярное произведение с ним
    evect = coords (vrtx_2, :) - coords (vrtx_1, :);
    nvect_0 = func_rot90 (evect);
    coords_shifted = coords - repmat (coords (vrtx_1, :), num_points, 1);
    tmp = coords_shifted * nvect_0';
    indices_left = find (0 > tmp);
    indices_rght = find (0 < tmp);
    %   Примитивная оптимизация быстродействия, работающая только в Matlab'е
    if numel (indices_left) > numel (indices_rght)
        tmp = indices_left;
        indices_left = indices_rght;
        indices_rght = tmp;
    end
    %   Для каждой точки одной части
    for l = 1 : numel (indices_left)
        %   Находим вектор, ортогональный прямой 1, проходящую через точку
        %   и первый конец отрезка
        nvect_1 = func_rot90 (coords_shifted (indices_left (l), :));
        %   Заслоняемые отрезком точки второй части лежат по ту же сторону
        %   от прямой 1, что и вектор отрезка
        flags = 0 < coords_shifted (indices_rght, :) * nvect_1' * (evect * nvect_1');
        %   Находим вектор, ортогональный прямой 2, проходящую через точку
        %   и второй конец отрезка
        nvect_2 = func_rot90 (coords (indices_left (l), :) - coords (vrtx_2, :));
        %   Заслоняемые отрезком точки второй части лежат по ту же сторону
        %   от прямой 2, что и вектор отрезка с обратным знаком
        tmp = coords - repmat (coords (vrtx_2, :), num_points, 1);
        flags = flags .* (0 > tmp (indices_rght, :) * nvect_2' * (evect * nvect_2'));
        indices_eclipsed = indices_rght (0 ~= flags);
        %   Отмечаем невидимые пары вершин
        visibility (indices_eclipsed, indices_left (l)) = 0;
        visibility (indices_left (l), indices_eclipsed) = 0;
    end
    num_vis_count = num_vis_count + numel (indices_left) * numel (indices_rght);
end
disp ([num_points num_pairs num_edges num_vis_count])

%subplot (111)
plot (coords (:, 1), coords (:, 2), 'k.')
hold on
for k = 1 : num_pairs
    if 1 == edges (k, 4)
        plot (coords (edges (k, 1 : 2), 1), coords (edges (k, 1 : 2), 2))
    end
end
hold off
for k = 1 : num_points
    text (0.02 + coords (k, 1), coords (k, 2), num2str (k), 'Color', 'r')
end
axis equal
xlim ([0 1])
ylim ([0 1])
grid on
title (num_points)
 

Файл test_contact_graph_pgen.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160335.html
%   B@R5uk, 2025, May

function coords = test_contact_graph_pgen (num_points, thresh_dist, thresh_area)

%   Массив для координат точек
coords = zeros (num_points, 2);
%   Оставшееся свободная площадь определяет вероятность добавления точки
free_area = 1;
%   Координаты первой точки
tmp = rand (1, 2);
%   Индекс точки
k = 0;
while 1
    k = k + 1;
    %   Добавляем новые координаты в список
    coords (k, :) = tmp;
    %   Уменьшаем оставшуюся площадь
    free_area = free_area - pi * thresh_dist ^ 2;
    %   Если вероятность добавления следующей точки слишком мала
    if thresh_area > free_area
        %   Заканчиваем, обрезая список координат
        coords = coords (1 : k, :);
        break
    end
    %   Если желаемое число точек достигнуто, то останавливаемся
    if num_points == k
        break
    end
    %   Генерируем координаты новой точки до тех пор, пока она не окажется
    %   на достаточном удалении от всех остальных уже добавленных
    while 1
        tmp = rand (1, 2);
        if thresh_dist < min (sqrt (sum ((coords (1 : k, :) - repmat (tmp, k, 1)) .^ 2, 2)))
            break
        end
    end
end

end
 

-- 07.05.2025, 15:54 --

Картинка с заполнением квадрата точками через золотое сечение:

(Оффтоп)

[B@R5uk]

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

Файл test_contact_graph_2.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160335.html
%   B@R5uk, 2025, May

clc
clearvars
format compact
%format long
format short

%   Число точек
num_points = 10;
%   Минимально допустимое расстояние между точками
thresh_dist = 0.02;
%   Порог отсечки на число попыток добавления очередной точки
thresh_area = 100;

%   Функция поворота двумерного вектора на 90 градусов против часовой стрелки
func_rot90 = @(v) v ([2 1]) .* [-1 1];
%   Функция дробной части числа
func_frac = @(x) x - floor (x);

%   Задание конфигурации точек
switch 1
    case 1
        %   Равномерное заполнение квадрата
        coords = test_contact_graph_pgen (num_points, thresh_dist, thresh_area);
        num_points = size (coords, 1);
    case 2
        %   Сегмент кольца (1/4) и равномерное заполнение угла
        coords = test_contact_graph_pgen (floor (0.8 * num_points), thresh_dist, thresh_area);
        coords = (0.6 + (sqrt (2) - 0.6) * coords (:, 1)) * [1 1] .* [cos(0.5 * pi * coords (:, 2)) sin(0.5 * pi * coords (:, 2))];
        coords = coords (1 > coords (:, 1) & 1 > coords (:, 2), :);
        coords = [0.2 * test_contact_graph_pgen(num_points - size (coords, 1), 2 * thresh_dist, thresh_area); coords];
        num_points = size (coords, 1);
    case 3
        %   Заполнение через иррациональный псевдогенератор случайности
        golden = (sqrt (5) + 1) / 2;
        coords = func_frac ((0.5 : num_points)' * golden);
        %coords = func_frac ((1 : num_points)' * sqrt (2));
        coords = [coords (0.5 : num_points)' / num_points];
end

%   Расстояние между точками
dists = repmat (permute (coords, [3, 1, 2]), num_points, 1);
dists = sqrt (sum ((permute (dists, [2, 1, 3]) - dists) .^ 2, 3));
dists = dists + diag (inf (num_points, 1));
%   Количество пар точек
num_pairs = num_points * (num_points - 1) / 2;
%   Таблица информации о потенциальных рёбрах: индексы вершин, длина, флаг наличия
edges = zeros (num_pairs, 4);
m = 1;
for k = 1 : num_points - 1
    for l = k + 1 : num_points
        edges (m, :) = [k, l, dists(k, l) 0];
        m = m + 1;
    end
end
%   Сортировка рёбер в порядке возрастания длины
edges = sortrows (edges, 3);
%   Таблица связности вершин рёбрами
connectivity = zeros (num_points);
%   Таблица видимости пар вершин (на пути нет ребра)
visibility = ones (num_points);
%   Количество рёбер
num_edges = 0;
%   Количество потенциально невидимых пар
num_vis_count = 0;
%   Перебираем все потенциальные рёбра
for k = 1 : num_pairs
    %   Индексы соединяемых вершин
    vrtx_1 = edges (k, 1);
    vrtx_2 = edges (k, 2);
    %   Если между вершинами нет видимости
    if 0 == visibility (vrtx_1, vrtx_2)
        %   Переход к следующему ребру
        continue
    end
    %   Увеличиваем счётчик рёбер
    num_edges = num_edges + 1;
    %   Добавляем индексы рёбер в таблицу связности вершин
    connectivity (vrtx_1, vrtx_2) = num_edges;
    connectivity (vrtx_2, vrtx_1) = num_edges;
    %   Отмечаем наличие ребра в таблице информации о рёбрах
    edges (k, 4) = 1;
    %   Разбиваем множество точек на две половины, расположенные по разные
    %   стороны прямой 0, содержащей ребро. Для этого находим вектор,
    %   ортогональный прямой 0, и считаем скалярное произведение с ним
    nvect_0 = func_rot90 (coords (vrtx_2, :) - coords (vrtx_1, :));
    tmp = (coords - repmat (coords (vrtx_1, :), num_points, 1)) * nvect_0';
    %   Индексы точек по разные стороны от прямой 0
    indices_left = find (0 > tmp);
    indices_rght = find (0 < tmp);
    %   Исключение идексов точек, задающих прямую, если они оказались в
    %   списке из-за ограниченности машинной точности
    indices_left = indices_left (vrtx_1 ~= indices_left & vrtx_2 ~= indices_left);
    indices_rght = indices_rght (vrtx_1 ~= indices_rght & vrtx_2 ~= indices_rght);
    %   Примитивная оптимизация быстродействия, работающая только в Matlab'е
    if numel (indices_left) > numel (indices_rght)
        tmp = indices_left;
        indices_left = indices_rght;
        indices_rght = tmp;
    end
    %   Для каждой точки одной части
    for l = 1 : numel (indices_left)
        %   Рассматриваем конечно порождённый конус (угол) с вершиной в
        %   точке и векторами направлений на концы добавляемого отрезка
        pvect = coords (indices_left (l), :);
        base = coords ([vrtx_1 vrtx_2], :) - repmat (pvect, 2, 1);
        %   Проверка регулярности базиса, то есть, что точка наблюдения
        %   не лежит на прямой с заграждающим отрезком
        [u s v] = svd (base);
        %   Если это так, то пропускаем точку (по хорошему,
        %   надо проверить, какие ещё точки лежат на прямой)
        if s (end) / s (1) < 2 * eps
            continue
        end
        %   Обратная матрица базиса
        base_inv = v * s ^ -1 * u';
        %   Координаты разложения по базису направлений на точки второй части
        tmp = (coords (indices_rght, :) - repmat (pvect, numel (indices_rght), 1)) * base_inv;
        %   Заслоняемые точки имеют только неотрицательные координаты
        indices_eclipsed = indices_rght (0 ~= prod (double (0 <= tmp), 2));
        %   Отмечаем невидимые пары вершин
        visibility (indices_eclipsed, indices_left (l)) = 0;
        visibility (indices_left (l), indices_eclipsed) = 0;
    end
    num_vis_count = num_vis_count + numel (indices_left) * numel (indices_rght);
end
%   Число точек, число пар точек, число рёбер, трудоёмкость алгоритма
disp ([num_points num_pairs num_edges num_vis_count])
%   Множители числа рёбер и трудоёмкости
%   Для плоскости число рёбер в пределе бесконечного равномерного
%   заполнения плоскости равно 3 на число точек. Трудоёмкость оказывается
%   почти пропорциональной кубу числа точек.
disp ([num_edges / num_points num_vis_count / num_points ^ 3])

%subplot (111)
%   Отображение точек
plot (coords (:, 1), coords (:, 2), 'k.')
hold on
%   Отображение рёбер
for k = 1 : num_pairs
    if 1 == edges (k, 4)
        plot (coords (edges (k, 1 : 2), 1), coords (edges (k, 1 : 2), 2))
    end
end
hold off
%   Отображение этикеток к точкам
for k = 1 : num_points
    text (0.02 + coords (k, 1), coords (k, 2), num2str (k), 'Color', 'r')
end
axis equal
xlim ([0 1])
ylim ([0 1])
grid on
title (num_points)
 

Файл test_contact_graph_pgen.m:

Используется синтаксис Matlab M

%   dxdy.ru/topic160335.html
%   B@R5uk, 2025, May

function coords = test_contact_graph_pgen (num_points, thresh_dist, thresh_fail)

%   Массив для координат точек
coords = zeros (num_points, 2);
%   Координаты первой точки
tmp = rand (1, 2);
%   Счётчик неудачных добавлений
num_fail = 0;
%   Индекс точки
k = 0;
%   Если найти новое положение для точки не удалось слишком много раз, то
%   скорее всего его не существует, поэтому останавливаемся
while thresh_fail > num_fail
    k = k + 1;
    %   Добавляем новые координаты в список
    coords (k, :) = tmp;
    %   Если желаемое число точек достигнуто, то останавливаемся
    if num_points == k
        break
    end
    %   Генерируем координаты новой точки до тех пор, пока она не окажется
    %   на достаточном удалении от всех остальных уже добавленных
    while thresh_fail > num_fail
        tmp = rand (1, 2);
        if thresh_dist < min (sqrt (sum ((coords (1 : k, :) - repmat (tmp, k, 1)) .^ 2, 2)))
            num_fail = 0;
            break
        end
        num_fail = num_fail + 1;
    end
end

end
 

(Оффтоп)

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

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

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

Выше я дал неверную оценку трудоёмкости алгоритма: я посчитал, что раз количество потенциальных рёбер имеет порядок квадрата точек (точнее, равно количеству пар точек), то и в конечном графе их количество тоже будет пропорционально квадрату. Эксперимент показывает, что это не так: на самом деле, на плоскости число рёбер пропорционально числу точек с коэффициентом 2...3. Откуда это число берётся и почему это так, станет очевидно, если рассмотреть предельный случай с плоскостью, полностью заполненной равноудалёнными точками. Другими словами, рассмотреть замощение плоскости равносторонними треугольниками, что эквивалентно задаче плотной упаковки кругов. Каждая вершина в этом замощении будет иметь валентность 6 (количество ближайших соседей), а поскольку каждое ребро соединяет 2 вершины, на каждую вершину в среднем придётся по 3 ребра. Отклонение (в меньшую сторону) среднего числа ребер на вершину от этого оптимального числа 3 в случае случайного заполнения квадрата точками можно рассматривать как некоторую меру неравномерности заполнения. В итоге получается, что алгоритм выше имеет кубическую сложность по числу точек, а не 4-й степени, как я предполагал.

По аналогии (не уверен, что это верно, но выглядит правдоподобно), количество рёбер на точку в трёхмерном пространстве будет связано с задачей плотной упаковки сфер. Каждая сфера будет иметь по 6 соседей в содержащем её слое и по 3 соседа в двух соседних слоях (вне зависимости от типа упаковки: ГЦК или ГПУ), то есть всего 12 соседей. Это означает, что оптимальным числом для среднего количество рёбер на вершину будет 6. Что будет происходить в четырёхмерном пространстве, а так же в пространствах высшей размерности я без понятия. Предположительно (опять же по аналогии с пространствами размерности 2 и 3) число рёбер будет тоже пропорционально числу вершин графа с неким целым коэффициентом. Известны точные решения для упаковок в пространствах размерности 8 и размерности 24 (решётки E8 и Лича, соответственно), но я не знаю, как из сказанного о них в вики можно получить число ближайших соседей. Если кто разбирается в вопросе и подскажет ответ, буду очень благодарен. В любом случае, трёхмерная реализация алгоритма тоже должна иметь кубическую сложность.

Небольшое уменьшение трудоёмкости алгоритма можно достигнуть, если перед проверкой "заслонённости" пары точек во внутреннем цикле добавить проверку того, что между ними уже пропала видимость за счёт ранее добавленного ребра/треугольника/симплекса. Поскольку одно сравнение проще, чем набор умножений, сложений и сравнений, то константа трудоёмкости заметно уменьшится, особенно для пространств большой размерности. Для увеличения эффективности алгоритма не на константу, а на порядки необходимо эксплуатировать тот факт, что пространство в процессе работы разбивается на симплексы, а из каждой точки, окружённой полным набором таких симплексов остальные точки не видны. Например, подход "разделяй и властвуй" будет пользоваться этим при объединении двух подмножеств: видимость будет достигаться только для граничных точек обоих множеств.

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

Вывод из всего выше сказанного для меня заключается в том, что для пространств более высокой размерности мне нужен другой, более простой алгоритм. Единственное строгое требование на него — это чтобы он в предельном случае равенства всех минимальных расстояний между точками давал граф из рёбер с этой минимальной длиной. Во всех остальных "неоптимальных" случаях он должен вести себя более-менее разумно, выбирая короткие рёбра со всех сторон вокруг каждой точки, но делая это максимально прозрачно и эффективно (что вполне может дать результаты, отличные от результатов алгоритма выше). Вопрос: как это формализовать?