Методы построения минимального эллипсоида

Кольцов Максим Алексеевич

Методы построения минимального эллипсоида

Кольцов Максим Алексеевич. Методы построения минимального эллипсоида. Научный руководитель профессор В. Н. Малозёмов. Направление математика, механика, кафедра исследования операций. Использовано 11 источников. В работе описано два известных алгоритма решения задачи о минимальном эллипсоиде. Использована теория двойственности по Лагранжу. На основе этих алгоритмов предложен и реализован комбинированный метод, достигающий хороших результатов на тестовых примерах. Также рассмотрен частный случай этой задачи, задача Сильвестра. Для неё получено решения средствами квадратичного программирования. Все алгоритмы реализованы в MATLAB и протестированы на модельных данных. Кольцов, М. А. Методы построения минимального эллипсоида: выпускная квалификационная работа бакалавра. / Кольцов Максим Алексеевич. — СПБ., 2016. — 33 с. — Библиогр.: с. 33.

Математика

Дипломы

Вуз: Санкт-Петербургский государственный университет (СПбГУ)

ID: 587d36575f1be77c40d58d28

UUID: cd0cc89f-4819-41a3-8c40-8c1761494ab0

Язык: Русский

Опубликовано: почти 8 лет назад

Просмотры: 27

Кольцов Максим Алексеевич

Источник: Санкт-Петербургский государственный университет

Комментировать 0

Рецензировать 0

Скачать - 529771 bytes

Поделиться работой

Санкт-Петербургский государственный университет Прикладная математика и информатика Исследование операций и принятие решений в задачах оптимизации, управления и экономики Кольцов Максим Алексеевич Методы построения минимального эллипсоида Бакалаврская работа Научный руководитель: д. ф.-м. н., профессор В. Н. Малозёмов Рецензент: к. ф.-м. н., научный сотрудник М. В. Долгополик Санкт-Петербург 2016

Saint Petersburg State University Applied Mathematics and Computer Science Operation Research and Decision Making in Optimisation, Control and Economics Problems Koltsov Maxim Alekseevich Methods of building minimal ellipsoid Bachelor’s Thesis Scientific Supervisor: Professor V. N. Malozemov Reviewer: Research fellow M. V. Dolgopolik Saint Petersburg 2016

3 Оглавление Введение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Глава 1. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Глава 2. Задача Сильвестра . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1. Постановка задачи . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.2. Общий вид задачи квадратичного программирования и двойственной к ней . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3. Двойственная задача для задачи Сильвестра и её решение . . . . . . . . 8 2.4. Решение задачи в среде MATLAB . . . . . . . . . . . . . . . . . . . . . . 9 Глава 3. Алгоритм Шора . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1. Оператор сжатия пространства . . . . . . . . . . . . . . . . . . . . . . . . 12 3.2. Итеративный алгоритм построения минимального эллипсоида . . . . . . 14 3.3. Практическое испытание алгоритма . . . . . . . . . . . . . . . . . . . . . 18 Глава 4. Алгоритм Хачияна . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.1. Вспомогательные сведения. . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2. Вывод двойственной задачи для задачи о минимальном эллипсоиде. . . 21 4.3. Алгоритм Хачияна. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 Глава 5. Комбинированный алгоритм . . . . . . . . . . . . . . . . . . . . . . 27 Глава 6. Сводное сравнение методов . . . . . . . . . . . . . . . . . . . . . . 30 Заключение . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Список литературы . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 Введение В данной работе рассматривается задача построения эллипсоида минимального объёма, содержащего заданное множество точек в R𝑛 . Как её частный случай также рассмотрена задача Сильвестра — задача построения минимального шара. Решение задачи Сильвестра получается с помощью двойственной задачи квадра тичного программирования. Оно реализовано в MATLAB. Продемонстрированы резуль таты работы программы на тестовых данных. Для задачи о минимальном эллипсоиде описано три алгоритма, дающих прибли жённое решение. Все три алгоритма используют вспомогательную задачу в простран стве R𝑛+1 , которая позволяет искать эллипсоид с фиксированным центром. Представлен вывод формул, по которым восстанавливается решение исходной задачи. Первый алгоритм, алгоритм Шора, использует оператор сжатия пространства. Оператор сжатия рассмотрен подробно и его работа проиллюстрирована на несколь ких примерах с шарами и эллипсами на плоскости. Вторым был рассмотрен алгоритм Хачияна. Этот алгоритм решает двойственную по Лагранжу задачу, по решению которой однозначно восстанавливается решение ис ходной. Подробно представлен вывод двойственной задачи и самого алгоритма. На основе алгоритма Хачияна разработан и описан комбинированный алгоритм, с помощью которого решается двойственная задача. Используется градиентный метод без точного линейного поиска и идея активного множества. Все три метода реализованы в MATLAB и протестированы на модельных задачах. Приведены сводные таблицы результатов, в которых представлены данные о точности вычислений и затраченном времени. Продемонстрировано превосходство комбинирован ного метода во всех тестовых случаях.

5 Глава 1 Постановка задачи Основной темой работы является задача о построении минимального эллипсоида: дано множество точек 𝑎𝑖 ∈ R𝑛 , 𝑖 ∈ 1 : 𝑚, требуется найти эллипсоид минимального объёма, который содержит все точки. Для существования решения будем предполагать, что афинная оболочка множества {𝑎𝑖 } совпадает с R𝑛 . Эллипсоид будем задавать таким образом: ℰ = {𝑥 ∈ R𝑛 | ⟨𝑀 (𝑥 − 𝑐), 𝑥 − 𝑐⟩ 6 1} , где 𝑀 — симметричная положительно-определённая матрица, 𝑐 — центр эллипсоида. Объём вычисляется по формуле vol(ℰ) = 𝜔𝑛 (det 𝑀 )−1/2 , 𝜋 𝑛/2 — объём единичного шара в 𝑛-мерном пространстве. Γ(𝑛/2 + 1) Пренебрегая множителем 𝜔𝑛 и переходя к логарифму целевой функции, задачу о где 𝜔𝑛 = минимальном эллипсоиде можно формально поставить как 𝑓 (𝑀, 𝑐) := − ln det 𝑀 → min ⟨𝑀 (𝑎𝑖 − 𝑐), 𝑎𝑖 − 𝑐⟩ 6 1, (1.1) 𝑖∈1:𝑚 𝑀 ≻ 0, где запись 𝑀 ≻ 0 означает положительную определённость матрицы 𝑀 . Основным пространством задачи является пространство симметричных положи тельно определённых матриц, на котором задана векторная норма — корень из суммы квадратов элементов матрицы. Таким образом, в задаче (1.1) имеется 𝑛(𝑛+1) 2 + 𝑛 неза висимых переменных. Отдельную сложность представляет то, что в задаче требуется найти и матрицу эллипсоида, и его центр. В работе будет представлен способ, позво ляющий искать только матрицу при фиксированном центре, решая задачу с (𝑛+1)(𝑛+2) 2 независимыми переменными. Также будет представлен способ перехода к задаче, в ко торой содержится 𝑚 независимых переменных. В работе будет рассмотрено несколько алгоритмов, решающих задачу (1.1) или эквивалентные ей.

6 Глава 2 Задача Сильвестра 2.1. Постановка задачи Рассмотрим предварительно задачу Сильвестра — по заданному набору точек 𝑎𝑖 ∈ R𝑛 , 𝑖 ∈ 𝑀 = 1 : 𝑚, найти наименьший по объёму шар, их содержащий. Чтобы формально поставить задачу, введем переменную 𝑥 ∈ R𝑛 — центр искомого шара. Оче видно, что для того, чтобы шар 𝐵R𝑛 (𝑥, 𝑟) содержал все заданные точки, необходимо, чтобы его радиус 𝑟 удовлетворял следующему условию: ∀𝑖 ∈ 𝑀 ‖𝑎𝑖 − 𝑥‖ 6 𝑟 То есть, в задаче Сильвестра требуется найти точку минимума функции max𝑖∈𝑀 {‖𝑎𝑖 − 𝑥‖} на R𝑛 , при этом радиус искомого минимального шара будет равен значению функ ции в этой точке. Заметим, что выражение под знаком максимума можно заменить на 1/2‖𝑎𝑖 −𝑥‖2 , что не меняет точку минимума. Теперь можно сказать так: решение задачи Сильвестра — это точка 𝑥* , являющаяся решением экстремальной задачи 𝜙(𝑥) := max{1/2‖𝑎𝑖 − 𝑥‖2 } → min𝑛 𝑖∈𝑀 𝑥∈R а радиус соответствующего минимального шара равен , (2.1) √︀ 2𝜙(𝑥* ). Доказательство суще ствования и единственности такого решения имеется в книге [1]. Задача 2.1 сводится к задаче квадратичного программирования. Раскроем по опре делению нормы: 𝜙(𝑥) = 1/2‖𝑥‖2 + max{−⟨𝑎𝑖 , 𝑥⟩ + 1/2‖𝑎𝑖 ‖2 } 𝑖∈𝑀 Теперь выражение под максимумом линейно по 𝑥, а значит сам максимум можно обо значить за 𝑡 и ввести набор ограничений. То есть, исходная задача эквивалентна задаче квадратичного программирования 1 ⟨𝐸𝑥, 𝑥⟩ + 𝑡 → min 2 1 ⟨𝑎𝑖 , 𝑥⟩ + 𝑡 > ‖𝑎𝑖 ‖2 , 𝑖 ∈ 𝑀 2

