Стивен Бойд: «В хорошей численной работе матрицы никто не обращает»

Stanford Online 7,2 тыс. 1 ч 18 мин 10 мин 23.03.2024
Главное

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

🧮 Демистификация линейной алгебры: как работают прямые методы решения уравнений 0:06

В рамках курса выпуклой оптимизации студентам предстоит за несколько недель пройти путь снизу вверх — от базовых матричных операций до верхнего уровня приложений. Итогом этого процесса станет самостоятельное написание собственного солвера для задач линейного программирования (LP). Стивен Бойд с иронией отмечает, что главная цель этого задания — демистифицировать процесс: если в будущем кто-то заявит, что создание LP-солвера является невероятно сложной задачей, студент сможет аргументированно ответить: «Au contraire, я это сделал».

Самым нижним уровнем численной оптимизации является линейная алгебра, знание идей которой критически необходимо любому инженеру. Для решения систем уравнений вида $Ax = b$ существуют итерационные методы, применяемые для гигантских систем с миллиардами переменных (например, при моделировании дифференциальных уравнений в частных производных). Однако в фокусе лекции находятся так называемые прямые методы, суть которых заключается в разложении (факторизации) матрицы $A$ на множители, каждый из которых легко обратить. При этом Стивен Бойд подчеркивает фундаментальное правило численного анализа: в хорошей инженерной практике реальное формирование обратной матрицы $A^{-1}$ происходит чрезвычайно редко. На самом деле математическая запись вида $A_1^{-1}b$ означает просто запуск функции решения системы уравнений с коэффициентами solve(A1, b).

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

Классическим примером выступает исключение Гаусса, или $LU$-разложение, при котором несингулярная матрица раскладывается на матрицу перестановок $P$, нижнюю треугольную матрицу $L$ и верхнюю треугольную матрицу $U$. Вычислительная стоимость самой факторизации составляет $O(n^3)$ флопс, в то время как последующее решение через прямую и обратную подстановку требует всего $O(n^2)$ флопс. На современном ноутбуке система уравнений размером $1000 \times 1000$ (содержащая миллион коэффициентов) решается всего за 100 миллисекунд. По словам профессора, удивительный для многих факт заключается в том, что решение 10 или даже 100 таких систем с разными правыми частями займет практически те же самые 100 миллисекунд, поскольку дорогостоящая факторизация $O(n^3)$ выполняется всего один раз. Под капотом современных программных пакетов, таких как библиотека CVXPY, эти свойства используются непрерывно.

🕸️ Разреженные и ленточные матрицы: обход NP-трудных задач ради миллисекундной скорости 8:49

Разреженные матрицы активно развиваются в вычислительной математике с 1960-х годов. Для них существует целая «параллельная вселенная» алгоритмов — от факторизации до сингулярного разложения ($SVD$) и поиска собственных значений, адаптированная под эффективное хранение только ненулевых элементов.

При факторизации разреженной матрицы $A$ обычно используются четыре фактора:

Математически для исключения Гаусса достаточно одной перестановки, но вторая вводится специально для борьбы с так называемым «заполнением» (fill-in). Заполнение происходит, когда в процессе разложения изначально нулевые элементы исходной матрицы превращаются в ненулевые, превращая разреженные матрицы $L$ и $U$ в плотные. Профессор Бойд утверждает, что работа с разреженной матрицей размером миллион на миллион строк на обычном ноутбуке не представляет проблемы, если в каждой строке содержится разумное количество (например, от 20 до 100) ненулевых элементов. Однако при неудачном выборе перестановок $P_1$ и $P_2$ матрица $L$ станет плотной, и операционная система просто откажет в выделении памяти для хранения $10^{12}$ элементов.

Задача поиска идеальных перестановок, минимизирующих заполнение, доказанно является NP-трудной. Тем не менее, лектор призывает не беспокоиться об этом, поскольку разработанные эффективные эвристики обеспечивают превосходный результат на практике.

Особым случаем разреженных структур являются ленточные матрицы шириной ленты $k$. Для них эвристические перестановки не требуются, а матрицы $L$ и $U$ гарантированно сохраняют исходную ширину ленты. Время факторизации ленточной матрицы падает до $O(nk^2)$, а обратная подстановка занимает всего $O(nk)$. Это позволяет решать систему из миллиона уравнений с шириной ленты 10 буквально за миллисекунды прямо на мобильном телефоне. По мнению Стивена Бойда, вся современная индустрия автоматического управления и обработки сигналов держится именно на этой математической особенности. Задача управления с 20 000 временных шагов и десятками актуаторов, кажущаяся стороннему наблюдателю вычислительным кошмаром, благодаря ленточной структуре щелкается современными гаджетами мгновенно.

📐 Факторизация Холецкого и симметрия: магия стреловидных структур 17:38

