Журнал вычислительной математики и математической физики, 2021, T. 61, № 7, стр. 1059-1069

1. ВВЕДЕНИЕ

В ранних работах Бориса Теодоровича Поляка (см., например, [1]) был предложен (суб)градиентный метод поиска экстремума в выпуклых задачах математического программирования. Большого внимания специалистов метод не привлек, что связано, очевидно, с некоторыми его “неудобными” особенностями. Во-первых, сходимость метода гарантировалась только при минимизации выпуклых функций, нарушение выпуклости могло приводить к нарушению релаксации на итерациях. Во-вторых, для работы метода требовалось априорное знание точного значения функции в экстремальных точках. Сформулированные требования и наличие развиваемых в те годы методов сопряженных градиентов Флетчера–Ривса и Полака–Поляка–Рибьера (см. [2]) привели к незаслуженному, на наш взгляд, снижению активности по исследованию и использованию предложенного алгоритма.

Новый всплеск интереса к методу появился лишь в последние годы, что связано, очевидно, с рассмотрением больших и сверхбольших задач оптимизации. В соответствии с современной классификацией (см., например, [3], [4]) задачи выпуклой оптимизации можно разделить по размерности на: “Small” – до 100 переменных, “Medium” – от ${{10}^{3}}$ до ${{10}^{4}}$, “Large” – от ${{10}^{5}}$ до ${{10}^{7}}$ и “Huge” – более ${{10}^{8}}$ переменных. В соответствии с доминирующей в настоящее время точкой зрения решать задачи выпуклой оптимизации размерности ${{10}^{8}}$ и более при отсутствии разреженности (при “плотных” векторах) с использованием современной вычислительной техники малореально.

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

В работе рассматриваются модификации метода Поляка, ориентированные на решение больших и свербольших задач оптимизации. На простейших примерах из семейств выпуклых функций сепарабельного и квазисепарабельного типа исследуются возможности численного решения гладких задач безусловной минимизации растущих размерностей. На примерах мультиэкстремальных задач оптимизации потенциалов атомно-молекулярных кластеров Морса (см. [5]) и Китинга (см. [6]) исследуются глобализующие свойства модификации метода Поляка. Для нелинейных задач оптимального управления оцениваются возможности решения аппроксимативных задач на сетках с уменьшающимся постоянным шагом. Рассматриваемые модификации метода реализованы на языке С (компиляторы BCC, GCC, Clang и ICC) под управлением операционных систем Linux, Mac OS X и Windows с использованием технологий распараллеливания OpenMP и MPI. Приводимые в работе вычислительные эксперименты выполнялись как на обычных персональных компьютерах, так и на кластерных системах МВС-100К Межведомственного Суперкомпьютерного Центра (МСЦ) РАН и Института прикладной математики им. М.В. Келдыша (ИПМ) РАН.

2. ОБЩАЯ СТРУКТУРА МЕТОДОВ ГРАДИЕНТНОГО ТИПА

Рассматривается стандартная задача конечномерной оптимизации:

Итерация классического алгоритма градиентного типа записывается в следующем традиционном виде:

где направление спуска ${{d}_{k}}$ выбирается как антиградиент минимизируемой функции $f(x)$ в точке ${{x}_{k}}$:

Размер шага ${{\alpha }_{k}}$ может выбираться разными способами, опирающимися на учет свойств целевой функции. Например, для дифференцируемых функций с известной константой Липшица оптимальным является выбор шага ${{\alpha }_{k}} = 1{\text{/}}L$. Однако такой способ не всегда применим для задач рассматриваемого класса в силу отсутствия информации о значении константы Липшица. Величина шага в направлении спуска может в общем случае вычисляться по формуле

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

3. МЕТОД Б.Т. ПОЛЯКА И ЕГО МОДИФИКАЦИЯ

Метод был предложен Б.Т. Поляком в [1] (и, по нашим сведениям, дальнейшего развития не получил):

где $f{\text{*}}$ – известное минимальное значение функции.

Для одномерного случая метод совпадает с методом Ньютона для решения уравнения $f(x) = f{\text{*}}$.

Поскольку в большинстве реальных постановок задач оптимизации $f{\text{*}}$ заранее не известно, предлагается ввести некоторую величину ${{\delta }_{f}} > 0$, определяющую прогнозируемое значение оптимизируемой функции на каждом шаге:

Сходимость метода следует из теоремы, доказанной в [1]:

Теорема 1. Пусть $f(x)$ является выпуклой, непрерывной, $Q$ выпукло и замкнуто, существует точка минимума $x{\kern 1pt} * \in Q$, $f(x{\kern 1pt} *) = f{\text{*}}$ и $\left\| {\nabla f(x)} \right\| \leqslant c$ на $S = \left\{ {x \in Q,\left\| {x - {{x}^{0}}} \right\| \leqslant \left\| {x{\kern 1pt} {\text{*}} - {{x}^{0}}} \right\|} \right\}.$ Тогда в методе (5) $f({{x}^{n}}) \to f{\text{*}}$, а ${{x}^{n}}$ слабо сходится к некоторой точке минимума.

