Полностью целочисленный разреженный симплекс-метод Текст научной статьи по специальности «Математика»

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

УДК 519.6

ПОЛНОСТЬЮ ЦЕЛОЧИСЛЕННЫЙ РАЗРЕЖЕННЫЙ СИМПЛЕКС-МЕТОД

© 2012 г. Н.Ю. Золотых, В.К. Кубарев

Нижегородский госуниверситет им. Н.И. Лобачевского nikolai. zolotykh@gmail. com

Поступила в редакцию 20.01.2012

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

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

Введение

Использование арифметики с плавающей запятой при решении задач линейного программирования с целыми или рациональными входными данными может оказаться неприемлемым. В этом случае используют точные варианты симплекс-метода, оперирующие только с целыми числами (или рациональными числами, представленными как отношение двух целых). В «наивной» реализации данного подхода для предотвращения экспоненциального роста промежуточных коэффициентов используется вычисление НОД и дальнейшее сокращение на НОД. Дж. Эдмондс [1] предложил вычислительную схему без многочисленных вычислений НОД, использующую вместо симплекс-таблицы введенную им целочисленную Q-матрицу. Данный подход гарантирует не более чем полиномиальный рост промежуточных коэффициентов. По существу, метод Эдмондса является адаптацией известного «полностью целочисленного» метода решения систем линейных уравнений, основанного на тождестве Сильвестра [2-7]. В [1] замечено, что эту схему можно применить и к модифицированному симплекс-методу (симплекс-методу с обратной матрицей). Реализация данного подхода описана в [8]. Заметим, что алгоритм, предложенный в [8], никак не использует возможную разреженность задачи, хотя именно такие задачи чаще всего встречаются на практике.

В настоящей работе предлагается вычислительная схема, основанная на симплекс-методе,

сохраняющая целочисленность промежуточных коэффициентов и предназначенная для больших разреженных задач. Вместо хранения присоединённой матрицы в явном виде (как это сделано в [8]), мы запоминаем только коэффициенты матричных преобразований, подобно тому, как это делается в симплекс-методе с мультипликативной формой обратной матрицы [9]. Вычислительный эксперимент на задачах из библиотеки NETLIB [10] показал преимущество нашего подхода.

1. Задача линейного программирования

Рассмотрим задачу линейного программирования (ЗЛП) в следующей форме:

z = cx ^ min , (1)

Ax = 0, (2)

l < x < u, (3)

где x - и-мерный вектор-столбец переменных, A

- m x n матрица, c - и-мерная вектор-строка, l, u

- и-мерные векторы-столбцы. Для компонент вектора и допустимы значения + да, а для компонент вектора l - значения - да.

Не нарушая общности, будем считать, что ранг матрицы A равен m. Пусть Ъ и Ж - упорядоченные множества, такие, что Ъ и Ж =

= {1,2, ,и}, и B = Аъ - невырожденная m x m

матрица, называемая базисной, а Ж = Аж - mx(n —

- m) матрица (здесь и далее через Aj обозначена подматрица матрицы A, составленная из столбцов с номерами из I). Матрица B называется базисной, а множество Ъ - базой.

Теперь систему ограничений (2) можно переписать в виде Вхъ + = 0, откуда получаем

Хъ = -В~1^к. Определим

и = съ В-1, d = ск - cъB~lN = ск - uN,

откуда

