3.3. PLU-разложение#

Разложение, представленное в предыдущем разделе, существует не всегда.

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

Применим разложение к следующей матрице.

A = [
     1  2  -1  9;
     1  2   1  3;
     5  1   8  7;
    -8  6   5  1;
]
L, U = lufact(A)
L

4×4 LowerTriangular{Float64, Matrix{Float64}}:
  1.0     ⋅      ⋅    ⋅ 
  1.0  NaN       ⋅    ⋅ 
  5.0  -Inf   NaN     ⋅ 
 -8.0   Inf   NaN    1.0

U

4×4 UpperTriangular{Float64, Matrix{Float64}}:
 1.0  2.0  -1.0    9.0
  ⋅   0.0   2.0   -6.0
  ⋅    ⋅   Inf   -Inf
  ⋅    ⋅     ⋅   NaN

Матрицы разложения имеют NaN и Inf значения, возникшие в результате деления на ноль. При этом сама матрица обратимая.

@show det(A);

det(A) = 1196.0

На первом шаге разложения всё в порядке

L = diagm(0 => ones(4))
U = zeros(4, 4)
A₁ = float(copy(A))
U[1, :] .= A₁[1, :]
L[:, 1] .= A₁[:, 1] ./ U[1, 1]
L

4×4 Matrix{Float64}:
  1.0  0.0  0.0  0.0
  1.0  1.0  0.0  0.0
  5.0  0.0  1.0  0.0
 -8.0  0.0  0.0  1.0

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

4×4 Matrix{Float64}:
 0.0   0.0   0.0    0.0
 0.0   0.0   2.0   -6.0
 0.0  -9.0  13.0  -38.0
 0.0  22.0  -3.0   73.0

Из вида матрицы A₂ понятно, что на втором шаге возникнет деление на ноль. А именно, при присваивании матрице U новой строки всё нормально, первое деление на ноль возникает при подсчёте столбца матрицы L, поскольку U[2, 2] == A₂[2, 2]

U[2, :] .= A₂[2, :]
L[:, 2] .= A₂[:, 2] ./ U[2, 2];
L

4×4 Matrix{Float64}:
  1.0  NaN   0.0  0.0
  1.0  NaN   0.0  0.0
  5.0  -Inf  1.0  0.0
 -8.0   Inf  0.0  1.0

Полученные NaN и Inf значения распроcтранятся в результате операции

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

4×4 Matrix{Float64}:
 NaN  NaN  NaN   NaN
 NaN  NaN  NaN   NaN
 NaN  NaN   Inf  -Inf
 NaN  NaN  -Inf   Inf

Чтобы подобной ситуации для обратимых матриц не возникало, достаточно на итерации выбирать главный элемент.

3.3.1. Выбор главного элемента#

Гланым элементом (англ., pivot) называют диагональный элемент матрицы U[k, k], который возникает в строке функции Функция 3.4

for k in 1:n-1
    U[k, :] .= Aₖ[k, :]
    L[:, k] .= Aₖ[:, k] ./ U[k, k]
    Aₖ .-= L[:, k] * U[k, :]'
end

Чтобы избежать делений на ноль, необходимо выбрать в качестве главного элемента ненулевой \(p\)-ый элемент \(k\)-ого столбца матрицы Aₖ. И затем уже работать со строкой Aₖ[p, :] в теле цикла, которую называют главной строкой.

Выбрать главный элемент можно разными способами, однако, численная стабильность разложения достигается при выборе элемента с максимальным абсолютным значением.

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

Рассмотрим ту же матрицу, что и в начале раздела.

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

Но на этот раз будем производить выбор главного элемента.

p = argmax(abs.(A₁[:, 1]))
@show p
@show A₁[p, 1];

p = 4
A₁[p, 1] = -8.0

Функция argmax возвращает индекс наибольшего элемента массива. В качестве массива используется столбец из абсолютных значений столбца матрицы A₁[:, 1].

Проведём обновление матриц разложения и убедимся, что на этот раз L не содержит неопределённых float-значений.

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

4×4 Matrix{Float64}:
 -0.125  0.0  0.0  0.0
 -0.125  1.0  0.0  0.0
 -0.625  0.0  1.0  0.0
  1.0    0.0  0.0  1.0

Завершим шаг

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

4×4 Matrix{Float64}:
 0.0  2.75  -0.375  9.125
 0.0  2.75   1.625  3.125
 0.0  4.75  11.125  7.625
 0.0  0.0    0.0    0.0

Поскольку главным элементом был четвёртый, то и A₂ имеет нулевую четвёртую строку.

Продолжим шаги дальше

@show p = argmax(abs.(A₂[:, 2]))
U[2, :] .= A₂[p, :]
L[:, 2] .= A₂[:, 2] ./ U[2, 2]
A₃ = A₂ - L[:, 2] * U[2, :]'

@show p = argmax(abs.(A₃[:, 3]))
U[3, :] .= A₃[p, :]
L[:, 3] .= A₃[:, 3] ./ U[3, 3]
A₄ = A₃ - L[:, 3] * U[3, :]'

@show p = argmax(abs.(A₄[:, 4]))
U[4, :] .= A₄[p, :]
L[:, 4] .= A₄[:, 4] ./ U[4, 4];

p = argmax(abs.(A₂[:, 2])) = 3
p = argmax(abs.(A₃[:, 3])) = 1
p = argmax(abs.(A₄[:, 4])) = 2

Мы действительно получили разложение \(\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

При этом матрица \(\mathbf{U}\) имеет верхнетреугольную структуру

U

4×4 Matrix{Float64}:
 -8.0  6.0    5.0       1.0
  0.0  4.75  11.125     7.625
  0.0  0.0   -6.81579   4.71053
  0.0  0.0    0.0      -4.61776

А матрица \(\mathbf{L}\) теперь не нижнетреугольная.

L

4×4 Matrix{Float64}:
 -0.125  0.578947   1.0       -0.0
 -0.125  0.578947   0.706564   1.0
 -0.625  1.0       -0.0       -0.0
  1.0    0.0       -0.0       -0.0

Чтобы восстановить структуру матрицы \(\mathbf{L}\), достаточно переставить в ней строки согласно последовательности выбора главных элементов. В нашем случае это были \(4\), \(3\), \(1\) и \(2\). В Julia это осуществляется с помощью следующего синтаксиса

L = L[[4, 3, 1, 2], :]

4×4 Matrix{Float64}:
  1.0    0.0       -0.0       -0.0
 -0.625  1.0       -0.0       -0.0
 -0.125  0.578947   1.0       -0.0
 -0.125  0.578947   0.706564   1.0

Voilà!

Утверждение 3.9

LU-разложение с выбором главной строки выполнимо, если производится над обратимой матрицей [DB17].

Несмотря на это утверждение, существуют примеры обратимых матриц, для которых численный алгоритм разложения всё-таки не работает, однако на практике такие матрицы встречаются редко.

3.3.2. PLU-разложение#

Теперь есть всё, чтобы построить LU-разложение с выбором главной строки. Такое разложение называют PLU-разложением (англ., row-pivoting LU factorization).

Определение 3.10 (PLU-разложение)

Пусть дана матрица \(\mathbf{A}\) размера \(n\times n\). PLU-разложением называется тройка из

  • нижнетреугольной матрицы \(\mathbf{L}\);

  • верхнетреугольной матрицы \(\mathbf{U}\);

  • и перестановки \(i_1, i_2, \ldots, i_n\) индексов \(1, \ldots, n\).

Эта тройка такова, что

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

где строки \(1, \ldots, n\) матрицы \(\hat{\mathbf{A}}\) это строки \(i_1, i_2, \ldots, i_n\) исходной матрицы \(\mathbf{A}\).

Итак, для реализации теперь необходимо возвращать не только матрицы разложения, но и вектор-перестановку.

Функция 3.11 (plufact)

LU-разложение с выбором главной строки

"PLU-разложение матрицы `A`. Возвращает `L`, `U` и вектор-перестановку."
function plufact(A::AbstractMatrix)
    n = size(A, 1)
    p = zeros(Int, n)
    
    U = float(similar(A))
    L = similar(U)
    Aₖ = float(copy(A))
    
    for k in 1:n-1
        p[k] = argmax(abs.(Aₖ[:, k]))
        U[k, :] .= Aₖ[p[k], :]
        L[:, k] .= Aₖ[:, k] ./ U[k, k]
        Aₖ .-= L[:, k] * U[k, :]'
    end
    
    p[n] = argmax(abs.(Aₖ[:, n]))
    U[n, n] = Aₖ[p[n], n]
    L[:, n] = Aₖ[:, n] / U[n, n]
    
    return LowerTriangular(L[p, :]), UpperTriangular(U), p
end
Демонстрация 3.12 (Пример решения системы)

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

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

Как и в случае LU, cначала выполним разложение

L, U, p = plufact(A);
p

4-element Vector{Int64}:
 3
 4
 2
 1

Однако, для решения системы нужно сначала сделать перестановку в векторе правой части \(\mathbf{b}\)

z = forwardsub(L, b[p])
x = backwardsub(U, z)

4-element Vector{Float64}:
  0.8767123287671232
  0.06849315068493142
  0.10958904109589052
 -0.38356164383561636

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

b - A*x

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

PLU-разложение является центральным алгоритмом для решения систем линейных уравнений общего вида. Его вычислительная сложность (число операций с плавающей точкой: флопс, flops) составляет \(\propto \frac{2}{3} n^3\), где \(n\) – размер матрицы системы.

В дальнейшем, упоминая LU-разложение, мы для краткости будем подразумевать PLU-разложение.

3.3.3. PLU в Julia#

В Julia PLU разложение называется просто LU разложением и доступно в модуле LinearAlgebra.

  • Функция lu(A) возвращает объект разложения F, содержащий

    • множители F.L и F.U

    • матрицу перестановок F.P

    • вектор перестановок F.p

  • Для решения системы уравнений можно пользоваться самим объектом разложения

  • Решение системы можно выполнять

    • оператором F \ b, тогда решение будет аллоцироваться

    • функцией ldiv!(x, F, b), тогда решение не будет аллоцироваться и запишется в x

julia> A = [5  1  0  9; 4  2 -1  4; 8 -1  4  1; 5  7  4  6];
julia> b = [1, 2, 7, 3];
julia> F = lu(A);

julia> F \ b - A \ b
4-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0

julia> norm(F.L * F.U - A[F.p, :])
5.570606194861277e-16

julia> F.p
4-element Vector{Int64}:
 3
 4
 2
 1

julia> ldiv!(x, F, b);
julia> norm(x - F \ b)
0.0

В чём разница между A \ b и F \ b? У оператора \ есть несколько методов. Наиболее общий из них принимает на вход AbstractMatrix и вектор AsbtractVector. Этот метод проверяет, квадратная ли система, и является ли она треугольной. Если условия не выполнены, то применяется LU разложение, и вызывается второй метод – F \ b. Таким образом, в первом случае производятся проверки и LU разложение, а во втором нет – разложение уже хранится в F и достаточно применить прямую/обратную подстановки.