Здесь функция $f(x)$ определена на некотором множестве $Q$. В случае неактивных границ и использования множеств простой структуры (например, параллелепипедного типа), приведенные результаты применимы к рассматриваемой постановке (1).

Алгоритм предлагаемой модификации метода Б.Т. Поляка строится следующим образом.

В качестве входных данных используются: начальная точка ${{x}_{0}}$, точность по норме градиента ${{\varepsilon }_{\nabla }} > 0$, точность по функции ${{\varepsilon }_{f}} > 0$, прогнозная оценка ${{\delta }_{f}} > 0$, коэффициент корректировки шага $0 < k < 1$, минимально допустимый шаг ${{\varepsilon }_{\lambda }} > 0$, $i = 0$ – номер итерации.

Algorithm 1.  Основной цикл алгоритма

 1: Вычисляется $f({{x}_{i}})$, $\nabla f({{x}_{i}})$, $\left\| {\nabla f({{x}_{i}})} \right\|$

 2: if $\left\| {\nabla f({{x}_{i}})} \right\| < {{\varepsilon }_{\nabla }}$ или $(f({{x}_{{i - 1}}}) - f({{x}_{i}})) < {{\varepsilon }_{f}}$ then

 3:  переходим на шаг 13

 4: end if

 5: Находится шаг $\lambda = \tfrac{{{{\delta }_{f}}}}{{{{{\left\| {\nabla f({{x}_{n}})} \right\|}}^{2}}}}$

 6: Вычисляется ${{x}_{{i + 1}}} = {{x}_{i}} - \lambda \cdot \nabla f({{x}_{i}})$

 7: if $f({{x}_{{i + 1}}}) < f({{x}_{i}})$, then

 8:  полагается $i = i + 1$ и выполняется переход на шаг 1

 9: end if

10: if $\lambda > {{\varepsilon }_{\lambda }}$ then

11:  $\lambda = \lambda \cdot k$ и выполняется переход к шагу 6

12: end if

13: Возвращается ${{x}_{i}}$ – как найденная точка минимума

Утверждение о сходимости предлагаемого метода для выпуклых задач с известным $f{\text{*}}$ и половинным делением шага ($k = 0.5$) является тривиальным следствием теоремы 1, поскольку последовательность модифицированного метода обязательно содержит в себе последовательность исходного, сходимость которого строго доказана.

Можно выделить следующие особенности предложенной модификации.

1. При соответствующей настройке алгоритмических параметров (${{\delta }_{f}}$) предложенный метод может делать “большие” шаги, что добавляет ему глобализующие свойства.

2. Метод имеет высокий потенциал параллелизма.

3. Алгоритм имеет простую конструкцию, что имеет большое значение при реализации метода на ряде высокопроизводительных архитектур, например, GPU.

Рассмотрим далее примеры использования предложенной модификации метода при решении экстремальных задач различных классов.

4. ВЫПУКЛАЯ ОПТИМИЗАЦИЯ. СЕПАРАБЕЛЬНАЯ ТЕСТОВАЯ ЗАДАЧА

Решение задач безусловной минимизации класса “Huge-Scale” естественным образом сопряжено с рядом сложностей, одна из которых – собственно, размерность. Увеличение числа переменных очевидным образом влечет за собой повышение требований к объему памяти доступной вычислительной системы. При рассмотрении задач с предельными размерностями возникают ситуации, когда требуемый объем памяти слишком велик для одного вычислительного узла, что влечет за собой необходимость использования распределенных суперкомпьютерных мощностей. Объем памяти, требуемый для задач различной размерности, приведен в табл. 1.

Можно увидеть, что с точки зрения потребления памяти задачи с размерностями до ${{10}^{9}}$ вполне возможно решать на современных вычислительных системах с общей памятью, задачи с большими размерностями требуют использования памяти суперкомпьютеров.

Рассмотрим сепарабельную функцию, имеющую следующий вид:

Очевидно, что данная функция имеет минимум в точке, где все компоненты вектора $x$ равны 0 и значение функции в этой точке также равно 0.

Расчеты проводились на суперкомпьютерах ИПМ РАН и МСЦ РАН, имеющих ограничение в 1 ГБ ОЗУ на процессорное ядро. Данное ограничение необходимо учитывать при выборе числа требуемых процессоров, чтобы не допустить “перерасхода” памяти на вычислительных узлах. Результаты расчетов для задач с различной размерностью приведены в табл. 2 и на фиг. 1.