г = (ск - съВГ1Щхк = (сж - иЩхх= dxк. Пусть множество номеров небазисных переменных Ж разбито на подмножества:

Ж = У и £ и "И и такие, что I) = - да и иу = + да для у є У, 1у Ф - да для у є £, и ф + да для у є "И, 1у = и для у є £. Базисным решением задачи (1)-(3), соответствующим базе Ъ, называется вектор х, такой, что

ХУ =

0, У є У,

1у, У є А

иу, У є и,

Л, У є ^

хъ = -В 1^х.

ХУ =

Базисное решение х называется (прямо) допустимым, если х удовлетворяет условиям (3). Очевидно, что для допустимости базисного решения достаточно потребовать, чтобы 1у < Ху < ц для у е Ъ. Базисное решение называется двойственно допустимым, если

dj = 0 для у е У, dj > 0 для у е £,

4 < 0 для у е и.

Хорошо известно, что если базисное решение одновременно прямо и двойственно допустимо, то оно оптимально.

2. Прямой модифицированный симплекс-метод

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

- это вычислительная схема, реализующая итерации симплекс-метода без явного вычисления всех элементов симплексной таблицы. Вместо этого используется некоторое представление матрицы Б1. Модифицированный симплекс-метод для общей задачи линейного программирования предложен Канторовичем и Залгалле-ром [11] и Данцигом с соавторами [12]. Метод учета двусторонних ограничений рассмотрен в [13]. В настоящем разделе мы описываем модифицированный симплекс-метод для задачи (1)-(3), следуя, в основном, [14].

Пусть известны Ъ = ^Р1,Р2,...^т^ и Ж = У и £ и и V, и $ - множества номеров базисных и небазисных переменных соответственно, определяющих некоторое начальное допустимое базисное решение задачи (1)-(3). Прямой модифицированный симплекс-метод решения задачи (1)-(3) включает следующие основные шаги.

Алгоритм 1. Модифицированный симплекс-метод

Шаг 1. Вычислить текущие значения небазисных и базисных переменных:

0 для у е У,

Л для у е £,

. „, хъ = -Б Шж.

для у е V.,

для у е £,

Шаг 2. Вычислить 4 = с к - съ Б~гМ.

Шаг 3. Найти 5 е Ж, такое, что 4, ^ 0 и

5 е У, или 4, < 0 и 5 е £, или 4, > 0 и 5 е и. Ес-

ли таких 5 нет, то СТОП - текущее базисное решение оптимально.

Шаг 4. Вычислить 5-й столбец ~ матрицы А: ~ = Ч ) = -Б.

Шаг 5. Вычислить 0г- (г = 1,2,.. ,,т):

а) если < 0, то

6, =

+ да,

мр, — ХР,

б) если ds > 0, то

+ да,

а. = 0,

Хр, мр,

~ . аі* <0.

I аь

Вычислить

0тт = тіп 0^ 0 = min{0mm,и— К}.

I

Если 0 = +да, то СТОП - минимум ЗЛП (1)-(3) не достигается.

Если 0 = 0р, то в Ъ = ^Р1,Р2,...,Р т) заменяем Рр на 5.

Если 0 = и.; - 4 < 0і (і = 1,2,.. ,,т), тогда в случае 5 є £ перемещаем 5 из £ в V., а в случае ; є V. перемещаем 5 из V. в £.

После выполнения указанных действий снова возвращаемся на шаг 1.

Каждое последовательное выполнение шагов 1-5 называется итерацией симплекс-метода.

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

3. Мультипликативная форма обратной матрицы

Один из способов представления обратной матрицы Б-1 - так называемая ее мультипликативная форма [9].

Рассмотрим процедуру замены столбца в матрице Б, проводимую на ^й итерации модифицированного симплекс-метода: 5-й столбец а5 матрицы А вводится на место р-го столбца матрицы Б. Нетрудно видеть, что в результате получаем матрицу Б' = БF, где

F = [е1 , ^ ..., ep-1, ~ , ep+1, ..., ет ], ~ = Б-1 ач , ег - г-й столбец единичной матрицы, откуда (Б')1 = ИБТ\ где

Н = ^ =[e1, ^ , ер+1,...,ет ] , (4)

, І Ф p,

a„.

i = P.

Алгоритм 2. Целочисленный модифицированный симплекс-метод

Шаг 1. Значения небазисных переменных вычисляем, как и в алгоритме 1. Вместо значений базисных переменных вычисляем текущие значения базисных переменных, умноженные на абсолютное значение определителя базисной матрицы: Дхв = -B+Nxk. Полученный вектор будет целочисленным.

Шаг 2. Вместо вектора d вычисляем Д^ = = Д^ск - csB N. Полученный вектор также будет целочисленным.

Шаг 3. Данный шаг такой же, как и в алгоритме 1, только вместо d используем Д-d.

Шаг 4. Вместо столбца as вычисляем Да = = (ais) = -B+as.

Шаг 5. Вычисляем Qi (i = 1,2,.. ,,m):

a) если ds < 0, то