7 2.2. Общий вид задачи квадратичного программирования и двойственной к ней Чтобы двигаться дальше, необходимо сначала рассмотреть понятия задачи квад ратичного программирование и двойственной задачи. Задача квадратичного програм мирования — это экстремальная задача вида 𝑄(𝑥) = 1 ⟨𝐷𝑥, 𝑥⟩ + ⟨𝑐, 𝑥⟩ → min 2 𝐴[𝑀1 , 𝑁 ] × 𝑥[𝑁 ] > 𝑏[𝑀1 ] 𝐴[𝑀2 , 𝑁 ] × 𝑥[𝑁 ] = 𝑏[𝑀2 ] 𝑥[𝑁1 ] > O[𝑁1 ] где 𝐷 = 𝐷[𝑁, 𝑁 ] — симметричная неотрицательно-определённая матрица. Как и в слу чае линейного программирования, можно ввести двойственную задачу и получить тео рему, связывающую её с исходной задачей — первую теорему двойственности. При этом двойственная задача выглядит так: 1 − ⟨𝐷𝑣, 𝑣⟩ + ⟨𝑏, 𝑢⟩ → max 2 −𝑣[𝑁 ] × 𝐷[𝑁, 𝑁1 ] + 𝑢[𝑀 ] × 𝐴[𝑀, 𝑁1 ] 6 𝑐[𝑁1 ] −𝑣[𝑁 ] × 𝐷[𝑁, 𝑁2 ] + 𝑢[𝑀 ] × 𝐴[𝑀, 𝑁2 ] = 𝑐[𝑁2 ] 𝑢[𝑀1 ] > O[𝑀1 ] Для удобства записи двойственной задачи можно использовать таблицу 𝑐 𝑣 −𝐷 − 12 𝐷𝑣 𝑢 𝐴 𝑏 Тогда целевая функция получается скалярным умножением первого и последнего столб ца таблицы, ограничения — первого и второго, а правые части ограничений стоят в пер вой строке, причём неравенствами являются ограничения, соответствующие знаковым ограничениям исходной задачи. Знаковые ограничения накладываются на двойствен ные переменные (компоненты вектора 𝑢), соответствующие исходным ограничениям неравенствам.

8 2.3. Двойственная задача для задачи Сильвестра и её решение Вернёмся к задаче Сильвестра. Из векторов-строк 𝑎𝑖 составим матрицу 𝐴0 и введём стандартные для квадратичной задачи обозначения: матрицу квадратичной формы 𝐷, вектор коэффициентов линейной формы 𝑐, матрицу ограничений 𝐴 и столбец правых частей 𝑏. Получаем следующее: ⎞ ⎛ 0 ⎟ ⎜ .. ⎟ ⎜ (︁ ⎜ 𝐸[𝑁, 𝑁 ] . ⎟ ⎟ 𝑐 = 0 ··· 𝐷=⎜ ⎟ ⎜ ⎜ 0⎟ ⎠ ⎝ 0 ··· 0 0 ⎛ 0 1 )︁ ⎜ ⎜ 𝐴 = ⎜𝐴0 ⎝ ⎞ 1 ⎟ .. ⎟ .⎟ ⎠ 1 1 𝑏[𝑖] = ‖𝑎𝑖 ‖2 2 Переменной в такой задаче является пара 𝑦 = (𝑥, 𝑡), а в двойственной вектор 𝑧 = ((𝑣, 𝑠), 𝑢). Чтобы записать двойственную задачу, воспользуемся таблицей 0 ··· 0 1 0 . 𝑣[𝑁 ] −𝐸[𝑁, 𝑁 ] .. − 12 𝑣[𝑁 ] 0 𝑠 𝑢 0 ··· 0 0 0 𝐴0 1 .. . 𝑏 1 Тогда по описанному общему алгоритму имеем 1 𝑔(𝑧) := − ⟨𝑣, 𝑣⟩ + ⟨𝑏, 𝑢⟩ → max 2 −𝑣 + 𝑢𝐴0 = 0 ∑︁ 𝑢[𝑖] = 1 𝑖∈𝑀 𝑢[𝑀 ] > O[𝑀 ] Видно, что здесь можно исключить вектор 𝑣 и получить задачу минимизации на стан дартном симплексе 1 ℎ(𝑢) := ‖𝑢𝐴0 ‖2 − ⟨𝑏, 𝑢⟩ → min 2 ∑︁ 𝑢[𝑖] = 1 𝑖∈𝑀 𝑢[𝑀 ] > O[𝑀 ]

9 Для последней задачи существует оптимальный план 𝑢* . В [1] доказано, что в таком случае вектор 𝑥* = 𝑢* 𝐴0 является решением задачи Сильвестра. Так как два оптималь ных плана пары двойственных задач имеют одно и то же значение целевой функции, то 𝜙(𝑥* ) = −ℎ(𝑢* ) (на последнем шаге целевая функция была умножена на −1, что бы перейти от задачи максимизации к задаче минимизации). Значит, радиус искомого √︀ минимального шара равен −2ℎ(𝑢* ). 2.4. Решение задачи в среде MATLAB Для удобства решения задачи Сильвестра была написана функция sylvester, при нимающая единственный аргумент — матрицу 𝐴, строки которой образованы исходным множеством точек. Алгоритм, реализованный в функции, следующий: 1. Вычислить матрицу 𝐷 = 𝐴𝐴𝑇 и вектор 𝑏, соответствующие двойственной задаче. 2. Вызвать встроенную функцию quadprog для решения двойственной задачи и по лучить её решение — вектор 𝑢* , и значение целевой функции 𝑓 на нём. 3. Вычислить радиус минимального шара по формуле 𝑟 = √ −2𝑓 . Написанная функция работает с пространствами любых размерностей, но для на глядности она была проверена на размерности 2. Для этого генерировались случайные множества точек, лежащих внутри и на границе круга (𝑥 − 4)2 + (𝑦 − 3)2 6 4 и для них решалась задача Сильвестра. Одно из таких решений с количеством точек 100 представлено на рис. 2.1. Стоит отметить, что уже при попытке решить задачу для 200 точек, мы получаем ошибку: Maximum number of iterations exceeded; increase options.MaxIter. To continue solving the problem with the current solution as the starting point, set x0 = x before calling quadprog. Это связано с тем, что по умолчанию количество итераций, которые делает функция quadprog, ограничено числом 200. Чтобы решить эту проблему, можно поднять огра ничение. При этом оказывается, что количество итераций зависит от количества точек примерно линейно (см рис. 2.2).

10 Sylvester problem solution Input points 5 Solution Resulting circle 4.5 4 3.5 3 2.5 2 1.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 6.5 Рис. 2.1. Решение задачи Сильвестра для 𝑚 = 100 Такой результат можно объяснить тем, что по умолчанию MATLAB выбирает для решения этой задачи алгоритм «active-set». Если же задать использование алгоритма «interior-point-convex», пригодного для решения выпуклых задач, то для 100 точек тре буется 6 итераций, а для 500 — всего 9. Также можно протестировать программу в размерности 3. Для этого было сге нерировано 100 случайных точек внутри и на границе стандартного куба [−1, 1]3 . На этом наборе точек программа выдает ожидаемый ответ — центр шара близок к 0, а его √ радиус — к 3.

11 220 200 180 Number of iterations 160 140 120 100 80 60 40 20 0 20 40 60 80 100 120 Number of points 140 160 180 200 Рис. 2.2. Зависимость числа итераций от числа точек

