CC BY
6
УДК 519.85 Вестник СПбГУ. Математика. Механика. Астрономия. 2022. Т. 9 (67). Вып. 1
MSC 90C20
Быстрый алгоритм решения квадратичной задачи о ранце
А. В. Плоткин
Санкт-Петербургский государственный университет,
Российская Федерация, 199034, Санкт-Петербург, Университетская наб., 7—9
Для цитирования: Плоткин А. В. Быстрый алгоритм решения квадратичной задачи о ранце // Вестник Санкт-Петербургского университета. Математика. Механика. Астрономия. 2022. Т. 9(67). Вып. 1. С. 76-84. https://doi.org/10.21638/spbu01.2022.108
В работе рассматривается задача квадратичного программирования со строго выпуклой сепарабельной целевой функцией, единственным линейным ограничением и двусторонними ограничениями на переменные. В англоязычной литературе эта задача называется Convex Knapsack Separable Quadratic Problem, или сокращенно — CKSQP. Нас интересует алгоритм решения задачи CKSQP с линейной сложностью. Посвященные этой теме работы содержат неточности в описании алгоритмов и неэффективные реализации. В данной работе удалось преодолеть имеющиеся трудности.
Ключевые слова: квадратичное программирование, задача о ранце, быстрые алгоритмы.
1. Постановка задачи. Рассмотрим следующую задачу квадратичного программирования :
7j{Hx,x) — (с,х) —> min,
n
aiXi = b, (1)
i=i
li ^ xi ^ ri, i G 1 : n.
Здесь b, cj, ai — вещественные параметры, li G R U {-то}, ri G R U {+то}, и H — диагональная матрица порядка n с положительными диагональными элементами hi,h,2,... ,hn. Целевая функция задачи (1) является строго выпуклой и сепарабельной. Через О обозначим множество планов, то есть векторов x = (xi, X2,..., xn), удовлетворяющих ограничениям задачи (1).
Данная задача рассматривалась, например, в диссертации [1], где имеется обширная библиография. В работе [2] приводится критический анализ существующих алгоритмов решения задачи (1).
В тексте статьи мы будем предполагать, что li, ri G R. Это позволит упростить изложение. При этом разработанное программное обеспечение будет учитывать случаи li = —то и ri = +то.
Во многих приложениях задача (1) встречается с единичной матрицей H, что соответствует проектированию точки c на множество О [3, 4]. Из других приложений можно выделить квадратичную задачу распределений ресурсов [5] и вычисление матрицы обновления в квазиньютоновском методе при наличии ограничений на якобиан [6].
© Санкт-Петербургский государственный университет, 2022
2. Упрощение постановки. Не умаляя общности, можно считать, что в задаче (1) выполнены условия li < ri и ai = 0 при всех i от 1 до п. Если ai = 0, то оптимальное значение xi находится путем решения тривиальной задачи:
^hixf — CiXi —> min, Ii < xi < ri.
Воспользовавшись предположением о том, что коэффициенты ai в определении множества Q отличны от нуля, выполним замену переменных:
Xi = Xi/ai + Ci/hi, i G 1: n. (2)
Получим задачу следующего вида:
Tj(Hx, х) —> min,
n
i=l
где
Xi = b, (3)
1i < Xi < ri, i G 1 : n,
2 и l„2-
H = diag(hi/a\,..., hn/an)
f , v"^ aici hi
i=1
li =
ai(li — ci/hi), если ai > 0, ai(ri — c-i/hi), если ai < 0,
{ai(ri — ci/hi), если ai > 0, ai(li — ci/hi), если ai < 0.
Задача (3) является частным случаем задачи (1). Преобразование (2) имеет линейную сложность и значительно упрощает дальнейшие вычисления. В связи с этим далее будем считать, что решаемая задача имеет следующий вид:
Q(x) := ^(Нх, х) —> min,
n
= b (4)
i= 1
li ^ Xi ^ ri, i G 1 : n.
3. Условия Куна — Таккера. Целевая функция задачи (4) непрерывна и строго выпукла, а значит критерием существования решения является непустота множества планов Q. В данном случае условие непустоты множества планов записывается следующим образом:
li < b п. (5)
li < b < у r
i= i i=i
Будем считать, что условие (5) выполнено, а значит множество планов О непусто.
Заметим, что множество П является выпуклым. Вместе со строгой выпуклостью целевой функции ф(ж) это гарантирует единственность решения задачи (4). Для его нахождения воспользуемся условиями Куна — Таккера:
= Л + Щ - V»,
«¿(ж» - /¿)=0, - Ж»)=0 ,
> г £ 1 : п, (6)
Щ > 0, V» > 0,
П
= 6. (7)
¿=1
4. Характеристика оптимального плана. Разберемся с условиями (6). Для этого нам потребуется вспомогательная лемма.
Лемма. Пусть выполнены соотношения
Лт = Л + и — V, (8)
и(т - /) = 0, V(г - т)=0, (9)
и > 0, V > 0, (10)
/ ^ т ^ г,
где Л, /, г — вещественные константы, причем с > 0 и / < г. Тогда необходимо, чтобы
{/, если Л/Л ^ /,
Л/Л, если / < Л/Л < г, (11)
г, если Л/Л ^ г.
Доказательство. Проведем анализ полного набора альтернатив для переменной т.
(а1) т = /. Согласно (9), V = 0. В силу (8) и (10) имеем Лт = Л + и ^ Л, так что
Л < Л/. (12)
(а2) / < т < г. Согласно (9), и = V = 0, так что в силу (8) имеем Л = Лт. В частности,
Л £ (Л/, Лг). (13)
(а3) т = г. Согласно (9), и = 0. В силу (8) и (10) имеем Лт = Л - V ^ Л, так что
Л > Лг. (14)
Отметим, что условия (12)-(14) представляют собой полный набор альтернатив для переменной Л £ К.
Переходим к доказательству формулы (11).
Если Л/Л ^ /, то Л £ (Л/, Лг) и Л ^ Л/ < Лг, поэтому альтернативы (а2) и (а3) исключаются. Остается одна возможность — т = /.
Если / < Л/Л < г, то Л > Л/ и Л < Лг, поэтому альтернативы (а1) и (а3) исключаются. Остается одна возможность — / < т < г. Согласно (а2), т = Л/Л.
Наконец, если А/к ^ г, то А 0 (к1,кг) и А ^ кг > к1, поэтому альтернативы (а1) и (а2) исключаются. Остается одна возможность — и = г.
Лемма доказана. □
Формула (11) определяет переменную и как функцию от А: и = и (А). На рис. 1 изображен график функции и (А).
П)' Ы
0 кг X
Рис. 1. График функции (Л). Функция и (А) допускает аналитическое представление (см., например, [7])
т(Х)=1 + ^Х-Ы)+-^Х-кг)+, (15)
где (щ)+ = тах{0,щ}.
5. Сведение к скалярному уравнению. Вернемся к условиям Куна — Таккера (6), (7).
Воспользуемся доказанной леммой при каждом 2 от 1 до п. Положив и = хг, к = кг, I = ¡г, г = гг, и = щ, V = V, получим представление хг как функций от А:
ж* (Л) = к + £7(А - Мг)+ - - г € 1 : п.
Каждая из функций хг (А) является непрерывной, неубывающей, имеет минимальное ¡г и максимальное гг значения.
Таким образом, условию (6) удовлетворяет множество
X = {(хх(А),...,хп(А)) | А е К} .
Добавив необходимость выполнения условия (7) на этом множестве, получим следующее уравнение:
(А) = Ь, А е К.
(16)
Свойства функций хг (А) вместе с условием (5) гарантируют существование решения данного уравнения. Подведем итоги.
Теорема. Пусть А* — решение уравнения (16), тогда х* = (х*, ...,хП), где х* = хх (А*),..., хП = хп (А*), является единственным решением задачи (4).
Далее нас будет интересовать быстрый алгоритм решения уравнения (16). Вестник СПбГУ. Математика. Механика. Астрономия. 2022. Т. 9(67). Вып. 1 79
6. Подготовка к построению алгоритма. Введем следующие обозначения:
= Ь
ai+n = ~TTi' Ai Zi,
Ai+n hiri?
> i £ 1 : n,
^i(A) = ai(A - Ai)+, i £ 1 : 2n,
2n
F (A) = £ &(A),
i=i
C = b li.
i=i
Тогда уравнение (16) можно переписать в следующем виде:
F (A) = C. (17)
При этом функция F(A) является неубывающей ломаной с узлами Ai,..., A2n.
Во избежание двусмысленности дадим определение медиане, которого будем придерживаться в дальнейшем.
Определение. Пусть задано конечное множество вещественных чисел A = {ai}"=1. Обозначим через b1, b2,..., bm последовательность, полученную в результате упорядочивания элементов множества A в порядке неубывания. Тогда элемент bj-m+i-j называется медианой множества А.
Напомним, что медиану множества A можно найти за линейное время [8, с. 250253].
Перейдем к построению алгоритма.
Условие существования решения гарантирует нам выполнение неравенства F( min Ai) < C < F( max Ai).
iei:2n iei:2n
Наш алгоритм будет вычислять значения функции F исключительно в узлах Ai, а в качестве результата будет возвращать узел A = max{Ai | F(Ai) ^ C}. Имея такой узел, несложно найти ближайший к нему узел справа (если такой имеется), вычислить значение функции в обоих узлах и выполнить обратную интерполяцию для нахождения A*. Такие дополнительные действия требуют линейного времени работы и не окажут влияния на итоговую вычислительную сложность алгоритма.
Для упрощения понимания сначала мы представим описание алгоритма, работающего за время O(n log n), а затем покажем, как добиться линейного времени работы.
7. Упрощенная версия — алгоритм 1. Опишем упрощенную версию алгоритма.
Начальный шаг. Положим /0 = {1, 2,..., 2n}.
Инвариант. Искомый узел Л содержится во множестве Лк, где Лк = {Ai}ie1fc.
Шаг алгоритма. Пусть имеется множество /k. Если множество /k содержит единственный элемент, то этот элемент является индексом искомого узла Л. Иначе, найдем медиану множества Ak и обозначим этот узел через Amed, где med — индекс в /к. Вычислим значение функции: Fmed = F(Amed). Если Fmed ^ C, то значение Л точно не меньше Amed. Если же Fmed > C, то Л точно меньше Amed.
Таким образом, положим
i{med} U {i £ 4 | Ai > Amed}, если Fmed < C, [{i £ /к | Ai < Amed}, если Fmed > C.
Полученное множество /k+i является подмножеством /k и имеет размер, не
- Г|iklо
превышающий Ч^1 .
Оценка сложности. Алгоритм делает O(log n) итераций. На каждой итерации алгоритма выполняется вычисление функции F за O(n) операций. Затраты на вычисление медианы суммарно составляют 0(2n) + 0(п) + О(^) + O(j) + ... + 0( 1) = O(n). Таким образом, итоговая сложность алгоритма равна O(n log n).
8. Быстрый алгоритм — алгоритм 2. Для достижения линейной сложности необходимо научиться вычислять значение целевой функции как можно быстрее. Добъемся того, чтобы вычисление функции F на k-й итерации требовало O(|/k|) операций.
Сделаем следующее наблюдение. Возьмем произвольный индекс io. Пусть в какой-то момент этот индекс вышел из последовательности множеств /: io £ /к, io £ /k+i. Тогда Ai0 ^ minAk+i или Ai0 ^ тахЛк+i по построению. Так как все дальнейшие вычисления функции F будут происходить только в узлах /к+1, то уже на этапе исключения индекса io можно однозначно описать вклад функции ^i0 (A) в F(A). Если Ai0 ^ minA,^, то ^i0 (A) = ai0 A — ai0 Ai0 для всех A £ /k+i. Если Ai0 > maxAk+i, то фг0 (A) = 0 для всех A £ /k+i.
Введем и будем поддерживать линейную функцию AkA + Bk, в которой будут «аккумулироваться» линейные функции узлов, вышедших из рассмотрения.
Начальный шаг. Положим /0 = {1, 2,..., 2n}, A0 =0, B0 = 0.
Инвариант 1. Искомый узел A содержится во множестве Ak, где Ak = {Ai}ie1fc.
Инвариант 2. Для всех узлов A из /k значение функции F может быть вычислено по формуле F(A) = У]^^ Фг(Л) + Ak A + Bk.
Шаг алгоритма. Пусть имеется множество /k. Если множество /k содержит единственный элемент, то этот элемент является индексом искомого узла Л. Иначе, найдем медиану множества Ak и обозначим этот узел через Amed, где med — индекс в /k. Вычислим значение функции: Fmed = ]Cie/fc ^i(Amed) + Ak Amed + Bk. Если Fmed ^ C, то значение A не меньше Amed. Если же Fmed > C, то Л точно меньше
Amed.
В случае Fmed ^ C положим
/k+i = {med} U {i £ /k | Ai > Amed},
Ak+i = Ak + ^ ai, ieifc\ifc+i
Bk+i = Bk — ^^ «iAi. ieifc\ifc+i
В случае Fmed > C положим
Ik+1 = {i G Ik | Aj < Amed},
Ak+1 = Ak, Bk+1 = Bk.
Оценка сложности. Теперь затраты на вычисление функций, пересчет коэффициентов и работу со множествами составляют 0(2n) + O(n) + ) + О(^) + • • • + O(1) = O(n). Таким образом, итоговая сложность алгоритма равна O(n).
9. Численные эксперименты. Проведем серию численных экспериментов.
Будем решать задачи вида (4) для разных n: 104,105,106,107. Коэффициенты h будем выбирать случайно из отрезка [0.1,1], коэффициенты I и r — из отрезка [-1,1]. Параметр b выберем случайно из отрезка [^П=12Г=1 r®] .
В нашем алгоритме для поиска медианы мы будем использовать алгоритм Qick-select [8, с. 245-249], осуществляющий поиск за O(n) в среднем, но показывающий хорошие результаты на практике.
Для каждого из рассматриваемых n были сгенерированы 100 тестовых задач. Вычисления производились на машине с процессором AMD Ryzen 5 3500U и 16 Гб ОЗУ. В табл. 1 приводится среднее время работы обеих версий алгоритма.
Таблица 1. Среднее время работы алгоритмов (мс) для решения задачи вида (4)
га Алгоритм 1: O(ralogra) Алгоритм 2: О(га)
104 2 1
10Б 16 11
106 105 71
107 1175 720
Далее было решено провести сравнение скорости работы полученной реализации со специализированной библиотекой алгоритмов CQKnPClass Project [9]. Один из представленных в библиотеке алгоритмов под названием DualCQKnP ориентирован на решение задач вида (1).
Сравнительный анализ производился с помощью тестовой процедуры MainRnd, использовавшейся авторами статьи [9].
Для n, равного 104, 105, 106, 107, были созданы 100 тестовых задач. Некоторые из них могли иметь пустое множество планов. В табл. 2 приводится среднее время работы сравниваемых алгоритмов.
Таблица 2. Среднее время работы алгоритмов (мс) для решения задачи вида (1)
га DualCQKnP Алгоритм 1: О (га 1с g га) Алгоритм 2: О (га)
104 1 1 1
10Б 20 11 4
106 303 144 43
107 4552 1588 486
Стоит отметить, что алгоритм DualCQKnP имеет вычислительную сложность O(n log n), и авторы ставили под сомнение возможность эффективной реализации алгоритма с линейным временем работы.
Реализация полученных алгоритмов и тестовые процедуры доступны в репози-тории [10].
Литература
1. Jeong J. Indefinite Knapsack Separable Quadratic Programming: Methods and Applications. Dr. Sci. thesis (2014).
2. Kiwiel K. Breakpoint searching algorithms for the continuous quadratic knapsack problem. Mathematical Programming 112, 473-491 (2008). https://doi.org/10.1007/s10107-006-0050-z
3. Dai Y. H., Fletcher R. New algorithms for singly linearly constrained quadratic programs subject to lower and upper bounds. Mathematical Programming 106, 403-421 (2006). https://doi.org/10.1007/s10107-005-0595-2
4. Nielsen S., Zenios S. Massively parallel algorithms for singly-constrained convex programs. ORSA J. Comput. 4, 166-181 (1992). https://doi.org/10.1287/ijoc.4.2.166
5. Patriksson M. A survey on the continuous nonlinear resource allocation problem. European J. Operational Research 185, 1-46 (2008). https://doi.org/10.1016/j.ejor.2006.12.006
6. Calamai P., More J. Quasi-Newton updates with bounds. SIAM J. Numer. Anal. 24, 1434-1441 (1987).
7. Малозёмов В. Н., Тамасян Г. Ш. Представления непрерывных кусочно-аффинных функций. В: Семинар «CNSA & NDO». Избранные доклады. 10 сентября 2019 г. (2019).
8. Кормен Т., Лейзерсон Ч., Ривест Р., Штайн К. Алгоритмы. Построение и анализ, пер. с англ. 3-е изд. Вильямс (2013).
9. Frangioni A., Gorgone E. A library for continuous convex separable quadratic knapsack problems. European J. Operational Research 229, 37-40 (2013). https://doi.org/10.1016/j.ejor.2013.02.038
10. GitHub. Доступно на: https://github.com/WhatElseIsThere/CKSQP (дата обращения: 30.11.2021).
Статья поступила в редакцию 10 июня 2021 г.;
доработана 16 июня 2021 г.; рекомендована к печати 2 сентября 2021 г.
Контактная информация:
Плоткин Артем Владимирович — аспирант; avplotkin@gmail.com
Fast algorithm for quadratic knapsack problem
A. V. Plotkin
St Petersburg State University, 7—9, Universitetskaya nab., St Petersburg, 199034, Russian Federation
For citation: Plotkin A. V. Fast algorithm for quadratic knapsack problem. Vestnik of Saint Petersburg University. Mathematics. Mechanics. Astronomy, 2022, vol. 9 (67), issue 1, pp. 76-84. https://doi.org/10.21638/spbu01.2022.108 (In Russian)
The paper considers a quadratic programming problem with a strictly convex separable objective function, a single linear constraint, and two-sided constraints on variables. This problem is commonly called the Convex Knapsack Separable Quadratic Problem, or CKSQP for short. We are interested in an algorithm for solving CKSQP with a linear time complexity. The papers devoted to this topic contain inaccuracies in the description of algorithms and ineffective implementations. In this work, the existing difficulties were overcome. Keywords: quadratic programming, knapsack problem, fast algorithms.
References
1. Jeong J. Indefinite Knapsack Separable Quadratic Programming: Methods and Applications. Dr. Sci. thesis (2014).
2. Kiwiel K. Breakpoint searching algorithms for the continuous quadratic knapsack problem. Mathematical Programming 112, 473-491 (2008). https://doi.org/10.1007/s10107-006-0050-z
3. Dai Y. H., Fletcher R. New algorithms for singly linearly constrained quadratic programs subject to lower and upper bounds. Mathematical Programming 106, 403-421 (2006). https://doi.org/10.1007/s10107-005-0595-2
4. Nielsen S., Zenios S. Massively parallel algorithms for singly-constrained convex programs. ORSA J. Comput. 4, 166-181 (1992). https://doi.org/10.1287/ijoc.4.2.166
5. Patriksson M. A survey on the continuous nonlinear resource allocation problem. European J. Operational Research 185, 1-46 (2008). https://doi.org/10.1016/j.ejor.2006.12.006
6. Calamai P., More J. Quasi-Newton updates with bounds. SIAM J. Numer. Anal. 24, 1434-1441 (1987).
7. Malozemov V. N., Tamasyan G.S. Representations of continuous piecewise affine functions. In: Seminar "CNSA & NDO". Selected reports. Sept. 10, 2019 (2019). (In Russian)
8. Cormen T. H., Leiserson C.E., Rivest R. L., Stein C. Introduction to Algorithms. 3rd ed. The MIT Press (2009). [Rus. ed.: Cormen T.H., Leiserson C.E., Rivest R. L., Stein C. Algoritmy. Postroenie i analiz. Williams Publ. (2013)].
9. Frangioni A., Gorgone E. A library for continuous convex separable quadratic knapsack problems. European J. Operational Research 229, 37-40 (2013). https://doi.org/10.1016/j.ejor.2013.02.038
10. GitHub. Available at: https://github.com/WhatElseIsThere/CKSQP (accessed: November 30, 2021).
Received: June 10, 2021 Revised: June 16, 2021 Accepted: September 2, 2021
Author's information:
Artyom V. Plotkin — avplotkin@gmail.com
