1. На главную
  2. Электронные версии
  3. Журнал вычислительной математики и математической физики
  4. T. 60, № 2, 2020

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

Локально полиномиальный метод решения систем линейных неравенств

Ю. Г. Евтушенко 1, *, А. А. Третьяков 1, 2, **

1 ВЦ ФИЦ ИУ РАН
119991 Москва, ул. Вавилова, 40, Россия

2 System Res. Inst., Polish Acad. Sie, Newelska 6, 01-447 Warsaw, and University of Siedlct, Faculty of Sciences
08-110 Siedlce, Poland

* E-mail: evt@ccas.ru
** E-mail: tret@ap.siedlce.pl

Поступила в редакцию 13.01.2017
После доработки 16.06.2017
Принята к публикации 17.11.2019

Полный текст (PDF)

Аннотация

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

Ключевые слова: система линейных неравенств, локально-полиномиальная сложность, сходимость, выпуклая функция.

Пусть $X \subset {{\mathbb{R}}^{n}}$ – непустое множество решений системы линейных неравенств

(1)
$Ax \leqslant b,$
где $A$ – матрица размерности $m \times n$, $n \geqslant m$, $x \in {{\mathbb{R}}^{n}}$, $b \in {{\mathbb{R}}^{m}}$. Ставится задача: найти произвольный элемент из множества $X$.

Для решения систем вида (1) существует большое количество численных методов. Укажем лишь на некоторые из них: прямые [1], [2] и итерационные [3]–[5]. Конечный алгоритм полиномиальной вычислительной сложности на данный момент пока не найден. Задача его нахождения включена в список Смейла [6] как одна из важнейших нерешенных проблем. В работах [7], [8] и [17] используется редукция (1) к задачам безусловной минимизации строго выпуклой кусочно-квадратичной функции, показано, что решение получается за конечное полиномиальное число вычислений, если начальная точка принадлежит достаточно малой окрестности единственного решения, либо мы попадаем в эту окрестность каким-либо методом. Однако к минимизируемой функции предъявляются весьма жесткие требования – она должна быть сильно выпуклой, собственные значения матрицы Гессе должны удовлетворять определенным соотношениям и т.д. Это приводит к существенным ограничениям на класс решаемых задач; требуется, например, единственность решения (1) и т.д. Идея метода, описанная в [7], [8], заключается в возможности получения информации о решаемой задаче на основе результатов расчетов в малой окрестности этого решения. Аналогичный подход был использован в работе [9] и [13] для построения 2-фактор-метода решения вырожденных систем нелинейных уравнений. Подобная идея применялась в [10], [11], [14], [15] для создания численных методов решения вырожденных задач нелинейного программирования. В работах [5], [16] построены точные методы решения задач квадратичного программирования на основе описываемой идеи и имеющие локально полиномиальную сложность. Отметим, что известный метод эллипсоидов, как показано в работе [12], в общем случае не является полиномиальным. В [4] предложен градиентно-проективный метод решения системы (1), который сходится за конечное число итераций и представляет собой объединение итерационного и прямого метода. В данной работе предлагается усовершенствованная схема подобного метода, более простая в своей реализации и обладающая локально полиномиальной оценкой сложности, т.е. если начальная точка принадлежит достаточно малой окрестности некоторой точки решения, то число вычислений, необходимых для получения решения, оценивается как полином от $m$ и $n$, т.е. $O(m \times {{n}^{3}})$. Причем используется идея проектирования на множество активных ограничений в окрестности решения.

Введем $m$ скалярных функций ${{f}^{i}}(x) = \left\langle {{{a}^{i}},x} \right\rangle - {{b}^{i}}$, где ${{a}^{i}} \in {{\mathbb{R}}^{n}}$ есть $i$-я строка матрицы $A$, $i \in D = \{ 1,2, \ldots ,m\} $. Шар радиуса $\varepsilon $ с центром в точке $x$ обозначим через ${{U}_{\varepsilon }}(x)$.

Для произвольных $x$, ${{x}_{*}}$, $i \in D$ выполнены неравенства Липшица

