LU-разложение: решаем много систем дёшево

Урок про LU-разложение — способ «запомнить» работу Гаусса, чтобы решать систему с новой правой частью почти мгновенно.

LU-разложение представляет матрицу как произведение нижнетреугольной $L$ (с единицами на диагонали) и верхнетреугольной $U$: $A = L \cdot U$. Тогда $Ax=b$ решается в два треугольных прохода.

Зачем отделять разложение от правой части

Представьте: матрица $A$ одна (например, описывает конструкцию), а нагрузок $b$ — сотня вариантов. Запускать Гаусс заново для каждой $b$ расточительно: вся дорогая работа ($O(n^3)$) зависит только от $A$, а не от $b$. LU-разложение делает эту работу один раз и сохраняет результат в виде $L$ и $U$. После этого каждая новая правая часть решается двумя дешёвыми треугольными проходами за $O(n^2)$.

Схема решения $LUx = b$: вводим $y = Ux$, тогда $Ly = b$. Сначала прямой подстановкой находим $y$ (из $L$, сверху вниз), затем обратной — $x$ (из $U$, снизу вверх).

def lu_разложение(A):
    n = len(A)
    L = [[0.0] * n for _ in range(n)]
    U = [[0.0] * n for _ in range(n)]
    for i in range(n):
        L[i][i] = 1.0
    for j in range(n):
        for i in range(j + 1):
            U[i][j] = A[i][j] - sum(L[i][k] * U[k][j] for k in range(i))
        for i in range(j + 1, n):
            L[i][j] = (A[i][j] - sum(L[i][k] * U[k][j] for k in range(j))) / U[j][j]
    return L, U

def реши_lu(L, U, b):
    n = len(b)
    y = [0.0] * n            # Ly = b, прямая подстановка
    for i in range(n):
        y[i] = b[i] - sum(L[i][k] * y[k] for k in range(i))
    x = [0.0] * n            # Ux = y, обратная подстановка
    for i in range(n - 1, -1, -1):
        x[i] = (y[i] - sum(U[i][k] * x[k] for k in range(i + 1, n))) / U[i][i]
    return x

A = [[4.0, 3.0], [6.0, 3.0]]
L, U = lu_разложение(A)
print("L =", L)
print("U =", U)
# Решим с ДВУМЯ разными правыми частями, не раскладывая заново
print("x для b=[10,12]:", реши_lu(L, U, [10.0, 12.0]))
print("x для b=[1, 0] :", реши_lu(L, U, [1.0, 0.0]))

Вывод:

L = [[1.0, 0.0], [1.5, 1.0]]
U = [[4.0, 3.0], [0.0, -1.5]]
x для b=[10,12]: [1.0, 2.0]
x для b=[1, 0] : [-0.5, 1.0]

Видно: разложили $A = LU$ один раз, а решили две системы. Третья, десятая, сотая правая часть обойдётся так же дёшево.

Где это особенно ценно

Обращение матрицы. $A^{-1}$ — это решение $AX = I$, то есть $n$ систем с правыми частями — столбцами единичной матрицы. LU делает это естественно.
Определитель. $\det(A) = \det(U)$ = произведение диагонали $U$ (с учётом знака перестановок) — почти бесплатно после разложения.
Итерационное уточнение. Имея LU, можно дёшево улучшать решение, решая систему для невязки.

Как работает под капотом

LU без выбора главного элемента (как выше) может споткнуться о нулевой $U[j][j]$ или потерять точность — поэтому на практике используют LU с частичным выбором, дающее $PA = LU$ с матрицей перестановок $P$. Тогда решают $Ly = Pb$, затем $Ux = y$. Для симметричных положительно определённых матриц есть более дешёвый и устойчивый вариант — разложение Холецкого $A = L \cdot L^{\mathsf{T}}$ (вдвое меньше работы, не нужен pivoting). В библиотеке всё это — одна строка: scipy.linalg.lu_factor + lu_solve.

# В библиотеке (SciPy) разложение и решение — две строки:
from scipy.linalg import lu_factor, lu_solve
lu, piv = lu_factor(A)            # один раз, O(n^3)
x = lu_solve((lu, piv), b)        # для каждой b, O(n^2)

Частые ошибки

Раскладывать заново для каждой $b$. Весь смысл LU — разложить один раз; не вызывайте факторизацию в цикле по правым частям.
Игнорировать выбор главного элемента. «Учебный» LU без pivoting неустойчив; для реальных данных берите $PA=LU$ или Холецкого.
Считать $A^{-1}$, чтобы потом умножать на $b$. Явное обращение дороже и менее точно, чем решение через LU; обратную матрицу почти никогда не нужно строить.

Итоги

LU-разложение: $A = L \cdot U$ — дорогая работа ($O(n^3)$) делается один раз и зависит только от $A$.
Каждую новую правую часть решают двумя треугольными проходами за $O(n^2)$.
Даёт почти бесплатно определитель и обращение; на практике — с выбором главного элемента ($PA=LU$) или Холецкого для SPD.
Главный сценарий выгоды — много систем с одной матрицей и разными $b$.