12 Глава 3 Алгоритм Шора Перейдём теперь к собственно задаче о минимальном эллипсоиде. В отличие от задачи Сильвестра, которую можно решить за конечное число шагов с помощью квад ратичного программирования, эту задачу удаётся решить лишь приближённо итератив ными методами. В этом разделе будет рассмотрен алгоритм Шора, строящий прибли жение к минимальному эллипсоиду. 3.1. Оператор сжатия пространства Алгоритм построения минимального эллипсоида использует оператор сжатия про странства, поэтому рассмотрим сначала его свойства. ОПРЕДЕЛЕНИЕ. Пусть 𝜉 ∈ R𝑛 , ‖𝜉‖ = 1 и 𝛾 ∈ (0, 1). Тогда оператор сжатия 𝑅𝛾 определяется формулой 𝑅𝛾 = 𝐸 − (1 − 𝛾)𝜉𝜉 𝑇 Здесь 𝜉𝜉 𝑇 — это матрица проектирования на прямую 𝑥 = 𝜆𝜉, 𝜆 ∈ R. Вычисле ние значения 𝑅𝛾 на конкретном векторе 𝑦 ∈ R𝑛 позволяет понять, как действует этот оператор: 𝑅𝛾 𝑦 = (𝐸 − (1 − 𝛾)𝜉𝜉 𝑇 )𝑦 = 𝑦 − (1 − 𝛾)𝜉𝜉 𝑇 𝑦 = 𝑦 − (1 − 𝛾) · ⟨𝜉, 𝑦⟩ · 𝜉 Из этого равенства и определения 𝑅𝛾 очевидны следующие свойства: 1. 𝑅𝛾 — симметричная матрица 2. 𝑅𝛾 𝜉 = 𝛾𝜉 3. 𝑅𝛾 𝑝 = 𝑝, если ⟨𝜉, 𝑝⟩ = 0 Значит, 𝑅𝛾 имеет собственное число 𝜆0 = 𝛾 кратности 1 с собственным вектором 𝜉 и соб ственное число 𝜆1 = 1 кратности 𝑛−1 с собственным подпространством, ортогональным 𝜉. Если 𝑦 = 𝛼𝜉 + 𝑝, где ⟨𝜉, 𝑝⟩ = 0, то 𝑅𝛾 𝑦 = 𝛾𝛼𝜉 + 𝑝. Кроме этого, матрица 𝑅𝛾 положительно определена при всех 𝛾 ∈ (0, 1).

