Журнал вычислительной математики и математической физики, 2019, T. 59, № 12, стр. 2086-2101

ВВЕДЕНИЕ

Работа посвящена теоретическому и экспериментальному изучению метода ньютоновского типа, использующего приближенные направления, получаемые с применением предобусловленного метода сопряженных градиентов, и является обобщением и уточнением результатов [1]. Ранее подобный метод успешно применялся в задачах оптимизации, возникающих при генерации численных сеток [2]–[4], см. также [5]. Здесь рассматривается применение данного метода для численного решения задач безусловной минимизации выпуклой кусочно-квадратичной функции. К таким задачам относятся нахождение нормального решения задачи линейного программирования (ЛП) [6], [7] и нахождение проекции заданной точки на множество решений задачи ЛП [8], [9]. С задачей ЛП тесно связано нахождение проекции точки на множество неотрицательных решений недоопределенной системы линейных уравнений [10], [11] и вычисление расстояния между двумя выпуклыми многогранниками [12].

Впервые О. Мангасариан показал эффективность применения обобщенного метода Ньютона для минимизации выпуклой кусочно-квадратичной функции. В своих работах [7], [13] он исходил из возможности простого вычисления обобщенной матрицы Гессе, предложенного в [11], так как кусочно-квадратичная функция непрерывно дифференцируема только один раз.

В разд. 1 задача нахождения неотрицательного решения недоопределенной системы линейных равенств рассматривается с точки зрения линейного программирования, что приводит к рассмотрению взаимно двойственных задач, одна из которых (прямая задача) является оштрафованной, а вторая (двойственная задача) – регуляризованной. При этом через решение задачи безусловной оптимизации кусочно-квадратичной функции легко вычисляется нормальное решение исходной линейной системы (проекция точки). В разд. 2 приведены некоторые вспомогательные технические результаты, связанные с кусочно-квадратичной функцией. В разд. 3 описывается метод ньютоновского типа для безусловной оптимизации кусочно-квадратичной функции. Раздел 4 содержит анализ сходимости метода ньютоновского типа с учетом специального правила остановки внутренних итераций предобусловленного метода сопряженных градиентов. В разд. 5 представлены численные результаты для различных модельных задач.

1. ЗАДАЧИ ОПТИМИЗАЦИИ И СИСТЕМА ЛИНЕЙНЫХ УРАВНЕНИЙ И НЕРАВЕНСТВ

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

Рассмотрим задачу нахождения неотрицательного решения системы линейных уравнений

Здесь матрица $A \in {{R}^{{m \times n}}}$, вектор $b \in {{R}^{m}}$ заданы и $m \leqslant n$. Через ${{0}_{i}}$ обозначен $i$-мерный нулевой вектор. Будем предполагать, что исходная система (1) разрешима. Оказывается, эту систему полезно рассмотреть как задачу линейного программирования (ЛП) в которой вектор целевой функции $c = {{0}_{n}}$. Двойственная к (2) имеет вид

Прямую задачу линейного программирования $P$ преобразуем следующими двумя способами: 1) с помощью квадратичного штрафа (преобразование $pen$), 2) используя квадратичную регуляризацию (преобразование $reg$). Можно продолжить преобразование полученных задач, применяя, соответственно, квадратичную регуляризацию и квадратичный штраф. Последовательность преобразований задачи $P$ можно представить схематически в виде ${{P}_{{pr}}}\;\mathop { \leftarrow \;}\limits^{reg} {{P}_{p}}\;\mathop { \leftarrow \;}\limits^{pen} P\;\mathop { \to \;}\limits^{reg} {{P}_{r}}\;\mathop { \to \;}\limits^{pen} {{P}_{{rp}}}$. Если применить квадратичный штраф к задаче $P$, то приходим к минимизации квадратичной функции на неотрицательном ортанте:

Вообще говоря, (3) есть задача наименьших квадратов для решения системы (1). В силу постоянства целевой функции в задаче $P$, которая всегда равна нулю, в задаче ${{P}_{p}}$ не требуется вводить коэффициент штрафа. Отметим, что квадратичная задача (3) всегда разрешима независимо от того, имеет ли решение исходная система (1).

