Быстрый приближенный алгоритм для задачи положительного линейного программирования Текст научной статьи по специальности «Математика»

Быстрый приближенный алгоритм для задачи положительного линейного программирования1

С. А. Фомин

Аннотация. Предлагается новый алгоритм нахождения е-оптимального решения задачи положительного линейного программирования, где все исходные данные неотрицательны — тах^а^А-с < Ь, х > 0}. Алгоритм имеет лучшую верх-

нию оценку (из известных алгоритмов для этой задачи) вычислительной сложно-/ N 1с^ \ ^

сти: О I-^ £ 1, где т х п — размер матрицы ограничении, а — количество

ненулевых элементов и имеет простую реализацию.

1. Введение

В этой работе будем рассматривать задачи положительного линейного программирования (ПЛП), также называемые задачами дробной упаковки, когда требуется найти вектор х = {х\,..., х„ ), такой что:

сх —> тах, Ах < Ь, > О,

где А — рациональная т х п матрица, ац > 0, с — рациональный вектор, С] > О, Ь — рациональный вектор, 6,; > 0, и задачи дробного покрытия (когда линейные ограничения есть ограничения вида Ах > Ь, причем все данные неотрицательны, и необходимо минимизировать целевую функцию).

Многие задачи оптимизации можно представить в виде задачи ПЛП, в частности задачи оптимизации потоков в сетях [1, 2]. Несмотря на существование полиномиальных алгоритмов решения задач линейного программирования [4, 5], рост размерности задачи, коммуникационные аспекты, возможность параллельной реализации подталкивают к поиску приближенных алгоритмов для решения задач ПЛП. В последние годы было разработано много ¿-приближенных ПЛП-алгоритмов [6, 7, 8, 1, 2, 9, 10].

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

1 Работа выполнена при поддержке РФФИ, проект 02-01-00713.

оценка вычислительной сложности О £ ^ , где N — число ненулевых

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

2. Постановка задачи ПЛП

Рассмотрим постановки задач ПЛП:

Задача 1. ПЛП типа упаковки

Имеется неотрицательная матрица ограничений Л: т строк, п столбцов, и- • О. положительные вектора с = (ci,...,cn), Cj > 0 и b = (bi,...,bm), bi >0.

п

CjXj —» max

3 =1 п

Vi Е ciijXj < bi

з=i

Vj Xj > 0.

Задача 2. ПЛП типа покрытия

Имеется неотрицательная матрица ограничений Л: т строк, п столбцов, п:; -П. положительные вектора с = (ci,...,cn), Cj > 0 и

Ъ = (bi,..., bm), bi > 0.

П

CjXj —> min

j=i

п

Vi Е ^ I j j х j bi

j=i

Vj Xj > 0.

Далее, везде n — это число переменных, m — число ограничений прямой задачи (не считая ограничений вида ж,; >0), а АГ — число ненулевых элементов матрицы ограничений. Мы будем использовать j для индексирования переменных прямой задачи и столбцов матрицы ограничений, и, где не будет явно указано обратное, будем считать, что j изменяется в диапазоне 1... п. Также мы будем использовать г для индексирования ограничений, и, соответственно, по умолчанию считаем, что г € [1... т].

Заметим, что задачи 1 и 2 можно представить в так называемой стандартной форме (см. [1]), где векторы стоимостей с и правых частей ограничений b состоят из единиц. Эта форма интересна тем, что одна матрица

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

Задача 3. Прямая задача ПЛП

Дана неотрицательная матрица ограничений Л: т строк, п столбцов, а .у > 0.

X = Ху —> шах

з

У г Л. = ^ ^ ' 1

з

^3 хі — 0-

Задача 4. Двойственная задача ПЛП

Дана неотрицательная матрица ограничений Л: т строк, п столбцов, а .у > 0.

V" = у.і —> П1ІП

г

\¡j ОС ^ ^ ац у і > 1

і

Уг уі > 0.