(2)
${{f}^{i}}(x) - \left\| {{{a}^{i}}} \right\|\left\| {x - {{x}_{*}}} \right\| \leqslant {{f}^{i}}({{x}_{*}}) \leqslant {{f}^{i}}(x) + \left\| {{{a}^{i}}} \right\|\left\| {x - {{x}_{*}}} \right\|.$

Определим индексные множества ${{J}_{0}}(x) = \{ i:i \in D,{{f}^{i}}(x) = 0\} $, ${{J}_{ - }}(x) = \{ i:i \in D,{{f}^{i}}(x) < 0\} $, ${{J}_{ + }}(x) = \{ i:i \in D,{{f}^{i}}(x) > 0\} $.

Для решения задачи (1) достаточно найти хотя бы одну точку ${{x}_{*}}$ такую, что множество ${{J}_{ + }}({{x}_{*}})$ пусто, либо, что то же самое, такую, что ${{J}_{0}}({{x}_{*}}) \cup {{J}_{ - }}({{x}_{*}}) = m$.

Лемма 1. Пусть ${{x}_{*}} \in X$, $i \in D$ и для любых $\varepsilon > 0$ существуют $x \in {{U}_{\varepsilon }}({{x}_{*}})$ такие, что ${{f}^{i}}(x) \geqslant 0$. Тогда $i \in {{J}_{0}}({{x}_{*}})$.

Доказательство. Из условия ${{x}_{*}} \in X$ следует, что либо ${{f}^{i}}({{x}_{*}}) = 0$, либо ${{f}^{i}}({{x}_{*}}) < 0$. В первом случае имеет место требуемое утверждение, во втором из (2) получим ${{f}^{i}}(x) \leqslant {{f}^{i}}({{x}_{*}}) + \left\| {x - {{x}_{*}}} \right\|\left\| {{{a}^{i}}} \right\|$ для любых $x \in {{U}_{\varepsilon }}({{x}_{*}})$. Отсюда при $0 < \varepsilon < - {{f}^{i}}({{x}_{*}}){\text{/}}\left\| {{{a}^{i}}} \right\|$, имеем ${{f}^{i}}(x) < 0$, что противоречит условию леммы о том, что на множестве ${{U}_{\varepsilon }}({{x}_{*}})$ есть точки, где ${{f}^{i}}(x) \geqslant 0$. Поэтому заключаем, что $i \in {{J}_{0}}({{x}_{*}})$.

Введем в рассмотрение функцию штрафа

(3)
$\varphi (x) = \frac{1}{2}{{\left\| {{{{(Ax - b)}}_{ + }}} \right\|}^{2}},$
где ${{(Ax - b)}_{ + }}$ – вектор, $i$-я компонента которого равна ${{(Ax - b)}^{i}}$, если ${{(Ax - b)}^{i}} > 0$, и равна нулю, если ${{(Ax - b)}^{i}} \leqslant 0$, $1 \leqslant i \leqslant m$. Очевидно, что

(4)
$X = \mathop {{\text{Argmin}}}\limits_{x \in {{\mathbb{R}}^{n}}} \,\varphi (x).$

Количество элементов, входящих в множество $Z$, обозначаем через $\left| Z \right|$ (его мощность). Для определенности считаем, что ${{J}_{ - }}({{x}_{*}}) = \{ 1, \ldots ,l\} $, ${{J}_{0}}({{x}_{*}}) = \{ l + 1, \ldots ,m\} $.

Если ранг матрицы $B$ размерности $m \times n$ равен $m$, то псевдообратной к $B$ будет $n \times m$ матрица ${{B}^{ + }} = {{B}^{ \top }}{{(B{{B}^{ \top }})}^{{ - 1}}}$. Квадратную матрицу $n \times n$ ортогонального проектирования на пространство строк матрицы $B$ определим как ${{({{B}^{ \top }})}^{\parallel }} = {{B}^{ \top }}{{(B{{B}^{ \top }})}^{{ - 1}}}B = {{B}^{ + }}B$, а проекцию на ортогональное дополнение ${{({{B}^{ \top }})}^{ \bot }} = J - {{({{B}^{ \top }})}^{\parallel }}$, где $J$ – единичная матрица $n \times n$, $n \geqslant m$.