13 𝑦 ⟨𝜉, 𝑦⟩ · 𝜉 (1 − 𝛾) ⟨𝜉, 𝑦⟩ · 𝜉 𝜉 𝑅𝛾 𝑦 O Рис. 3.1. Оператор сжатия с 𝛾 = 0.2 ЛЕММА 1. Оператор 𝑅𝛾 обратим и обратный к нему оператор задаётся формулой 𝑅𝛾−1 = 𝐸 + 1−𝛾 𝑇 𝜉𝜉 𝛾 Доказательство. Вычислим произведение 𝑅𝛾 · 𝑅𝛾−1 : (︂ )︂ (︀ )︀ 1−𝛾 𝑇 𝑇 𝐸 − (1 − 𝛾)𝜉𝜉 𝜉𝜉 𝐸+ = 𝛾 (1 − 𝛾)2 𝑇 𝑇 1−𝛾 𝑇 𝜉𝜉 − (1 − 𝛾)𝜉𝜉 𝑇 − 𝜉𝜉 𝜉𝜉 = =𝐸+ ⏟ ⏞ 𝛾 𝛾 =1 )︀ 1 − 𝛾 (︀ 𝑇 =𝐸+ 𝜉𝜉 − 𝛾𝜉𝜉 𝑇 − (1 − 𝛾)𝜉𝜉 𝑇 = 𝐸 𝛾 К обратному оператору применимы те же рассуждения о собственных числах и собственных векторах, в частности, вектор 𝜉 собственный с собственным числом 𝛾1 . Теперь поймём, как оператор сжатия действует на шары и эллипсоиды в R𝑛 . Пусть, для начала, имеется шар с центром в точке 𝑐 и радиусом 𝑟, который задан уравнением ⟨𝑥 − 𝑐, 𝑥 − 𝑐⟩ 6 𝑟2 Подействуем на него оператором 𝑅𝛾 , получив новые точки 𝑦 = 𝑅𝛾 𝑥 и новый центр 𝑎 = 𝑅𝛾 𝑐. Так как 𝑅𝛾 обратим, то 𝑥 = 𝑅𝛾−1 𝑦 и 𝑐 = 𝑅𝛾−1 𝑎. Подставим эти выражения в уравнение шара ⟨𝑥 − 𝑐, 𝑥 − 𝑐⟩ 6 𝑟2 ⟨︀ −1 ⟩︀ 𝑅𝛾 (𝑦 − 𝑎), 𝑅𝛾−1 (𝑦 − 𝑎) 6 𝑟2 ⟨︀ −2 ⟩︀ (𝑅𝛾 )(𝑦 − 𝑎), 𝑦 − 𝑎 6 𝑟2 ⟨ −2 ⟩ 𝑅𝛾 (𝑦 − 𝑎), 𝑦 − 𝑎 6 1 𝑟2

14 Получилось уравнение эллипсоида с центром в точке 𝑎 = 𝑅𝛾 𝑐 и матрицей 𝑀 := 𝑅𝛾−2 . 𝑟2 Как было установлено выше, матрица 𝑀 имеет собственный вектор 𝜉, которому соот ветствует собственное число имеют собственные числа 1 . 𝑟2 1 , 𝑟2 𝛾 2 а остальные собственные векторы ортогональны 𝜉 и Таким образом, оператор сжатия превращает шар радиуса 𝑟 в эллипсоид, одной из полуосей которого является 𝜉 с длиной 𝑟𝛾, а длины остальных полуосей равны 𝑟 — шар «сплющивается» в направлении 𝜉. Разобравшись с шаром, подействуем теперь оператором 𝑅𝛾 на эллипсоид, задан ный уравнением ⟨𝑀 (𝑥 − 𝑐), 𝑥 − 𝑐⟩ 6 1 Положим аналогично 𝑥 = 𝑅𝛾−1 𝑦, 𝑐 = 𝑅𝛾−1 𝑎 и получим новое уравнение эллипсоида в виде ⟨︀ −1 ⟩︀ 𝑅𝛾 𝑀 𝑅𝛾−1 (𝑦 − 𝑎), 𝑦 − 𝑎 6 1 Наглядное представление о работе оператора сжатия даёт рис. 3.2. 3.2. Итеративный алгоритм построения минимального эллипсоида Для задачи о минимальном эллипсоиде известен итеративный алгоритм, строящий последовательные приближения к ответу. Этот алгоритм, описанный Н.З. Шором в статье [2], основан на достаточно простых идеях. Введём вспомогательную задачу: вложим точек 𝑎𝑖 в гиперплоскость ⎛ множество ⎞ 𝑎𝑖 𝑥𝑛+1 = 1 пространства R𝑛+1 . Обозначим 𝑞𝑖 = ⎝ ⎠. В новом пространстве будем искать 1 минимальный эллипсоид с фиксированным центром в нуле. В статье [2] утверждается, что сечение такого эллипсоида плоскостью 𝑥𝑛+1 = 1 и будет решением исходной задачи (см рис. 3.3). ̃︀ порядка 𝑛 + 1: Поставим формально вспомогательную задачу в терминах матрицы 𝐵 ̃︀ := − ln det 𝐵 ̃︀ → min 𝑓 (𝐵) ⟨ ⟩ ̃︀ 𝑖 , 𝑞𝑖 6 1, 𝑖 ∈ 1 : 𝑚 𝐵𝑞 ̃︀ ≻ 0. 𝐵 Рассмотрим вопрос построения сечения подробнее. Итак, пусть {︁ ⟨ ⟩ }︁ ˜ 𝑥˜, 𝑥˜ 6 1 ℰ̃︀ = 𝑥˜ ∈ R𝑛+1 | 𝐵 (3.1)

15 1 1 0.8 0.8 0.6 0.6 0.4 0.4 0.2 0.2 0 0 −0.2 −0.2 −0.4 −0.4 −0.6 −0.6 −0.8 −0.8 −1 −1 −0.5 0 0.5 −1 −1 1 −0.5 0 0.5 1 (a) Простейший случай — сжатие единичного круга в направлении (1, 1) с 𝛾 = 0.2 3 3 2.5 2.5 2 2 1.5 1.5 1 1 0.5 0.5 0 −1 0 −0.5 0 0.5 1 1.5 2 (b) Сжатие круга в направлении (1, 1) с 𝛾 = 0.2 0 0.5 1 1.5 2 (c) Сжатие эллипса в направле нии (0, 1) с 𝛾 = 0.4 Рис. 3.2. Действие оператора сжатия на множество, выделенное красной линией. Результат выделен зелёным цветом

16 𝑥2 1 𝑥1 Рис. 3.3. Сечение эллипса с центром в нуле прямой 𝑥2 = 1 — минимальный эллипсоид, содержащий точки 𝑞𝑖 , 𝑖 ∈ 1 : 𝑚. Возьмём две различные точки 𝑞𝑖0 и 𝑞𝑖1 . Отметим, что ⟨ ⟩ ⟨ ⟩ ⟨ ⟩ ⟨ ⟩ ̃︀ 𝑖0 + 𝑞𝑖1 ), 𝑞𝑖0 + 𝑞𝑖1 + 𝐵(𝑞 ̃︀ 𝑖0 − 𝑞𝑖1 ), 𝑞𝑖0 − 𝑞𝑖1 = 2 𝐵𝑞 ̃︀ 𝑖0 , 𝑞𝑖0 + 2 𝐵𝑞 ̃︀ 𝑖1 , 𝑞𝑖1 𝐵(𝑞 Обозначим ⎛ ⎞ 𝑥0 1 𝑥˜0 = (𝑞𝑖0 + 𝑞𝑖1 ) = ⎝ ⎠ , 2 1 (3.2) ⎛ ⎞ 𝑥1 1 𝑥˜1 = (𝑞𝑖0 − 𝑞𝑖1 ) = ⎝ ⎠ 2 0 Поделив равенство (3.2) на 4, получим ⟩ ⟨ ⟩ ⟨ ⟩ ⟨ ⟩ ⟨ ̃︀ 𝑖0 , 𝑞𝑖0 + 1 𝐵𝑞 ̃︀ 𝑖1 , 𝑞𝑖1 − 𝐵 ̃︀ 𝑥˜1 , 𝐵 ˜ 𝑥˜1 ̃︀ 𝑥˜0 , 𝑥˜0 = 1 𝐵𝑞 𝐵 2 2 По определению точки 𝑞𝑖0 и 𝑞𝑖1 принадлежат ℰ̃︀ и 𝑥˜1 ̸= O, поэтому ⟨ ⟩ ̃︀ 𝑥˜0 , 𝑥˜0 < 1 𝐵 (3.3) ˜ в виде Представим матрицу 𝐵 ⎞ 𝐵 𝑏 ˜=⎝ ⎠ 𝐵 𝑏𝑇 ˆ𝑏 ⎛ где 𝐵 — главный минор порядка 𝑛. Очевидно, что 𝐵 симметрична и положительно определена (по критерию Сильвестра). Аккуратно раскроем скобки в (3.3): ⎞ ⎛ ⎞⟩ ⟨⎛ ⟨ ⟩ 𝐵𝑥0 + 𝑏 𝑥 ̃︀ 𝑥˜0 , 𝑥˜0 = ⎝ ⎠ , ⎝ 0 ⎠ = ⟨𝐵𝑥0 + 𝑏, 𝑥0 ⟩ + ⟨𝑏, 𝑥0 ⟩ + ˆ𝑏 = 𝐵 ⟨𝑏, 𝑥0 ⟩ + ˆ𝑏 1 ⟨︀ ⟩︀ ⟨︀ ⟩︀ = 𝐵(𝑥0 + 𝐵 −1 𝑏), 𝑥0 + 𝐵 −1 𝑏 − 𝐵(𝑥0 + 𝐵 −1 𝑏), 𝐵 −1 𝑏 + ⟨𝑏, 𝑥0 ⟩ + ˆ𝑏 = ⟨︀ ⟩︀ ⟨︀ ⟩︀ = 𝐵(𝑥0 + 𝐵 −1 𝑏), 𝑥0 + 𝐵 −1 𝑏 − 𝐵 −1 𝑏, 𝑏 + ˆ𝑏 < 1 ЛЕММА 2. Справедливы неравенства ⟨︀ ⟩︀ 0 < − 𝐵 −1 𝑏, 𝑏 + ˆ𝑏 < 1 (3.4)

17 ⎛ Доказательство. Пусть, во-первых, 𝑥˜ = ⎝ ⟨ −𝐵 −1 𝑏 1 ⎞ ˜ положительно ⎠ ̸= O. Матрица 𝐵 ⟩ ˜ 𝑥˜, 𝑥˜ > 0. Рассуждения, аналогичные (3.4), приводят к равен определена, так что 𝐵 ⟨ ⟩ ˜ 𝑥˜, 𝑥˜ = − ⟨𝐵 −1 𝑏, 𝑏⟩ + ˆ𝑏, что доказывает левое неравенство. ству 𝐵 Для доказательство правого неравенства заметим, что в неравенстве (3.4) ⟨𝐵(𝑥 + 𝐵 −1 𝑏), 𝑥 + 𝐵 −1 𝑏⟩ > 0 в силу положительной определённости 𝐵. Тогда имеем ⟨︀ ⟩︀ ⟨︀ ⟩︀ ⟨︀ ⟩︀ 1 > 𝐵(𝑥 + 𝐵 −1 𝑏), 𝑥 + 𝐵 −1 𝑏 − 𝐵 −1 𝑏, 𝑏 + ˆ𝑏 > − 𝐵 −1 𝑏, 𝑏 + ˆ𝑏 Неравенство доказано. Теперь последнее неравенство из (3.4) поделим на положительное число 1+⟨𝐵 −1 𝑏, 𝑏⟩− ˆ𝑏, введём обозначения 𝑐 := −𝐵 −1 𝑏, ⟨︀ ⟩︀ 𝑀 := 𝐵/(1 + 𝐵 −1 𝑏, 𝑏 − ˆ𝑏) (3.5) и окончательно получим ⟨𝑀 (𝑥 − 𝑐), 𝑥 − 𝑐⟩ 6 1 — уравнение минимального эллипсоида в R𝑛 Алгоритм решения задачи (3.1) устроен следующим образом. Текущее множество (𝑘) точек 𝑞𝑖 хранится в виде столбцов матрицы 𝑋𝑘 , текущая матрица последовательно сти операторов сжатия обозначается 𝐴𝑘 . Изначально 𝑋0 состоит из исходных точек 𝑞𝑖 ∈ R𝑛+1 , а 𝐴0 = 𝐸𝑛+1 . Один шаг работы алгоритма состоит из нескольких простых действий: (𝑘) ∙ Найти среди векторов 𝑞𝑖 (𝑘) (𝑘) вектор 𝑞𝑖𝑘 с максимальной нормой 𝜏𝑘+1 = ‖𝑞𝑖𝑘 ‖ (𝑘) ∙ Вычислить единичный вектор направления сжатия 𝜉𝑘+1 𝑞 = 𝑖0 𝜏𝑘+1 ∙ Построить оператор сжатия 𝑅𝑘+1 вдоль вектора 𝜉𝑘+1 с коэффициентом 𝛼𝑘+1 (𝑘+1) ∙ Вычислить новые точки 𝑞𝑖 , умножив матрицу 𝑅𝑘+1 на матрицу 𝑋𝑘 , и сохранить результат в 𝑋𝑘+1 ∙ Добавить матрицу 𝑅𝑘+1 к последовательности операторов сжатия, умножив её на матрицу 𝐴𝑘 , и сохранить результат в 𝐴𝑘+1

18 Коэффициенты 𝛼𝑘 в соответствии со статьей [2] выбираются из условий 𝛼𝑘 = 1 − 𝛽𝑘 , 𝛽𝑘 > 0, 𝛽𝑘 −−−→ 0, 𝑘→∞ ∞ ∑︁ 𝛽𝑘 = ∞ 𝑘=1 Пусть вычисления остановлены после шага 𝑚. Тогда из описания алгоритма сле дует, что 𝑋𝑚 = 𝑅𝑚 · 𝑅𝑚−1 · · · 𝑅1 · 𝑋0 = 𝐴𝑚 · 𝑋0 (𝑚) Также известно, что точки 𝑞𝑖 лежат внутри сферы с радиусом 𝜏𝑚 и центром в начале координат. С учётом этих двух фактов можно получить матрицу искомого эллипсоида: ⟨ ⟩ (𝑚) (𝑚) 2 𝑞𝑖 , 𝑞𝑖 6 𝜏𝑚 2 ⟨𝐴𝑚 𝑞𝑖 , 𝐴𝑚 𝑞𝑖 ⟩ 6 𝜏𝑚 ⟩︀ ⟨︀ 𝑇 2 𝐴𝑚 𝐴𝑚 𝑞𝑖 , 𝑞𝑖 6 𝜏𝑚 ⟨ 𝑇 ⟩ 𝐴𝑚 𝐴𝑚 𝑞𝑖 , 𝑞𝑖 6 1 2 𝜏𝑚 ̃︀ := 𝐴𝑇𝑚 𝐴𝑚 /𝜏 2 симметрична, положительно определена и является (при Матрица 𝐵 𝑚 ближённым) решением вспомогательной задачи о минимальном эллипсоиде. Оконча тельное решение основной задачи получается по формулам (3.5). 3.3. Практическое испытание алгоритма Алгоритм был реализован в среде MATLAB и испытан на тестовых данных. Вычис ления производились в два этапа: пока нормы 𝜏𝑘 меняются больше, чем на 𝜀1 , использо вался постоянный коэффициент 𝛽𝑘 = 0.2. На втором этапе использовались убывающие коэффициенты 𝛽𝑘 = 1/(𝑘 + 1), пока нормы 𝜏𝑘 меняются больше, чем на 𝜀2 . В качестве тестовых данных было сгенерировано множество случайных точек внутри эллипсоида с центром в точке (1, 2) с полуосями 2, 1, который вытянут в направлении (−1, 1). Резуль таты вычисления при различных 𝜀1 и 𝜀2 приведены в таблице 3.1. В столбце «точность» содержится относительная ошибка 𝛿, то есть ⃒ ⃒ ⃒(det 𝑀 )−1/2 − (det 𝑀 * )−1/2 ⃒ 𝛿= , (det 𝑀 * )−1/2 где 𝑀 * — известная матрица минимального эллипсоида. Как видно из таблицы, при 𝜀1 = 10−6 и 𝜀2 = 10−10 достигается наилучшее время работы, но наихудшая точность. Подобрать параметры 𝜀1 , 𝜀2 так, чтобы обеспечить и высокую точность, и низкое время работы, не удалось.

19 4 3.5 3 2.5 2 1.5 1 0.5 0 −1 −0.5 0 0.5 1 1.5 2 2.5 3 Рис. 3.4. Тестовое множество точек. Синим отмечена граница минимального эллипсоида Таблица 3.1. Результаты алгоритма Шора при разных 𝜀1 и 𝜀2 . 𝜀1 𝜀2 10−5 точность число итераций время (с) 10−10 7 · 10−5 21990 2.4 10−4 10−9 8 · 10−5 20591 2.3 10−6 10−10 4 · 10−4 4121 0.4 10−4 10−10 2 · 10−5 100000 21.2 Кроме того, в таблице 3.2 приведены результаты вычисления при 𝜀1 = 10−6 , 𝜀2 = 10−10 и различных комбинациях 𝑛 и 𝑚. Таблица 3.2. Результаты алгоритма Шора при разных 𝑛 и 𝑚. 𝑛 𝑚 число итераций точность время (с) 2 104 4121 4 · 10−4 0.2 2 504 4118 5 · 10−4 0.2 5 510 16729 8 · 10−4 1.1 10 1020 3533 3 · 10−3 0.6

20 Глава 4 Алгоритм Хачияна Как видно из таблиц, алгоритм Шора не достигает высоких результатов по точно сти полученного эллипсоида. Но, несмотря на это, в нём содержится удачная идея — решение задачи в расширенном пространстве (3.1). В этой главе будет рассматриваться ещё один алгоритм решения расширенной задачи. Решение исходной задачи будет получаться по формулам (3.5). 4.1. Вспомогательные сведения. Пусть 𝑞 — вектор, а 𝑀 — симметричная матрица соответствующей размерности. Найдём другое представление для 𝑞 | 𝑀 𝑞, используя индексную технику: ⟨ 𝑛 ⟩ 𝑛 ∑︁ ∑︁ 𝑞[𝑘] ⟨𝑀 [1 : 𝑛, 𝑘], 𝑞⟩ = 𝑞 | 𝑀 𝑞 = ⟨𝑀 𝑞, 𝑞⟩ = 𝑀 [1 : 𝑛, 𝑘]𝑞[𝑘], 𝑞 = 𝑘=1 𝑘=1 = = 𝑛 ∑︁ 𝑘=1 𝑛 ∑︁ (︃ 𝑞[𝑘] 𝑛 ∑︁ )︃ 𝑀 [𝑗, 𝑘]𝑞[𝑗] 𝑗=1 = 𝑛 𝑛 ∑︁ ∑︁ 𝑀 [𝑗, 𝑘] (𝑞[𝑗]𝑞[𝑘]) = (4.1) 𝑘=1 𝑗=1 𝑀 [𝑗, 𝑘] (𝑞𝑞 | ) [𝑗, 𝑘] = tr(𝑀 · 𝑞𝑞 | ). 𝑘=1 Нам понадобится ещё одно равенство, называемое «лемма об определителе матри цы» в англоязычной литературе [3]. Пусть 𝐴 — обратимая матрица, 𝑢, 𝑣 – векторы. Тогда det(𝐴 + 𝑢𝑣 | ) = (1 + 𝑣 | 𝐴−1 𝑢) · det 𝐴. (4.2) Прежде чем проверить это равенство, заметим, что det(𝐴 + 𝑢𝑣 | ) = det 𝐴 · det(𝐸 + (𝐴−1 𝑢)𝑣 | ). Следовательно, достаточно проверить только случай 𝐴 = 𝐸. В этом случае необходимое равенство следует из свойств определителя и равенства ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞ 𝐸 0 𝐸 + 𝑢𝑣 | 𝑢 𝐸 0 𝐸 𝑢 ⎝ ⎠·⎝ ⎠·⎝ ⎠=⎝ ⎠. | | | 𝑣 1 0 1 −𝑣 1 0 1+𝑣 𝑢 Действительно, определитель матрицы справа равен 1 + 𝑣 | 𝑢, а определитель средней матрицы в левой части равен det(𝐸 +𝑢𝑣 | ). Кроме того, определители остальных матриц равны 1. Значит, равенство доказано.

21 Рассмотрим теперь общую задачу выпуклого программирования с ограничениями: 𝑓 (𝑥) → min ℎ𝑖 (𝑥) 6 0, 𝑖 ∈ 1 : 𝑚, 𝑥 ∈ 𝐷, где 𝐷 – открытое выпуклое множество. Введём неотрицательные числа 𝜆𝑖 , соответ ствующие ограничениям задачи. Тогда функцией Лагранжа [4] этой задачи называется функция ℒ(𝑥, 𝜆) := 𝑓 (𝑥) + 𝑚 ∑︁ 𝜆𝑖 ℎ𝑖 (𝑥), 𝑖=1 а числа 𝜆𝑖 называются множителями Лагранжа. Используя функцию Лагранжа, можно записать двойственную по Лагранжу зада чу: 𝑔(𝜆) := inf ℒ(𝑥, 𝜆) → max 𝑥∈𝐷 𝜆𝑖 > 0, 𝑖 ∈ 1 : 𝑚. Для любой пары векторов 𝑥 и 𝜆, удовлетворяющих прямой и двойственной задачам соответственно, верно 𝑓 (𝑥) > 𝑔(𝜆). Отсюда, если 𝑥* и 𝜆* — решения пары двойственных задач, то 𝑓 (𝑥* ) > 𝑔(𝜆* ). 4.2. Вывод двойственной задачи для задачи о минимальном эллипсоиде. Попробуем применить двойственность Лагранжа к задаче (3.1). Введём множители ⟨ ⟩ ̃︀ 𝑖 , 𝑞𝑖 − 1 6 0, а ограничение Лагранжа 𝑝𝑖 > 0, соответствующие ограничениям 𝐵𝑞 𝑀 ≻ 0 будем считать нефункциональным и включим в область определения ℒ. Тогда функция Лагранжа записывается как ̃︀ 𝑝) := − ln det 𝐵 ̃︀ + ℒ(𝐵, 𝑚 ∑︁ (︀ )︀ ̃︀ 𝑖 − 1 . 𝑝𝑖 𝑞𝑖| 𝐵𝑞 𝑖=1 Или, с учётом формулы (4.1), ̃︀ 𝑝) := − ln det 𝐵 ̃︀ + ℒ(𝐵, 𝑚 ∑︁ 𝑖=1 ̃︀ − 𝑝𝑖 tr(𝑞𝑖 𝑞𝑖| · 𝐵) 𝑚 ∑︁ 𝑝𝑖 . 𝑖=1 {︁ }︁ ̃︀ 𝑝)|𝐵 ̃︀ ≻ 0, 𝑝 ∈ R𝑚 . Здесь и далее будем считать, что dom ℒ = (𝐵,

22 По определению целевая функция двойственной к (3.1) задачи равна [︃ ]︃ 𝑚 𝑚 ∑︁ ∑︁ ̃︀ 𝑝) = inf − ln det 𝐵 ̃︀ + ̃︀ − 𝑔(𝑝) := inf ℒ(𝐵, 𝑝𝑖 tr(𝑞𝑖 𝑞 | · 𝐵) 𝑝𝑖 . ̃︀ 𝐵≻0 𝑖 ̃︀ 𝐵≻0 𝑖=1 𝑖=1 Для того, чтобы упростить это выражение, найдём инфимум аналитически. Известно ̃︀ и tr(𝑞𝑖 𝑞 | · 𝐵) ̃︀ выпуклы по 𝐵 ̃︀ на открытом множестве (см. [4]), что функции − ln det 𝐵 𝑖 положительно определённых матриц. Значит, инфимум достигается в точке, где произ ̃︀ выражения под ним равна нулю. Возьмём эту производную: водная по 𝐵 𝑚 ∑︁ ̃︀ 𝑝) 𝜕ℒ(𝐵, ̃︀ −1 + ̃︀ −1 + 𝑄𝑃 𝑄| , = −𝐵 𝑝𝑖 𝑞𝑖 𝑞𝑖| = −𝐵 ̃︀ 𝜕𝐵 𝑖=1 где 𝑄 — матрица, составленная из точек 𝑞𝑖 , а 𝑃 = diag(𝑝). Для дальнейших рассуждений необходимо, чтобы матрица 𝑄𝑃 𝑄| была положительно определена. Проверим, так ли это. Зафиксируем произвольный ненулевой вектор 𝑥 ∈ R𝑛+1 . По определению, для по ложительной определённости 𝑄𝑃 𝑄| необходимо, чтобы ⟨𝑄𝑃 𝑄| 𝑥, 𝑥⟩ > 0. Распишем ска лярное произведение: ⟨ 𝑚 ∑︁ ⟩ 𝑝𝑖 𝑞𝑖 𝑞𝑖| · 𝑥, 𝑥 = 𝑖=0 𝑚 ∑︁ 𝑝𝑖 ⟨𝑞𝑖 𝑞𝑖| · 𝑥, 𝑥⟩ = 𝑖=1 𝑚 ∑︁ 𝑝𝑖 ⟨𝑞𝑖 , 𝑥⟩2 𝑖=1 Так как среди точек 𝑞𝑖 содержится аффинно независимое множество, то вектор 𝑥 не может быть ортогонален сразу всем векторам 𝑞𝑖 . Значит, если для всех 𝑖 𝑝𝑖 > 0, то матрица 𝑄𝑃 𝑄| заведомо положительно определена. На самом деле, достаточно, чтобы числа 𝑝𝑖 были положительны на индексах 𝑖, соответствующих афинно независимому множеству точек 𝑞𝑖 . Будем теперь рассматривать только такие 𝑝, что 𝑄𝑃 𝑄| ≻ 0. В следующей части будет описан алгоритм, на каждом шаге которого выполнено условие 𝑝 > 0. ̃︀ 𝑝) 𝜕ℒ(𝐵, = 0 при фиксирован ̃︀ 𝜕𝐵 ̃︀ = (𝑄𝑃 𝑄| )−1 . Подставим это значение в ℒ(𝐵, ̃︀ 𝑝) и получим ном 𝑝 является матрица 𝐵 Таким образом, единственным решением уравнения выражение для целевой функции двойственной задачи 𝑔(𝑝): 𝑔(𝑝) = ln det 𝑄𝑃 𝑄| + 𝑚 ∑︁ 𝑚 [︀ ]︀ ∑︁ 𝑝𝑖 tr 𝑞𝑖| 𝑞𝑖 · (𝑄𝑃 𝑄| )−1 − 𝑝𝑖 . 𝑖=1 𝑖=1 Обозначим 𝐴 = 𝑄𝑃 𝑄| и распишем выражение под первой суммой: (︃ 𝑚 )︃ 𝑚 ∑︁ ∑︁ [︀ | ]︀ 𝑝𝑖 tr 𝑞𝑖 𝑞𝑖 · (𝑄𝑃 𝑄| )−1 = tr 𝑝𝑖 𝑞𝑖 𝑞𝑖| 𝐴−1 = tr 𝐴𝐴−1 = 𝑛 + 1 𝑖=1 𝑖=1 (4.3)

23 В итоге имеем | 𝑔(𝑝) = ln det 𝑄𝑃 𝑄 − 𝑚 ∑︁ 𝑝𝑖 + 𝑛 + 1 (4.4) 𝑖=1 Значит, двойственная задача имеет вид | 𝑔(𝑝) := ln det 𝑄𝑃 𝑄 − 𝑚 ∑︁ 𝑝𝑖 + 𝑛 + 1 → max (4.5) 𝑖=1 𝑝𝑖 > 0, 𝑖 ∈ 1 : 𝑚. Чтобы ещё упростить вид целевой функции, запишем для этой задачи условия Куна-Таккера. Чтобы сделать это, посчитаем сначала производную ln det 𝑄𝑃 𝑄| , акку ратно применив правило цепочки. (︂ )︂ (︀ )︀ 𝜕 ln det 𝑄𝑃 𝑄| 𝜕𝑄𝑃 𝑄| 𝜕 | | | −1 ln det 𝑄𝑃 𝑄 = tr = tr (𝑄𝑃 𝑄 ) · 𝑞 𝑞 = 𝑞𝑖| (𝑄𝑃 𝑄| )−1 𝑞𝑖 · 𝑖 𝑖 𝜕𝑝𝑖 𝜕𝑄𝑃 𝑄| 𝜕𝑝𝑖 Следовательно, условия Куна-Таккера для (4.5) с множителями 𝜆𝑖 , соответствующими ограничениям 𝑝𝑖 > 0, выглядят как 𝑞𝑖| (𝑄𝑃 𝑄| )−1 𝑞𝑖 − 1 + 𝜆𝑖 = 0 𝜆𝑖 𝑝𝑖 = 0 𝜆𝑖 > 0 Избавимся от переменных 𝜆𝑖 : 𝑞𝑖| (𝑄𝑃 𝑄| )−1 𝑞𝑖 6 1 (︀ )︀ 𝑝𝑖 1 − 𝑞𝑖| (𝑄𝑃 𝑄| )−1 𝑞𝑖 = 0 (4.6) (4.7) Просуммировав (4.7) по 𝑖 и воспользовавшись ещё раз формулой (4.3) получаем 𝑚 ∑︁ 𝑝𝑖 = 𝑛 + 1. 𝑖=1 Теперь двойственную задачу можно записать как 𝑔(𝑝) := ln det 𝑄𝑃 𝑄| → max 1| 𝑝 = 𝑛 + 1 𝑝 > 0. С помощью замены переменных 𝑢 := 𝑝 1 , 𝑈 := diag(𝑢) = 𝑃 получаем итого 𝑛+1 𝑛+1 вый вид двойственной задачи 𝑔(𝑢) ̂︀ := ln det 𝑄𝑈 𝑄| → max 1| 𝑢 = 1 𝑢 > 0. (4.8)

24 Таким образом, как и в задаче Сильвестра, мы свели нахождение минимального эллипсоида к максимизации функции на стандартном симплексе в R𝑚 . Пусть найдено решение двойственной задачи 𝑢* . Из того, что в точке 𝑢* дости гается максимум функции 𝑔, ̂︀ следует, что 𝑔(𝑢 ̂︀ * ) > −∞, а значит 𝑄𝑈 * 𝑄| положитель ̃︀ * = (𝑛 + 1) · (𝑄𝑈 * 𝑄| )−1 . В силу (4.6) эта но определена. Положим 𝑝* = (𝑛 + 1)𝑢* , 𝐵 матрица будет удовлетворять ограничениям задачи (3.1). Очевидно, что выполняет ̃︀ * ) = 𝑔(𝑝* ). В силу неравенства между целевыми функциями пары ся равенство 𝑓 (𝐵 ̃︀ удовлетворяющей ограничениям, верно двойственных задач, для любой матрицы 𝐵, ̃︀ > 𝑔(𝑝* ) = 𝑓 (𝐵 ̃︀ * ), следовательно, 𝐵 ̃︀ * — решение задачи (3.1). 𝑓 (𝐵) 4.3. Алгоритм Хачияна. В статье [5] описан простой алгоритм, решающий двойственную задачу (4.8). При ведём его вывод. Во-первых, заметим, что формула (4.6) для изменённой задачи (4.8) имеет вид 𝑞𝑖| (𝑄𝑈 𝑄| )−1 𝑞𝑖 6 𝑛 + 1 Как было вычислено выше, 𝑔̂︀𝑗 (𝑢) := 𝜕 𝑔(𝑢) ̂︀ = 𝑞𝑗| (𝑄𝑈 𝑄| )−1 𝑞𝑗 . Значит, если 𝑢* — реше 𝜕𝑢𝑗 ние, то 𝑔̂︀𝑗 (𝑢* ) 6 𝑛 + 1, ∀𝑗 ∈ 1 : 𝑚. (4.9) Пусть на 𝑘-м шаге имеется вектор 𝑢𝑘 . Найдём, с какой ошибкой 𝑢𝑘 удовлетворяет условию (4.9). То есть, найдём такое минимальное число 𝜀, что 𝑔̂︀𝑗 (𝑢𝑘 ) 6 (1 + 𝜀)(𝑛 + 1), ∀𝑗 ∈ 1 : 𝑚. Пусть 𝑅 ⊂ 1 : 𝑚 — множество индексов, на которых неравенство обращается в равен ство. Тогда для 𝑗 ̸∈ 𝑅 𝑔̂︀𝑗 (𝑢𝑘 ) < (1 + 𝜀)(𝑛 + 1) = 𝑔̂︀𝑟 (𝑢𝑘 ), ∀𝑟 ∈ 𝑅. Выберем какой-нибудь 𝑟 ∈ 𝑅. Из формулы видно, что 𝑔̂︀𝑟 (𝑢𝑘 ) = max {̂︀ 𝑔𝑗 (𝑢𝑘 )|𝑗 ∈ 1 : 𝑚} и 𝜀= 𝑔̂︀𝑟 (𝑢𝑘 ) − (𝑛 + 1) . 𝑛+1 Хачиян предлагает следующее: в качестве 𝑢𝑘+1 выбрать точку из отрезка [𝑢𝑘 , 𝑒𝑟 ], где 𝑒𝑟 — 𝑟-й орт R𝑚 , максимизирующую 𝑔(𝑢 ̂︀ 𝑘+1 ). То есть, выбрать число 𝛼 ∈ (0, 1) и поло жить 𝑢𝑘+1 = (1 − 𝛼)𝑢𝑘 + 𝛼𝑒𝑟 .

25 Покажем, как найти такое 𝛼. Для этого распишем 𝑔(𝑢 ̂︀ 𝑘+1 ): ]︀ [︀ ]︀ [︀ 𝑔(𝑢 ̂︀ 𝑘+1 ) = ln det 𝑄 ((1 − 𝛼)𝑈𝑘 + 𝛼𝐸𝑟 ) 𝑄| = ln det (1 − 𝛼)𝑄𝑈𝑘 𝑄| + 𝛼𝑞𝑟 𝑞𝑟| = [︂ (︂ )︂]︂ 𝛼 𝑛+1 | | = ln (1 − 𝛼) det 𝑄𝑈𝑘 𝑄 + 𝑞𝑟 𝑞 . 1−𝛼 𝑟 Применим формулу (4.2): )︂ (︂ )︂ (︂ 𝛼 | 𝛼 | −1 | | det 𝑄𝑈𝑘 𝑄 + 𝑞𝑟 𝑞 = 1 + 𝑞 (𝑄𝑈𝑘 𝑄 ) 𝑞𝑟 det 𝑄𝑈𝑘 𝑄| . 1−𝛼 𝑟 1−𝛼 𝑟 Таким образом, ]︂ )︂ 𝛼 | | −1 | 𝑔(𝑢 ̂︀ 𝑘+1 ) = ln (1 − 𝛼) 𝑞 (𝑄𝑈𝑘 𝑄 ) 𝑞𝑟 det 𝑄𝑈𝑘 𝑄 = 1−𝛼 𝑟 [︀ (︀ )︀ ]︀ = ln (1 − 𝛼)𝑛 1 + 𝛼(̂︀ 𝑔𝑟 (𝑢𝑘 ) − 1) det 𝑄𝑈𝑘 𝑄| = (︀ )︀ = 𝑛 ln(1 − 𝛼) + ln 1 + 𝛼(̂︀ 𝑔𝑟 (𝑢𝑘 ) − 1) + ln det 𝑄𝑈𝑘 𝑄| . [︂ 𝑛+1 (︂ 1+ Возьмём производную по 𝛼: d̂︀ 𝑔(𝑢𝑘+1 ) 𝑛 𝑔̂︀𝑟 (𝑢𝑘 ) − 1 =− + . d𝛼 1 − 𝛼 1 + 𝛼(̂︀ 𝑔𝑟 (𝑢𝑘 ) − 1) Приравнивая нулю и решая относительно 𝛼, получаем единственное решение 𝛼= 𝑔̂︀𝑟 (𝑢𝑘 ) − (𝑛 + 1) . (𝑛 + 1)(̂︀ 𝑔𝑟 (𝑢𝑘 ) − 1) Так как 𝑔̂︀𝑟 (𝑢𝑘 ) = (1 + 𝜀)(𝑛 + 1) > 𝑛 + 1 > 1, то 𝛼 > 0. Кроме того, 𝛼 < 1 ⇔ 𝑔̂︀𝑟 (𝑢𝑘 ) − 𝑛 − 1 < 𝑛 · 𝑔̂︀𝑟 (𝑢𝑘 ) − 𝑛 + 𝑔̂︀𝑟 (𝑢𝑘 ) − 1 ⇔ 𝑔̂︀𝑟 (𝑢𝑘 ) > 0. Значит, 𝛼 действительно лежит в (0, 1). Так как 𝑔(𝑢 ̂︀ 𝑘+1 ) выпукла по 𝛼, 𝑢𝑘+1 — искомая точка максимума на отрезке (0, 1). В качестве начального приближения 𝑢0 можно взять вектор (1, . . . , 1)/𝑚. Так как 𝑢0 > 0, то 𝑄𝑈0 𝑄| ≻ 0 и целевая функция 𝑔(𝑢 ̂︀ 0 ) корректно определена. Кроме того, на каждом шаге 𝛼𝑘 < 1, значит если 𝑢𝑘 > 0, то и 𝑢𝑘+1 > 0. Таким образом, алгоритм корректный. Опишем ещё раз последовательность действий на 𝑘-м шаге алгоритма: ∙ Вычислить ∇̂︀ 𝑔(𝑢𝑘 ) ∙ Найти максимальную компоненту градиента 𝑔̂︀𝑟 (𝑢𝑘 )

26 ∙ Вычислить 𝜀 = 𝑔̂︀𝑟 (𝑢𝑘 ) − (𝑛 + 1) 𝑛+1 ∙ Если 𝜀 меньше требуемой погрешности 𝜀0 , закончить вычисления. В противном 𝑔̂︀𝑟 (𝑢𝑘 ) − (𝑛 + 1) случае, вычислить 𝛼 := (𝑛 + 1)(̂︀ 𝑔𝑟 (𝑢𝑘 ) − 1) ∙ Положить 𝑢𝑘+1 :=(1 − 𝛼)𝑢𝑘 + 𝛼𝑒𝑟 и перейти на следующий шаг Кроме того, заметим, что на каждом шаге матрица 𝑄𝑈 𝑄| умножается на 1 − 𝛼 и к ней прибавляется одноранговая матрица 𝛼𝑞𝑟 𝑞𝑟| . В таком случае нет необходимости каж дый раз заново вычислять обратную матрицу, а можно применить формулу обновления (см. [5]): | −1 (𝑄𝑈𝑘+1 𝑄 ) (︂ = 1+ 𝜀 𝑛 · (1 + 𝜀) )︂ (𝑄𝑈𝑘 𝑄| )−1 − 𝜀 𝑏𝑘 𝑏|𝑘 , 2 𝑛 · (1 + 𝜀) где 𝑏𝑘 :=(𝑄𝑈𝑘 𝑄| )−1 𝑞𝑟 . Были проведены численные эксперименты при 𝜀0 = 10−5 и различных размерно стях основного пространства 𝑛 и различном количестве двойственных переменных 𝑚 (то есть количестве точек 𝑞𝑖 ). Результаты приведены в таблице 4.1. Таблица 4.1. Численные результаты алгоритма Хачияна 𝑛 𝑚 число итераций точность время (с) 2 104 199993 7 · 10−6 17.2 2 504 199994 7 · 10−6 78.8 5 510 499990 1.4 · 10−5 261.3 Если в алгоритме критерий останова по 𝜀 (по точности выполнения ограничений) заменить на останов по внутренней сходимости, то число итераций и время работы сократится, но относительная точность не изменится. Отметим, что даже в случае с четырьмя точками в R2 , одна из которых лежит внутри треугольника, образованного тремя другими, алгоритм делает порядка 80 тысяч шагов и работает около 4 секунд.

27 Глава 5 Комбинированный алгоритм Как было продемонстрировано в предыдущей главе, применение одного только ал горитма Хачияна не приносит желаемых результатов из-за низкой точности и большого времени работы. Поэтому были предприняты попытки изменить алгоритм так, чтобы он работал быстрее и точнее. В ходе численных экспериментов была обнаружена следующая закономерность — в начале своей работы алгоритм делает достаточно большие шаги 𝛼, быстро получая эллипсоид, близкий к искомому. Затем же алгоритм начинает делать очень короткие шаги, с 𝛼 ≈ 10−5 , постепенно увеличивая эллипсоид и приближаясь к решению. Это связано с тем, что на каждом шаге алгоритм производит изменение в направлении толь ко одной координаты, как бы «подтягивая» эллипсоид к точке, которая дальше всего от него находится. Таким образом, если текущий эллипсоид уже достаточно близок к оптимальному, то алгоритму приходится поочерёдно немного подтягивать его к точкам, лежащим вне. На основе этого наблюдения и статьи [6] возникла идея улучшения алгоритма. Во-первых, в 1948 году Джон доказал, что минимальный эллипсоид определяется не более (𝑛2 + 3𝑛)/2 точками, где 𝑛 — размерность пространства. То есть, в нашем случае (︀ )︀ только 𝑙 = (𝑛 + 1)2 + 3(𝑛 + 1) /2 точек из исходных 𝑚 определяют искомый эллипсоид. Кроме того, из формулы (4.6) и её эквивалента для вектора 𝑢 видно, что те точки 𝑎𝑖 , для которых 𝑔̂︀𝑖 (𝑢) > 𝑛 + 1, лежат на границе или вне эллипсоида, задаваемого 𝑢. Значит, если их количество не превосходит 𝑙, можно предположить, что только эти точки и задают эллипсоид, который мы ищем. Будем дальше работать уже только с ними. Для корректности задания целевой функции на этих точках необходимо потребовать, чтобы они так же содержали аффинно независимое множество. Таким образом, в «активное» множество всегда будет входить не более 𝑙 и не менее 𝑛 + 1 точек. Во-вторых, раз увеличение эллипсоида в направлении только одной точки приво дит к очень маленьким шагам, разумно попытаться увеличивать его сразу в нескольких направлениях. В первую очередь для этого был испытан обычный градиентный метод наискорейшего подъёма. В процедуре линейного поиска использовано условие Армихо. В качестве направления подъёма выбирается направление градиента, спроецированное

28 на стандартный симплекс по алгоритму, описанному в [7]. Кроме того, линейный поиск останавливается при достижении какой-либо компонентой вектора нуля. Таким образом обеспечивается выполнение ограничений двойственной задачи на каждом шаге гради ентного метода. Критерий останова используется тот же, что и в алгоритме Хачияна. На рисунке 5.1 наглядно представлена работа такого комбинированного подхода, а подробные результаты его применения приведены в таблице 5.1. Здесь столбцы «алго ритм Хачияна» и «градиентный метод» содержат количество шагов соответствующих алгоритмов, в столбце «активные точки» содержится число точек в активном множе стве после остановки алгоритма Хачияна, а столбец «точность» имеет тот же смысл, что и в предыдущей таблице. Для вычислений использовалась погрешность 𝜀0 = 10−5 . Таблица 5.1. Численные результаты комбинированного алгоритма 𝑛 𝑚 алгоритм Хачияна активные точки градиентный метод точность время (с) 2 104 18 7 40 2 · 10−9 0.05 2 504 123 9 87 1.5 · 10−8 0.18 5 510 92 26 75 1.5 · 10−7 0.2 10 1020 80 74 217 3.6 · 10−6 0.9 2.5 2.5 2 2 1.5 1.5 1 1 0.5 0.5 0 0 -0.5 -0.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -1 -0.5 0 0.5 1 1.5 2 2.5 3 (a) Тестовое множество и известный ми (b) Активное множество и построенный нимальный эллипсоид (зелёный) алгоритмом Хачияна эллипсоид (синий) Рис. 5.1. Иллюстрация работы комбинированного алгоритма при 𝑛 = 2, 𝑚 = 104. К сожалению, в некоторых случаях алгоритм зацикливается или выдает неверный ответ. Такое происходит, например, при 𝑛 = 30, 𝑚 = 560 — в этом случае градиентный метод начиная с 1890-го шага начинает делать циклические шаги, не имея возможности никак остановиться. При этом достигается относительная точность только 1 · 10−4 . Это

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

30 Глава 6 Сводное сравнение методов Приведём ещё раз таблицы сравнения трёх методов: алгоритма Шора, алгоритма Хачияна с остановом по внутренней сходимости и комбинированного метода. Таблица 6.1. Алгоритм Шора 𝑛 𝑚 число итераций точность время (с) 2 104 4121 4 · 10−4 0.2 2 504 4118 5 · 10−4 0.2 5 510 16729 8 · 10−4 1.1 10 1020 3533 3 · 10−3 0.6 Таблица 6.2. Алгоритма Хачияна 𝑛 𝑚 число итераций точность время (с) 2 104 86594 1.8 · 10−5 6.2 2 504 86595 1.7 · 10−6 36.5 5 510 94854 8 · 10−5 56.2 10 1020 97451 2.8 · 10−3 343.5 Таблица 6.3. Комбинированный алгоритм 𝑛 𝑚 число итераций точность время (с) 2 104 40 2 · 10−9 0.05 2 504 87 1.5 · 10−8 0.18 5 510 75 1.5 · 10−7 0.2 10 1020 217 3.6 · 10−6 0.9 Также все три алгоритма были запущены на множестве, минимальным эллипсои дом которого является круг, и их результаты сравнены с результатом программы для решения задачи Сильвестра из первой главы. Данные приведены в таблицах 6.4 и 6.5.

31 Таблица 6.4. Сравнение алгоритмов — точность 𝑛 𝑚 задача Сильвестра алгоритм Шора алгоритм Хачияна комбинированный 2 104 1.4 · 10−11 4 · 10−4 1.7 · 10−5 7.5 · 10−9 2 504 1.8 · 10−4 3.7 · 10−4 1.7 · 10−5 1.7 · 10−9 5 510 6.6 · 10−5 6.6 · 10−4 7.9 · 10−5 2.8 · 10−7 10 1020 1.8 · 10−5 1.7 · 10−3 2.8 · 10−4 9 · 10−7 Таблица 6.5. Сравнение алгоритмов — время работы (в секундах) 𝑛 𝑚 задача Сильвестра алгоритм Шора алгоритм Хачияна комбинированный 2 104 0.2 0.3 6.5 0.3 2 504 1.4 0.5 51 0.3 5 510 1.4 1.4 100 0.2 10 1020 8.2 0.7 398 1.4

32 Заключение В рамках данной работы проведён анализ решения задачи о минимальном эллип соиде с применением теории двойственности по Лагранжу. Рассмотрен также частный случай — задача Сильвестра, для которой получено решение с помощью двойственно сти в квадратичном программировании. Одним из результатов работы является разработка комбинированного алгоритма, использующего идею повышения размерности из алгоритма Шора, идею двойственно сти по Лагранжу из алгоритма Хачияна и градиентный метод без точного линейного поиска. Все алгоритмы реализованы в MATLAB и проверены на множестве тестовых при меров. Приведены сводные таблицы точности и времени работы всех методов, выполне но сравнение их производительности на задаче Сильвестра, для которой имеется точное решение с помощью квадратичного программирования. Предварительные результаты были опубликованы в [8, 9, 10, 11].

33 Список литературы 1. Гавурин М. К., Малозёмов В. Н. Экстремальные задачи с линейными ограничения ми. Л.: Изд-во ЛГУ, 1984. С. 176. 2. Шор Н. З., Стеценко С. И. Алгоритм последовательного сжатия пространства для построения описанного эллипсоида минимального объёма // Исследование методов решения экстремальных задач, ИК АН УССР, Киев. 1990. С. 25–29. 3. Brookes M. The Matrix Reference Manual. URL: http://www.ee.imperial.ac.uk/ hp/staff/dmb/matrix/intro.html. 4. Boyd S., Vandenberghe L. Convex Optimization. New York, NY, USA: Cambridge University Press, 2004. 5. Khachiyan L. G. Rounding of Polytopes in the Real Number Model of Computation // Mathematics of Operations Research. 1996. Vol. 21, no. 2. P. 307–320. 6. Sun P., Freund R. M. Computation of Minimum-Volume Covering Ellipsoids // Opera tions Research. 2004. Vol. 52, no. 5. P. 690–706. 7. Малозёмов В. Н. Проектирование точки на подпространство и на стандартный симплекс. Семинар «DHA & CAGD». Избранные доклады. 28.02.2013. URL: http: //dha.spb.ru/reps13.shtml#0228. 8. Кольцов М. А. Решение задачи Сильвестра в MATLAB. Семинар «CNSA & NDO». Избранные доклады. 26.02.2015. URL: http://www.apmath.spbu.ru/cnsa/reps15. shtml#0226. 9. Кольцов М. А. Построение минимального эллипсоида. Семинар «CNSA & NDO». Избранные доклады. 14.05.2015. URL: http://www.apmath.spbu.ru/cnsa/reps15. shtml#0514. 10. Кольцов М. А. Построение минимального эллипсоида: алгоритм Хачияна. Семинар «CNSA & NDO». Избранные доклады. 21.04.2016. URL: http://www.apmath.spbu. ru/cnsa/reps16.shtml#0421. 11. Кольцов М. А. Построение минимального эллипсоида // Процессы управления и устойчивость. Т. 3 (19). 2016. (принята в печать).

Рецензии:

Авторизуйтесь, чтобы добавить рецензию

- у работы пока нет рецензий -

Отзывы:

Авторизуйтесь, чтобы оставить отзыв