Линейная алгебра
Матричные уравнения: Сильвестр, Ляпунов, и стабильность нейросетей
Квадрокоптер DJI удерживает горизонт при порывах ветра благодаря LQR-регулятору: каждые 10 мс внутри контроллера решается матричное уравнение Риккати. Фильтр Калмана в GPS, стабильность LSTM при обучении - это уравнения, где неизвестное само является матрицей.
- Управление: LQR-регулятор дрона решает уравнение Риккати AX + XAᵀ - XBB тXᵀ + Q = 0
- Фильтр Калмана: матрица ковариации ошибки - уравнение Ляпунова на каждом шаге
- RNN стабильность: уравнение Ляпунова AP + PAᵀ = -Q определяет затухание/взрыв градиентов
- Domain adaptation в ML: subspace alignment через уравнение Сильвестра
- Системное управление: размещение полюсов замкнутой системы - уравнение Сильвестра
**Квадрокоптер DJI держит горизонт при порывах ветра - потому что внутри микроконтроллера за 10 миллисекунд решается матричное уравнение.** Фильтр Калмана в GPS-навигаторе, LQR-регулятор промышленного манипулятора, анализ устойчивости рекуррентной нейросети - всё это уравнения вида AX + XB = C, где неизвестное само является матрицей. Обычное уравнение Ax = b - частный случай, где B = 0 и X - вектор.
**О чём урок**: три семейства матричных уравнений в порядке возрастания сложности. AX = B - просто несколько систем одновременно. Уравнение Сильвестра AX - XB = C - ключ к domain adaptation в ML. Уравнение Ляпунова - критерий стабильности динамических систем и RNN. Каждое имеет scipy-решатель в одну строку.
AX = B: несколько систем одновременно
Уравнение AX = B, где A квадратная и X, B - матрицы, читается как набор отдельных систем: каждый столбец B порождает свою правую часть, каждый столбец X - своё решение. Но решать их совместно эффективнее: LU-факторизация A считается один раз.
**Когда A прямоугольная**: если A имеет размер m×n (m ≠ n), то вместо solve используют `np.linalg.lstsq(A, B)` - решение в смысле наименьших квадратов. Именно это происходит при нормальном уравнении линейной регрессии: X = (AᵀA)⁻¹AᵀB, что эквивалентно lstsq.
Уравнение Сильвестра: AX - XB = C
**Уравнение Сильвестра** AX - XB = C - первый пример, где X не вынести за скобку обычным образом: X стоит справа от A и слева от B одновременно. Единственное решение существует тогда и только тогда, когда A и B не имеют общих собственных значений.
Сильвестр в ML: domain adaptation
**Domain adaptation** - задача переноса модели между доменами (обучение на ImageNet, инференс на медицинских снимках). Один из подходов - **subspace alignment**: найти матрицу X, которая выравнивает подпространство признаков исходного домена A с целевым B. Это буквально уравнение Сильвестра.
Уравнение Ляпунова: стабильность систем
Система дифференциальных уравнений dx/dt = Ax описывает динамику - от электрической цепи до скрытого состояния RNN. Вопрос стабильности: затухает ли x(t) к нулю при t → ∞? Ответ - через **уравнение Ляпунова**.
Ляпунов в ML: стабильность RNN
Рекуррентная нейросеть: h[t+1] = tanh(W·h[t] + U·x[t]). Если пренебречь tanh (линеаризация), получаем h[t+1] = W·h[t]. Система устойчива тогда и только тогда, когда все |λᵢ(W)| < 1. **Взрывающийся градиент** в RNN - это нарушение устойчивости: |λ_max(W)| > 1. Уравнение Ляпунова позволяет проверить это аналитически.
**Gradient clipping в PyTorch** (`torch.nn.utils.clip_grad_norm_`) - эвристическая защита от взрывающегося градиента. Уравнение Ляпунова даёт теоретически обоснованный ответ: если спектральный радиус W < 1, градиенты не взрываются. Именно поэтому в LSTM добавили gates - они динамически управляют спектральным радиусом.
Уравнение Риккати и оптимальное управление
**Алгебраическое уравнение Риккати (ARE)** - нелинейное по X (квадратичный член XBR⁻¹BᵀX), что делает его значительно сложнее Ляпунова. Решение ARE даёт оптимальный линейный регулятор (LQR), минимизирующий квадратичный функционал качества.
**Фильтр Калмана** решает двойственную задачу: вместо управления - оценка состояния. Уравнение Риккати там тоже присутствует, но для матрицы ковариации ошибки. Навигация GPS, инерциальные системы в смартфонах, трекинг объектов в self-driving - всё это фильтр Калмана, то есть уравнение Риккати в реальном времени.
Сравнение трёх уравнений
**Ловушка при численном решении Ляпунова**: `solve_continuous_lyapunov(A, Q)` решает AX + XAᵀ = Q (без минуса). Чтобы решить AX + XAᵀ = -Q, нужно передать `-Q`. Это частая ошибка - неправильный знак даёт решение, не удовлетворяющее исходному уравнению, но без ошибки в scipy.
Практика: уравнение Сильвестра
AX = B: несколько систем одновременно
Уравнение $AX = B$, где $X$ и $B$ - матрицы, это просто $k$ независимых систем $Ax_i = b_i$ (по столбцам). Решение: $X = A^{-1}B$ или через LU-разложение с несколькими правыми частями. Если $A$ не квадратная или вырожденная, используют псевдообратную: $X = A^+ B$.
**Уравнение Сильвестра**: $AX - XB = C$. Левая часть - линейное отображение на пространстве матриц. Единственное решение существует тогда и только тогда, когда $A$ и $B$ не имеют общих собственных значений. В scipy: `scipy.linalg.solve_sylvester(A, B, C)`.
Калибровка нескольких камер
AX = B на практике
Стерео-система из 4 камер: для каждой нужно найти матрицу трансформации $X_i \in \mathbb{R}^{4 \times 4}$. Всё это - $AX = B$ с матричными правыми частями. LU-факторизация $A$ один раз, затем 4 back-substitution - экономия 3x по сравнению с 4 независимыми решениями.
При каком условии уравнение Сильвестра $AX - XB = C$ имеет единственное решение?
Уравнение Сильвестра $AX - XB = C$ задаёт линейное отображение $\mathcal{L}(X) = AX - XB$. Оно обратимо (единственное решение для любого $C$) тогда и только тогда, когда $\lambda_i(A) \neq \mu_j(B)$ для всех $i, j$.
Уравнение Ляпунова и стабильность динамических систем
Уравнение Ляпунова: $AP + PA^\top = -Q$, где $Q$ - симметричная положительно определённая матрица. Если оно имеет решение $P \succ 0$, то динамическая система $\dot{x} = Ax$ асимптотически устойчива (все собственные значения $A$ имеют отрицательную вещественную часть).
**Связь с RNN**: рекуррентная нейросеть $h_t = Wh_{t-1} + \ldots$ устойчива, если матрица весов $W$ удовлетворяет уравнению Ляпунова с $P \succ 0$. Это даёт альтернативный взгляд на проблему затухающих/взрывных градиентов: спектральный радиус $\rho(W) < 1$ - достаточное условие.
Фильтр Калмана и уравнение Риккати
Управление и оценивание
Фильтр Калмана: оптимальная оценка состояния системы $\dot{x} = Ax + Bu + w$ по наблюдениям $y = Cx + v$. Матрица ковариации ошибки $P$ удовлетворяет уравнению Риккати (обобщение Ляпунова). В GPS-навигаторе это решается за микросекунды - специализированные алгоритмы (DARE: Discrete Algebraic Riccati Equation) из scipy.
Уравнение Ляпунова $AP + PA^\top = -Q$ ($Q \succ 0$) имеет решение $P \succ 0$. Что это говорит о системе $\dot{x} = Ax$?
Существование $P \succ 0$ при $Q \succ 0$ эквивалентно тому, что $V(x) = x^\top P x$ - функция Ляпунова: $\dot{V} = x^\top(AP + PA^\top)x = -x^\top Q x < 0$. Все собственные значения $A$ имеют $\text{Re}(\lambda_i) < 0$.
AX = B, A - (n×n), X и B - (n×k) Экивалентно k отдельным системам: A·x₁ = b₁ A·x₂ = b₂ ... A·xₖ = bₖ Наивно: k × O(n³) операций Через LU: O(n³) + k × O(n²) - Факторизуем PA = LU один раз - k раз решаем L(Ux) = Pb за O(n²) каждый При k=100, n=1000: 10¹¹ vs 10⁹+10⁸ = 100x ускорение
AX - XB = C Если бы X были числа (1×1 матрицы): ax - xb = c -> x(a-b) = c -> x = c/(a-b) Не работает когда a = b (делитель 0) Для матриц: нельзя вынести X за скобку, потому что: AX ≠ XA (матрицы не коммутируют) AX - XB (X умножается на A слева и на B справа) Вместо этого - оператор Кронекера: vec(X) - вектор из столбцов X, длина mn (I ⊗ A - Bᵀ ⊗ I) · vec(X) = vec(C) Это линейная система mn×mn - правильная постановка, но алгоритм Бартелса-Стюарта (Шур + обратная подстановка) решает её за O(m³ + n³) вместо O(m³n³)
УРАВНЕНИЕ: AX + XAᵀ = -Q Где Q - произвольная положительно определённая матрица (обычно I). ТЕОРЕМА ЛЯПУНОВА: Если Q > 0 (положительно определена), то: Существует X > 0 <-> A устойчива (Re(λᵢ) < 0 для всех λ) ФИЗИЧЕСКИЙ СМЫСЛ: V(x) = xᵀX x - функция Ляпунова (аналог "энергии") dV/dt = xᵀ(AᵀX + XA)x = -xᵀQx < 0 Энергия убывает вдоль любой траектории -> система устойчива Для дискретной системы x[t+1] = Ax[t]: AXAᵀ - X = -Q (дискретное уравнение Ляпунова) X > 0 <-> все |λᵢ| < 1
УРАВНЕНИЕ РИККАТИ (непрерывное): AᵀX + XA - XBR⁻¹BᵀX + Q = 0 Где: A - матрица динамики системы (dx/dt = Ax + Bu) B - матрица управления Q >= 0 - штраф за отклонение состояния R > 0 - штраф за управление ОПТИМАЛЬНЫЙ РЕГУЛЯТОР: u*(t) = -R⁻¹BᵀX · x(t) (линейная обратная связь) МИНИМИЗИРУЕТ: J = ∫₀^∞ (xᵀQx + uᵀRu) dt Отличие от Ляпунова: Ляпунов: линейный по X -> одно решение Риккати: квадратичный по X -> несколько решений, нас интересует стабилизирующее X > 0
| Уравнение | Вид | Сложность | Применение |
|---|---|---|---|
| AX = B | Линейное | O(n³), np.linalg.solve | Несколько правых частей, meta-learning |
| Сильвестр: AX - XB = C | Линейное по X | O(m³+n³), solve_sylvester | Domain adaptation, декаплирование систем |
| Ляпунов: AX + XAᵀ = -Q | Линейное по X (Сильвестр с B=A) | O(n³), solve_continuous_lyapunov | Стабильность RNN, H2-норма, Kalman |
| Риккати: AᵀX+XA-XBR⁻¹BᵀX+Q=0 | Квадратичное по X | O(n³) через Гамильтониан | LQR, H∞ управление, оптимальный Kalman |
- **LQR-регулятор (дроны, роботы)**: ARE -> оптимальный линейный регулятор
- **Фильтр Калмана (навигация)**: Riccati -> оценка состояния + ковариации ошибки
- **Анализ стабильности RNN**: Lyapunov -> проверка спектрального радиуса, H2-норма
- **Domain adaptation (ML)**: Sylvester -> выравнивание подпространств признаков
- **H∞ управление (авиация, энергетика)**: Модифицированная ARE -> робастный регулятор
Упражнения
- При каком условии уравнение Сильвестра AX - XB = C имеет единственное решение? Что происходит, когда это условие нарушается? — Условие: λ(A) ∩ λ(B) = ∅ - спектры не пересекаются; При общем λ: оператор X → AX - XB вырожден (имеет ненулевое ядро); Если решение существует - их бесконечно много (свободный параметр); Это аналог det(A) = 0 для обычной системы Ax = b
- Как уравнение Ляпунова связано с проблемой взрывающегося градиента в RNN? — Дискретная RNN h[t+1] = Wh[t] устойчива тогда и только тогда, когда |λ_i(W)| < 1; Дискретное уравнение Ляпунова WXWᵀ - X = -Q имеет X > 0 ↔ |λ_max| < 1; Взрывающийся градиент = |λ_max| > 1 = нет решения X > 0; Gradient clipping - эвристика; правильное решение - инициализировать W с |λ_max| < 1 или использовать LSTM/GRU
- Почему уравнение Риккати сложнее уравнения Ляпунова, хотя оба имеют scipy-решатель в одну строку? — Ляпунов: линейное по X -> одно решение, вычисляется напрямую через преобразование Кронекера; Риккати: содержит XBR⁻¹BᵀX -> нелинейное, квадратичное -> несколько решений; Для ARE нужно стабилизирующее решение X > 0, которое выбирается через структуру Гамильтониановой матрицы; scipy.linalg.solve_continuous_are внутри вычисляет инвариантное подпространство 2n×2n матрицы
Что унести из урока
- **AX = B** - несколько систем с общей матрицей A: np.linalg.solve принимает матричную правую часть, LU считается один раз
- **Сильвестр AX - XB = C**: единственное решение при непересекающихся спектрах; scipy.linalg.solve_sylvester; применение - domain adaptation
- **Ляпунов AX + XAᵀ = -Q**: X > 0 эквивалентно устойчивости A (Re(λ) < 0); ключ к анализу RNN и расчёту H2-нормы
- **Взрывающийся градиент в RNN** = |λ_max(W)| > 1 = нарушение условия Ляпунова; LSTM решает это через gate механизм
- **Риккати**: квадратичное по X, несколько решений, нужно стабилизирующее; даёт оптимальный LQR регулятор и фильтр Калмана
- **Фильтр Калмана = Riccati в реальном времени**: GPS, INS, трекинг объектов - всё решает ARE на каждом шаге
- **Знак важен**: solve_continuous_lyapunov(A, Q) решает AX + XAᵀ = Q, для -Q нужно передавать -Q явно
Связи
Матричные уравнения на пересечении управления, оптимизации и теории устойчивости
- Разреженная линейная алгебра — Для больших систем управления уравнения Ляпунова и Риккати решаются итерационными sparse методами
- Собственные векторы и диагонализация — Устойчивость определяется собственными значениями A; алгоритм Бартелса-Стюарта использует разложение Шура
- Матричные функции — Матричная экспонента eᴬᵗ - решение dx/dt = Ax; стабильность системы связана с уравнением Ляпунова через интегральное представление