Если к задаче $P$ применить квадратичную регуляризацию (регуляризация по Тихонову [15]), то получим задачу строго выпуклого квадратичного программирования при ограничениях

В этой задаче также не требуется вводить коэффициент регуляризиции.

Задачу минимизации на неотрицательном ортанте ${{P}_{p}}$ можно регуляризовать, и тогда приходим к следующей строго выпуклой квадратичной задаче:

где $\varepsilon > 0$ представляет собой коэффициент регуляризации. Решение $x(\varepsilon )$ задачи (5) единственно при любом $\varepsilon > 0$. При $\varepsilon \to 0$ решение $x(\varepsilon )$ сходится к нормальному решению системы (1) ([15]), т.е. к решению задачи (4).

Если в задаче ${{P}_{r}}$ снять ограничения с помощью квадратичного штрафа, то получаем следующую задачу:

Очевидно, что задачи ${{P}_{{rp}}}$ и ${{P}_{{pr}}}$ эквивалентны в смысле совпадения их решений при любом положительном $\varepsilon $.

Выпишем теперь двойственные задачи к приведенным выше оптимизационным задачам. Начнем с задачи ${{P}_{p}}$, двойственной для которой является задача строго вогнутого квадратичного программирования:

Хорошо известно (см., например, [16]), что задачу безусловной оптимизации можно легко свести к задаче с ограничениями с помощью дополнительных переменных и естественно возникающих ограничений. Для такой задачи уже можно ввести функцию Лагранжа и стандартным образом получить двойственную задачу. Задачи ${{P}_{p}}$ и ${{D}_{r}}$ являются взаимно двойственными. По решению ${{x}_{ * }}$ задачи ${{P}_{p}}$ легко вычисляется решение задачи ${{D}_{r}}$ по формуле ${{u}_{ * }} = b - A{{x}_{ * }}$. Следует отметить, что квадратичные задачи ${{P}_{p}}$ и ${{D}_{r}}$ всегда разрешимы вне зависимости от разрешимости задач линейного программирования $P$ и $D$. Задача ${{D}_{r}}$ получается применением квадратичной регуляризации к задаче $D$.

Двойственной к ${{P}_{r}}$ является следующая оштрафованная задача $D$

Здесь и ниже ${{a}_{ + }}$ обозначает вектор $a$, у которого все отрицательные компоненты заменены на нули. По решению ${{u}_{ * }}$ задачи ${{D}_{p}}$ находится решение ${{x}_{ * }} = {{({{A}^{{\text{т}}}}{{u}_{ * }})}_{ + }}$ задачи ${{P}_{r}}$, т.е. нормальное решение исходной линейной системы (1).

Приведем двойственную задачу для ${{P}_{{rp}}}$

которая является задачей безусловной максимизации строго вогнутой кусочно-квадратичной функции. При любом $\varepsilon > 0$ по решению $u(\varepsilon )$ задачи ${{D}_{{pr}}}$ вычисляется решение $x(\varepsilon ) = {{({{A}^{{\text{т}}}}u(\varepsilon ))}_{ + }}$ задачи ${{P}_{{rp}}}$. При $\varepsilon \to 0$ решение $x(\varepsilon )$ сходится к решению задачи ${{P}_{r}}$, т.е. к нормальному решению исходной системы (1). Итак, оштрафованную задачу ${{P}_{r}}$, т.е. задачу минимизации ${{P}_{{rp}}}$ на положительном ортанте с $n$ переменными, заменили задачей безусловной максимизации ${{D}_{{pr}}}$ другой размерности с $m$ переменными.

В свою очередь решение $u(\varepsilon )$ задачи ${{D}_{{pr}}}$ находится через решение $x(\varepsilon )$ задачи ${{P}_{{rp}}}$ минимизации строго выпуклой кусочно-квадратичной функции на неотрицательном ортанте по простой формуле $u(\varepsilon ) = \tfrac{1}{\varepsilon }(b - Ax(\varepsilon ))$.

Аналогично получается двойственная задача для ${{P}_{{pr}}}$:

При любом $\varepsilon > 0$ решения взаимно двойственных задач ${{P}_{{pr}}}$ и ${{D}_{{rp}}}$ выражаются друг через друга. По решению $u(\varepsilon )$ задачи ${{D}_{{rp}}}$ легко находится решение $x(\varepsilon ) = \tfrac{1}{\varepsilon }{{({{A}^{{\text{т}}}}u(\varepsilon ))}_{ + }}$ задачи ${{P}_{{pr}}}$ и из решения $x(\varepsilon )$ задачи ${{P}_{{pr}}}$ имеем решение $u(\varepsilon ) = b - Ax(\varepsilon )$ задачи ${{D}_{{rp}}}$.

Приведенную связь оштрафованных и регуляризованных задач с помощью двойственности можно представить в виде следующей схемы:

В эту же схему укладываются и задачи нахождения проекции заданной точки $\hat {x}$ на множество решений системы (1). Регуляризованная задача $P$ будет иметь вид

взаимно двойственной к ней будет задача При этом через ее решение ${{u}_{ * }}$ находится решение ${{\hat {x}}_{ * }} = {{(\hat {x} + {{A}^{{\text{т}}}}{{u}_{ * }})}_{ + }}$ задачи ${{P}_{r}}$ (6), т.е проекция точки $\hat {x}$ на множество решений системы (1). Задача ${{D}_{p}}$ (7) есть оштрафованная задача $D$, в которой ограничения сдвинуты на вектор $\hat {x}$. Отметим, что задача ${{D}_{p}}$ (7) безусловной максимизации от $m$ переменных имеет в общем случае неединственное решение.

Если оштрафовать задачу ${{P}_{r}}$ (6), то приходим к задаче минимизации строго выпуклой квадратичной функции от $n$ переменных на неотрицательном ортанте:

взаимно двойственная к которой имеет вид Задача (9) является задачей безусловной максимизации строго вогнутой кусочно-квадратичной функцией от $m$ переменных. Эта задача получается регуляризацией задачи ${{D}_{p}}$ (7). При любом $\varepsilon > 0$ решения $x(\varepsilon )$ и $u(\varepsilon )$ соответственно задач (8) и (9) выражаются по простым формулам $x(\varepsilon ) = {{(\hat {x} + {{A}^{{\text{т}}}}u(\varepsilon ))}_{ + }}$ $u(\varepsilon ) = \tfrac{1}{\varepsilon }(b - Ax(\varepsilon ))$.

Далее остановимся на рассмотрении задачи (7), которую перепишем в следующем виде безусловной минимизации кусочно-квадратичной функции

Задача (10) является взаимно двойственной к задаче нахождения проекции заданного вектора $\hat {x}$ на множество неотрицательных решений недоопределенной системы линейных уравнений:

Решение ${{x}_{ * }}$ этой задачи выражается через решение ${{u}_{ * }}$ задачи (10) по простой формуле ${{x}_{ * }} = {{(\hat {x} + {{A}^{{\text{т}}}}{{u}_{ * }})}_{ + }}$.

Итак, введем в рассмотрение выпуклую дифференцируемую кусочно-квадратичную функцию $\varphi (u):{{R}^{m}} \to {{R}^{1}}$:

Ее градиент $g(u) = {\text{grad}}\varphi (u)$ имеет вид а обобщенную матрицу Гессе [14] определим как Связь между $H(u)$, $\varphi (u)$ и $g(u)$ будет установлена ниже в Замечании 1.

2. ТЕЙЛОРОВСКОЕ РАЗЛОЖЕНИЕ ФУНКЦИИ $( \cdot )_{ + }^{2}$

Следующий результат представляет собой специальный случай разложения Тейлора с остаточным членом в интегральной форме.

Лемма 1. Для любых действительных скалярных величин $\eta $ и $\zeta $ справедливо

