Журнал вычислительной математики и математической физики, 2020, T. 60, № 5, стр. 815-827

1. ВВЕДЕНИЕ

Пусть задано натуральное $N \geqslant 2$. Задача Фекете порядка $N$ формулируется следующим образом [1]. Найти множество точек ${\mathbf{x}} = {\text{\{ }}{{x}_{1}},\; \ldots ,\;{{x}_{N}}{\text{\} }} \in {{\mathbb{R}}^{{3N}}}$ с ограничением для каждой точки ${{x}_{i}}$ на единичной сфере $S_{i}^{2} = {{S}^{2}} = \{ x \in {{\mathbb{R}}^{3}}:{{\left\| x \right\|}_{2}} = 1\} $ такое, что оно минимизирует потенциальную энергию системы:

При $s = 1$ получается кулоновский потенциал и соответственно задача Томсона. Последняя была поставлена Дж.Дж. Томсоном еще в 1904 г. [2], однако, до сих пор остается нерешенной для произвольного $N$. Задача Фекете привлекла особое внимание математиков, потому что она появляется в списке Смейла [3] – нерешенных задач для 21-го века. Кроме того, задача имеет применения в химии и биологии [4]–[7], а также в силу своей сложности стала эталоном для тестирования численных методов оптимизации [8]. Стоит отметить, что в современных приложениях, например, в моделировании макромолекул [9] также применяются численные методы для минимизации потенциала. Конечно, часто потенциалы имеют более сложный вид и решается задача безусловной оптимизации, однако эти задачи имеют структуру, похожую на задачу Фекете в том смысле, что потенциал является суммой некоторых симметричных функций $K:\mathcal{X} \times \mathcal{X} \to \mathbb{R}$, что позволяет эффективно пересчитывать значение потенциала при смещении одной точки. В контексте вышеупомянутых приложений стоит обозначить, что целью работы является именно локальная минимзация функции (1), и в дальнейшем под ${\mathbf{x}}{\text{*}}$ понимается локальный минимум, ближайший к начальному положению ${{{\mathbf{x}}}^{0}}$.

Учитывая широкое применение, задача Фекете была хорошо изучена [10]–[12]. Известны аналитические решения при $s = 1$ для $N = 2{\kern 1pt} - {\kern 1pt} 6$ и $12$. Для различных $s$ и больших $N$ получены решения с помощью численных методов [13]. В данной статье мы остановимся на обсуждении метода проективного покоординатного градиентного спуска. Рассматриваемый метод использует быстрый счет функции и метод дробления шага. Также предлагается алгоритм второго порядка для решения данной задачи.

Определение точек Фекете есть невыпуклая задача оптимизации с нелинейными ограничениями. Следует отметить, что здесь можно легко избавиться от ограничений, параметризуя каждую точку на сфере двумя углами [14], [15], и решать задачу нелинейной оптимизации без ограничений от $2N$ переменных вместо $3N$. Однако на практике применение градиентных методов в таких координатах существенно замедляет счет оптимизируемой функции и градиента из-за относительно медленного вычисления тригонометрических функций. Поэтому остановимся подробнее на решении задачи в исходных координатах. Можно применять градиентные методы, рассматривая зависимость энергии от всех $3N$ переменных, и на каждой итерации перемещать все точки одновременно. Сравнительный анализ в разд. 5 показывает, что покоординатный метод оказывается вполне конкурентоспособным. Более того, появляется возможность использовать покоординатную версию метода второго порядка, что оказывается затруднительным для полноградиентной версии (так как это требует подсчета и обращения матрицы размера $3N \times 3N$). Поэтому интересно рассматривать проективный покоординатный спуск в исходных координатах. Еще одним аргументом в пользу изучения покоординатного спуска на примере задачи вида (1) является то, что впоследствии это дает возможность эффективно распараллелить вычисления. Если удается эффективно организовать приближенный счет функции и его градиента по ближайшим соседям, то можно на каждой итерации случайно разбивать множество точек так, чтобы в одну группу попадали близкие точки, и проводить блочно-покоординатный спуск параллельно для каждой группы [9]. В отличие от [16], в описанном ниже методе координата точки сразу после шага проектируется на поверхность сферы, а не раскладывается на нормальную и касательную составляющие.

2. МЕТОД ПЕРВОГО ПОРЯДКА

Итак, мы минимизируем функцию $3N$ переменных ${\mathbf{x}} = ({{x}_{1}},\; \ldots ,\;{{x}_{N}})$

Зафиксируем положения $N - 1$ точек (без $k$-й) на единичной сфере. Обозначим через

Градиент и гессиан (3) можно вычислить аналитически, при $s > 0$ и $y: = {{x}_{i}}$:

Оценим норму гессиана (8). В силу эквивалентности норм в конечномерных пространствах ${{l}_{p}}$ и ${{l}_{q}}$ (размерность пространства $N - 1$) при $s \geqslant 1$ найдется константа $C = C(s,N)$:

Тогда если обозначить $Q = \left\{ {{\mathbf{x}} \in \prod\nolimits_{i = 1}^N {S_{i}^{2}} :U({\mathbf{x}}) \leqslant U({{{\mathbf{x}}}^{0}})} \right\}$, то на таком множестве будет выполнено

Несложно показать, что минимизацией квадратичного приближения (6):

Утверждение 1. Пусть $f({\kern 1pt} \cdot {\kern 1pt} ):{{\mathbb{R}}^{n}} \to \mathbb{R}$ дифференцируема в точке $x \in {{\mathbb{R}}^{n}}:\left\| x \right\| = 1$ и $\left\| {f{\kern 1pt} {\text{'}}(x)} \right\| \ne \left| {(f{\kern 1pt} {\text{'}}(x),x)} \right|$. Тогда существует длина шага ${{t}_{0}} > 0$, что

Доказательство 1. Обозначим $y = y(x,t): = \tfrac{{x - tf{\kern 1pt} {\text{'}}(x)}}{{\left\| {x - tf{\kern 1pt} {\text{'}}(x)} \right\|}}$. В силу дифференцируемости $f({\kern 1pt} \cdot {\kern 1pt} )$:

Тогда $o\left( {\left\| {y - x} \right\|} \right) = o(t)$. Следовательно, найдется ${{t}_{2}} > 0$: $\forall t \in (0,{{t}_{2}})$

Определение 1. Точку ${\mathbf{x}} \in \prod\nolimits_{i = 1}^N \,S_{i}^{2}$ будем называть стационарной точкой функции $U( \cdot )$, если

В недавней работе [17] получены оценки сходимости для задачи вида (6) при некоторых предположениях. В частности, при условии липшицевости градиента на сфере метод (12) генерирует невозрастающую по функции последовательность точек. В нашем случае применим метод дробления шага, который заключается в том, чтобы уменьшать шаг ${{t}_{0}}$ вдвое, при невыполнении условия (13).

Утверждение 2. Пусть $L = C{{(U({{{\mathbf{x}}}^{0}}))}^{{\tfrac{s}{{s + 2}}}}}$ для некоторого начального положения ${{{\mathbf{x}}}^{0}}$, причем $s \geqslant 1$. Тогда последовательность точек ${\text{\{ }}{{{\mathbf{x}}}^{n}}{\text{\} }}_{{n = 1}}^{\infty }$, полученная из (12) (например, последовательным обновлением компонент ${{x}_{k}}$) методом дробления шага, монотонно $U({{{\mathbf{x}}}^{{n + 1}}}) \leqslant U({{{\mathbf{x}}}^{n}})$  $\forall n \in \mathbb{N}$ сходится к стационарной точке функции $U( \cdot )$.

Отметим, что для функции $U({\kern 1pt} \cdot {\kern 1pt} )$, условие (10) не выполнено на всем множестве $\prod\nolimits_{i = 1}^N \,S_{i}^{2}$, однако, ясно, что если (9) выполнено для некоторого ${\mathbf{x}}$, то после шага (12) данная оценка не ухудшится и (10) будет справедливо на траектории.

Итак, опишем общую структуру метода.

1. Выбрать начальное приближение ${{{\mathbf{x}}}^{0}}$.

2. Пока не достигнут критерий останова на $n$-й итерации.

2.1. Выбрать индекс $k$-й точки, которую будем изменять.

2.2. Выбрать шаг ${{t}_{0}} = 1{\text{/}}L$.

2.3. Вычислить градиент $\nabla {{U}_{k}}(x_{k}^{n},{\mathbf{\alpha }}_{k}^{n}))$ относительно $k$-й точки.

2.4. Сделать шаг (12).

2.5. Если выполнено условие (13), обновить $k$-ю компоненту ${{{\mathbf{x}}}^{n}}$.

Иначе дробить шаг ${{t}_{0}}: = {{t}_{0}}{\text{/}}2$ и перейти к п. 2.4.

3. МЕТОД ВТОРОГО ПОРЯДКА

Идея метода Ньютона и других методов второго порядка заключается в более точном приближении оптимизируемой функции квадратичной функцией вместо параболоида вращения (11), что ведет к большему выигрышу по функции на каждой итерации. Однако несмотря на быструю сходимость таких алгоритмов по количеству итераций, их главным недостатком остается сложность вычисления и обращения матрицы вторых производных. В связи с чем активно развиваются различные модификации метода Ньютона, основанные на приближенном подсчете гессиана (квазиньютоновские методы) [18], [19], а также различные покоординатные версии [20], [21]. Здесь мы остановимся на покоординатной версии метода Ньютона для оптимизации с ограничением на сфере. Разложим функцию ${{U}_{k}}(z,{{\alpha }_{k}}):{{\mathbb{R}}^{3}} \to \mathbb{R}$, зависящую только от координаты одной точки, в ряд Тейлора до второго порядка в окрестности точки ${{x}_{k}}$:

Решение данной подзадачи опирается на следующий результат.

Утверждение 3. Если ${{\lambda }_{1}} < {{\lambda }_{2}} \leqslant \; \ldots \; \leqslant {{\lambda }_{m}}$ – собственные значения матрицы ${\mathbf{A}}$, а ${{\phi }_{1}},\; \ldots ,\;{{\phi }_{m}}$ – соответствующие ортонормированные собственные векторы и ${{\beta }_{i}} = {{{\mathbf{b}}}^{{\text{т}}}}(x){{\phi }_{i}}$, $i = 1,\; \ldots ,\;m$, причем ${{\beta }_{1}} \ne 0$. Тогда $\phi = \sum\nolimits_{i = 1}^m {{{c}_{i}}} {{\phi }_{i}}$ будет решением задачи (15), если ${{c}_{i}} = \tfrac{{{{\beta }_{i}}}}{{{{\lambda }_{i}} + \mu }}$, $i = 1,\; \ldots ,\;m$, $\mu > - {{\lambda }_{1}}$, выбирается из условия

Это позволяет построить эффективный метод нахождения решения задачи (15), используя метод деления отрезка пополам. Для этого достаточно ограничить $\mu $ из утверждения (3). Поскольку $0 < {{\lambda }_{1}} + \mu < {{\lambda }_{i}} + \mu $, $1 \leqslant i \leqslant m$, то

Итак, опишем общую структуру проективного покоординатного метода Ньютона.

1. Выбрать начальное приближение ${{{\mathbf{x}}}^{0}}$.

2. Пока не достигнут критерий останова на $n$-й итерации.

2.1. Выбрать индекс $k$ точки, которую будем изменять.

2.2. Вычислить градиент $\nabla {{U}_{k}}(x_{k}^{n},{\mathbf{\alpha }}_{k}^{n}))$ и гессиан ${{\nabla }^{2}}{{U}_{k}}(x_{k}^{n},{\mathbf{\alpha }}_{k}^{n}))$.

2.3. Сделать шаг (15), обновив $k$-ю компоненту ${\mathbf{x}}$.

4. ВЫБОР НАЧАЛЬНОГО ПРИБЛИЖЕНИЯ, БЫСТРЫЙ СЧЕТ ФУНКЦИИ И ГРАДИЕНТА

5. РЕЗУЛЬТАТЫ ВЫЧИСЛИТЕЛЬНЫХ ЭКСПЕРИМЕНТОВ

Для ясности дальнейшего анализа результатов, следуя [26], введем понятие “эпоха” (epoch) – $N$ последовательных итераций покоординатного спуска, где $N$ – количество точек. Итерацией по-прежнему будем называть шаг вида (12). При сравнении с полноградиентным методом естественно сопоставить один шаг полноградиентного метода с одной эпохой покоординатного спуска.

В разд. 4 предложен способ быстрого счета функции и градиента. Время выполнения метода первого порядка при фиксированном количестве итераций определяется линейными членами по количеству точек $N$. Мы исследовали зависимость реального времени выполнения метода (от количества точек $N$), запущенного на современном ноутбуке (процессор 2.8 GHz, оперативная память 16 ГБ). Как видно из фиг. 1, гипотеза о линейности оправдана.

Фиг. 1.

Линейное время выполнения.

В предположении условия (10) было показано свойство релаксации. Интересно экспериментально выяснить скорость сходимости метода. Мы ожидаем получить сходимость по функции $U({\mathbf{x}})$ со скоростью геометрической прогрессии с параметрами $a$ и $b$:

Таким образом, величина $a$ приближенно равна знаменателю геометрической прогрессии, с которой убывает отклонение функции от минимального, приближенно равного $\tfrac{b}{{1 - a}}$. Для подбора параметров $a$ и $b$ используется нелинейный метод наименьших квадратов. График на фиг. 2 показывает скорость сходимости метода относительно функции (1) при $s = 1$ в режиме $epoch \leqslant N$. Это соответствует так называемому режиму “large scale” оптимизации, так как для данного покоординатного метода одна итерация оказывается в $N$ раз дешевле подсчета полного градиента. Параметры МНК и минимально найденное значение функции $U({{{\mathbf{x}}}^{n}})$ для различных $N$ указаны в табл. 1. В четвертом столбце рассчитано предельное значение аппроксимации $\tfrac{b}{{1 - a}}$, к которому ассимптотически стремится модель (17), (18). Сравнивая столбцы $U({{{\mathbf{x}}}^{n}})$ и $U({\mathbf{x}}*)$ табл. 1, можно заметить, что алгоритм достаточно близко приближается к минимально известным значениям. Различие составляет порядка 10^–3%.