Из индексного множества ${{J}_{ + }}(x)$ выделим подмножество ${{\bar {J}}_{ + }}(x)$ – совокупность индексов, соответствующих максимальному набору линейно независимых строк матрицы $A$ с номерами из ${{J}_{ + }}(x)$. Пусть

(5)
$M(\bar {x}) = \{ z \in {{\mathbb{R}}^{n}}:{{f}^{i}}(z) = 0,i \in {{J}_{ + }}(\bar {x})\} .$
Обозначим через $\bar {A}(\bar {x})$ матрицу, состоящую из строк, соответствующих индексному множеству ${{\bar {J}}_{ + }}(\bar {x})$. Аналогично определим вектор $\bar {b}$. Введем оператор проектирования точек $\bar {x} \in {{\mathbb{R}}^{n}}$ на множество $M(\bar {x})$:
(6)
$z(\bar {x}) = {{P}_{{M(\bar {x})}}}(\bar {x}) = \mathop {\left( {{{{\bar {A}}}^{ \top }}(\bar {x})} \right)}\nolimits^ \bot \bar {x} + {{\bar {A}}^{ + }}(\bar {x})\bar {b},$
где ${{\bar {b}}^{i}} = {{b}^{i}}$, $i \in {{\bar {J}}_{ + }}(x)$. Вычисление указанной проекции предполагает выделение максимальной системы линейно независимых строк матрицы $A$ со строками ${{a}_{j}}$, $j \in {{J}_{ + }}(\bar {x})$. Данный процесс легко описывается как последовательное добавление строк данной матрицы, начиная с произвольной ненулевой строки. При этом новая строка ${{a}_{j}}$ добавляется в том случае, если она не представима в виде линейной комбинации уже добавленных строк. Эта проверка осуществима при помощи стандартных программ линейной алгебры, см., например, https://topepo.github.io/caret/pre-processing.html#indep. Аналогично вычисление псевдообратной матрицы, см. https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/ginv.html.

Опишем алгоритм нахождения решения системы (1) по информации, полученной в текущей точке $\bar {x}$.

АЛГОРИТМ

1. В текущей точке $\bar {x}$ определяем набор индексов ${{J}_{ + }}(\bar {x})$. Если множество ${{J}_{ + }}(\bar {x})$ – пусто, то $\bar {x}$ – решение системы (1).

2. По формуле (6) определим проекцию $z(\bar {x})$ и проверим выполнение неравенств

(7)
${{f}^{i}}\left( {z(\bar {x})} \right) \leqslant 0,\quad \quad i = 1, \ldots ,m.$

3. Если неравенство (7) выполнено для всех $i \in D$, то $z(\bar {x})$ – решение системы (1).

Если для некоторых $i \in D{\backslash }{{J}_{ + }}(\bar {x})$ неравенство (7) нарушается, то выбираем индекс ${{i}_{0}}$ такой, что

(8)
${{i}_{0}} = \mathop {{\text{argmin}}}\limits_{i \notin {{J}_{ + }}(\bar {x})} \,\left| {{{f}^{i}}\left( {z(\bar {x})} \right)} \right|.$
(Минимум может достигаться на нескольких индексах ${{i}_{0}} \notin {{J}_{ + }}(\bar {x})$. В этом случае берем любой из них.)

4. Добавим ${{i}_{0}}$ к множеству ${{J}_{ + }}(\bar {x})$

(9)
${{J}_{ + }}(\bar {x}): = {{J}_{ + }}(\bar {x}) \cup \{ {{i}_{0}}\} $
и перейдем к 2.

Теоpема 1. При достаточно малых $\varepsilon > 0$ и любых $\bar {x} \in {{U}_{\varepsilon }}({{x}_{*}})$. Алгоритм определяет решение ${{z}_{*}} = z(\bar {x})$ системы (1) за число вычислений порядка $O(m \times {{n}^{3}})$.

