LU-разложение матрицы алгоритм: схема и пример

LU-разложение матрицы - это представление квадратной матрицы $A$ в виде произведения двух треугольных множителей: нижней треугольной $L$ и верхней треугольной $U$, то есть $A = LU$. Эта факторизация лежит в основе большинства прямых методов решения систем линейных уравнений и фактически является матричной записью обычного исключения Гаусса. Зная алгоритм LU-разложения матрицы, можно один раз разложить $A$, а затем дёшево решать систему $Ax = b$ для множества разных правых частей, считать определитель и обращать матрицу. Ниже разберём схему построения $L$ и $U$, роль перестановок и пошаговый пример.

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

Пусть дана невырожденная матрица $A \in \mathbb{R}^{n \times n}$. Её LU-разложение - это пара матриц

$$L = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ l_{21} & 1 & \cdots & 0 \\ \vdots & & \ddots & \vdots \\ l_{n1} & l_{n2} & \cdots & 1 \end{pmatrix}, \qquad U = \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ \vdots & & \ddots & \vdots \\ 0 & 0 & \cdots & u_{nn} \end{pmatrix},$$

таких, что $A = LU$. Матрица $L$ имеет единицы на главной диагонали (это так называемое разложение Дулиттла), а $U$ совпадает с матрицей, которую мы получаем в конце прямого хода метода Гаусса. Множители $l_{ij}$ - это в точности коэффициенты, на которые умножали строки при исключении. Существование разложения без перестановок гарантируется, если все ведущие (угловые) миноры матрицы отличны от нуля.

Дальше - интерактивный построитель: введите свою матрицу, и модель выполнит LU-разложение по шагам, покажет $L$, $U$ и (при необходимости) матрицу перестановок $P$.

Алгоритм LU-разложения матрицы

В основе лежит исключение по Гауссу. Идём по столбцам $k = 1, 2, \ldots, n-1$ и для каждой строки $i > k$ вычисляем множитель

$$l_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}},$$

после чего обнуляем поддиагональные элементы, обновляя оставшуюся часть матрицы:

$$a_{ij}^{(k+1)} = a_{ij}^{(k)} - l_{ik}\, a_{kj}^{(k)}, \qquad j = k, k+1, \ldots, n.$$

Получаемые множители $l_{ik}$ складываются в нижнюю треугольную $L$, а итоговая верхняя треугольная матрица и есть $U$. Главное наблюдение алгоритма: обнулённые места под диагональю можно «переиспользовать» для хранения $l_{ik}$, поэтому на практике $L$ и $U$ хранят прямо на месте исходной матрицы $A$, не тратя дополнительной памяти. Диагональный элемент $a_{kk}^{(k)}$, на который делим, называется ведущим (pivot).

Пример LU-разложения

Разложим матрицу

$$A = \begin{pmatrix} 2 & 1 & 1 \\ 4 & 3 & 3 \\ 8 & 7 & 9 \end{pmatrix}.$$

Первый столбец: ведущий элемент $a_{11} = 2$. Множители $l_{21} = 4/2 = 2$ и $l_{31} = 8/2 = 4$. Вычитаем из второй строки удвоенную первую, из третьей - учетверённую:

$$\begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 3 & 5 \end{pmatrix}.$$

Второй столбец: ведущий $a_{22} = 1$, множитель $l_{32} = 3/1 = 3$. Вычитаем из третьей строки утроенную вторую - получаем $U$. Собираем результат:

$$L = \begin{pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 3 & 1 \end{pmatrix}, \qquad U = \begin{pmatrix} 2 & 1 & 1 \\ 0 & 1 & 1 \\ 0 & 0 & 2 \end{pmatrix}.$$

Легко проверить, что $LU = A$. Определитель сразу читается с диагонали $U$: $\det A = 2 \cdot 1 \cdot 2 = 4$.

Решение СЛАУ через LU

Главная ценность алгоритма LU-разложения матрицы - быстрое решение системы $Ax = b$. Подставив $A = LU$, получаем $LUx = b$. Вводим вспомогательный вектор $y = Ux$ и решаем задачу в два прохода:

$$Ly = b \quad \text{(прямая подстановка)}, \qquad Ux = y \quad \text{(обратная подстановка)}.$$

Так как $L$ нижняя треугольная, $y$ находится сверху вниз:

$$y_i = b_i - \sum_{j=1}^{i-1} l_{ij}\, y_j.$$

Затем из верхней треугольной $U$ находим $x$ снизу вверх:

$$x_i = \frac{1}{u_{ii}}\left(y_i - \sum_{j=i+1}^{n} u_{ij}\, x_j\right).$$

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

Выбор ведущего элемента и матрица перестановок

Если на очередном шаге ведущий элемент $a_{kk}^{(k)}$ оказывается нулём, делить нельзя - разложение $A = LU$ без перестановок не существует. Даже когда pivot просто мал по модулю, деление на него раздувает погрешности округления. Лекарство - частичный выбор ведущего элемента (partial pivoting): на шаге $k$ среди элементов столбца ниже диагонали ищут максимальный по модулю и переставляют соответствующую строку наверх.

Перестановки строк фиксируются в матрице перестановок $P$, и разложение принимает вид

$$PA = LU.$$

Тогда система решается как $L y = Pb$, $U x = y$. Частичный выбор делает алгоритм численно устойчивым практически для всех задач, встречающихся на практике, поэтому именно $PA = LU$ реализовано в библиотеках вроде LAPACK и в функциях lu пакетов NumPy/SciPy. Близкая по духу идея устойчивости разбирается в заметке про число обусловленности матрицы: чем хуже обусловлена $A$, тем важнее аккуратный выбор pivot.

Определитель и обращение через LU

Поскольку определитель произведения равен произведению определителей, а $\det L = 1$ (единичная диагональ), определитель исходной матрицы получается почти бесплатно:

$$\det A = \det L \cdot \det U = \prod_{i=1}^{n} u_{ii}.$$

При наличии перестановок добавляется знак: $\det A = (-1)^{s}\prod_{i} u_{ii}$, где $s$ - число выполненных перестановок строк. Обратную матрицу $A^{-1}$ тоже находят через LU: решают $n$ систем $A x = e_k$ с единичными ортами $e_k$ в правых частях, используя одно и то же разложение. Это заметно дешевле, чем считать $\det A$ по определению или обращать матрицу методом присоединённой матрицы.

Когда LU не подходит

LU-разложение универсально, но не всегда оптимально. Для симметричной положительно определённой матрицы выгоднее разложение Холецкого $A = LL^{\top}$ - оно вдвое экономнее и не требует выбора ведущего элемента. Для больших разреженных систем прямое LU-разложение «забивает» нули (fill-in) и съедает память, поэтому переходят к итерационным методам - например, к методу Якоби или методу сопряжённых градиентов. Для прямоугольных и плохо обусловленных задач вместо LU берут QR-разложение или SVD. То есть LU - это рабочая лошадка для плотных квадратных систем среднего размера.

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

• Деление на нулевой или малый pivot без перестановок. Если на диагонали оказался ноль, нужен частичный выбор и переход к $PA = LU$; иначе алгоритм просто упадёт или выдаст мусор.
• Путаница, где $L$, а где $U$. $L$ - нижняя с единицами на диагонали, $U$ - верхняя; множители $l_{ik}$ записывают именно в $L$, а не в $U$.
• Забытый знак определителя. При нечётном числе перестановок строк определитель меняет знак - множитель $(-1)^{s}$ терять нельзя.
• Повторное разложение для каждой правой части. Смысл LU в том, чтобы разложить $A$ один раз; на каждый новый $b$ нужны только две дешёвые подстановки.
• Применение к разреженным матрицам без учёта fill-in. Для больших разреженных систем плотное LU неэффективно - там уместнее итерационные методы.

FAQ

Чем LU-разложение отличается от метода Гаусса? По сути это одно и то же исключение, но записанное как факторизация $A = LU$. Метод Гаусса решает систему за один проход, а LU отделяет «дорогую» часть (разложение) от «дешёвой» (подстановки), позволяя переиспользовать $L$ и $U$ для разных $b$.

Всегда ли существует LU-разложение? Без перестановок - только если все угловые миноры ненулевые. С частичным выбором ведущего элемента разложение $PA = LU$ существует для любой невырожденной матрицы.

Какова сложность алгоритма? Само разложение требует около $\tfrac{2}{3}n^3$ арифметических операций, а каждая последующая прямая/обратная подстановка - порядка $n^2$.

Коротко

LU-разложение матрицы представляет $A$ как произведение нижней и верхней треугольных множителей $L$ и $U$ и является матричной формой метода Гаусса. Алгоритм строит $L$ из множителей исключения, а $U$ - из верхней треугольной части; при нулевых или малых ведущих элементах добавляют частичный выбор и матрицу перестановок $P$, получая $PA = LU$. Готовое разложение позволяет дёшево решать $Ax = b$ для многих правых частей через прямую и обратную подстановку, мгновенно считать определитель как произведение диагонали $U$ и обращать матрицу.