Для симметричных положительно определенных матриц применяется факторизация Холецкого, позволяющая представить матрицу в виде произведения $A = LL^T$. Если потребовать, чтобы диагональные элементы матрицы $L$ были строго положительными, такое разложение становится абсолютно уникальным. Стоимость факторизации Холецкого составляет $1/3 n^3$ флопс, что в два раза быстрее общего случая $LU$-разложения благодаря исключению избыточных симметричных вычислений. Бойд призывает не переоценивать точные коэффициенты в оценках флопс, поскольку реальная скорость на современных процессорах зависит от оптимизации микрокода, кэширования и многопоточности, за разработку которых инженерам чипов стоит быть глубоко благодарными. В выпуклой оптимизации системы положительно определенных уравнений лежат в самом сердце большинства алгоритмов.

В разреженном варианте Холецкого выбор правильной перестановки строк и столбцов имеет решающее значение. Профессор приводит пример так называемой «матрицы-стрелы» (arrow matrix), у которой ненулевые элементы расположены только на главной диагонали, а также в первой строке и первом столбце. Если запустить алгоритм Холецкого для такой структуры напрямую, матрица $L$ получится абсолютно плотной, заполнив все пустые места. Но если переставить первую и последнюю строки и столбцы местами, структура превратится в «перевернутую стрелу», и факторизация Холецкого пройдет без единого лишнего ненулевого элемента (zero fill-in). Время решения в таком случае сократится до линейного — $O(n)$.

Процесс поиска оптимальной перестановки называется символьной фазой решения (или фазой перестановки). На этом этапе процессор выполняет исключительно графовые вычисления без использования вещественных чисел. Лектор делится практическим инсайдом: в бортовых системах ракет или реактивных двигателей, где алгоритм оптимизации должен циклически отрабатывать на частоте 100 Гц, инженеры могут «потратить целый вечер» на оффлайн-поиск идеальной перестановки для фиксированной структуры графа, чтобы затем проводить вычисления в реальном времени с максимальной скоростью.

Для систем, которые являются симметричными, но не обязательно положительно определенными, применяется факторизация вида $LDL^T$, где матрица $D$ является блочно-диагональной с блоками размера $1 \times 1$ или $2 \times 2$, которые также легко поддаются обращению.

🧩 Дополнение Шура и блочное исключение: когда размер имеет значение 30:51

Концепция дополнения Шура (Schur complement) универсальна и под разными именами встречается практически во всех инженерных и научных дисциплинах:

Если разбить матрицу $A$ на блоки $A_{11}, A_{12}, A_{21}, A_{22}$ и conformally разделить векторы $x$ и $b$ на две соответствующие части, то при условии обратимости блока $A_{11}$ исходную систему можно свести к уравнению относительно блока $x_2$. Матрица этого уравнения $S = A_{22} - A_{21}A_{11}^{-1}A_{12}$ и называется дополнение Шура.

Общий алгоритм блочного исключения включает в себя:

  1. Факторизацию блока $A_{11}$ и вычисление промежуточных матриц через вызовы solve для каждого столбца $A_{12}$.
  2. Формирование матрицы дополнения Шура $S$.
  3. Решение результирующей системы $Sx_2 = \tilde{b}$ через факторизацию и обратную подстановку.
  4. Нахождение оставшейся части вектора $x_1$.

Стивен Бойд честно признает: для обычной плотной матрицы этот громоздкий метод на уровне флоп-каунта не дает абсолютно никакого выигрыша по сравнению со стандартным исключением Гаусса. Однако блочное исключение становится прорывным, когда блок $A_{11}$ устроен просто — например, является диагональной матрицей, чей фактор-кост равен нулю. В этом случае сложность вычисления становится строго линейной относительно размерности $n_1$. Аналогично, если структура представляет собой ленточную матрицу с добавлением нескольких плотных строк и столбцов (блочно-стреловидная структура), алгоритм блочного исключения позволяет решать систему за линейное время, избавляя инженера от необходимости загружать гигантские плотные аналоги в память.

📈 Лемма об обращении матриц и финансовые факторные модели 42:59

Уравнение вида $(A + BC)x = b$ часто называют низкоранговым возмущением или обновлением матрицы $A$. Математический инструмент для его аналитического решения известен как формула Шермана — Моррисона — Вудбери, или лемма об обращении матриц. Профессор Бойд шутит, что его близкий друг по используемому термину может безошибочно определить географию обучения собеседника — учился ли он в Оксфорде, Кембридже, Гарварде или Массачусетском технологическом институте.

Широкое практическое применение эта лемма нашла в количественных финансах (quantitative finance), где ковариационная матрица доходностей активов моделируется как диагональная матрица плюс матрица низкого ранга. В таких факторных моделях (factor models) ковариация тысяч акций описывается через небольшое число общих движущих сил рынка. Диагональная часть $A$ отражает так называемый идиосинкразический риск (idiosyncratic risk) — индивидуальные колебания акций компаний вроде Google, Coca-Cola или Chevron, не объясняемые общим трендом рынка. Вторая часть, имеющая малую ширину (например, 100–200 факторов на 80 000–100 000 активов), задает влияние глобальных факторов. Благодаря лемме об обращении матриц, сложнейшие портфельные расчеты и бэк-тесты инвестиционных фондов выполняются на компьютерах практически мгновенно.

