Симплекс-метод

Аккорд можно записать как 88 чисел — по одному на клавишу рояля: единица, если клавиша нажата, ноль, если нет. Пусть прибор не видит клавиш, а умеет только складывать их громкости с какими-то знаками: «первая плюс вторая минус третья…» — и выдаёт 24 таких суммы. Уравнений 24, неизвестных 88: решений бесконечно много. Но если среди всех неотрицательных решений взять то, у которого сумма чисел наименьшая, на всех 38 различных аккордах прелюдии Шопена до минор получается ровно тот аккорд, что звучал. Эта задача — линейная, а решает её симплекс-метод.

Кружки на клавишах — настоящий аккорд, оранжевые столбики — что нашёл симплекс-метод, серые — самое короткое решение, которое дают наименьшие квадраты. Серые размазаны по шести десяткам клавиш и ни разу не совпадают с аккордом. Сдвиньте число сумм: при 16 точно находится меньше половины аккордов, при 24 — все. Кривая под клавиатурой — доля аккордов выбранной пьесы, найденных точно, в зависимости от числа сумм: у разреженной прелюдии Баха, где звучат две-три ноты, все аккорды находятся уже при 18 суммах, у плотных аккордов Шопена — при 22.

Идея

Линейная задача в стандартной форме: найти неотрицательные x₁, …, x_n с A·x = b — это m равенств, — при которых c·x наименьшее. Любые неравенства сводятся к этому виду: «≤» становится равенством с лишней неотрицательной переменной. Допустимые x образуют многогранник, а линейная функция на многограннике, если её минимум существует, достигает его в вершине.

Вершина — это решение, у которого ненулевых переменных не больше m, а их столбцы в A линейно независимы; такой набор столбцов называется базисом. Симплекс-метод стоит в вершине и спрашивает про каждый столбец вне базиса: если начать увеличивать эту переменную, оставаясь на многограннике, функция будет расти или падать? Если какая-то переменная обещает уменьшение, она входит в базис, а выходит та, что первой упрётся в ноль, — это шаг по ребру в соседнюю вершину. Если никакая не обещает, текущая вершина оптимальна.

Отсюда и аккорды. Вершина многогранника A·x = b, x ≥ 0 имеет не больше 24 ненулевых клавиш, а сумма x наименьшая там, где «громкость» не размазана, а собрана в немногие клавиши. Если аккорд достаточно разрежен, а знаки в суммах случайны, такая вершина единственна, и это сам аккорд. Это доказали Эммануэль Кандес, Джастин Ромберг и Теренс Тао и, независимо, Дэвид Донохо в 2006 году; здесь мы этого не доказываем, а только проверяем на нотах. Наименьшие квадраты ищут самое короткое x, а у самого короткого вектора энергия распределена по всем клавишам понемногу.

Шаг метода

Пусть в базисе столбцы B. Переменные вне базиса равны нулю, а базисные из равенств: x_B = B⁻¹b. Если увеличить внебазисную переменную x_j на t, базисные изменятся на −t·B⁻¹A_j, а функция — на t·d_j, где

d_j = c_j − c_B·B⁻¹A_j

— приведённая стоимость. Отрицательная — переменная выгодна. Тогда t растёт, пока какая-то базисная не упрётся в ноль (наименьшее отношение), и они меняются местами. В таблице симплекс-метода всё это — одна операция метода Гаусса: строка с ведущим элементом делится на него, из остальных строк он вычитается.

Остаётся найти первую вершину. Для этого к каждому равенству добавляют свою «искусственную» неотрицательную переменную: при x = 0 искусственные равны правым частям (строки с отрицательной правой частью сначала умножают на −1), и это вершина. Фаза 1 тем же методом уменьшает сумму искусственных; если она дошла до нуля, найдена вершина исходной задачи. Фаза 2 уменьшает настоящую функцию.

Это первый аккорд прелюдии, до минор, — шесть клавиш по 24 суммам. Сначала в базисе 24 искусственные переменные, и за 38 шагов фаза 1 вытесняет их клавишами; сумма искусственных падает до нуля, а сумма x в найденной вершине — 10,24. Фаза 2 за четыре шага опускает её до шести — ровно столько нот в аккорде, и это уже ответ. Но метод об этом не знает: ещё 174 шага он меняет базис, не сдвигая точку, пока все приведённые стоимости не станут неотрицательными. Оптимальная вершина вырожденная — из 24 базисных переменных ненулевых только шесть, — и базисов у неё много.

Вырожденность лечат «шевелением»: правые части b сдвигают на случайные числа порядка 10⁻⁷. У сдвинутой задачи вершины невырожденные, и каждый шаг действительно сдвигает точку. Найденный базис потом применяют к настоящему b — обратная матрица базиса хранится в столбцах искусственных переменных. Переключите сцену: с шевелением тот же аккорд находится за 71 шаг вместо 216. Так поступают и настоящие решатели.

Код

