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

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

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

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

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

Схема решения 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·Lᵀ (вдвое меньше работы, не нужен 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·U — дорогая работа (O(n³)) делается один раз и зависит только от A.
Каждую новую правую часть решают двумя треугольными проходами за O(n²).
Даёт почти бесплатно определитель и обращение; на практике — с выбором главного элемента (PA=LU) или Холецкого для SPD.
Главный сценарий выгоды — много систем с одной матрицей и разными b.