Доказательство. Пусть $f(\xi ) = \tfrac{1}{2}{{({{\xi }_{ + }})}^{2}}$, тогда $f{\kern 1pt} {\text{'}}{\kern 1pt} (\xi ) = {{\xi }_{ + }}$ и $f{\kern 1pt} {\text{''}}(\xi ) = {\text{sign}}({{\xi }_{ + }})$ (без потери общности будем полагать $f{\text{''}}(0) = 0$) [14]. Подставляя эти равенства в тейлоровское разложение

приходим к нужному результату.

Лемма 2. Для любых $n$-мерных векторов $y$ и $z$ справедливо

где

Доказательство. Полагая $\eta = {{y}_{j}}$, $\zeta = {{z}_{j}}$ в (14) и суммируя по всем $j = 1, \ldots ,\;n$, приходим к нужной формуле. Заметим, что использование множителя 2 в подынтегральном выражении приводит к оценке $\left\| {{\text{Diag}}(d)} \right\| \leqslant 1$.

Лемма 3. Функция (11) и ее градиент (12) удовлетворяют равенству

где

Доказательство. Полагая $y = \hat {x} + {{A}^{{\text{т}}}}u$ и $z = {{A}^{{\text{т}}}}q$ в (15) и (16), легко получаем требуемый результат (с учетом приведения линейных членов, включающих b в левой части (17)).

Замечание 1. Как видно из (18), если условие

справедливо для любого $0 \leqslant \vartheta \leqslant 1$, то (17) приводится к простейшему виду: где обобщенная матрица Гессе $H(u)$ задана формулой (13). Это объясняет ключевую роль $H(u)$ в рассмотренном ниже методе ньютоновского типа. Заметим, что достаточным условием для выполнения (19) является т.е. если некоторые компоненты приращения $q$ относительно малы, то $\varphi $ в точности квадратична (20) в соответствующей окрестности $u$.

3. МЕТОД НЬЮТОНОВСКОГО ТИПА ДЛЯ ДВОЙСТВЕННОЙ ЗАДАЧИ (10)

Из условия (19) и его следствия (20) вытекает, что можно численно минимизировать $\phi $ методом ньютоновского типа ${{u}_{{k + 1}}} = {{u}_{k}} - {{d}_{k}}$, где ${{d}_{k}} = H{{({{u}_{k}})}^{{ - 1}}}g({{u}_{k}})$. Заметим, что согласно (20) точка минимума ${{u}_{ * }} = {{u}_{{k + 1}}}$ получается за один шаг, если компоненты ${{d}_{k}}$ достаточно малы и удовлетворяют (21), где $u = {{u}_{k}}$ and $q = - {{d}_{k}}$. Однако начальные ${{u}_{k}}$ могут быть еще далеки от решения и поэтому возможны только постепенные улучшения точности приближений. Во-первых, должен использоваться коэффициент ${{\alpha }_{k}}$ для гарантии монотонной сходимости (в смысле уменьшения $\varphi ({{u}_{k}})$ при увеличении $k$). Во-вторых, вместо $H({{u}_{k}})$ должна использоваться некоторая ее аппроксимация ${{M}_{k}}$, являющаяся обратимой матрицей с не слишком большой нормой обратной $\left\| {M_{k}^{{ - 1}}} \right\|$. Поэтому предлагается следующий прототип метода:

где Выбор параметров $0 < {{\alpha }_{k}} \leqslant 1$ и $0 \leqslant \delta \ll 1$ существенен для лучшей сходимости. Кроме того, на начальных этапах итераций наиболее эффективной стратегией является использование приближенных ньютоновских направлений ${{d}_{k}} \approx M_{k}^{{ - 1}}g({{u}_{k}})$, которые могут быть получены предобусловленным методом сопряженных градиентов, примененным к решению уравнения ${{M}_{k}}{{d}_{k}} = g({{u}_{k}})$. Как будет показано ниже, достаточно использовать любой вектор ${{d}_{k}}$, который удовлетворяет условиям с $0 < {{\vartheta }_{k}} < 1$, достаточно хорошо отделенным от нуля. При любом предобусловливании метода сопряженных градиентов такие аппроксимации удовлетворяют (23) (см. п. 4.3 ниже). Для выбора параметра $0 < \alpha \leqslant 1$ используется правило Армийо где выбирается максимально возможное значение $0 < \alpha \leqslant 1$, при котором выполняется условие (24). Соответствующий алгоритм ньютоновского типа можно представить в виде:

АЛГОРИТМ 1

 Input:

 $A \in {{R}^{{m \times n}}}$, $b \in {{R}^{m}}$, $\hat {x} \in {{R}^{m}}$;

 Initialization:

 $\delta = {{10}^{{ - 6}}}$, $\varepsilon = {{10}^{{ - 12}}}$, $\tau = {{10}^{{ - 15}}}$;

 ${{k}_{{max}}} = 2000$, ${{l}_{{max}}} = 10$, ${{u}_{0}} = {{0}_{m}}$;

 Iterations:

 for $k = 0,1,\; \ldots ,\;{{k}_{{max}}} - 1$:

  ${{x}_{k}} = {{(\hat {x} + {{A}^{{\text{т}}}}{{u}_{k}})}_{ + }}$

  ${{\varphi }_{k}} = \tfrac{1}{2}{{\left\| {{{x}_{k}}} \right\|}^{2}} - {{b}^{{\text{т}}}}{{u}_{k}}$

  ${{g}_{k}} = A{{x}_{k}} - b$

  if $(\left\| {{{g}_{k}}} \right\| \leqslant \varepsilon \left\| b \right\|)$ return ${\text{\{ }}{{x}_{k}},{{u}_{k}},{{g}_{k}}{\text{\} }}$

  find ${{d}_{k}} \in {{R}^{m}}$ such that

   $d_{k}^{{\text{т}}}{{g}_{k}} = d_{k}^{{\text{т}}}{{M}_{k}}{{d}_{k}} = \vartheta _{k}^{2}{{g}_{k}}M_{k}^{{ - 1}}{{g}_{k}}$,

   where ${{M}_{k}} = A{\text{Diag}}({\text{sign}}({{x}_{k}})){{A}^{{\text{т}}}} + \delta {\text{Diag}}(A{{A}^{{\text{т}}}})$

  ${{\alpha }^{{(0)}}} = 1$

  $u_{k}^{{(0)}} = {{u}_{k}} - {{d}_{k}}$

  for $l = 0,1,\; \ldots ,\;{{l}_{{max}}} - 1$:

   $x_{k}^{{(l)}} = {{(\hat {x} + {{A}^{{\text{т}}}}u_{k}^{{(l)}})}_{ + }}$

   $\varphi _{k}^{{(l)}} = \tfrac{1}{2}{{\left\| {x_{k}^{{(l)}}} \right\|}^{2}} - {{b}^{{\text{т}}}}u_{k}^{{(l)}}$

   $\zeta _{k}^{{(l)}} = \left( {\tfrac{1}{2}{{\alpha }^{{(l)}}}d_{k}^{{\text{т}}}{{g}_{k}} + \varphi _{k}^{{(l)}}} \right) - {{\varphi }_{k}}$

   if ($\zeta _{k}^{{(l)}} > \tau \left| {{{\varphi }_{k}}} \right|$) then

    ${{\alpha }^{{(l + 1)}}} = {{\alpha }^{{(l)}}}{\text{/}}2$

    $u_{k}^{{(l + 1)}} = {{u}_{k}} - {{\alpha }^{{(l + 1)}}}{{d}_{k}}$

   else

    $u_{k}^{{(l + 1)}} = u_{k}^{{(l)}}$

    go to NEXT

   end if

  end for

  NEXT: ${{u}_{{k + 1}}} = u_{k}^{{(l + 1)}}$

 end for

Далее исследуем свойства сходимости этого алгоритма.

4. АНАЛИЗ СХОДИМОСТИ МЕТОДА НЬЮТОНОВСКОГО ТИПА

Оказывается, что анализ сходимости Алгоритма 1 можно выполнить тем же образом, что и в работах [4] и [17] (где, однако, рассматривается нелинейная задача наименьших квадратов для гладких функций). Для полноты изложения и совместимости обозначений мы воспроизводим здесь основные результаты из [4].

АЛГОРИТМ 2

 ${{r}^{{(0)}}} = - g,$ ${{d}^{{(0)}}} = {{s}^{{( - 1)}}} = {{t}^{{( - 1)}}} = 0,$ ${{\zeta }^{{( - 1)}}} = 0;$

 $i = 0,1,\; \ldots ,\;i{{t}_{{max}}}:$

  ${{w}^{{(i)}}} = C{{r}^{{(i)}}},$

  ${{z}^{{(i)}}} = M{{w}^{{(i)}}},$

  ${{\gamma }^{{(i)}}} = {{({{r}^{{(i)}}})}^{{\text{т}}}}{{w}^{{(i)}}},$  ${{\xi }^{{(i)}}} = {{({{w}^{{(i)}}})}^{{\text{т}}}}{{z}^{{(i)}}},$ ${{\eta }^{{(i - 1)}}} = {{({{s}^{{(i - 1)}}})}^{{\text{т}}}}{{t}^{{(i - 1)}}},$

  ${{\zeta }^{{(i)}}} = {{\zeta }^{{(i - 1)}}} + {{\eta }^{{(i - 1)}}},$

  if $((\varepsilon _{{{\text{CG}}}}^{{ - 1}} + i){{\eta }^{{(i - 1)}}} \leqslant {{\zeta }^{{(i)}}})$ or $({{\gamma }^{{(i)}}} \leqslant \varepsilon _{{{\text{CG}}}}^{2}{{\gamma }^{{(0)}}})$ return ${\text{\{ }}{{d}^{{(i)}}}{\text{\} }}$;

  if ($k = 0$) then

    ${{\alpha }^{{(i)}}} = - {{\gamma }^{{(i)}}}{\text{/}}{{\xi }^{{(i)}}},\quad {{\beta }^{{(i)}}} = 0;$

  else

    ${{\delta }^{{(i)}}} = {{\gamma }^{{(i)}}}{\text{/}}({{\xi }^{{(i)}}}{{\eta }^{{(i - 1)}}} - {{({{\gamma }^{{(i)}}})}^{2}})$,

    ${{\alpha }^{{(i)}}} = - {{\eta }^{{(i - 1)}}}{{\delta }^{{(i)}}}$, ${{\beta }^{{(i)}}} = {{\gamma }^{{(i)}}}{{\delta }^{{(i)}}}$;

  end if

  ${{t}^{{(i)}}} = {{z}^{{(i)}}}{{\alpha }^{{(i)}}} + {{t}^{{(i - 1)}}}{{\beta }^{{(i)}}},$ ${{r}^{{(i + 1)}}} = {{r}^{{(i)}}} + {{t}^{{(i)}}},$

  ${{s}^{{(i)}}} = {{w}^{{(i)}}}{{\alpha }^{{(i)}}} + {{s}^{{(i - 1)}}}{{\beta }^{{(i)}}},$ ${{d}^{{(i + 1)}}} = {{d}^{{(i)}}} + {{s}^{{(i)}}}.$

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

Несмотря на некоторое увеличение затрат рабочей памяти и количества векторных операций алгоритма 2 по сравнению со стандартным (см. следующий раздел), он может допускать более эффективную параллельную реализацию операций скалярного произведения при большом числе процессов MPI и значительной латентности операций обмена данными. Так, на каждой итерации алгоритма 2 достаточно применения одной операции сборки/рассылки MPI$\_$AllReduce($ * $,$ * $,3,…), в то время как стандартный алгоритм требует двух таких операций MPI$\_$AllReduce($ * $,$ * $,1,…). Аналогичные алгебраически эквивалентные переформулировки алгоритма предобусловленного метода сопряженных градиентов можно найти в работе [22] и цитированных в ней источниках.

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

Ниже рассмотрены два семейства тестовых задач, которые могут быть решены сведением к безусловной минимизации кусочно-квадратичных функций. Первое семейство — нахождение проекции заданной точки на множество неотрицательных решений системы линейных уравнений (см. также [10]), второе представлено задачей вычисления расстояния между двумя выпуклыми многогранниками, рассмотренной в работе [12]. Последняя задача важна, например, в таких приложениях, как робототехника и компьютерная анимация.