Предположим, что начальная базисная матрица - единичная, Е. К концу і-й итерации матрицу В1 можно представить в виде:

В 1 = нкнк_ 1 ...Н 2 НЕ, (5)

где Ні - матрица вида (4) (і = 1,2,.,к). Представление (5) известно как мультипликативная форма обратной матрицы [9]. Заметим, что для хранения элементов матрицы Ні достаточно п ячеек памяти.

Таким образом, часть симплексной таблицы, представляющую матрицу коэффициентов, на (к + 1)-м шаге можно представить как А(к+1) =(В^к+1)) -1А = НшА(к).

Можно заметить, что умножение матрицы А® на матрицу Нк+1 слева соответствует шагу исключения метода Жордана-Г аусса с ведущим элементом а{р1 = ар5, а коэффициенты вектора

П соответствуют коэффициентам, с которыми ведущая строка складывается с остальными строками матрицы А®.

4. Модифицированный симплекс-метод над кольцом целых чисел

Предположим, что все входные данные задачи (1)-(3) целочисленные. Тогда, вообще говоря, симплекс-метод требует арифметики рациональных чисел. Следуя [1, 8], покажем, как изменить алгоритм 1, чтобы все вычисления производились над кольцом целых чисел.

Вместо В1 основную роль теперь будет играть матрица В+ =[П^е^|[(В)|[1]]В-1, с точностью

*

до знака совпадающая с матрицей В , присоединенной к В. Обозначим А =[П^е^|[(В)|П].

Є, =

Д • /р - Д • хр

б) если ds > 0, то

0, =

+ да,

Д • xpi — Д • lpi

Д • хр — Д • мр

а. < 0;

Вычисляем 0min = min 0;, 0 = min{0min,u - 4}.

i

Остальные действия шага 5 такие же, как и в алгоритме 1.

Заметим, что при вычислении и сравнении 0 нет необходимости выполнять непосредственное деление на ~ .

Используя различные способы представления матрицы В+, можно получать различные вычислительные схемы симплекс-метода. В следующем разделе мы покажем, как адаптировать мультипликативную форму обратной матрицы к целочисленному случаю.

5. Мультипликативная форма присоединённой матрицы

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

Ч.Л. Доджсона [2] (где использовался для вычисления определителя). В работе [3] он приоб-

ретает современный вид. Метод несколько раз переоткрывался [4-7]. В одном из его вариантов (адаптация метода исключения Жордана-Гаусса) вычисления на к-й итерации над элементами матрицы выполняются по формуле (для простоты мы предполагаем, что ведущий элемент находится на диагонали)

,(k—1)„( k-1)

'a,

kk

- a,

(k-1) (k-1)

kj

a

ik

(k-1)

k —i.k —I

,(k—1)

i = k.

(6)

Пользуясь тождеством Сильвестра (см., например, [15]), можно показать, что результат деления в (6) - целое число.

В [1] введено понятие ^-матрицы. Опираясь на него, введем мультипликативную форму присоединенной матрицы. Пусть Ъ =

= (Р1,Р2,...,Р т) - некоторая база в А. Тогда Q-матрицей Q(A,Ъ) называется т*и матрица с элементами цгу = det Аъ-рг+у, где матрица Аъ-рг+у получена из Аъ = Б заменой рг-го столбца нау-й. Отметим следующие свойства матрицы Q(A,Ъ):

• Q(A,Ъ) = Б .А, где Б - матрица, присоединённая к Б;

• если А - целочисленная, то матрица Q(A,Ъ) тоже целочисленная;

• если Ъ' получается из Ъ = ^Р1,Р2,...,Р^ заменой Рр на 5, то переход от матрицы Q = = Q(A,Ъ) к Q' = Q(A,Ъ') осуществляется с помощью одной итерации метода Доджсона с ведущим элементом др5:

9ц =

4j4pS — 4pj4is

det B

qj,

i ^ P, i = P,

причем det Б' = др„ где Б' = Аъ'.

Таким образом, справедливо

Q(A,Ъ') = (Б)А = WkQ(A,Ъ) = WkБ*А, где Wк - матрица преобразования, проведённого по правилам метода (6):

5k 5k 5k

5 1’ ? 2 2 P—1

k—1 5k—1 5k—1

5k 5k

Ц s , ^ep+1,•••,

5k —1 5k —1

Vis =

(7)

q,s

1,

г = P,

5к - ведущий элемент на к-й итерации.

Пусть матрица А содержит единичную подматрицу, тогда её можно взять в качестве начальной базисной подматрицы. После проведения к исключений получим целочисленную матрицу Q(A, Ък) = Wk. W2W1A, где Wi - матрица вида (7). Итак, Q(A,Ъk) = Б*А, где

B= Wb..W2W1. (8)

Представление (8) назовем мультипликативной формой присоединенной матрицы.

Заметим, что если при выполнении преобразования Доджсона мы будем умножать ведущую строку на знак ведущего элемента, то в итоге получим матрицу B+ = |det(B)|B 2, которая с точностью до знака совпадает с присоединённой матрицей B .

6. Экспериментальные результаты

В таблице «Экспериментальные данные» приведены результаты решения задач целочисленного линейного программирования из библиотеки NETLIB [4] модифицированным симплекс-методом с мультипликативным представлением обратной матрицы, использующим рациональную арифметику, и предложенным здесь симплекс-методом с мультипликативной формой присоединенной матрицы, использующим целочисленную арифметику.

Экспериментальная программа написана на языке программирования C++. Все вычисления точные. Для представления целых чисел использовалась библиотека GMP [16]. Для представления рациональных чисел использовался шаблонный тип рациональных чисел из библиотеки Arageli [17] на основе целых чисел из библиотеки GMP. Создание и ведение проекта, редактирование и компиляция кода производились с помощью интегрированной среды разработки программного обеспечения Visual Studio 9.0 [18].

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

Тем не менее, как видно из таблицы 1, на некоторых задачах рациональный симплекс-метод оказался существенно быстрее целочисленного, что могло произойти в результате двух причин:

• Изменение исходных данных в результате преобразования задачи к целочисленному виду, привело к тому, что в процессе выполнения алгоритма выбирались другие промежуточные базы, что изменило количество итераций для решения задачи. Также выбор других промежуточных баз мог сказаться на величине промежуточных коэффициентов (задачи bore3d, Israel, ship08s).

• Хотя рост коэффициентов ограничен для целочисленного симплекс-метода, он может быть значительно больше, чем для рационального, так как отсутствуют операции деления на НОД (задачи Finnis, scsd6, ship04s, ship08s).

Таблица

Экспериментальные данные

Имя задачи Размерность задачи m x n Количество ненулевых элементов Количество итераций Время выполнения задачи, мс Максимальное количество бит на представление целого числа после чтения задачи Максимальное количество бит на представление целого числа при решении задачи

Q Z Q Z Q Z Q Z

Adlitle 57x97 465 138 106 672 125.8 17 27 386 785

Afiro 27x32 88 18 18 9.75 2.24 12 12 63 208

Agg 489x163 2451 153 163 3760.8 2397.2 23 32 458 2525

agg2 517x302 4515 162 159 1482.97 1289.3 27 36 339 4212

agg3 517x302 4531 166 160 1989.76 1428.6 27 36 514 4451

Beaconfd 174x262 3476 111 111 499.7 235.4 14 17 385 1747

boeingl 351x384 3865 1116 1114 119847 39590.1 31 32 655 3911

boeing2 167x143 1396 295 272 5421.9 1178.55 29 34 402 852

bore3d 234x315 1525 183 258 4273.7 4879.6 25 25 809 4060

Finnis 498x614 2714 1008 1125 150479 868037 31 50 859 14241

fitld 25x1026 14430 1269 1477 44690 4442.8 11 12 543 1385

Israel 174x142 2269 345 446 30072 36067 20 24 426 2261

Lofti 154x308 1086 328 376 3122.5 3078.8 8 34 192 1446

sc50a 51x48 131 49 49 61.72 11.87 9 9 238 739

sc50b 51x48 119 51 51 68.39 13.36 8 9 74 136

sc205 206x203 552 226 231 7417.3 961.6 9 9 358 321

scagr25 472x500 2029 917 915 132933 127294 18 18 564 1863

scfxml 331x457 2612 550 581 49364 33628.8 17 19 910 4803

Scorpion 389x358 1708 386 385 3946.3 3387.1 20 20 208 4212

scsd6 148x1350 5666 477 477 96315 120582 29 29 1564 11660

sctapl 301x480 2052 363 36 10380 5171.7 7 7 479 1394

sharelb 118x225 1182 379 391 34842 8699.6 25 25 1243 2906

ship04s 403x1458 5810 487 434 5500.2 11183 24 26 171 13416

ship08s 779x2387 9501 706 786 16035 116400 24 27 184 21594

Seba 516x1028 4874 611 598 22493 4500.2 17 19 48 132

Shell 537x1775 4900 757 757 10275.5 2356 19 19 188 1409

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

1. Edmonds J., Maurras J.F. Notes sur les Q-matrices d’Edmonds // Recherche operationelle. 1997. V. 31. № 2. P. 203-209.

2. Dodgson C.L. Condensation of determinants // Proceedings of Royal Society of London. 1856. V. 15. P. 150-155.

3. Waugh F.V., Dwyer P.S. Compact computation of the inverse of a matrix // Annals of Mathematical Statistics. 1945. V. 16. P. 259-271.

4. Edmonds J. Systems of distinct representatives and linear algebra // Journal of Research of the National Bureau of Standards. 1967. V. 71B. № 4. P. 241-245.

5. Bareiss E.H. Sylvester's identity and multistep integer-preserving Gaussian elimination // Mathematics of Computation. 1968. V. 22. № 103. P. 565-578.

6. Bareiss E.H. Computational solutions of matrix problems over an integral domain // J. Inst. Maths Ap-plics. 1972. V. 10. P. 68-104.

7. Малашонок Г.И. Решение системы линейных уравнений в целостном кольце // Журнал вычислительной математики и математической физики. 1983. Т. 23. № 6. С. 1497-1500.

8. Azulay D., Pique J.F. A revised simplex method with integer Q-matrices // ACM Transactions on Mathematical Software. 2001. V. 27. № 3. P. 350-360.

9. Dantzig G.B., Orchard-Hays W. Alternate algorithm for the revised simplex method using product form of the inverse. The RAND Corporation Research Memorandum RM-1268. 1953.

10. netlib. URL: http://netlib.org/

11. Канторович Л.В., Залгаллер В.А. Расчет рационального раскроя промышленных материалов. М.: Госстатиздат, 1951.

12. Dantzig G.B., Orden A., Wolfe P. The generalized simplex method. The RAND Corporation Research Memorandum RM-1264. 1954.

13. Dantzig G.B. Upper bounds, secondary constraints and block triangularity in linear programming // Econometrica. 1955. V. 23. № 2. P. 174-183.

14. Makhorin A. GNU Linear Programming Kit. Version 2.1. Implementation of the Revised Simplex Method. 2001. URL: http://www.ime.unicamp.br/~mo-retti/ms428/1sem2006/simplex.pdf

15. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1966.

16. GMP. URL: http://gmplib.org/

17. Arageli. URL: http://www.arageli.org/

18. Microsoft Visual Studio 9.0. URL: http://msdn. microsoft.com/en-us/vstudio/ff431703

ALL-INTEGER SPARSE SIMPLEX METHOD N. Yu.Zolotykh, V.K. Kubarev

A variant of the modified simplex method where all operations are performed over the ring of integers is proposed. The suggested method is intended to solve large sparse linear programming problems. The results of experiments show the viability of the suggested method.

Keywords: modified simplex method, sparse simplex method, integer simplex method, integer Q-matrix, multiplicative form of complementary matrix, Arageli.