Сразу заметим, что ПЛП-задача 1, легко преобразуется в ПЛП-задачу 3, посредством преобразования матрицы ограничений:

(і \

а>з = -ГГ- (!)

Очевидно такое преобразование не повлияет на значение целевой функции. А решения х* задачи 3 и х* задачи 1 соотносятся очевидным образом:

X*

хз = —■ сз

Ниже (см. параграф 2.1) мы покажем, что без потери общности мы можем считать, что для коэффициентов матрицы ограничений входной задачи ПЛП выполняется атіп < а.у < атах, где = 0((^)С0П8‘),

I = тіп(то, п).

В случае, если имеется экспоненциальный разброс в коэффициентах, то можно провести округление, в результате которого, для данного т, а'тах = О(-Ц-), а значение решение округленной задачи 3 будет не меньше

0>тъп 'Г

(1 — г), умноженного на значение решения неокругленной задачи.

2.1. Округление

Обозначим а* = шах,; а.у, а* = тих,- а*, I = ппп(то, п).

Пусть апшх = . Тогда для всех а* > атах, зафиксируем Xj = 0.

Заметим, что значение оптимального решения исходной задачи X* > ^¡г, что следует из допустимости решения Хк = -¿г, где к = ащпшг^а*.

Для потерь из-за округления АХ с одной стороны АХ < ——, с другой АХ < т . Получаем:

0"тах

АХ < —— < -^ < -X*. (2)

&тах 2(2* 2

Теперь, пусть атЛп = Для всех ау < атЛп, положим ау = а„„:п. Заметим также, что X* = У* < ^¡г, что следует из допустимости двойственного решения V?' € Б у, = -¿г, где 5'-минимальное множество ограничений, в котором участвуют все переменные: |5| < то., |5| < п —> |5| < I.

Оценим максимальное увеличение левых частей округленных ограничений, при подстановке оптимального решения неокругленной задачи 3:

та* I т

ААтаг, ^ X • а п!/п ^ , ~ ^ •

21 а* 2

Отсюда для значения округленного оптимума будет

X' > 7—Т- (3)

1 2

Из (2) и (3) имеем, что задача, с коэффициентами, усеченными по заданным ашах и аШгП, имеет решение допустимое для исходной задачи 3, и

для его значения выполняется

X' > 1 ~Т^2Х* > (1 - т)Х*.

- 1 + т/2 - 1 '

3. Обзор алгоритмов для приближенного решения задачи ПЛП

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

Определение 2. Алгоритм с мультипликативной точностью (1+ь) называется £—оптимальным (подразумевается, что £ близко к нулю).

Термин є-оптимальное решение используется для обозначения допустимого решения со значением целевой функции, отличающемся от оптимума не более, чем в (1 + є) раз.

В [7] был описан t-оптимальный алгоритм для сопряженных задач 1 и

2 состоящий из 0( т”/£) итераций сложности О (піп).

Затем в [1] был описан приближенный алгоритм который, что для заданных £ > О И 0 < Г < 1п( ^р2~ ) находил допустимые решения X и у, для сопряженных задач 3 и 4 соответственно, что Xj > г_^^'Еу • При г « є алгоритм можно было рассматривать, как t-оптимальный.

Алгоритм также состоял из итераций сложности 0(тп), причем в ста,-

I , і .. / і. \' іу 2(' пі) ' т п /;- ) , m2

тье [1] число итерации оценено, как (J( w—'rggv ), гДе 1 = что

при ?’ « t эквивалентно 0( 1о° Е ^). К сожалению, в анализе времени выполнения, приведенном в [1] нами обнаружена ошибка, и согласно нашим

.. ґ / . і 1 / ¡j (' ■ m j2 11 j_' (' ■ m n /;-) ,

оценкам, правильная оценка числа итерации будет (J( ———' г2°2 ------—) или

і 3 / ТГ1TL \

(при г « t) 0( °ё 4 £ ), что асимптотически эквивалентно числу итераций

алгоритма из [7].

В работе [2] был представлен t-оптимальный алгоритм решения задачи плп, для которого доказана верхняя оценка числа итераций 0(р- log(m)), каждая итерация сложности 0(N), (напомним, что N обозначает число ненулевых элементов в матрице ограничений задачи ПЛП). Заметим, что при фиксированном є алгоритмы из работ [7, 1] имеют полилогарифмиче-ские верхние оценки числа итераций, а алгоритм из [2] имеет полиномиальную верхнюю оценку, причем существуют входные данные, на которых эта оценка достигается.

В работе [9] был представлен t-оптимальный алгоритм для сопряжен-

1 О » ҐЛ! ІО 22(піп/є)\

ных задач 1 и 2 временной сложности (J(mn є2 )•

В статье [10] рассмотрен более общий класс задач — приближенное решение смешанной задачи типа покрытия и упаковки2:

Задача 5. Смешанная ПЛП задача типа покрытия и упаковки:

Даны, неотрицательные матрицы Р,С, вектора р,с, е Є (0,1), необходимо найти допустимый вектор х > 0, для которого

Рх < (1 + е)р

Сх > с.

В [10] представлен алгоритм для решения задачи 5 с вычислительной

2Approximate Mixed Packing and Covering

сложностью 0(т Однако, хотя теоретически задача 1 сводится к

задаче 5, представленный в [10] алгоритм сведения достаточно сложен.

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

г* ( N log ^ \

асимптотической оценкой сложности — U I —g2 g I.

Уі

X

Y

Y

Y

Axj

\i

Amax

aj

Ot-min

Переменные задачи 3 Переменные задачи 4 Значение целевой функции задачи 3 Значение целевой функции задачи 4 Верхняя граница для У

Верхняя оценка для значения оптимума задач 3 и 4 Дискретные приросты переменных задачи 3 Значение левых частей ограничений задачи 3

Максимальное значение Л^. -----допустимое решение задачи 3.

Значение левых частей ограничений задачи 4

Минимальное значение а ,

— допустимое решение задачи 4.

Параметры сходимости. Необходимое условие: £ + [i + < є.

Вход: т,п, е Є (0,1) неотрицательная т х п матрица А. Выход: х, у — е-оптимальные решения задач 3 и 4.

£ = At'!— |> X = (1 + 0(1 + А4)

(l+e) log ?71+g\

Y = e 1+6-x

Vi Д.т, <- « = 5 J a. max,-3 г an

Vi xj <- 0, X <- 0 {X

Vi Уі <-1, Y m {Vi

Vj (Xj <- E іаіїУі &тгп minj OLj

— „А,-

repeat

for all j : — >

J (У.----

xj <- xj + Ахз

Ej aijxj

Vi уі <- eXi,

(1+aO

do

Vi Аг

Vi y, V? o. j Y <- -

a

end for until (y^-

return "T

At7

■ <— maxj Ai Ег Vі

.. <— Inin ay

X

E,

> _1_ or у > У)

Алгоритм 1: Упрощенная идея алгоритма PLPAPX

4. Новый алгоритм для приближенного решения задачи ПЛП

Поясним основные идеи предлагаемого алгоритма (см. алгоритмы 1,2). Алгоритм является усовершенствованием алгоритма из [9], но его идея алгоритма близка к идеям алгоритмов из [7, 1, 2], содержащих разные формы применения релаксаций Лагранжа, и основана на последовательном увеличении переменных прямой задачи-упаковки в зависимости от переменных двойственной задачи-покрытия, соответствующих ограничениям прямой задачи, т.е. алгоритм одновременно ищет решение прямой и двойственной задачи.

Сначала поясним основную идею алгоритма (см. алгоритм 1). Получив на вход матрицу ограничений ПЛП-задачи, алгоритм инициализирует параметры (/л, £, Аху), определяющие его поведение, и обнуляет значения переменных задачи-упаковки х.

Далее, основная работа алгоритма состоит из последовательного увеличения переменных х, и перерасчета вектор-переменных А, у, а (в указанном порядке). На основе значений а вновь выбираются переменные для увеличения, и цикл повторяется.

Выполнение алгоритма прерывает достижение оптимальности. Действительно, если найдены решения прямой и двойственной задачи, отношение значений которых меньше требуемого мультипликативного фактора (1 + е), то, т.к. оптимумы прямой и двойственной задачи совпадают, имеем

> ——У* = ~^—х*.

Лтах (1+е) 1 + £ 1 + £

Либо, в случае если У > У, оптимальность гарантируется леммой 2.

Теперь рассмотрим алгоритм 2. Основная проблема алгоритма 1, это неоптимальный перерасчет зависящих от х переменных, при увеличении каждого компонента Xj. Зависимые переменные можно подразделить на две группы. В первую войдут переменные А, Хтах, у, X, У, которых можно перерасчитать при увеличении компоненты х^, за время 0(\1^\). Тем самым, вклад сложности этих операций в общую вычислительную сложность алгоритма можно оценить как 0(|Х/|) на число увеличений х^ (см. лемму 4).

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

Поэтому, глобальный пересчет а, имеющий сложность в худшем случае

О(ЛГ), происходит не более двух раз внутри внешнего цикла алгоритма, и его вклад в сложность алгоритма можно оценить как О (Ж) на число внешних циклов (см. лемму 3).

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

Напомним, что если проводилось округление с параметром т (см. параграф 2.1) то могла появится мультипликативная погрешность не превышающая (1 — т) (см. параграф 2.1). Поэтому, если входной параметр точности е, параметр округления т, то алгоритм 2 следует запускать с параметром точности е = £ — т — £т, чтобы компенсировать погрешность округления. Тогда мы получим, что мультипликативная точность, с учетом округления:

1 — т 1 — т 1

1 + 6 ( 1 — 1~)( 1 + ^ ) 1 + Ь

Выбор т < £, т.е. распределение допустимой погрешности между процедурой округления и самим алгоритмом 2 предоставляется на усмотрение реализатора. Например, если известно, что матрицы не нуждаются в округлении, то можно положить е = £.

Кроме этого при реализации можно выбрать параметры £ и /х, при условии, что £, ¡1 > 0 и ц+^+ц/х < е. И все же, хотя выбор конкретных значений параметров £, /л и т может оказать существенное влияние на эффективность реализации, все оценки сложности и ¿-оптимальность алгоритма доказаны для произвольных значений параметров £, /л и т, удовлетворяющих вышеупомянутым условиям.

Ниже мы приведем оценки вычислительной сложности нового алгоритма и доказательство оптимальности.

Сначала небольшая вспомогательная лемма. Из определения алгоритма

2 непосредственно следует:

Лемма 1. При увеличении любого хVI ДА,; < £ и Зг ДА,; = £.

Теперь оптимальность:

Лемма 2. Если при выполнении алгоритма выполнилось условие У > У , то найдено е-оптимальное решение.

Доказательство. Выведем оценку зависимости У от X. Назовем шагом, увеличение компоненты х в алгоритме. Пусть Ь номер некоторого шага,

Аз

С^тгп

СХ-тгп

{г : ац >0} индексы строк с ненулевым коэффициентом для столбца у Нижняя оценка для тіп^ а,^. Нужна для выбора увеличиваемых х^. Верхняя оценка для тіп^ а,у. Нужна для выбора пересчитываемых .

Вход: т,п, е £ (0,1) неотрицательная т х п матрица Л.

Выход: х, у — €-оптимальные решения задач 3 и 4 соответственно.

£ = а4'!—|> х = (1 + 0(1 + м)

(1 + e) logm + gx

У = е 1+*-\

V? Ах3 ^ 4r =

V/ Xj <— О Vi А і <— О,

Vi yi <— 1,

Vj o,j <— 'YLi aij>

repeat

for all j : — >

J a} —

a3 Eiai]Vi

Y

maxj a^j

Amax ^ max* A^

Y m

Qtmin * minj otj

Y

(і+д)

do

while >

Y

OLj — (1 + m)

; + Аж 7,

and Y < Y do

X <- X + Аж7

{ПрирОСТ Xj}

Уг{еАХі ~ 1),

for all i G /, do

AA8 <

Ath <■ ai aj + °,ij Aj/,

end for

Otmin — Olj • end while end for

Y? G {j : (Xj < amin(l + m)} aj <-until ( Y^m°J < (1 + e) or y > y) return -v-25—.—2—•

А і + AA і ■ Уі + Аг/,,

Arr

У -

r <— max(A„ ■ У + Аг/г

■ Ei aij'yiі <>n

Y =

Алгоритм 2: Алгоритм PLPAPX

а )— индекс увеличившейся на этом шаге компоненты х.

уЕ = 1>^ = ХУ~1еАА<-

г г

Согласно лемме 1 ДА,; < £ < 1. Используя неравенство ех < 1 + (1 + х)х верное для 0 < х < 1, получаем

¥<¥ь < (1 + (1+^)ЛА0

= + (1 +

= ¥ь-г + (1 + ^у^ац

= Vе-1 + {1 + 0Ах^а1^-1.

С другой стороны, для любого шага Ь и выбранного на нем индекса ) должно выполнятся

у£—1 у у* аЕ 1 1 + /л 1 + Ц 1 + Ц (4)

(напомним, что X* = У* — значение оптимального решения), или

„1.-1 > У'~Н1+!>) > X• (5)

Получаем по индукции:

¥ь < Vе-1 ^ ^ + ^1 + А^') ^ г (1 + 6)(1 + ^) д . < 1е х* 3 < У е X* Х (1+€)(1+м) у = те х* ^ .

Логарифмируя последнее выражение, получаем:

X* ^ (1 + 0(1 +М) Л X ~ log(Уi/то) log(Уi/то)

Оценим \%гах. Во-первых, заметим, что на шаге Ь — 1:

1 = ^ОёУ^1 < г

Во-вторых, учитывая лемму 1, имеем

Ата.ж < А+ £ < 1<^ У + £.

36

Итак, оценим отношение оптимума X* к значению решения \г~'-

•У 'А(, „ Х*(1оёУ + 0 X - X

(!+£)(!+ И)

< (iogy + o

< (iogy + o

\og(YL/т) X

log(y) — log m 1 + e) log m + + £ (1 + e - )

-\(1 + e)(logm + Q =

.\(bgm + 0 1 +£j'

□

4.1. Оценка вычислительной сложности

Сначала оценим вычислительную сложность пересчета а^, обозначим ее Со,.

Лемма 3. Са = 0Г°^ )

Доказательство. Оценим число внешних циклов. Очевидно, что значение с каждой итерацией внешнего цикла должно увеличиваться не менее чем в 1 + ц раз. Учитывая, что aEvin < amaxYL и a°min > amin, а верхняя оценки сложности пересчета — 0(N), получаем

C« = 0(iVlog1+/1^)=o(ivilog^)

\ min / \ г1 min J

= О (N-log \ () f

У ^ 11 III /п J \ ^ ^

Теперь оценим вычислительную сложность остальных операций, обозначим ее Сх.

Лемма 4.

/ log^i

Сх = О N

Доказательство. Сначала оценим максимальное число увеличений каждой из переменных Xj. Заметим, что для каждого ) существует по крайней мере один А,;, который будет увеличиваться на £, при каждом увеличении Xj. При этом сложность пересчета переменных после увеличения х^ будет 0(1з). Таким образом,

С, = £0(7;)^ = о = О ( N

□

Из лемм 3 и 4 получаем оценку общей сложности алгоритма.

Теорема 1. Вычислительная сложность алгоритма равна О (^Nl°°e2E ) .

5. Вычислительные эксперименты

Ниже мы приведем результаты ряда вычислительных экспериментов (см. таблицу 1).

Тестовыми данными были задачи линейного программирования, где векторы стоимости и ограничений были заполнены единицами: с = (1,...,1), b = (1,..., 1), а 0/1-матрица ограничений порождалась случайным размещением заданного количества ненулевых элементов (единиц).

Время решения тестовых задач (в секундах) приводится в сравнении с временем точного решения этих задач алгоритмами симплекс-метода и метода внутренней точки, реализованных в пакете GLPK [11].

По результатам [12] регулярного тестирования различных пакетов линейной оптимизации, пакет GLPK является лидером среди некоммерческих пакетов, и несильно уступает ведущим коммерческим программам для линейной оптимизации.

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

Информация о процессоре используемом в тестировании:

Model : AMD Athlon(tm) XP 2000+

Co-Processor (FPU) : Built-in

Строк Столбцов Ненулевых Симплекс метод Метод ВТ PLPAPX е = 0.1

1000 1000 00200000 207.153 418.235 1.297

1000 1000 00300000 218.901 495.268 1.297

1000 1000 00400000 311.504 591.932 0.953

1000 1000 00500000 336.478 662.848 0.844

1000 4000 00400000 432.328 1097.92 19.508

1000 4000 00800000 896.733 1820.3 14.303

1000 4000 01200000 1380.57 2515.66 14.752

1000 4000 01600000 1962.68 3068.66 16.438

1000 4000 02000000 2060.53 3870.04 16.203

1000 7000 00700000 599.666 1961.27 38.157

1000 7000 01400000 1176.04 4070.44 37.594

1000 7000 02100000 1827.29 6173.1 52.401

1000 7000 02800000 2063.51 7697.57 60.343

1000 7000 03500000 3157.05 10450 58.719

4000 1000 00400000 869.033 * 6.677

Таблица 1: Результаты вычислительных экспериментов

Speed : 1.67GHz

Performance Rating : PR2427 (estimated)

Литература

[1] Yair Bartal, John W. Byers, and Danny Raz. Global optimization using local information with application to flow control. In IEEE FOGS, pages 303-311, 1997.

[2] N. Garg and J. Ivoenemann. Faster and simpler algorithms for multicommodity flow and other fractional packing problems. In IEEE FOGS, pages 300-309, 1998.

[3] Luca Trevisan and Fatos Xhafa. The parallel complexity of positive linear programming. Parallel Processing Letters, 8(4):527-533, 1998.

[4] L.G. Ivhachiyan. A polynomial-time algorithm for linear programming. Soviet Math. Dokl., 20( 1): 191-194, 1979.

[5] N. Ivarmarkar. A new polynomial-time algorithm for linear programming. Combinatorial, 4:373-396, 1984.

[6] S. Plotkin, D. Shmoys, and E. Tardos. Fast approximation algorithms for fractional packing and covering problems. In Proc. 32nd IEEE FOCS 91, pages 495-504, 1991.

[7] M. Luby and N. Nisan. A parallel approximation algorithm for positive linear programming. In Proc. of 25th ACM STOC, pages 448-457, 1993.

[8] M.D.Grigoriadis and L.G.Ivhachiyan. Fast approximation schemes for convex programs with many blocks and coupling constraints. SIAM J. Optimization, 4:86-107, 1994.

[9] C.A. Фомин. Новый приближенный алгоритм для решения задачи положительного линейного программирования. Дискретный анализ и исследование операций. Серия 2., 8(2), 2001.

[10] Neal Е. Young. Sequential and parallel algorithms for mixed packing and covering. In IEEE Symposium on Foundations of Computer Science, pages 538-546, 2001.

[11] Andrew Makhorin. GLPIv Reference Manual.

http://www.gnu.org/software/glpk/glpk.html, 2003.

[12] Hans Mittelmann. Benchmarks for Optimization Software.

Department of Mathematics and Statistics, Arizona State University, http://plato.la.asu.edu/bench.html, 2003.