С технической точки зрения для вывода этой формулы применяется метод «разукрупнения» (unelimination) — искусственное введение новых переменных для расширения системы уравнений, например, искусственное раздувание системы с размера $1000 \times 1000$ до $10000 \times 10000$. Профессор предостерегает студентов от любительского стремления всегда уменьшать размер систем: увеличение размерности с одновременным сохранением жесткой разреженной структуры в линейной алгебре работает гораздо эффективнее. Более того, Бойд заявляет, что развитие современных разреженных солверов во многом сделало ненужными классические академические лекции прошлых десятилетий (например, сложные изящные аналитические выводы фильтров Калмана): сегодня инженеру достаточно просто составить гигантскую разреженную систему уравнений и передать её готовому программному пакету.

🏃 Безусловная минимизация: спуск по гладким функциям 55:08

Поднимаясь на следующий этаж иерархии оптимизации, лектор вводит тему безусловной минимизации гладких функций $f$, предполагая их дважды непрерывно дифференцируемыми. Сюда не относятся популярные задачи оптимизации $L_1$-нормы (например, лассо или робастная аппроксимация), поскольку они содержат изломы и не являются дифференцируемыми. В общем виде аналитических решений для таких задач не существует, поэтому применяются итерационные методы, генерирующие последовательность точек, сходящуюся к оптимальному значению по величине функции или по уменьшению модуля градиента. Здесь, в отличие от прямых методов линейной алгебры, критически важен критерий остановки вычислений.

Для математического доказательства сходимости алгоритмов вводятся условия замкнутости подмножеств (надграфика функции) и предположение о существовании минимальной кривизны. Минимальная кривизна означает, что минимальное собственное значение гессиана функции ограничено снизу некоторым строго положительным числом $m > 0$. По мнению Бойда, в реальной инженерной практике создатели алгоритмов чаще всего понятия не имеют, чему равно это число $m$. Редким исключением являются задачи мультиномиальной логистической регрессии с квадратичной регуляризацией (ридж-регрессия), где сам параметр регуляризации напрямую задает гарантированную минимальную кривизну $m$. Знание параметра $m$ позволяет строго ограничить расстояние текущей точки до истинного оптимума через норму градиента:

$$f(x) - p^* \le \frac{1}{2m}|\nabla f(x)|^2_2$$

Данная формула служит строгим математическим обоснованием для остановки итераций, как только градиент становится достаточно малым.

📉 Методы спуска и парадокс линейного поиска: от теории к зигзагам градиента 1:02:19

Общая схема методов спуска заключается в поиске направления шага $\Delta x$ (поискового направления) и подборе положительного скаляра $t$ (размера или длины шага). Скалярное произведение выбранного направления с градиентом должно быть строго отрицательным — только тогда значение функции гарантированно пойдет вниз.

Существует два основных подхода к выбору величины шага $t$:

Стивен Бойд раскрывает главный парадокс линейного поиска: интуитивное обывательское представление о том, что точный поиск на линии всегда эффективнее грубого возвратного алгоритма, абсолютно ошибочно. На практике точный поиск тратит слишком много ресурсов и не дает преимуществ, а зачастую работает даже хуже. Простого эмпирического backtracking-алгоритма с параметром $\alpha = 0.01$ (согласие на малейшее улучшение) более чем достаточно для быстрой работы.

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

Если же линии уровня сильно вытянуты (высокое число обусловленности, большая разница в кривизне по разным осям), антиградиент перестает указывать на центр, и алгоритм начинает метаться из стороны в сторону. Этот неприятный феномен в оптимизации называют эффектом Маратоса или просто «зигзагами» (zig-zagging). Профессор переводит это явление на язык смежных дисциплин: радиотехники назвали бы это необходимостью пропустить поисковый вектор через фильтр низких частот, а физики предложили бы добавить в алгоритм «терм вязкости» или импульс (momentum), чтобы сгладить поперечные колебания и ускорить движение к цели. В чистом виде градиентный спуск демонстрирует лишь слабую линейную сходимость (прямая линия на полулогарифмическом графике), требуя сотен итераций для снижения ошибки даже на простых двумерных неквадратичных функциях.

💬 Цитаты

«В хорошей численной работе формирование обратной матрицы происходит невероятно редко.»

Стивен Бойд 02:51

«Развитие современных разреженных солверов сделало ненужными множество академических лекций прошлых десятилетий.»

Стивен Бойд 52:52
👥 Спикер
🔗 Упомянутые сайты и проекты
📖 Термины
Факторизация матрицы
Разложение исходной матрицы на произведение нескольких матриц специального (например, треугольного) вида для упрощения последующих вычислений.
Дополнение Шура
Матричный блок, возникающий в процессе исключения части переменных из системы линейных уравнений.
Эффект Маратоса (зигзаги)
Феномен в оптимизации, при котором градиентный алгоритм совершает множество осциллирующих шагов из стороны в сторону вместо прямого движения к минимуму.
📊 Цифры
⚖️ Другая сторона
Математика и физика Стивен Бойд Выпуклая оптимизация Факторизация Холецкого Дополнение Шура Градиентный спуск