3.2. LU-разложение

На практике оказывается, что проще работать не с самой матрицей, а с представлением её в виде произведения нескольких матриц, которое называется разложением. В данном разделе выводится алгоритм решения системы уравнений, представляющий исходную матрицу в виде произведения \(\mathbf{A}=\mathbf{L}\mathbf{U}\), после чего пользующийся прямой и обратной подстановкой.

3.2.1. Внешнее произведение

Вывод алгоритма разложения основывается на внешнем произведении векторов (outer product). Для двух векторов \(\mathbf{u} \in \real^m\), \(\mathbf{v} \in \real^n\) внешнее произведение является матрицей размера \(m \times n\)

\[\begin{split}\mathbf{u}\mathbf{v}^\top = \begin{bmatrix} u_1 v_1 & u_1 v_2 & \cdots & u_1 v_n \\ u_2 v_1 & u_2 v_2 & \cdots & u_2 v_n \\ \vdots & \vdots & & \vdots \\ u_m v_n & u_m v_n & \cdots & u_m v_n \end{bmatrix}.\end{split}\]

Другими словами, столбец \(\mathbf{u}\) умножается на строку \(\mathbf{v}\).

Произведение матриц \(\mathbf{A}\mathbf{B}\) является суммой внешних произведений столбцов \(\mathbf{a}_{k}\) матрицы \(\mathbf{A}\) на строки \(\mathbf{b}^{\top}_{k}\) матрицы \(\mathbf{B}\), т.е.

(3.1)\[\mathbf{A}\mathbf{B} = \sum_{k=1}^n \mathbf{a}_k \mathbf{b}_k^\top.\]

Рассмотрим теперь произведение нижнетреугольной матрицы на верхнетреугольную \(\mathbf{L}\mathbf{U}\). Обозначим \(\mathbf{l}_k\) столбцы матрицы \(\mathbf{L}\), а \(\mathbf{u}_k^\top\) – строки матрицы \(\mathbf{U}\). Тогда произведение запишется в виде (3.1)

\[\mathbf{L}\mathbf{U} = \sum_{k=1}^n \mathbf{l}_k \mathbf{u}_k^\top.\]

Первая строка этого произведения имеет вид (здесь \(\mathbf{e}_1^\top\) – первая строка единичной матрицы, домножение на неё слева выделяет первую строку)

(3.2)\[\mathbf{e}_1^\top \sum_{k=1}^n \mathbf{l}_k \mathbf{u}_k^\top = \sum_{k=1}^n (\mathbf{e}_1^\top \mathbf{l}_k) \mathbf{u}_k^\top = L_{11} \mathbf{u}_1^\top,\]

поскольку только первый элемент столбца \(\mathbf{l}_k\) ненулевой только у \(\mathbf{l}_1\). Посмотрите на схему ниже

     L               U

| □ ⋅ ⋅ ⋅ |     | □ □ □ □ |
| □ □ ⋅ ⋅ |     | ⋅ □ □ □ |
| □ □ □ ⋅ |     | ⋅ ⋅ □ □ |
| □ □ □ □ |     | ⋅ ⋅ ⋅ □ |

Первый же столбец \(\mathbf{L}\mathbf{U}\) имеет вид

(3.3)\[\Bigg( \sum_{k=1}^n \mathbf{l}_k \mathbf{u}_k^\top \Bigg) \mathbf{e}_1 = \sum_{k=1}^n \mathbf{l}_k (\mathbf{u}_k^\top \mathbf{e}_1) = U_{11} \mathbf{l}_1.\]

Воспользуемся свойствами (3.2) и (3.3) для вывода LU-разложения.

3.2.2. Разложение

Определение 3.2 (LU-разложение)

Для квадратной матрицы \(\mathbf{A}\) размера \(n \times n\) необходимо найти

\[ \mathbf{A} = \mathbf{LU}, \]

где \(\mathbf{L}\) – нижнетреугольная матрица, а \(\mathbf{U}\) – верхнетреугольная матрица.

LU-разложение сводится к нахождению \(n^2 + n\) неизвестных.

Пусть диагональные элементы \(\mathbf{L}\) единичны \(L_{11} = L_{22} = \cdots = L_{nn} = 1\). Теперь, зная \(L_{11} = 1\), применим (3.2), и получим

\[\mathbf{a}_1^\top = \mathbf{e}_1^\top \mathbf{A} = \mathbf{e}_1^\top (\mathbf{LU}) = \mathbf{u}_1^\top,\]

т.е. первая строка матрицы \(\mathbf{U}\) равняется первой строке матрицы \(\mathbf{A}\), \(\mathbf{u}_1^\top = \mathbf{a}_1^\top\).

На этом шаге нам стало известно значение \(U_{11}\). Воспользуемся теперь (3.3)

\[\mathbf{a}_1 = \mathbf{A} \mathbf{e}_1 = (\mathbf{LU}) \mathbf{e}_1 = U_{11} \mathbf{l}_1,\]

откуда следует выражение для первого столбца матрицы \(\mathbf{L}\)

\[\mathbf{l}_1 = \frac{\mathbf{a}_1}{U_{11}}.\]

Теперь вычтем из исходной матрицы внешнее произведение

\[\mathbf{A}_2 = \mathbf{A} - \mathbf{l}_1 \mathbf{u}_1^\top = \sum_{k=2}^n \mathbf{l}_k \mathbf{u}_k^\top.\]

Полученная матрица \(\mathbf{A}_2\) имеет нулевую первую строку и нулевой первый столбец.

Таким образом, задача свелась к предыдущей, но, фактически, на квадратной матрице меньшего размера. Повторяя эту процедуру дальше, получим остальные строки и столбцы разложения. Ниже приведена схема вычислений и демонстрация на небольшой матрице.

  A₁ = A             A₂              A₃              A₄

| □ □ □ □ |     | ⋅ ⋅ ⋅ ⋅ |     | ⋅ ⋅ ⋅ ⋅ |     | ⋅ ⋅ ⋅ ⋅ |
| □ □ □ □ |     | ⋅ □ □ □ |     | ⋅ ⋅ ⋅ ⋅ |     | ⋅ ⋅ ⋅ ⋅ |
| □ □ □ □ |     | ⋅ □ □ □ |     | ⋅ ⋅ □ □ |     | ⋅ ⋅ ⋅ ⋅ |
| □ □ □ □ |     | ⋅ □ □ □ |     | ⋅ ⋅ □ □ |     | ⋅ ⋅ ⋅ □ |

Демонстрация 3.3

Рассмотрим пример работы алгоритма на матрице размера \(4 \times 4\). Исходную матрицу обозначим A₁, нижнетреугольную L проинициализируем единичной матрицей, а верхнетреугольную U нулями.

A₁ = Float64[
   5  1  0  9;
   4  2 -1  4;
   8 -1  4  1;
   5  7  4  6;
]
L = diagm(0 => ones(4))
U = zeros(4, 4);

Первый шаг алгоритма основывается на (3.2): вычисляем первую строку U

U[1, :] .= A₁[1, :]
U

4×4 Matrix{Float64}:
 5.0  1.0  0.0  9.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

Теперь вычислим первый столбец L (3.3)

L[:, 1] .= A₁[:, 1] ./ U[1, 1]
L

4×4 Matrix{Float64}:
 1.0  0.0  0.0  0.0
 0.8  1.0  0.0  0.0
 1.6  0.0  1.0  0.0
 1.0  0.0  0.0  1.0

Вычтем из A₁ внешнее произведение первого столбца L и первой строки U

A₂ = A₁ - L[:, 1] * U[1, :]'

4×4 Matrix{Float64}:
 0.0   0.0   0.0    0.0
 0.0   1.2  -1.0   -3.2
 0.0  -2.6   4.0  -13.4
 0.0   6.0   4.0   -3.0

\(\mathbf{A}_2 = \mathbf{l}_2 \mathbf{u}_2^\top + \mathbf{l}_3 \mathbf{u}_3^\top + \mathbf{l}_4 \mathbf{u}_4^\top\). Так, задача свелась к предыдущей, на матрице размера \(3 \times 3\).

U[2, :] .= A₂[2, :]
L[:, 2] .= A₂[:, 2] ./ U[2, 2]
A₃ = A₂ - L[:, 2] * U[2, :]'

4×4 Matrix{Float64}:
 0.0  0.0  0.0        0.0
 0.0  0.0  0.0        0.0
 0.0  0.0  1.83333  -20.3333
 0.0  0.0  9.0       13.0

Предпоследний шаг

U[3, :] .= A₃[3, :]
L[:, 3] .= A₃[:, 3] ./ U[3, 3]
A₄ = A₃ - L[:, 3] * U[3, :]'

4×4 Matrix{Float64}:
 0.0  0.0  0.0    0.0
 0.0  0.0  0.0    0.0
 0.0  0.0  0.0    0.0
 0.0  0.0  0.0  112.818

И последний

U[4, 4] = A₄[4, 4];

Теперь известно разложение матрицы

U

4×4 Matrix{Float64}:
 5.0  1.0   0.0        9.0
 0.0  1.2  -1.0       -3.2
 0.0  0.0   1.83333  -20.3333
 0.0  0.0   0.0      112.818

L

4×4 Matrix{Float64}:
 1.0   0.0      0.0      0.0
 0.8   1.0      0.0      0.0
 1.6  -2.16667  1.0      0.0
 1.0   5.0      4.90909  1.0

Проверим, что \(\mathbf{A} = \mathbf{LU}\). Из-за погрешность вычислений с плавающей точкой, возможны неточности.

A₁ - L * U

4×4 Matrix{Float64}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

Следуя примеру выше несложно написать реализацию.

Функция 3.4 (lufact)

LU-разложение (нестабильное)

"Нестабильное LU-разложение квадратной матрицы `A`. Возвращает `L`, `U`."
function lufact(A::AbstractMatrix)
    n = size(A, 1)
    L = diagm(0 => ones(n))
    U = zeros(n, n)
    Aₖ = float(copy(A))
    
    for k in 1:n-1
        U[k, :] .= Aₖ[k, :]
        L[:, k] .= Aₖ[:, k] ./ U[k, k]
        Aₖ .-= L[:, k] * U[k, :]'
    end
    U[n, n] = Aₖ[n, n]
    return LowerTriangular(L), UpperTriangular(U)
end

Примечание

В строке Aₖ = float(copy(A)) используется копирование входной матрицы, поскольку в Julia массивы передаются в функцию по ссылке.

В теле цикла использован броадкаст операций для избежания выделений памяти. Например, операция Aₖ[:, k] / U[k, k] выделяет вектор Aₖ[:, k], а затем выделяет новый вектор в результате деления. В случае с броадкастом операций такие выделения памяти не происходят.

3.2.3. Решение системы уравнений

Алгоритм 3.5 (Решение системы уравнений через LU-разложение)

Пусть дана система уравнений \(\mathbf{A x} = \mathbf{b}\). Алгоритм решения следующий

1. Разложить матрицу \(\mathbf{A} = \mathbf{LU}\);

2. Решить \(\mathbf{L z} = \mathbf{b}\), используя прямую подставновку;

3. Решить \(\mathbf{U x} = \mathbf{z}\), используя обратную подстановку.

Решение системы уравнений через LU-разложение имеет преимущество над методом Гаусса. Так, в методе Гаусса составляется матрица \([\mathbf{A}\:\mathbf{b}]\), где левая часть приводится к треугольному виду. В результате, для нескольких систем с разными правыми частями приходится проводить процедуру заново. В LU-разложении достаточно один раз разложить исходную матрицу, а затем воспользоваться шагами 2 и 3 алгоритма выше для разных правых частей \(\mathbf{b}\).

Демонстрация 3.6 (Пример решения)

Рассмотрим решение системы

A = [5  1  0  9; 4  2 -1  4; 8 -1  4  1; 5  7  4  6]
b = [1, 2, 7, 3];

Сначала выполним разложение, а затем последовательно воспользуемся прямой и обратной подстановками

L, U = lufact(A)
z = forwardsub(L, b)
x = backwardsub(U, z)

4-element Vector{Float64}:
  0.876712328767123
  0.06849315068493178
  0.10958904109589054
 -0.38356164383561636

Невязка при этом составила

b - A*x

4-element Vector{Float64}:
 -4.440892098500626e-16
  2.220446049250313e-16
  8.881784197001252e-16
 -1.7763568394002505e-15