Максимальная размерность тестовой сепарабельной задачи, которую удалось решить – ${{10}^{{11}}}$ переменных. На фиг. 1 хорошо видна линейная зависимость времени расчета от числа переменных, что говорит о хорошей масштабируемости предложенной модификации алгоритма.

5. ВЫПУКЛАЯ ОПТИМИЗАЦИЯ. НЕСЕПАРАБЕЛЬНАЯ ТЕСТОВАЯ ЗАДАЧА

Рассмотрим тестовую несепарабельную функцию вида

При реализации варианта алгоритма для этого класса задач потребовалось применение технологии обмена сообщениями MPI для кластерных вычислительных систем. Произведенные оценки показали, что обмен данными между вычислительными узлами при нахождении градиента функции составляют незначительную величину от общего времени расчетов. Основным лимитирующим фактором для алгоритма все так же является объем доступной памяти.

Результаты расчетов для задач с различной размерностью приведены в табл. 3 и на фиг. 1. Произведенные вычислительные эксперименты показали, что для несепарабельных задач также имеет место линейная зависимость временных затрат от размерности. В частности, решение “максимальной” задачи с ${{10}^{{11}}}$ переменными с использованием 1000 процессоров потребовало всего 6 мин.

6. НЕВЫПУКЛАЯ ОПТИМИЗАЦИЯ. ПОИСК МОЛЕКУЛЯРНЫХ КЛАСТЕРОВ С МИНИМАЛЬНОЙ ЭНЕРГИЕЙ

Рассматривается потенциальная функция Морса (см. [5])

Задача минимизации данного потенциала является популярной “мультиэкстремальной” задачей глобальной оптимизации с астрономически растущим от размерности числом локальных экстремумов. Данная функция относится к классу несепарабельных, и оптимальное значение $f{\text{*}}$ неизвестно.

Традиционно для решения задач такого типа применяются специализированные методы невыпуклой оптимизации. Предложенная модификация алгоритма, дополненная механизмом случайного мультстарта, позволила успешно найти глобальные решения для кластеров, содержащих до 80 атомов. При этом проявились нестандартные “глобализующие” свойства изначально локального метода, аппроксимирующего случайное распределение проб в окрестности рекордного значения. Примеры проекций траекторий спусков приведены на фиг. 2.

7. МИНИМИЗАЦИЯ ПОТЕНЦИАЛА КИТИНГА

Рассмотрим задачу минимизации потенциала Китинга (см. [6]), моделирующего межатомное взаимодействие в кристаллических структурах “кремний–германий”:

С применением предложенной модификации метода авторами производилась минимизация потенциала Китинга для кристалла с 1 029 000 переменными (343 000 атомов в кристаллической решетке). Проводилось сравнение результатов расчетов с результатами вычислений, полученных с помощью метода сопряженных градиентов (вариант Полак–Поляк–Рибьер) и “метода Коши” (модификация метода наискорейшего спуска). Результаты расчетов приведены в табл. 4. На фиг. 3 показано поведение (сходимость) сравниваемых алгоритмов оптимизации.

Проводилось исследование поведения модифицированного алгоритма при различных значениях параметра ${{\delta }_{f}}$ на задаче минимизации потенциала Китинга для кристалла с 24 000 переменными (8000 атомов), результаты показаны на фиг. 4.

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

8. ЗАДАЧА ОПТИМАЛЬНОГО УПРАВЛЕНИЯ

Традиционная постановка простейшей задачи оптимального управления (так называемая задача со свободным правым концом) заключается в поиске минимума терминального функционала

определенного на траекториях управляемой системы дифференциальных уравнений (дифференциального включения) заданной на интервале времени $t \in [{{t}_{0}},{{t}_{1}}]$. Начальный фазовый вектор ${{x}_{0}} = x({{t}_{0}})$ предполагается фиксированным, на управления наложены ограничения $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{u} \leqslant u(t) \leqslant \bar {u}$, где $\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{u} $, $\bar {u}$ – фиксированные векторы. Рассматриваемая задача занимает центральное место в теории оптимального управления, к такой постановке сводится большое количество задач более сложных классов (с интегральными функционалами, с терминальными и фазовыми ограничениями, с нефиксированным временем и др.). Для задачи в такой постановке получен значительный объем глубоких теоретических результатов и предложено множество вычислительных методов. Однако практика численного решения задач такого класса, особенно прикладных, несмотря на усилия множества специалистов по теории управления, позволяет сделать вывод о доминировании методов, основанных на редукции к задачам конечномерной оптимизации. После введения сетки разбиения интервала времени, например, с равномерным шагом $h$, оказывается возможным сформулировать задачу конечномерной оптимизации, аппроксимирующую задачу в исходной постановке. В аппроксимативной задаче динамика системы описывается рекуррентными соотношениями $x(t + h) = {{f}_{h}}(x(t),u(t),t,h)$, все остальные условия сохраняются. Аппроксимативные задачи, очевидно, представляют собой специальный класс задач математического программирования с условиями типа равенств и неравенств. Если зафиксировать разностную схему, по которой выполняется дискретизация (Эйлера, Рунге–Кутты или другую), то можно воспользоваться предложенным в теории оптимального управления методом оценки градиента функционала, основанным на применении формализма сопряженных переменных. Получаемые по такой схеме алгоритмы оказываются весьма близкими к алгоритмам, которые генерируются с помощью методик быстрого автоматического дифференцирования (см., например, [7]). Вычисление градиента функционала по любой из указанных схем требует решения всего двух задач Коши. Таким образом, задача оптимального управления со свободным правым концом может быть рассмотрена как задача безусловной минимизации с “дешевым”, алгоритмически оцениваемым градиентом, но вычисляемым с некоторой погрешностью, зависящей от размеров шага дискретизации $h$. Размер шага дискретизации непрерывной задачи задает число переменных в аппроксимативной задаче, равное, грубо говоря, $N = ({{t}_{1}} - {{t}_{0}}){\text{/}}h$.

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

