LU-разложение: метод Гаусса в разобранном виде

**scipy.linalg.lu_factor решает задачу, от которой метод Гаусса падает в обморок.** Рассмотрим задачу: 1000 разных правых частей для одной матрицы 500x500 - например, разные граничные условия в симуляции конечных элементов. Метод Гаусса: O(n^3) * 1000 операций. LU: O(n^3) один раз + O(n^2) * 1000 раз. Разница - в тысячи раз. LU-разложение - это метод Гаусса, записанный в разобранном виде.

**О чём урок**: A = LU - это не магия факторизации, а просто аккуратная запись шагов исключения Гаусса. L хранит множители (что мы вычитали), U - результат (что получилось). Но эта запись открывает главное преимущество: матрицу раскладывают один раз, а потом подставляют сколько угодно правых частей за O(n^2) каждый раз.

Структура: что такое L и U

A = LU, где L - нижнетреугольная (Lower) с единицами на диагонали, U - верхнетреугольная (Upper). Размер обеих - n x n.

Откуда берутся L и U: шаг за шагом

Метод Гаусса исключает переменные, вычитая строки с множителями. LU просто **запоминает эти множители** в матрицу L.

Решение системы: два треугольных шага

Как только A = LU, система Ax = b превращается в два простых шага. Треугольные системы решаются подстановкой за O(n^2) - без всякого исключения.

**numpy.linalg.solve** внутри делает то же самое: вызывает LAPACK функцию dgesv, которая сначала выполняет PLU разложение, потом два шага подстановки. Если нужно решать много систем с одной матрицей - lu_factor/lu_solve явно эффективнее.

PLU: перестановки строк для устойчивости

Простое LU ломается, если на диагонали встречается ноль (деление на ноль). И даже при маленьком диагональном элементе - численная нестабильность из-за ошибок округления. Решение: **частичный выбор главного элемента** (partial pivoting) - перестановки строк так, чтобы наибольший элемент оказался на диагонали.

Cholesky: LU для симметричных положительно определённых матриц

Если матрица симметричная и положительно определённая (SPD) - существует более эффективное разложение Холецкого. Используется в Gaussian Processes, Kalman Filter, MCMC с матрицами ковариации.

**Sparse LU в конечно-элементном анализе**: при симуляции физических процессов (FEM/CFD) матрица системы разреженная (sparse) - большинство элементов нули. Специализированные sparse LU solvers (UMFPACK, SuperLU, PARDISO) используют структуру разреженности чтобы хранить и обрабатывать только ненулевые элементы. Матрицы 10^6 x 10^6 решаются на обычном сервере.

Сложность и сравнение методов

LU-разложение в реальных системах

От numpy до физических симуляций Components: numpy.linalg.solve / scipy.linalg.lu_factor, Gaussian Process Regression (sklearn, GPyTorch), FEM / CFD симуляции, Kalman Filter, Нормальное уравнение в регрессии

Практика: LU-разложение

- Почему для решения 1000 систем Ax = b с одной матрицей A лучше использовать lu_factor/lu_solve вместо 1000 вызовов np.linalg.solve? - Когда использовать Cholesky вместо LU? Какие условия нужны? - Что такое partial pivoting и зачем он нужен? Что происходит без него?

L = | 1 0 0 | U = | u11 u12 u13 | | l21 1 0 | | 0 u22 u23 | | l31 l32 1 | | 0 0 u33 | L: единицы на диагонали, ниже - множители исключения U: верхнетреугольная - результат метода Гаусса Свойства: det(L) = 1 (произведение диагональных = 1*1*...*1) det(A) = det(L) * det(U) = 1 * u11*u22*...*unn = произведение диагональных элементов U

A = | 2 1 1 | | 4 3 3 | | 8 7 9 | Шаг 1: обнуляем первый столбец l21 = 4/2 = 2: R2 = R2 - 2*R1 l31 = 8/2 = 4: R3 = R3 - 4*R1 -> | 2 1 1 | | 0 1 1 | | 0 3 5 | Шаг 2: обнуляем второй столбец l32 = 3/1 = 3: R3 = R3 - 3*R2 -> U = | 2 1 1 | | 0 1 1 | | 0 0 2 | Множители идут в L: L = | 1 0 0 | | 2 1 0 | | 4 3 1 | Проверка: L * U = | 1*2+0+0 1*1+0+0 1*1+0+0 | | 2 1 1 | | 2*2+1*0 2*1+1*1 2*1+1*1 | = | 4 3 3 | = A ✓ | 4*2+3*0 4*1+3*1 4*1+3*1+1*2 | | 8 7 9 |

Ax = b LUx = b Обозначим y = Ux, тогда Ly = b ШАГ 1 - прямая подстановка (forward substitution): решаем Ly = b L нижнетреугольная -> y находим сверху вниз за O(n^2) y1 = b1 y2 = b2 - l21*y1 y3 = b3 - l31*y1 - l32*y2 ШАГ 2 - обратная подстановка (back substitution): решаем Ux = y U верхнетреугольная -> x находим снизу вверх за O(n^2) x3 = y3 / u33 x2 = (y2 - u23*x3) / u22 x1 = (y1 - u12*x2 - u13*x3) / u11 Итого: LU один раз O(n^3) + каждая новая b за O(n^2)

ПРОБЛЕМА без pivoting: A = | 0 1 | Шаг 1: нужно разделить на a11 = 0 | 2 3 | Деление на ноль! Алгоритм падает. С partial pivoting: Найти строку с max |элемент| в первом столбце: строка 2 (|2| > |0|) Поменять строки: P*A = | 2 3 | | 0 1 | Теперь можно продолжать нормально. P*A = L*U, где P - матрица перестановок Частичный pivoting гарантирует: |l_ij| <= 1 для всех i,j Это означает числовую устойчивость. det(A) = (-1)^s * u11 * u22 * ... * unn где s - число перестановок строк

Для SPD матрицы A: A = L * L^T L - нижнетреугольная с положительными диагональными элементами Преимущества: Операций: n^3/3 (вдвое меньше чем n^3/1.5 для LU) Память: хранить только L (вдвое меньше L + U) Устойчивость: всегда стабильно без pivoting (если A SPD) Когда встречается в ML: - Gaussian Process: K = L*L^T, затем L^{-1} y для предсказания - Multivariate Normal: sampling через L*z, z ~ N(0, I) - Kalman Filter: ковариационные матрицы - всегда SPD - Adam optimizer: second-moment estimate (диагональная, но принцип тот же)