Доказательство. Идея доказательства состоит в том, что в соответствии с леммой 1 при $\bar {x}$ достаточно близких к ${{x}_{*}}$ ограничения, которые нарушаются, т.е. ${{f}^{i}}(\bar {x}) > 0$, должны соответствовать активным ограничениям в точке ${{x}_{*}}$, т.е. ${{f}^{i}}({{x}_{*}}) = 0$, ${{J}_{ + }}(\bar {x}) \in \left\{ {\ell + 1, \ldots ,m} \right\}$. Нулевые ограничения в точке $\bar {x}$, т.е. ${{f}^{j}}(\bar {x}) = 0$, опять по лемме 1 соответствуют ${{f}^{j}}({{x}_{*}}) = 0$. Остаются только ограничения ${{f}^{i}}(\bar {x}) < 0$, среди которых могут быть соответственно активные в точке ${{x}_{*}}$. Пусть это $\left\{ {\ell + 1, \ldots ,p} \right\}$, т.е. ${{f}^{j}}({{x}_{*}}) = 0$ для $\ell + 1 \leqslant j \leqslant p$. А для остальных ${{f}^{j}}({{x}_{*}}) \leqslant - c$, $1 \leqslant j \leqslant \ell $, и $0 < c$ – некоторая независимая константа. Поэтому ${{f}^{j}}(\bar {x}) \leqslant - c{\text{/}}2$, $1 \leqslant j \leqslant \ell $, а остальные ограничения будут меньше нуля в точке $\bar {x}$, но заведомо ${{f}^{j}}(\bar {x}) > - c{\text{/}}2$, при $\ell + 1 \leqslant j \leqslant p$. Эти ограничения занумеруем по убыванию $\left| {{{f}^{i}}(\bar {x})} \right|$, $i \in \left\{ {1, \ldots ,p} \right\}$ и будем выбирать последовательно индексы

${{i}_{0}} = \mathop {{\text{argmin}}}\limits_{1 \leqslant i \leqslant p} \,\left| {{{f}^{i}}\left( {\bar {x}} \right)} \right|,$
т.е. последовательно выделять в точке ${{x}_{*}}$ подозреваемые на активные ограничения и добавлять их в ${{J}_{ + }}(\bar {x}): = {{J}_{ + }}(\bar {x}) \cup \{ {{i}_{0}}\} $. Таким образом, на некотором шаге будет ${{J}_{ + }}(\bar {x}) = {{J}_{0}}({{x}_{*}})$, так как ограничения неактивные даже не будут проверяться в силу того, что мы получим решение, как только ${{J}_{ + }}(\bar {x}) = {{J}_{0}}({{x}_{*}})$.

На основании вышеизложенного доказательство можно записать в следующем виде.

В соответствии с леммой 1 в множестве индексов ${{J}_{ + }}(\bar {x}) \in \{ \ell + 1, \ldots ,m\} $ и среди ограничений ${{f}^{i}}(\bar {x}) < 0$, $i \notin {{J}_{ + }}(\bar {x})$, могут быть такие, что $i \in {{J}_{0}}({{x}_{*}})$. Пусть это для определенности $i \in \{ \ell + 1, \ldots ,p\} $, т.е. ${{f}^{j}}(\bar {x}) < 0$ для $\ell + 1 \leqslant j \leqslant p$, а для остальных индексов имеем ${{f}^{j}}(\bar {x}) \leqslant - \tfrac{c}{2}$, $1 \leqslant j \leqslant \ell $, $c > 0$, т.е. получаем заведомо допустимые ограничения для неактивных $i \in D$. Очевидно, что при $\varepsilon > 0$ достаточно малых,

(10)
$\left| {{{f}^{i}}(\bar {x})} \right| < \left| {{{f}^{j}}(\bar {x})} \right|,\quad j \in \{ 1, \ldots ,\ell \} ,\quad i \in \{ \ell + 1, \ldots ,p\} ,$
так как $\bar {x} \in {{U}_{\varepsilon }}({{x}_{*}})$ и ${{f}^{i}}({{x}_{*}}) = 0$, $i = \ell + 1, \ldots ,p$. При этом в точке $z(\bar {x}) = {{P}_{{M(\bar {x})}}}(\bar {x})$ некоторые из ограничений могут нарушаться и нужно изменять ${{J}_{ + }}(\bar {x})$ в соответствии с (9).