Динамика модели описывается системой (см., например, [8])

определенной на интервале времени $t \in [0,\;5]$ с начальными условиями ${{x}_{1}}(0) = 5$, ${{x}_{2}}(0) = 0$. Требуется минимизировать функционал $I(u) = x_{1}^{2}(5) + x_{2}^{2}(5)$ при ограничениях на управления $\left| {u(t)} \right| \leqslant 1$. В задаче имеет место два локальных экстремума, глобально оптимальное значение функционала $I(u*) = 11.908013767$, полученное с применением метода интегрирования DOPRI8 (8-й порядок точности), достигается на релейном управлении:

Для решения задач Коши использовался метод интегрирования из семейства Рунге–Кутты 2-го порядка (так называемый метод “Эйлера с итерациями”). Расчет прерывался при выполнении условия $I({{u}^{{k - 1}}}) - I({{u}^{k}}) < {{10}^{{ - 6}}}$. Все задачи решались, начиная с одинаковых начальных приближений ${{u}^{0}}(t) \equiv - 0.5$. В этой задаче, очевидно, известна гарантированная, но неоптимальная, нижняя оценка значения функционала, равная $0$.

Эксперимент выполнялся с использованием стандартного программного комплекса (ПК) OPTCON-A (см. [9], [10]) в состав которого включена предложенная модификация алгоритма. Результаты эксперимента приведены в табл. 5.

Максимальной размерностью задач, решенных при использовании ПК OPTCON" A, оказалась размерность 819 201, что связано, очевидно, с излишним объемом внутренней памяти комплекса, используемой для работы других алгоритмов. Процессорное время решения при размерностях задач от 101 до 3201 переменной было менее 0.5 с и оценено в 0 с, поскольку используемый датчик времени позволяет учитывать его с точностью до 1. Начиная с размерности 401 переменных, во всех задачах получены значения функционала, с точностью до четырех знаков совпадающие с известным оптимальным значением. Быстрый рост процессорного времени с ростом размерности задач связан с неэффективностью работы процедур аппроксимации управления общего назначения, не позволяющих учитывать специфические особенности алгоритма.

Второй эксперимент проводился со специальной программной реализацией алгоритма, проверенной с использованием нескольких популярных компиляторов на нескольких различных компьютерах. В эксперименте участвовали следующие пары “процессор–компилятор”.

1. Intel Core i7-2640M; Linux; icc-15.0.1, gcc-4.8.2, gcc-4.9.1, clang-3.4.2, clang-3.5.0.

2. Intel Core i7-2640M; Mac OS X; gcc-4.8.3, clang-3.4.2.

3. Intel Core i5 4260U; Mac OS X; gcc-4.8.3, clang-600.0.56 (based on LLVM 3.5svn).

4. AMD Opteron 6220; Linux; icc-14.0.0.

Результаты выполненных расчетов приведены в табл. 6. Наилучшую производительность показала версия алгоритма, реализованная с использованием компилятора icc.

Проведенные эксперименты показали принципиальную возможность решения аппроксимативных задач оптимального управления с размерностями, превышающими ${{10}^{8}}$ даже без применения параллельных вычислительных технологий.

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

Рассмотренная в работе модификация метода Поляка позволила создать вычислительные технологии, обладающие способностью эффективно решать экстремальные задачи размерностей “Huge-Scale” из различных классов. На основе численных экспериментов выявлены ограничения по размерности решаемых “плотных” задач, которые оказались значительно большими, чем предполагалось теоретически (см. [3]). Метод Б.Т. Поляка 1969 г. с применением простейших модификаций может рассматриваться в качестве удобного и высокомасштабируемого инструмента для численного решения “больших” задач конечномерной и бесконечномерной оптимизации.