Фиг. 2.

Скорость сходимости, $N = 972$, $s = 1.$

В разд. 4 также описывается возможность выбора точки для движения за время $O(N)$ на каждой итерации. На фиг. 3 изображены графики сходимости для различных стратегий при $s = 0$. Оказывается, что на начальном этапе стратегия выбора точки с максимальной проекцией градиента (max_grad) достаточно сильно обгоняет циклический перебор (cyclic). Стратегия (max_grad_prob), которая заключается в том, что точки выбираются пропорционально квадратам проекций, также оказывается достаточно хорошей. Вид графика оказывается типичным и при других $s$.

Фиг. 3.

Сравнение способов выбора точек для выполнения шага, $N = 500$, s = 0.

На фиг. 5 и фиг. 6 изображены графики сходимости на последних 100 из 5000 итераций для $N = 30$ для методов первого и второго порядка соответственно. Начальные положения точек генерировались методом Монте-Карло. Можно заметить, что функция (1) для $s = 1$ является многоэкстремальной уже при $N = 30$. Интересно также, что метод второго порядка намного чаще достигает глобального минимума, а метод первого порядка застревает в точках, достаточно далеких от оптимума. Прерывистыми линиями выделены графики при генерации точек методом generalized spiral.

Фиг. 4.

Сравнение полноградиентного метода с покоординатными методами с градиентным и ньютоновским шагами, N = 470, s = 1.

Фиг. 5.

Локальные минимумы для $N = 30$, $s = 1$. Метод 1 порядка.

Фиг. 6.

Локальные минимумы для $N = 30$, $s = 1$. Метод 2 порядка.

В разд. 3 предложен покоординатный аналог метода Ньютона с ограничением на сфере. Ожидается, что скорость сходимости по количеству итераций будет достаточно высокой, чтобы оправдать решение более сложной подзадачи на каждой итерации. Были проведены численные эксперименты с целью сравнить два предложенных метода в режиме $epoch \leqslant N$. На фиг. 4 изображены графики сходимости для данных двух методов и полноградиентного метода при $s = 1$. Оказывается, что время работы алгоритмов proj_newt и proj_grad в среднем отличается примерно в два раза, при этом из графиков видно, что в рассмотренном режиме proj_newt обгоняет proj_grad больше, чем в два раза по функции. То есть количество эпох, которое требуется методу proj_newt, чтобы достичь значения функции, которое достиг proj_grad за $N$ эпох меньше, чем $N{\text{/}}2$. Также можно заметить, что оба метода значительно обгоняют полноградиентный метод при рассмотрении сходимости по функции. В табл. 2 и 3 приведены результаты работы двух методов. В табл. 3 критерием останова, кроме количества итераций, задано изменение по функции на последних двух эпохах (если оно оказывается меньше ${{10}^{{ - 12}}}$, то алгоритм останавливается). Можно заметить, что метод 2 порядка (proj_newt) значительно чаще достигает меньших минимумов, при этом ему часто для этого требуется меньшее число итераций.

Расположение точек на сфере независимо от вида задачи (1) наглядно представляeтся с помощью сферических диаграмм Вороного. Данные диаграммы представляют собой разбиение поверхности сферы на области (многоугольники), в которых каждая точка является ближайшей к данной точке (электрону). Большинство областей – шестиугольники, многоугольники с иным числом сторон выделены цветом. Интуитивно можно ожидать, что чем меньше таких “неправильных” многоугольников, тем лучше прооптимизирована целевая функция. На фиг. 7–10 представлены диаграммы Вороного, в каждой паре слева – диаграмма при начальной генерации точек методом generalized spiral, справа – диаграмма после оптимизации с функцией (1) при $s = 1$. Отсюда можно заметить, что выбранное начальное приближение располагает точки достаточно равномерно, структура диаграмм после оптимизации согласуeтся с результатами из литературы [11], [27].

Фиг. 7.

Диаграмма Вороного до оптимизации, N = 100.

Фиг. 8.

Диаграмма Вороного после оптимизации, N = 100.

Фиг. 9.

Диаграмма Вороного до оптимизации, N = 972.

Фиг. 10.

Диаграмма Вороного после оптимизации, N = 972.

6. ЗАКЛЮЧЕНИЕ

Предложены два проективных покоординатных метода для поиска локального минимума в задаче Фекете. Алгоритмы используют быстрый счет функции, градиента и выбор наилучшей компоненты для оптимизации. Наши численные эксперименты показывают превосходство предложенных методов перед полноградиентным методом при количестве точек $N < 1000$. Можно ожидать, что предложенные методы окажутся предпочтительными для задач больших размерностей.