Отметим, что при вычислении точки $z(\bar {x})$, система ${{f}^{i}}(z) = 0$, $i \in \mathop {\bar {J}}\nolimits_ + (\bar {x})$ и система ${{f}^{i}}(z) = 0$, $i \in {{J}_{ + }}(\bar {x})$ будут совместны, поскольку совместны системы

${{f}^{i}}(z) = 0,\quad i \in \mathop {\bar {J}}\nolimits_ + ({{x}_{*}})\quad {\text{и}}\quad {{f}^{i}}(z) = 0,\quad i \in {{J}_{ + }}({{x}_{*}}),$
а точка $\bar {x}$ принадлежит достаточно малой окрестности ${{x}_{*}}$.

Изменение множества ${{J}_{ + }}(\bar {x})$ происходит за счет добавления ограничений с индексами из множества $\{ \ell + 1, \ldots ,p\} $. Считаем, что множество ${{J}_{0}}({{x}_{*}})$ не пусто, т.к. в противном случае ${{J}_{ + }}(\bar {x}) = 0$ и все ${{f}^{i}}(\bar {x}) < 0$. Тогда число изменений множества ${{J}_{ + }}(\bar {x})$ не больше, чем $p - \ell - 1 \leqslant m$. Поэтому за число проектирований не более, чем $m$, мы получим проекцию ${{z}_{*}} = z(\bar {x}) = {{P}_{{M(\bar {x})}}}(\bar {x})$. Число операций, требуемых для проектирования точки $x$ на $M(\bar {x})$, будет порядка $O({{n}^{3}})$, а число изменений множества ${{J}_{ + }}(\bar {x})$ – не более $m$. Поэтому общее число итераций для отыскания ${{z}_{*}}$ будет порядка $O(m \times {{n}^{3}})$.

Остается описать способ нахождения точки $\bar {x}$ из достаточно малой окрестности ${{U}_{\varepsilon }}({{x}_{ * }})$ некоторого фиксированного решения ${{x}_{*}} \in X$. Для этого используем простейший градиентный метод