Листинг решает каждую задачу дважды — как есть и с шевелением. Входящий столбец — с самой отрицательной приведённой стоимостью, правило Данцига. После шага, который не сдвинул точку, — первый подходящий по номеру, а при равных отношениях выходит переменная с меньшим номером: это правило Роберта Бленда (1977), и с ним метод не зацикливается на вырожденных вершинах. Доказательство этого длинное, и мы его не приводим. Суммы считаются обычным сложением по порядку, а не встроенной sum: начиная с Python 3.12 она складывает числа с плавающей точкой с поправкой на ошибку округления, как суммирование Кэхэна, и вырожденные шаги — сравнения с допуском — в Python и Go пошли бы по-разному.

Сколько это стоит

Шаг — один пересчёт таблицы, O(m·n). Сколько шагов, заранее не сказать. В 1972 году Виктор Кли и Джордж Минти построили слегка сплющенный n-мерный куб, на котором правило Данцига обходит все 2ⁿ вершин. Но на практике шагов обычно порядка числа строк, и объяснение этому дали Дэниел Спилман и Шан-Хуа Тэн в 2001 году: если чуть-чуть случайно пошевелить данные, среднее число шагов становится многочленным. На аккордах шагов в несколько раз больше строк — из-за вырожденности: при 24 суммах в среднем 139 шагов на аккорд, с шевелением — 62.

Перебор всех вершин — всех наборов из m столбцов — растёт как число сочетаний: при двадцати неизвестных в Python это больше полсекунды, и каждые четыре неизвестных добавляют множитель больше десяти. Симплекс-метод с шевелением решает задачу с 256 неизвестными за две с половиной секунды в Python и за сотые доли секунды в Go. Дальше учебная таблица начинает буксовать: при 384 неизвестных та же программа на Go сделала 57 674 шага — ошибки округления копятся в таблице от шага к шагу, и сдвиги порядка 10⁻⁷ перестают отличать вершины друг от друга. Настоящие решатели не хранят таблицу целиком, а держат разложение базисной матрицы и время от времени пересчитывают его заново. Многочленные по времени алгоритмы для линейных задач тоже есть: метод эллипсоидов Леонида Хачияна (1979) и методы внутренней точки Нарендры Кармаркара (1984). Современные решатели держат и симплекс-метод, и методы внутренней точки.

Где легко ошибиться

• Зацикливание. В вырожденной вершине шаг может не сдвинуть точку, и с правилом Данцига метод способен ходить по кругу между базисами одной вершины. Правило Бленда или случайный выбор среди подходящих столбцов это исключают.
• Допуски. Сравнения «приведённая стоимость меньше нуля» и «отношение меньше» в числах с плавающей точкой делаются с допуском; без него шум округления выбирает столбцы, которые на самом деле не улучшают.
• Фаза 1 после нуля. Когда сумма искусственных дошла до нуля, фазу 1 надо останавливать сразу: дальше её приведённые стоимости — шум округления, и метод может тысячи шагов переставлять базис ради ничего.
• Искусственные в базисе после фазы 1. Если сумма искусственных дошла до нуля, но какая-то осталась в базисе на нулевом уровне, перед фазой 2 её нужно заменить настоящим столбцом с ненулевым элементом в той же строке — иначе фаза 2 может поднять её обратно.
• Неотрицательность. Симплекс-метод работает с x ≥ 0. Свободную переменную записывают разностью двух неотрицательных. Для аккордов неотрицательность — подарок: громкость клавиши не бывает отрицательной, и это помогает восстановлению.
• Наименьшие квадраты вместо наименьшей суммы. Самый короткий вектор почти никогда не бывает разреженным; восстановить разреженный сигнал помогает именно сумма модулей, а не сумма квадратов.

Где работает

Решатели. HiGHS (его вызывает linprog в SciPy), GLPK, CPLEX и Gurobi решают задачи с миллионами переменных; внутри — двойственный симплекс-метод, который хорошо продолжает работу после добавления ограничения, и методы внутренней точки.

Целые числа. В методе ветвей и границ для задач с целыми переменными оценка ветви — решение той же задачи без требования целых чисел, и считают его симплекс-методом, сотни тысяч раз за одно решение.

Планирование. Перевозки, смены экипажей и расписания, смеси сырья на заводах — задачи, ради которых линейное программирование и придумали. Поток наименьшей стоимости — тоже линейная задача, и у неё есть свой, сетевой симплекс-метод.

Сжатые измерения. Восстановление разреженного сигнала по малому числу измерений минимизацией суммы модулей — basis pursuit (Скотт Чен, Дэвид Донохо, Майкл Сондерс, 1998) — стало основой ускоренной съёмки в магнитно-резонансной томографии: снимков меньше, картинка та же.

История

Леонид Канторович описал задачи линейного программирования и способ их решать в 1939 году; Джордж Данциг в 1947 году предложил симплекс-метод — этой истории посвящена глава о системах в курсе математики. Виктор Кли и Джордж Минти в 1972 году показали, что метод может быть экспоненциально медленным, Леонид Хачиян в 1979 году — что линейную задачу можно решить за многочленное время, Нарендра Кармаркар в 1984 году — что многочленный метод бывает и быстрым на практике. Роберт Бленд в 1977 году нашёл простое правило без зацикливания, а Дэниел Спилман и Шан-Хуа Тэн в 2001 году объяснили, почему на практике шагов мало.

Задачи

Задача о смеси, ради которой линейное программирование и появилось: самый дешёвый сплав из имеющегося лома.