Максиминное расположение точек на сфере

[мат-ламер]

B@R5uk в сообщении #1682004 писал(а):

Возможно, это и делается в статье китайцев
, что я приводил ранее. Надо уже в конце-концов прочитать её и вникнуть в суть дела там.

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

мат-ламер в сообщении #1681275 писал(а):

Думаю попробовать такой подход.

У меня руки до такого подхода пока не дошли. Китайцы считают, что этот метод неудовлетворительный

[B@R5uk]

Интересно, есть ли какая-нибудь обзорная инфа по тому, какие методы и на сколько удачно применяются для нахождения локальных минимумов задачи Таммеса? Пока без относительно поиска глобального. Чтобы если и пытаться придумать что-нибудь новое, то хотя бы знать, что уже опробовано и на сколько высока планка для перепрыгивания. Китайцы, судя по параграфу:

Цитата:

The subproblem dened by Eqs. (5) and (6) is still a constrained optimization problem which is not easy to handle by popular local optimization methods like the LBFGS method. To perform the local optimization, we convert the problem further to an unconstrained optimization problem by using the spherical coordinate transformation of points on S²

тоже используют что-то стороннее и не заморачиваются особо. И дальше, на странице 13 они упоминают алгоритм LBFGS и дают ссылку на статью On the limited memory BFGS method for large scale optimization авторов Dong C. Liu и Jorge Nocedal. Непонятно, правда, они самостоятельно его пишут или пользуются чей-то сторонней библиотекой, как, например, Laszlo Hars пользовался LD_SLSQP из библиотеки NLopt в своей работе.

[zgemm]

B@R5uk в сообщении #1682176 писал(а):

И дальше, на странице 13 они упоминают алгоритм LBFGS и дают ссылку на статью On the limited memory BFGS method for large scale optimization авторов Dong C. Liu и Jorge Nocedal. Непонятно, правда, они самостоятельно его пишут или пользуются чей-то сторонней библиотекой, как, например, Laszlo Hars пользовался LD_SLSQP из библиотеки NLopt в своей работе.

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

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

[B@R5uk]

Вот есть теорема (лемма?) Фейеш-Тота, которая даёт верхнюю границу для оптимальной длины (угла) задачи Таммеса:

$$\cos\beta\ge\frac{1}{2}\left(\ctg^2\left(\frac{n}{n-2}\frac{\pi}{6}\right)-1\right)$$

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

Самому в голову приходит только идея о расположении на сфере на равных угловых расстояниях порядка $\sqrt{n}$ окружностей в параллельных плоскостях и упаковка точек на эти окружности с расстоянием между точками не более расстояния между окружностями. Но это больше получается примитивный оптимизационный алгоритм с обращением сложной суммы с округлениями, чем удобосчитаемая формула для быстрой оценки. Что-то как-то даже не верится, что этот вопрос никем раньше не рассматривался.

-- 19.04.2025, 13:38 --

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

[Sender]

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

[мат-ламер]

B@R5uk
Из асимптотического поведения конкретно для этой задачи вряд ли возможно получить достаточно точные и близкие друг к другу оценки. Дело в том, что если вы построите график зависимости оцениваемой величины (косинус угла) от размерности задачи, то на нём будут периодические всплески для некоторых значений $n$ , для которых будет наблюдаться повышенная симметрия. Поэтому между оценкой сверху и оценкой снизу будет зазор. Однако, может вам удастся тут самому что-то оценить? Но тема эта трудная.

[B@R5uk]

Edith Mooers в статье Tammes's Problem в четвёртом параграфе пишет:

Цитата:

This problem is quite a bit harder than it would appear to be at first glance. Induction seems out of the question since, in general, the best set of points for n bears no relation to the best set of points for (n+1). Although good upper and lower bounds exist for the distances between points, actual solutions vary greatly as to how close they come to reaching these bounds. As noted above, the problem has been solved only for very low values of n. If there exists a good technique for showing that a given vertex set is optimal for large values of n, no one has yet come across it.

Выделил упоминание жирным. В дальнейшем в статье он раскрывает вопрос оценки сверху, приводя теорему Фейеш-Тота, а про оценку с низу даже не заикается. Щас вот закинул себе Thomas Ericson, Victor Zinoviev - Codes on Euclidean spheres (Elsevier, 2001) почитать ещё разок, может там что написано по этому поводу, в первый раз не присматривался к этому вопросу.

Самостоятельное решение, я уже предлагал выше. Зададимся неким сферическим расстоянием d. Зенитный угол в пи радиан делится этим расстоянием на m параллелей, равноудалённых на расстояние d:

$$m=\left\lfloor\frac{\pi}{d}\right\rfloor$$

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

$$\theta(m,k)=\frac{\pi-(m+1)d}{2}+kd,\qquad\qquad 1\le k\le m$$

а радиусы будут

$$r(m,k)=\cos\theta(m,k)=\cos\left(\frac{\pi-(m+1)d}{2}+kd\right)$$

Теперь бы было бы не плохо прикинуть количество точек, которые можно разместить на каждой из параллелей, поделив на угловое расстояние d их длину:

$$l(m,k)=2\pi r(m,k)$$

Однако, расчёт в лоб даст неправильный результат, потому что угловое сферическое расстояние между точками на параллели, отличной от экватора, короче дуги параллели между этими точками. Надо внести поправку в расстояние d для каждой параллели:

$$\left(\begin{matrix}r\\0\\\sqrt{1-r^2}\end{matrix}\right)^T\left(\begin{matrix}r\cos d'\\r\sin d'\\\sqrt{1-r^2}\end{matrix}\right)=r^2\cos d'+1-r^2=\cos d$$

Откуда поправленное расстояние вдоль каждой параллели

$$d'(m,k)=\arccos\left(1-\frac{1-\cos d}{r^2(m,k)}\right)$$

Количество точек n, умещающихся на сфере с расстоянием не менее d друг от друга, теперь можно прикинуть как сумму количества точек, умещающихся на каждой параллели:

$$n\ge\sum\limits_{k=1}^{m}\left\lfloor\frac{l(m,k)}{d'(m,k)}\right\rfloor$$

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

Не представляю как это можно было бы сделать в этом подходе, не сломав моск. Наверное, надо придумать другой подход.

[SergeyGubanov]

B@R5uk в сообщении #1681694 писал(а):

$$\theta=\arccos\left(1-\frac{2n-1}{N}\right),\qquad\varphi=\pi\left(\sqrt{5}\pm1\right)n,\qquad 1\le n\le N$$

$$\begin{cases}x=\sin\theta\cos\varphi\\y=\sin\theta\sin\varphi\\z=\cos\theta\end{cases}$$

Забавно. Попробовал другие иррациональные числа. Визуально $\sqrt{6}$ примерно так же хорош как и $\sqrt{5}$. С несколькими другими наугад взятыми числами равномерного распределения не получается, а получаются красивые узоры.

Правильно ли я понимаю, что для трёхмерной сферы в четырёхмерном пространстве это можно пытаться обобщить как-то вот так:

$$\begin{cases} x_1 = \sin\theta_n \cos\left( \pi \left( \sqrt{5} - 1\right) n \right) \\ x_2 = \sin\theta_n \sin\left( \pi \left( \sqrt{5} - 1\right) n \right) \\ x_3 = \cos\theta_n \cos\left( \pi \left( \sqrt{6} - 1\right) n \right) \\ x_4 = \cos\theta_n \sin\left( \pi \left( \sqrt{6} - 1\right) n \right) \end{cases}$$

Или я что-то не правильно понял в этом методе?

[SergeyGubanov]

(Оффтоп)

Еще порисовал. Корень из 6 не годится, корень из 2 лучше:

$$\begin{cases} x_1 = \sin\theta_n \cos\left( \pi \left( \sqrt{5} - 1\right) n \right) \\ x_2 = \sin\theta_n \sin\left( \pi \left( \sqrt{5} - 1\right) n \right) \\ x_3 = \cos\theta_n \cos\left( \pi \sqrt{2} n \right) \\ x_4 = \cos\theta_n \sin\left( \pi \sqrt{2} n \right) \end{cases}$$

[B@R5uk]

SergeyGubanov в сообщении #1684919 писал(а):

Правильно ли я понимаю, что для трёхмерной сферы в четырёхмерном пространстве это можно пытаться обобщить как-то вот так... Или я что-то не правильно понял в этом методе?

По ссылке на стэковерфлоу в ответе начинающимся "The golden spiral method" есть объяснение откуда берутся формулы для θ в плоском (раздел Spacing the points radially) и сферическом (раздел Now doing the sunflower on a sphere) случаях. Основная идея: воспользоваться формулами теории вероятности для преобразования равномерного на квадрате $[0,1]^2$ распределения в равномерное на круге или сфере распределение. Это нужно, чтобы спираль с точками равномерно заполняла двумерные области: круг и сферу. При этом плотность распределение от угла φ не зависит (является равномерным по нему), поэтому для этого угла нужно всего лишь подобрать такой коэффициент шага, чтобы точки при набрасывании на спираль в целом не выстраивались в ряды. По какой-то волшебной причине, этим коэффициентом является золотое сечение. В начале ответа по ссылке приводится попытка объяснить, почему так получается (типа "золотое сечение является наиболее иррациональным числом из всех возможных"), но я на объяснение не купился, поэтому пока остаюсь в неведении о причинах. Поскольку поверхностная плотность распределения точек от угла φ не зависит, то его равномерность определяется только формулой для угла θ.

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

Как вы вообще визуализируете получающийся нелинейный четырёхмерный объект, чтобы можно было оценивать расстояния между точками и убеждаться в их приблизительной одинаковости, а так же в том, что точки не укладываются на линии?

[B@R5uk]

B@R5uk в сообщении #1684958 писал(а):

Как действовать в случае трёхмерного множества точек я не представляю.

ОК, вру. Идея появилась. Будем действовать по аналогии с двумерным случаем. Но сначала надо понять, что в нём на самом деле происходит. У нас имеется квадрат $[0,1]^2$, который по двум формулам заполняется точками. Вдоль одного направления (соответствующего углу φ) точки в квадрате размещаются с фиксированным шагом, при этом при выходе за квадрат координата отбрасывается назад. Фактически это даже не квадрат, а внешняя сторона цилиндра: по углу φ он цикличен. При отображении в круг, нижняя граница цилиндра стягивается в точку, и вся внешняя сторона сплющивается в круг. При этом заполнение по второй координате выбрано так, чтобы в результате этого сплющивания плотность распределения точек получается равномерной.

В случае заполнения сферы обе границы — верхняя и нижняя — внешней стороны цилиндра стягиваются в точку (каждая — в свою), которые образуют полюса, а сам цилиндр надувается до сферы. Формула для угла θ гарантирует, что после этого преобразования распределение точек на сфере будет равномерным.

Если действовать по аналогии, то для трёхмерной сферы, вложенной в четырёхмерное пространство, нам придётся в качестве начальной области взять куб $[0,1]^3$ и зациклить его по двум координатам. Одна из этих зацикленных координат (пусть это будет угол φ для удобства) будет всё так же меняться линейно с номером точки, а другая (а так же третья, не зацикленная) — нелинейно. Формулы для них придётся рассчитать по теории вероятности. При этом надо будет ввести два (а не один как раньше) коэффициента, чтобы гарантировать разброс точек. Они будут добавляться после нелинейности (для второй координаты), но до зацикливания (для первой и второй координаты). Возможно, второй коэффициент должен быть просто квадратом золотого сечения.

То есть должно получаться что-то в этом духе:

$$\begin{cases}x=\sin\theta_2\sin\theta_1\cos\varphi\\y=\sin\theta_2\sin\theta_1\sin\varphi\\z=\sin\theta_2\cos\theta_1\\v=\cos\theta_2\end{cases}$$

$$\theta_2=f_2\left(\frac{2n-1}{2N}\right),\qquad\theta_1=\pi\left\{\frac{\gamma}{\pi}\left[\sqrt{N}\right]f_1\left(\frac{2n-1}{2N}\right)\right\},\qquad\varphi=2\pi\gamma^2n,\qquad 1\le n\le N$$

$$\gamma=\frac{\sqrt{5}+1}{2}$$

Фигурные скобки в формуле для $\theta_1$ означают взятие дробной части числа, а квадратные — округление до ближайшего целого. Округлённый квадратный корень из числа точек N даёт периодичность средней координаты значительно большую левой координаты (у которой она равна 1) и значительно меньшую правой координаты (у которой она порядка N). Это должно позволить периодически заполнить весь куб точками. Наличие золотого сечения в формулах, по идее, должно предотвратить выстраивание точек в ряды. Монотонные функции f с индексом рассчитываются по теории вероятности таким образом, чтобы распределение

$$\theta_2=f_2\left(\xi_2\right),\qquad\theta_1=f_1\left(\xi_1\right),\qquad\varphi=2\pi\xi_0,\qquad\xi_k\sim U[0,1],\qquad k=0,1,2$$

давало равномерное распределение на гиперсфере. Формулы сам не пробовал, результат не гарантирован.

SergeyGubanov, вопрос визуализации, если вы что-то такое действительно проделываете, очень актуален.

[SergeyGubanov]

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

Что касается меня, то я лишь убедился в том, что для двумерной сферы $\sqrt{2}$ визуально даёт картинку не хуже чем $\sqrt{5} - 1$ (несколько других наугад выбранных чисел давали визуально легко наблюдаемые регулярные узоры):

R1 = Table[{Sin[ArcCos[1 - (2*n - 1)/nmax]]*Cos[Pi*(Sqrt[5] - 1)*n], Sin[ArcCos[1 - (2*n - 1)/nmax]]*Sin[Pi*(Sqrt[5] - 1)*n], Cos[ArcCos[1 - (2*n - 1)/nmax]]}, {n, 1, nmax}];
R2 = Table[{Sin[ArcCos[1 - (2*n - 1)/nmax]]*Cos[Pi*Sqrt[2]*n], Sin[ArcCos[1 - (2*n - 1)/nmax]]*Sin[Pi*Sqrt[2]*n], Cos[ArcCos[1 - (2*n - 1)/nmax]]}, {n, 1, nmax}]; 
ListPointPlot3D[R1, PlotTheme -> "Detailed", PlotRange -> All, BoxRatios -> Automatic]
ListPointPlot3D[R2, PlotTheme -> "Detailed", PlotRange -> All, BoxRatios -> Automatic]

И в том, что псевдо-случайные числа порождаемые через $\sqrt{2}$ не очень сильно коррелируют с псевдо-случайными числами порождаемыми через $\sqrt{5} - 1$, это означает, что можно будет двумерную поверхность покрыть более-менее (об идеальной равномерности речи, к сожалению, не идёт):

L = Table[{ArcCos[Cos[Pi*(Sqrt[5] - 1)*n]], ArcCos[Cos[Pi*Sqrt[2]*n]]}, {n, 100000, 101000}]; 
ListPlot[L, AspectRatio -> Automatic]

[B@R5uk]

SergeyGubanov в сообщении #1684973 писал(а):

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

Жаль, было бы здорово найти что-нибудь работоспособное на будущее.

SergeyGubanov в сообщении #1684973 писал(а):

...я лишь убедился в том, что для двумерной сферы $\sqrt{2}$ визуально даёт картинку не хуже чем $\sqrt{5} - 1$... И в том, что псевдо-случайные числа порождаемые через $\sqrt{2}$ не очень сильно коррелируют с псевдо-случайными числами порождаемыми через $\sqrt{5} - 1$, это означает, что можно будет двумерную поверхность покрыть более-менее...

Да, вы правы. Спасибо. Я тут поигрался с равномерным заполнением куба, проверяя результат тремя проекциями на координатные плоскости и рассматривая 3d-картинку. Ваш подход работает значительно лучше того, что я выше предлагал с округлённым корнем. Возведение золотого сечения в квадрат — это вообще просто увеличение числа на единицу, и поэтому никак не изменяет картинку псевдослучайного генератора. Поэтому в трех измерениях без корня из двух или другого подобного основания для генератора не обойтись. Интересно, есть ли какая-то теория, которая объясняет, почему эти иррациональности работают хорошо?

(Оффтоп)

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

Программа с моими экспериментами:

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

clc
clearvars
format compact
%format long
format short

num_points = 258;

frac = @(x) x - floor (x);

g = 0.5 * (sqrt (5) + 1);
h = sqrt (2);

k1 = g * num_points;
k2 = h * num_points;
%k2 = h * round (sqrt (num_points));

nn = 1 : num_points;
nn = (nn - 0.5) / num_points;

coords = [
    nn
    frac(k1 * nn)
    frac(k2 * nn)
];
coords = coords';

subplot (221)
plot (coords (:, 1), coords (:, 2), 'k.')
axis square
title ('OXY projection')
xlabel ('X')
ylabel ('Y')

subplot (222)
plot (coords (:, 3), coords (:, 2), 'k.')
axis square
title ('OZY projection')
xlabel ('Z')
ylabel ('Y')

subplot (223)
plot (coords (:, 1), coords (:, 3), 'k.')
axis square
title ('OXZ projection')
xlabel ('X')
ylabel ('Z')

subplot (224)
plot3 (coords (:, 1), coords (:, 2), coords (:, 3), 'k.')
axis vis3d
xlabel ('X')
ylabel ('Y')
zlabel ('Z')
 

[dgwuqtj]

B@R5uk в сообщении #1684987 писал(а):

Интересно, есть ли какая-то теория, которая объясняет, почему эти иррациональности работают хорошо?

По идее это диофантова аппроксимация, раздел теории чисел. В двумерном случае надо брать иррациональность, у которой при разложении в цепную дробь как можно меньшие коэффициенты, то есть лучше всего $\frac {\sqrt 5 + 1} 2 = 1 + \frac 1 {1 + \frac 1{1 + \ldots}}$.

[B@R5uk]

dgwuqtj, спасибо! А в трёхмерном что лучше? А то, вон, на картинке под катом точки в параллельные плоскости расслаиваются.