(11)
${{x}_{{k + 1}}} = {{x}_{k}} - \alpha \varphi {\text{'}}({{x}_{k}}),$
где градиент $\varphi {\text{'}}(x) = {{A}^{ \top }}{{(Ax - b)}_{ + }}$ удовлетворяет условию Липшица с константой $L = \left\| {{{A}^{ \top }}A} \right\|$. Справедлива следующая

Теоpема 2. Пусть ${{x}_{0}} \in {{\mathbb{R}}^{n}}$ и последовательность $\{ {{x}_{k}}\} $ строится по схеме (11) с постоянным шагом $\alpha = 1{\text{/}}2L$. Тогда ${{x}_{k}} \to {{x}_{*}} \in X$ при $k \to \infty $ и $\left\| {{{x}_{{k + 1}}} - y} \right\| < \left\| {{{x}_{k}} - y} \right\|$ для любого $y \in X$ и $k = 0,1, \ldots $ .

Доказательство можно найти, например, в [1].

Таким образом, метод (11) реализует последовательность, которая сходится к некоторой точке ${{x}_{*}} \in X$, т.е. заведомо для каждого $\varepsilon > 0$ достаточно малого, начиная с некоторого $\bar {k} = k(\varepsilon )$, все $\{ {{x}_{k}}\} \in {{U}_{\varepsilon }}({{x}_{*}})$, $k \geqslant \bar {k}$. Это означает, что на некотором шаге с номером $\bar {k}$ выполнены условия теоремы 1 и получается решение системы (1).

Опишем окончательно метод решения системы (1).

Теоpема 3. Выполняются следующие действия.

1. Полагаем $k = 0$ и берем в качестве ${{x}_{0}}$ любую точку из ${{\mathbb{R}}^{n}}$.

2. Определяем итерацию ${{x}_{{k + 1}}} = {{x}_{k}} - \alpha \varphi {\text{'}}({{x}_{k}})$, $k = k + 1$.

3. Вычисляем по алгоритму точку $z({{x}_{k}})$.

4. Проверяем: будет ли $z({{x}_{k}})$ решением системы (1):

− если нет, то переходим на п. 2.

− если да, то найдено конечное $\bar {k}$ такое, что $z({{x}_{{\bar {k}}}})$ – решение системы (1).

Доказательство. Поскольку последовательность $\{ {{x}_{k}}\} $ сходится к фиксированной точке ${{x}_{*}} \in X$, то на некотором шаге $\bar {k}$ будут выполнены условия теоремы 1 и получится решение ${{z}_{*}} = {{P}_{{M({{x}_{{\bar {k}}}})}}}({{x}_{{\bar {k}}}}) \in X$.

Замечание. Если в алгоритме нарушается условие совместности системы (5), то считаем, что алгоритм не дает решения и переходим к п. 2.

Отметим, что требование непустоты множества $X$, введенное выше, можно опустить. Тогда описанный метод находит псевдорешение системы (1).

Список литературы

  1. Карманов В.Г. Математическое программирование. М.: Физматлит, 2000.

  2. Голиков А.И., Евтушенко Ю.Г. Теоремы об альтернативах и их применение в численных методах // Журн. вычисл. мат. и матем. физ. 2003. Т. 43. № 3. С. 354–375.

  3. Evtushenko Y.G., Golikov A.I. New perspective on the theorems of alternative // High Performance Algorithms and Software for Nonlinear Optimization. Springer US, 2003. P. 227–241.

  4. Tret’yakov A.A. A finite-termination gradient projection method for solving systems of linear inequalities // Russian Journal of Numerical Analysis and Mathematical Modelling. 2010. V. 25(3). P. 279–288.

  5. Tret’yakov A., Tyrtyshnikov E. Exact differentiable penalty for a problem of quadratic programming with the use of a gradient-projective method // Russian Journal of Numerical Analysis and Mathematical Modelling, 2015. V. 30 (2). P. 121–128.

  6. Smale S. Mathematical problems for the next century, Mathematics: Frontiers and perspectives, 271–294 // American Mathematical Society, Providence, RI. 2000.

  7. Mangasarian O.L. A finite Newton method for classification problems. Technical Report 01-11, Data Mining Institute, Computer Sciences Department, University of Wisconsin, Madison, Wisconsin, 2001. C. 01–11.

  8. Mangasarian O.L. A Newton method for linear programming // J. Optimizat. Theory and Applicat. 2004. T. 121. № 1. C. 1–18.

  9. Белаш К.Н., Третьяков А.А. Методы решения вырожденных задач // Ж. вычисл. матем. и матем. физ. 1988. Т. 28. № 7. С. 1097–1102.

  10. Brezhneva O.A., Tret’yakov A.A. The p-factor-Lagrange methods for degenerate nonlinear programming // Numerical Functional Analysis and Optimization. 2007. T. 28. № 9–10. C. 1051–1086.

  11. Брежнева О.А., Евтушенко Ю.Г., Третьяков А.А. $2$-фактор метод модифицированной функции Лагранжа для решения вырожденных задач условной оптимизации // Доклад. АН. 2006. Т. 408. № 4. С. 439–442.

  12. Goffin J.L. On the non-polynomiality of the relaxation method for systems of linear inequalities // Mathematical Programming. 1982. T. 22. № 1. C. 93–103.

  13. Facchinei F., Fischer A., Kanzow C. On the accurate identification of active constraints // SIAM Journal on Optimization. 1998. T. 9. № 1. C. 14–32.

  14. Szczepanik E., Tret’yakov A. p-factor methods for nonregular inequality-constrained optimization problems // Nonlinear Analysis. 2008. T. 69. № 12. C. 4241–4251.

  15. Wright S.J. An algorithm for degenerate nonlinear programming with rapid local convergence // SIAM Journal on Optimization. 2005. T. 15. № 3. C. 673–696.

  16. Tret’yakov A., Tyrtyshnikov E.E. A finite gradient-projective solver for a quadratic programming problem // Russian Journal of Numerical Analysis and Mathematical Modelling. 2013. T. 28. № 3. C. 289–300.

  17. Han S. P. Least-Squares Solution of Linear Inequalities. WISCONSIN UNIV-MADISON MATHEMATICS RESEARCH CENTER, 1980. № MRC-TSR-2141.

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал