$$
\begin{array}{cllll|l}
1 & c_{1,2} & \dots & c_{1,n-1} & c_{1,n} & y_1 \\
& a^{(1)}_{2,2} & \dots & a^{(1)}_{2,n-1} & a^{(1)}_{2,n} & b^{(1)}_2 \\
& \vdots & \vdots & \vdots & \vdots & \vdots \\
& a^{(1)}_{n-1, 2} & \dots & a^{(1)}_{n-1, n-1} & a^{(1)}_{n-1, n} & b^{(1)}_{n-1} \\
& a^{(1)}_{n, 2} & \dots & a^{(1)}_{n, n-1} & a^{(1)}_{n,n} &b^{(1)}_n
\end{array}
$$
В каждом столбце обнуляем часть под диагональю с использованием элементарного преобразования «вычесть строку с коэффициентом» (метод применим только в случае det Ak ≠ 0). Сложность — $2/3 n^3 + O(n^2) $
$$ \mathbf{a}_k = \dfrac{\mathbf{a}_k}{a_{k,k}}, \qquad y_k = \dfrac{b_k}{a_{k,k}}; \qquad \mathbf{a}_i = \mathbf{a}_i - \dfrac{a_{i,k}}{a_{k,k}} \mathbf{a}_k, \qquad b_i = b_i - \dfrac{a_{i,k}}{a_{k,k}} y_k $$
$$
\begin{array}{ccllll|l}
& \bbox[#afa,3pt]{i} & \bbox[#afa,3pt]{⟶} \\
& 1 & c_{1,2} & \dots & c_{1,n-1} & c_{1,n} & y_1 \\
\bbox[#ff0,3pt] j & & a^{(1)}_{2,2} & \dots & a^{(1)}_{2,n-1} & a^{(1)}_{2,n} & b^{(1)}_2 \\
\bbox[#ff0,3pt] ↓& & \vdots & \vdots & \vdots & \vdots & \vdots \\
& & a^{(1)}_{n-1, 2} & \dots & a^{(1)}_{n-1, n-1} & a^{(1)}_{n-1, n} & b^{(1)}_{n-1} \\
& & a^{(1)}_{n, 2} & \dots & a^{(1)}_{n, n-1} & a^{(1)}_{n,n} &b^{(1)}_n \\
& & \bbox[#faa,3pt]k & \bbox[#faa,3pt]⟶
\end{array}
$$
for i := 0…(n-1) for k := (i+1)…(n-1) A[i, k] := A[i, k] / A[i, i] A[i, i] := 1 for j := (i+1)…(n-1) for k := (i+1)…(n-1) A[j, k] := A[j, k] - A[i, k] * A[j, i] for k := 0…(m-1) # размер правой части (1 или n) B[j, k] := B[j, k] - B[i, k] * A[j, i]
На самом деле, предыдущие преобразования могут быть представлены в форме умножения на матрицы особого вида.
for k := (i+1)…(n-1) A[i, k] := A[i, k] / A[i, i] A[i, i] := 1
$A'^{(i)} = D_i · A^{(i-1)} $ , где $D_i = \mathrm{diag} [ 1, 1, 1, …, {(a^{(i-1)}_{i,i})^{-1}}, 1, …, 1] $ (единичная с обратным к элементу $a_{i,i} $ на $i $-м месте)
for k := (i+1)…(n-1) A[j, k] := A[j, k] - A[i, k] * A[j, i]
$A^{(i)} = L_i · A'^{(i)} $, где $L_i $: (единичная со столбцом обратных по сложению элементов под $i $-м диагональным элементом)
$$
\begin{pmatrix}
1 \\
& \ddots & \\
& & 1 \\
& & -a^{(i-1)}_{i+1, i} & 1 \\
& & \vdots & & \ddots \\
& & -a^{(i-1)}_{n, i} & & & 1
\end{pmatrix}
$$
Таким образом, применение прямого хода метода Гаусса можно представить как:
Ln · Dn · … · L1 · D1 · A = U
где $U $ — верхнетреугольная матрица с единицами на диагонали.
Тогда $A = LU $, где $L = (L_n · D_n · … · L_1 · D_1)^{-1} $ — нижнетреугольная матрица (т.к. нижнетреугольные матрицы замкнуты по умножению).
С использованием такого разложения можно далее решать линейные системы с правой частью $A $ и произвольной левой частью за $O(n^2) $ операций:
L · U · x = b
Формулы построены таким образом, что обе матрицы $L, U $ можно (и нужно для выполнения требований) хранить на месте исходной матрицы $A $ (единичная диагональ $U $ не хранится).
$$
\begin{pmatrix}
{l_{1,1}} & u_{1,2} & \bbox[#faf,3px]{u_{1, 3}} & \bbox[#ffa,3px]{u_{1, 4}} & … & & & u_{1,n} \\
{l_{2, 1}} & {l_{2, 2}} & \bbox[#faf,3px]{u_{2, 3}} & \bbox[#ffa,3px]{u_{2, 4}} & … & & & u_{2, n} \\
\bbox[#aff,3px]{l_{3, 1}} & \bbox[#aff,3px]{l_{3, 2}} & \bbox[#faa,3px]{l} & \bbox[#afa,3px]{u} & \bbox[#afa,3px]⟶ & & & a_{2, n} \\
\bbox[#faa,3px]↪︎ & \vdots & \vdots & & \ddots & & & \vdots \\
l_{n, 1} & l_{n, 2} & a_{n,3} & & … & & & a_{n, n}
\end{pmatrix}
$$
$$
\begin{array}{ccccc|c}
c_{1,1} & c_{1,2} & \dots & c_{1,n-1} & {c_{1,n}} & y_1 \\
& \bbox[#ffa,3px]{c_{2,2}} & \bbox[#faf,3px]\dots & \bbox[#faf,3px]{c_{2,n-1}} & \bbox[#faf,3px]{c_{2,n}} & \bbox[#fea,3px]{y_2} \\
& &\ddots & \vdots & {\vdots} & \bbox[#aef,3px]{\vdots} \\
& & & c_{n-1,n-1} & {c_{n-1, n}} & \bbox[#aef,3px]{y_{n-1}} \\
& & & & c_{n,n} & \bbox[#aef,3px]{y_n}
\end{array}
$$
$ x_n = y_n \big/ c_{n,n} \qquad \qquad \bbox[#fea,3px]{x_i} = \left( \bbox[#fea,3px]{y_i }- \sum_{j=i+1}^n \bbox[#faf,3px]{c_{i,j}} · \bbox[#aef,3px]{x_j} \right) \big/ \bbox[#ffa,3px]{c_{i,i}} $
Сложность — $\dfrac{n·(n-1)}{2} ≈ O(n^2)$
Во внутреннем цикле придерживаемся упорядоченного доступа к памяти.
for i:=(n-1) … 0
for j:= (i+1) … n
b[i] := b[i] - A[i, j] · b[j]
b[i] := b[i] / A[i,i]
for i := (n-1) … 0
for j := (i+1) … n
for k := 0 … (n-1)
B[i, k] := B[i, k] - A[i, j] * B[j, k]
for k := 0 … (n-1)
B[i, k] := B[i, k] / A[i, i]
$$
\begin{array}{cllll|l}
1 & \bbox[#afa,5px]{0} & \dots & c_{1,n-1} & c_{1,n} & y_1 \\
& 1 & \dots & a^{(2)}_{2,n-1} & a^{(2)}_{2,n} & b^{(2)}_2 \\
& & \vdots & \vdots & \vdots & \vdots \\
& & \dots & a^{(2)}_{n-1, n-1} & a^{(2)}_{n-1, n} & b^{(2)}_{n-1} \\
& & \dots & a^{(2)}_{n, n-1} & a^{(2)}_{n,n} &b^{(2)}_n
\end{array}
$$
Единственное отличие — обнуляются элементы столбца не только под, но и над диагональю.
for i := 0…(n-1) MultiplyRow(A[i, i…(n-1)], 1 / A[i, i]) for j := 0…(i-1), (i+1)…(n-1) A[j, i…(n-1)] := SubtractRow(A[j, i…(n-1)], A[i, i…(n-1)], A[j, i]) B[j, 0…(n-1)] := SubtractRow(B[j, 0…(n-1)], B[i, 0…(n-1)], A[j, i])
В чистом виде метод Гаусса и разновидности (LU-разложение, метод Жордана) неприменим к матрицам, у которых есть нулевой главный минор
Например, метод не сработает если a1, 1 = 0
Для исправления ситуации можно использовать элементарные преобразования вида «перестановка строк» или «перестановка столбцов» — переставим на нужное место строку (столбец), где нужный элемент отличен от нуля.
Как выбирать какую именно строку или столбец переставлять?
По итогам первого шага (деление строки на диагональный элемент) погрешность изменяется обратно пропорционально модулю диагонального элемента:
$$ ε^{(1)}_{i,j} = ε_{i,j} - ε_{1,j} · \dfrac{a_{i,1}}{a_{1,1}} $$
⇒ выгоднее выбирать наибольший по модулю элемент
Переставляем строки — «с выбором главного элемента по столбцу»
Переставляем столбцы — «с выбором главного элемента по строке»
И строки и столбцы — «с выбором главного элемента по всей матрице»
Асимптотика — по столбцу или по строке сохраняется ($2 \big/ 3 n^3 + O(n^2)$), по всей матрице — n3 + O(n2)). Перестановка элементов в памяти требует лишние O(n2) операций.
Возможен вариант (используется в реализациях с распределенной памятью), при котором перестановка не выполняется, а используется массив индексов перестановок.
Все методы с выбором главного элемента применимы к любым невырожденным матрицам.
Преобразования правой части:
Если мы делаем перестановку столбцов, нужно запоминать все транспозиции (τi — транспозиция $i ↔ $ <номер столбца, в котором находится главный элемент $i$-го шага>):
σ = τ1 · τ2 · … · τn
Для получения переменных в правильной последовательности нужно делать обратные транспозиции, начиная с последней (τn−1 · τn − 1−1 · …) — фактически, перемешать строки правой части согласно обратной перестановке столбцов левой.
Rows := [0, 1, 2, …, (n-1)]
Cols := [0, 1, 2, …, (n-1)]
Get(A, i, j) := A[Rows[i], Cols[j]]
for i := 0…(n-1)
max_col, max_row :=i, max_abs := |Get(A, i, i)|
for j := i…(n-1), k:= i…(n-1)
if |Get(A, j, k)| > max_abs
max_abs, max_col, max_row := |Get(A, j, k)|, j, k
swap(Rows, i, j)
swap(Cols, j, k)
MultiplyRow(Get(A, i, i…(n-1)), 1 / Get(A, i, i))
for j := 0…(i-1), (i+1)…(n-1)
B[j, 0…(n-1)] := SubtractRow(Get(B, j, 0…(n-1)), Get(B, i, 0…(n-1)), Get(A, j, i))
A[j, i…(n-1)] := SubtractRow(Get(A, j, i…(n-1)), Get(A, i, i…(n-1)), Get(A, j, i))
for i:=0…(n-1)
SwapRows(B, i, Cols[i])