Методы последовательного исключения
Зачастую конечным результатом применения тех или иных сложных численных методов решения дифференциальных или иных уравнений является система линейных алгебраических уравнений (СЛАУ), матричная форма которой имеет вид:
\begin{equation}
\begin{bmatrix}
a_{11} & \dots & a_{1n} \\
\vdots & & \vdots \\
a_{n1} & \dots & a_{nn}
\end{bmatrix}
\begin{bmatrix}
x_{1} \\
\vdots \\
x_{n}
\end{bmatrix}
=
\begin{bmatrix}
b_{1} \\
\vdots \\
b_{n}
\end{bmatrix}
\end{equation}
или кратко:
\begin{equation}
\label{eq:matrix_eq}
\bm{A} \bm{x} = \bm{b},
\end{equation}
где $\bm{A} \in \mathbb{R}^{n \times n}$, $\bm{x} \in \mathbb{R}^{n}$ и $\bm{b} \in \mathbb{R}^{n}$.
Так как в сложных численных методах решение СЛАУ после некоторых простых манипуляций дает полное численное решение задачи, очевидна необходимость нахождения оптимальных методов решения СЛАУ. Глобально большинство методов решения СЛАУ можно разделить на две группы:
Прямые методы позволяют найти точное решение СЛАУ \eqref{eq:matrix_eq}, в то время как итерационные методы рассчитывают такое $\bm{x}^{(k)}$, что $\bm{A} \bm{x}^{(k)} - \bm{b} \to \bm{0}$ при $k \to \infty$. Прямые методы обычно используются для СЛАУ малой размерности, в то время как для СЛАУ большой размерности используют итерационные методы.
Разнообразие как прямых, так и итерационных методов связано с разнообразием особенностей матриц, составляющих СЛАУ. Так, например, если матрица является положительно определенной, то из класса прямых методов логично выбрать разложение Холецкого.
Мы начнем рассмотрение методов численной линейной алгебры с прямых методов.
Метод Гаусса (метод последовательного исключения)
Метод Гаусса является базовым методом решения СЛАУ и состоит в приведении матрицы $\bm{A}$ к треугольному виду с помощью элементарных преобразований, после чего решение $\bm{x}$ находится сравнительно легко. Рассмотрим следующую расширенную матрицу:
\begin{equation}
\tilde{\bm{A}} = [\bm{A}, \bm{b}]
=
\left[\begin{array}{cccc|c}
a_{11} & a_{12} & \dots & a_{1n} & b_1 \\
a_{21} & a_{22} & \dots & a_{2n} & b_2 \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
a_{n1} & a_{n2} & \dots & a_{nn} & b_n
\end{array}\right].
\end{equation}
Для того, чтобы привести матрицу $\bm{A}$ к треугольному виду, необходимо последовательно обнулять элементы, находящиеся под главной диагональю. Так, если из строк $i = 2, 3, \dots, n$ отнять первую строку, домноженную на $\frac{a_{i1}}{a_{11}}$, то все элементы, расположенные под $a_{11}$ будут равны нулю. Мы обозначим это элементарное преобразование следующим образом:
\begin{equation}
(i) \longrightarrow (i) - \frac{a_{i1}}{a_{11}} (1), \quad i = 2, \dots, n
\end{equation}
где $(i)$ обозначает $i$-ю строку. Тогда, применив операцию к расширенной матрицу, имеем:
\begin{equation}
\tilde{\bm{A}}
=
\left[\begin{array}{cccc|c}
a_{11} & a_{12} & \dots & a_{1n} & b_1 \\
0 & a_{22}^{(1)} & \dots & a_{2n}^{(1)} & b_2^{(1)} \\
\vdots & \vdots & \vdots & \vdots & \vdots \\
0 & a_{n2}^{(1)} & \dots & a_{nn}^{(1)} & b_n^{(1)}
\end{array}\right],
\end{equation}
где $a_{ij}^{(1)} = a_{ij} - \frac{a_{i1}}{a_{11}} a_{1j}$ и $b_{i}^{(1)} = b_{i} - \frac{a_{i1}}{a_{11}} b_{1}$. Аналогично, для обнуления элементов, находящихся под диагональным элементом $a_{22}^{(1)}$, нам необходимо применить следующее элементарное преобразование:
\begin{equation}
(i) \longrightarrow (i) - \frac{a_{i2}}{a_{22}^{(1)}} (2), \quad i = 3, \dots, n
\end{equation}
Применяя подобные преобразования каскадно вплоть до последнего диагольного элемента, мы получаем верхнюю треугольную матрицу для $\bm{A}$:
\begin{equation}
\tilde{\bm{A}}
=
\left[\begin{array}{cccc|c}
a_{11} & a_{12} & \dots & \dots & a_{1n} & b_1 \\
0 & a_{22}^{(1)} & \dots & \dots & a_{2n}^{(1)} & b_2^{(1)} \\
\vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\
0 & 0 & \dots & a_{n-1, n-1}^{(n-2)} & a_{n-1, n}^{(n-2)} & b_n^{(n-2)} \\
0 & 0 & \dots & 0 & a_{nn}^{(n-1)} & b_n^{(n-1)}
\end{array}\right].
\end{equation}
Этот этап метода Гаусса называется
прямым ходом. Очевидно, что тот же метод можно использовать для получения нижней треугольной матрицы. Можно заметить, что решение СЛАУ, в которой матрица имеет треугольную форму, легко находится с помощью так называемого
обратного хода метода Гаусса:
\begin{align*}
x_n &= \frac{b_n^{(n-1)}}{a_{nn}^{(n-1)}}, \\
x_{n-1} &= \frac{b_{n-1}^{(n-2)} - a_{n-1, n}^{(n-2)} x_n}{a_{n-1, n-1}^{(n-2)}}, \\
&\dots \\
x_{1} &= \frac{b_1 - \sum_{i=2}^{n} a_{1i} x_i}{a_{11}}.
\end{align*}
Рассчитаем количество операций, необходимых для нахождения решения СЛАУ методом Гаусса. Так как вычислительное время, потребное для операций умножения и деления, сильно больше, чем вычислительное время сложения и вычитания, мы сконцентрируемся на расчете числа умножений и делений. В первую очередь заметим, что на $k$-й итерации из всего $(n-1)$ итераций прямого хода требуется $(n-k)$ делений для вычисления множителя $m_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$. Затем требуется для $(n-k)$ строк требуется произвести $(n-k+1)$ перемножений, т.е. суммарно $(n-k)(n-k+1)$. Таким образом, на $k$-й итерации мы имеем $n-k + (n-k)(n-k+1) = (n-k)(n-k+2)$ умножений и делений. Суммируя по $k = 1, \dots, n-1$, получаем общее число умножений и делений в прямом ходе метода Гаусса:
\begin{equation}
\begin{split}
\sum_{k=1}^{n-1} (n-k)(n-k+2) &= \sum_{k=1}^{n-1} (n-k)^2 + 2 \sum_{k=1}^{n-1} (n-k) \\
&= \sum_{k=1}^{n-1} k^2 + 2 \sum_{k=1}^{n-1} k \\
&= \frac{(n-1)n(2n-1)}{6} + 2 \frac{(n-1)n}{2} \\
&= \frac{2n^3 + 3n^2 -5n}{6}.
\end{split}
\end{equation}
В обратном ходе метода Гаусса для нахождения $x_k$ требуется $n-k$ перемножений и одно деление, что в сумме дает:
\begin{equation}
\begin{split}
\sum_{k=1}^{n} (n-k + 1) &= \sum_{k=1}^{n} k \\
&= \frac{n(n+1)}{2} \\
&= \frac{n^2 + n}{2}.
\end{split}
\end{equation}
Наконец, суммарно требуется следующее число операций для нахождения решения методом Гаусса:
\begin{equation}
\begin{split}
\frac{2n^3 + 3n^2 -5n}{6} + \frac{n^2 + n}{2} = O(n^3).
\end{split}
\end{equation}
Метод Гаусса с выбором главного элемента
В случае, когда на $k$-й итерации алгоритма диагональный элемент оказывается равен нулю, мы не имеем возможности рассчитать множитель $m_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}$. Эта проблема решается соответствующим переставлением строк. Если элемент $a_{kk}^{(k)}$ мал по сравнению с $a_{ik}^{(k)}$, множитель $m_{ik}$ будет иметь большое значение, что приводит к усилению погрешности округления в результате ряда последующих умножений на $m_{ik}$. Более того, деление на малый элемент $a_{kk}^{(k)}$ происходит также во время обратного года Гаусса и усиливает накопленные погрешности в числителе. Простейший способ избежать подобного поведения состоит в перестановке строк матрицы так, что диагональным элементом становится наибольший по модулю элемент $k$-го столбца:
\begin{equation}
|a_{kk}^{(k)}| = \max_{k \leq i \leq n} |a_{ik}^{(k)}|.
\end{equation}
Такой подход называется
частичным выбором главного элемента. Логичным развитием является перестановка строк и столбцов так, что диагональным элементом становится наибольший по модулю элемент $|a_{ij}^{(k)}|$, где $k \leq i \leq n$ и $k \leq j \leq n$:
\begin{equation}
|a_{kk}^{(k)}| = \max_{k \leq i, j \leq n} |a_{ij}^{(k)}|.
\end{equation}
Такой подход именуется
полным выбором главного элемента. Легко убедиться, что полный выбор главного элемента требует суммарно $O(n^3)$ сравнений и таким образом рекомендуется к использованию только тогда, когда решаемая задача особенно чувствительна к погрешности округления.
LU-разложение
Предположим, что матрица $\bm{A}$ может быть разложена в нижнюю и верхнюю треугольные матрицы:
\begin{equation}
\bm{A} = \bm{L}\bm{U}.
\end{equation}
Такое разложение называется
LU-разложение. В таком случае система линейных уравнений $\bm{A} \bm{x} = \bm{b}$ решается в два шага. Обозначив вектор $\bm{y}$ как $\bm{y} = \bm{U} \bm{x}$, мы находим его как решение уравнения:
\begin{equation}
\bm{L} \bm{y} = \bm{b}.
\end{equation}
Вектор $\bm{x}$ тогда является решением уравнения:
\begin{equation}
\bm{U} \bm{x} = \bm{y}.
\end{equation}
Каждый из шагов требует $O(n^2)$ операций. Таким образом, один раз разложив матрицу $\bm{A}$ в нижнюю и верхнюю треугольные матрицы, мы можем находить решение $\bm{x}$ для различных $\bm{b}$ в $O(n^2)$ шагов. LU-разложение связано с прямым ходом метода Гаусса и соответственно требует $O(n^3)$ операций.
Получим явные выражения для матриц $\bm{L}$ и $\bm{U}$. Пусть перед $k$-й итерацией прямого хода метода Гаусса матрица $\bm{A}$ находится в форме $\bm{A}^{(k)}$. Тогда $k$-я итерация эквивалентна домножению матрицы $\bm{A}^{(k)}$ и вектора $\bm{b}^{(k)}$ слева на матрицу $\bm{M}^{(k)}$, имеющую вид:
\begin{equation}
\bm{M}^{(k)}
=
\begin{bmatrix}
1 & 0 & \dots & \dots & \dots & \dots & \dots & 0 \\
0 & \ddots & \ddots & & & & & \vdots \\
\vdots & \ddots & \ddots & \ddots & & & & \vdots \\
\vdots & & 0 & \ddots & \ddots & & & \vdots \\
\vdots & & \vdots & -m_{k+1, k} & \ddots & \ddots & & \vdots \\
\vdots & & \vdots & \vdots & 0 & \ddots & \ddots & \vdots \\
\vdots & & \vdots & \vdots & \vdots & \ddots & \ddots & 0 \\
0 & \dots & 0 & -m_{nk} & 0 & \dots & 0 & 1
\end{bmatrix}.
\end{equation}
Верхняя треугольная матрица, получающаяся в результате обратного хода метода Гаусса, тогда выражается следующим образом:
\begin{equation}
\bm{A}^{(n)} = \bm{M}^{(n-1)} \cdots \bm{M}^{(1)} \bm{A}.
\end{equation}
В контексте LU-разложения мы таким образом получаем матрицу $\bm{U} = \bm{A}^{(n)}$:
\begin{equation}
\bm{U} = \bm{M}^{(n-1)} \cdots \bm{M}^{(1)} \bm{A}.
\end{equation}
Матрица $\bm{L}$ в таком случае может быть получена следующим образом:
\begin{equation}
\begin{split}
&\bm{A} = \left(\bm{M}^{(n-1)} \cdots \bm{M}^{(1)}\right)^{-1} \bm{U} \\
\Longrightarrow &\quad \bm{A} = \left(\bm{M}^{(1)}\right)^{-1} \cdots \left(\bm{M}^{(n-1)}\right)^{-1} \bm{U} \\
\Longrightarrow &\quad \bm{L} = \left(\bm{M}^{(1)}\right)^{-1} \cdots \left(\bm{M}^{(n-1)}\right)^{-1}.
\end{split}
\end{equation}
Явные выражения для обратных матриц могут быть получены, если мы заметим, что операция, обратная $(i) \longrightarrow (i) - \frac{a_{ik}}{a_{kk}^{(1)}} (k), \quad i = k+1, \dots, n$, есть операция $(i) \longrightarrow (i) + \frac{a_{ik}}{a_{kk}^{(1)}} (k), \quad i = k+1, \dots, n$. Тогда обратная матрица $\left(\bm{M}^{(k)}\right)^{-1}$ имеет вид:
\begin{equation}
\left(\bm{M}^{(k)}\right)^{-1}
=
\bm{L}^{(k)}
=
\begin{bmatrix}
1 & 0 & \dots & \dots & \dots & \dots & \dots & 0 \\
0 & \ddots & \ddots & & & & & \vdots \\
\vdots & \ddots & \ddots & \ddots & & & & \vdots \\
\vdots & & 0 & \ddots & \ddots & & & \vdots \\
\vdots & & \vdots & m_{k+1, k} & \ddots & \ddots & & \vdots \\
\vdots & & \vdots & \vdots & 0 & \ddots & \ddots & \vdots \\
\vdots & & \vdots & \vdots & \vdots & \ddots & \ddots & 0 \\
0 & \dots & 0 & m_{nk} & 0 & \dots & 0 & 1
\end{bmatrix}.
\end{equation}
Легко убедиться, что перемножение таких матриц дает:
\begin{equation}
\bm{L}
=
\bm{L}^{(1)} \cdots \bm{L}^{(n-1)}
=
\begin{bmatrix}
1 & 0 & \dots & 0 \\
m_{21} & 1 & \ddots & \vdots \\
\vdots & \ddots & \ddots & 0 \\
m_{n1} & \dots & m_{n, n-1} & 1
\end{bmatrix}.
\end{equation}
В данном выводе LU-разложения мы не предполагали перестановок строк, которые зачастую являются необходимыми. Подобные перестановки можно учесть, используя
матрицу перестановки. Матрица перестановки $\bm{P}$ образуется с помощью перестановки строк единичной матрицы $\bm{E}$. Несложно проверить, что при умножении такой матрицы слева на матрицу $\bm{A}$ происходит та же перестановка строк, что использовалась при построении матрицы $\bm{P}$ из единичной матрицы. При умножении матрицы перестановки справа на матрицу $\bm{A}$ происходит перестановка столбцов. Пусть мы имеем матрицу необходимых перестановок $\bm{P}$. Тогда мы можем найти LU-разложение для матрицы:
\begin{equation}
\bm{P} \bm{A} = \bm{L} \bm{U},
\end{equation}
и найти решение матричного уравнения:
\begin{equation}
\bm{P} \bm{A} \bm{x} = \bm{P} \bm{b}.
\end{equation}
Матрицы с диагональным преобладанием
Матрица $\bm{A}$ обладает свойством диагольного преобладания, если
\begin{equation}
|a_{ii}| \geq \sum_{\substack{j = 1\\ i \neq j}}^n |a_{ij}|, \quad i = 1, \dots, n.
\end{equation}
Для матрицы со строгим диагональным преобладанием верны соответствующие строгие неравенства:
\begin{equation}
|a_{ii}| > \sum_{\substack{j = 1\\ i \neq j}}^n |a_{ij}|, \quad i = 1, \dots, n.
\end{equation}
Матрицы со строгим диагональным преобладанием примечательны тем, что они всегда являются невырожденными. Более того, метод Гаусса может быть использован для них без перестановок строк и столбцов и является вычислительно устойчивым.
Матрица $\bm{A}$ со строгим диагональным преобладанием является невырожденной.
Рассмотрим доказательство от обратного. Пусть матрица $\bm{A}$ является вырожденной. Это означает, что уравнение
\begin{equation}
\bm{A} \bm{x} = \bm{0}
\end{equation}
имеет нетривиальное решение. Пусть $k$ – индекс, для которого верно
\begin{equation}
x_k = \max_{1 \leq j \leq n} |x_j|.
\end{equation}
Тогда имеем:
\begin{equation}
\begin{split}
&\sum_{1 \leq j \leq n} a_{kj} x_j = 0 \\
\Longrightarrow \quad &a_{kk} x_k = -\sum_{\substack{1 \leq j \leq n \\ k \neq j}} a_{kj} x_j.
\end{split}
\end{equation}
В соответствии с неравенством треугольника имеем:
\begin{equation}
|a_{kk}| |x_k| \leq \sum_{\substack{1 \leq j \leq n \\ k \neq j}} |a_{kj}| |x_j|,
\end{equation}
что дает
\begin{equation}
\begin{split}
|a_{kk}| &\leq \sum_{\substack{1 \leq j \leq n \\ k \neq j}} |a_{kj}| \frac{|x_j|}{|x_k|}, \\
&\leq \sum_{\substack{1 \leq j \leq n \\ k \neq j}} |a_{kj}|,
\end{split}
\end{equation}
что противоречит свойству строгого диагонального преобладания.
Положительно определенные матрицы
Матрица $\bm{A}$ называется
положительно определенной, если она симметричная, и верным является неравенство $\bm{x}^T \bm{A} \bm{x} > 0$ для любого вектора $\bm{x} \neq \bm{0}$ подходящей размерности. Легко проверить, что выражение $\bm{x}^T \bm{A} \bm{x}$ имеет вид:
\begin{equation}
\bm{x}^T \bm{A} \bm{x} = \sum_{i = 1}^n \sum_{j = 1}^n a_{ij} x_i x_j.
\end{equation}
Положительно определенные матрицы рядом полезных свойств, которые мы докажем в следующей теореме.
Пусть $\bm{A}$ – положительно определенная матрица. Тогда верно следующее:
- существует $\bm{A}^{-1}$;
- $a_{ii} > 0, i = 1, \dots, n$;
- $\max_{1\leq i \leq n} |a_{ii}| \geq \max_{1\leq k, j\leq n} |a_{kj}|$;
- $a_{ii} a_{jj} > a_{ij}^2, i \neq j$.
Первый пункт доказывается через определение положительно определенной матрицы. Если $\bm{A} \bm{x} = \bm{0}$, то верно и $\bm{x}^T \bm{A} \bm{x} = 0$. По определению это возможно только в том случае, когда $\bm{x} = \bm{0}$. Соответственно $\bm{A} \bm{x} = \bm{0}$ имеет только тривиальное решение, что эквивалентно существованию обратной матрицы для $\bm{A}$.
Второй пункт доказывается через введение ненулевого вектора $\bm{x^{(i)}}$ для $i = 1, \dots, n$:
\begin{equation}
\begin{cases}
x^{(i)}_{i} = 1, \\
x^{(i)}_{j} = 0, \text{~при~} j \neq i.
\end{cases}
\end{equation}
По определению имеем:
\begin{equation}
\begin{split}
&\left(\bm{x^{(i)}}\right)^T \bm{A} \bm{x^{(i)}} > 0 \\
\Longrightarrow &\left(\bm{x^{(i)}}\right)^T \bm{A} \bm{x^{(i)}} = \left(\bm{x^{(i)}}\right)^T \begin{bmatrix}a_{1i} \\ a_{2i} \\ \vdots \\ a_{ni}\end{bmatrix} = a_{ii} > 0.
\end{split}
\end{equation}
Схожим образом доказывается и третий пункт. Определим $\bm{z^{(j, k)}}$ как:
\begin{equation}
\begin{cases}
z^{(j, k)}_{j} = 1, \\
z^{(j, k)}_{k} = -1, \\
z^{(j, k)}_{k} = 0, \text{~при~} i \neq j, i \neq k.
\end{cases}
\end{equation}
Тогда имеем:
\begin{equation}
\left(\bm{z^{(j, k)}}\right)^T \bm{A} \bm{z^{(j, k)}} = \left(\bm{z^{(j, k)}}\right)^T \begin{bmatrix}a_{1j} - a_{1k} \\ a_{2j} - a_{2k} \\ \vdots \\ a_{nj} - a_{nk}\end{bmatrix} = a_{jj} - a_{jk} - (a_{kj} - a_{kk}) > 0.
\end{equation}
В силу симметричности матрицы получаем:
\begin{equation}
a_{kj} < \frac{a_{jj} + a_{kk}}{2}.
\end{equation}
Теперь определим $\bm{x^{(j, k)}}$ как:
\begin{equation}
\begin{cases}
x^{(j, k)}_{j} = 1, \\
x^{(j, k)}_{k} = 1, \\
x^{(j, k)}_{k} = 0, \text{~при~} i \neq j, i \neq k.
\end{cases}
\end{equation}
Аналогичным образом получаем:
\begin{equation}
\left(\bm{x^{(j, k)}}\right)^T \bm{A} \bm{x^{(j, k)}} = \left(\bm{x^{(j, k)}}\right)^T \begin{bmatrix}a_{1j} + a_{1k} \\ a_{2j} + a_{2k} \\ \vdots \\ a_{nj} + a_{nk}\end{bmatrix} = a_{jj} + a_{jk} + a_{kj} + a_{kk} > 0,
\end{equation}
из чего следует:
\begin{equation}
-a_{kj} < \frac{a_{jj} + a_{kk}}{2}.
\end{equation}
Иными словами, мы имеем:
\begin{equation}
|a_{kj}| < \frac{a_{jj} + a_{kk}}{2}.
\end{equation}
Несложно заметить, что:
\begin{equation}
\frac{a_{jj} + a_{kk}}{2} \leq \max_{1 \leq i \leq n} |a_{ii}|,
\end{equation}
что в результате дает:
\begin{equation}
\max_{1 \leq i \leq n} |a_{ii}| \geq \max_{1 \leq k, j \leq n} |a_{kj}|.
\end{equation}
Последний пункт доказывается через модификацию вектора $\bm{x^{(j, k)}}$:
\begin{equation}
\begin{cases}
x^{(j, k)}_{j} = \alpha, \\
x^{(j, k)}_{k} = 1, \\
x^{(j, k)}_{k} = 0, \text{~при~} i \neq j, i \neq k.
\end{cases}
\end{equation}
Тогда получаем:
\begin{equation}
\left(\bm{x^{(j, k)}}\right)^T \bm{A} \bm{x^{(j, k)}} = \left(\bm{x^{(j, k)}}\right)^T \begin{bmatrix}\alpha a_{1j} + a_{1k} \\ \alpha a_{2j} + a_{2k} \\ \vdots \\ \alpha a_{nj} + a_{nk}\end{bmatrix} = \alpha (\alpha a_{jj} + a_{jk}) + \alpha a_{kj} + a_{kk} > 0,
\end{equation}
После упрощений имеем:
\begin{equation}
a_{jj} \alpha^2 + 2 a_{kj} \alpha + a_{kk} > 0.
\end{equation}
Это неравенство выполняется, когда дискриминант меньше нуля:
\begin{equation}
D = 4 a_{kj}^2 - 4 a_{kk} a_{jj} < 0,
\end{equation}
из чего следует:
\begin{equation}
a_{kk} a_{jj} > a_{kj}^2.
\end{equation}
Особенностью положительно определенных матриц является тот факт, что метод Гаусса для них является устойчивым с вычислительной точки зрения и не требует перестановки строк. Более того, для положительно определенных матриц существует два специальных вида разложений – $LDL^T$ и $LL^T$ разложения. Разложение вида $LDL^T$ существует также для симметричных матриц. Для случая положительно определенных матриц мы сформулируем эти разложения в виде двух теорем без доказательства.
Матрица $\bm{A}$ является положительно определенной тогда и только тогда, когда существует разложение $\bm{A} = \bm{L} \bm{D} \bm{L}^T$, где $\bm{L}$ – нижняя треугольная матрица с единицами на диагонали, $\bm{D}$ – диагональная матрица с положительными элементами.
В качестве примера $LDL^T$-разложения рассмотрим положительно определенную матрицу размерности $3\times 3$:
\begin{equation}
\begin{split}
\begin{bmatrix}
a_{11} & a_{21} & a_{31} \\
a_{21} & a_{22} & a_{32} \\
a_{31} & a_{32} & a_{33} \\
\end{bmatrix}
&=
\begin{bmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1 \\
\end{bmatrix}
\begin{bmatrix}
d_1 & 0 & 0 \\
0 & d_2 & 0 \\
0 & 0 & d_3 \\
\end{bmatrix}
\begin{bmatrix}
1 & l_{21} & l_{31} \\
0 & 1 & l_{32} \\
0 & 0 & 1 \\
\end{bmatrix}
\\
&=
\begin{bmatrix}
1 & 0 & 0 \\
l_{21} & 1 & 0 \\
l_{31} & l_{32} & 1 \\
\end{bmatrix}
\begin{bmatrix}
d_1 & d_1 l_{21} & d_1 l_{31} \\
0 & d_2 & d_2 l_{32} \\
0 & 0 & d_3 \\
\end{bmatrix}
\\
&=
\begin{bmatrix}
d_1 & d_1 l_{21} & d_1 l_{31} \\
d_1 l_{21} & d_1 l_{21}^2 + d_2 & d_1 l_{21} l_{31} + d_2 l_{32} \\
d_1 l_{31} & d_1 l_{21} l_{31} + d_2 l_{32} & d_1 l_{31}^2 + d_2 l_{32}^2 + d_3 \\
\end{bmatrix}
\end{split}.
\end{equation}
Тогда неизвестные элемента матриц вычисляются следующим образом:
\begin{align}
d_1 &= a_{11}, \\
l_{21} &= \frac{a_{21}}{d_1}, \\
l_{31} &= \frac{a_{31}}{d_1}, \\
d_2 &= a_{22} - d_1 l_{21}^2, \\
l_{32} &= \frac{1}{d_2}\left(a_{32} - d_1 l_{21} l_{31}\right), \\
d_3 &= a_{33} - d_1 l_{31}^2 - d_2 l_{32}^2.
\end{align}
Матрица $\bm{A}$ является положительно определенной тогда и только тогда, когда существует разложение $\bm{A} = \bm{L} \bm{L}^T$, называемое разложением Холецкого, где $\bm{L}$ – нижняя треугольная матрица с ненулевыми элементами на диагонали.
Из $LDL^T$-разложения можно легко получить разложение Холецкого, если разложить $\bm{D}$ как $\bm{D} = \sqrt{\bm{D}} \sqrt{\bm{D}}$. Тогда имеем:
\begin{equation}
\bm{L} \bm{D} \bm{L}^T = \bm{L} \sqrt{\bm{D}} \sqrt{\bm{D}} \bm{L}^T = \left(\bm{L} \sqrt{\bm{D}}\right) \left(\bm{L} \sqrt{\bm{D}}\right)^T = \tilde{\bm{L}} \tilde{\bm{L}}^T.
\end{equation}
В очередной раз воспользуемся положительно определенной матрицей размерности $3\times 3$ для демонстрации разложения Холецкого:
\begin{equation}
\begin{split}
\begin{bmatrix}
a_{11} & a_{21} & a_{31} \\
a_{21} & a_{22} & a_{32} \\
a_{31} & a_{32} & a_{33} \\
\end{bmatrix}
&=
\begin{bmatrix}
l_{11} & 0 & 0 \\
l_{21} & l_{22} & 0 \\
l_{31} & l_{32} & l_{33} \\
\end{bmatrix}
\begin{bmatrix}
l_{11} & l_{21} & l_{31} \\
0 & l_{22} & l_{32} \\
0 & 0 & l_{33} \\
\end{bmatrix}
\\
&=
\begin{bmatrix}
l_{11}^2 & l_{11} l_{21} & l_{11} l_{31} \\
l_{11} l_{21} & l_{21}^2 + l_{22}^2 & l_{21} l_{31} + l_{22} l_{32} \\
l_{11} l_{31} & l_{21} l_{31} + l_{22} l_{32} & l_{31}^2 + l_{32}^2 + l_{33}^2 \\
\end{bmatrix}
\end{split}.
\end{equation}
Неизвестные элементы матрицы $\bm{L}$ тогда вычисляются следующим образом:
\begin{align}
l_{11} &= \sqrt{a_{11}}, \\
l_{21} &= \frac{a_{21}}{l_{11}}, \\
l_{31} &= \frac{a_{31}}{l_{11}}, \\
l_{22} &= \sqrt{a_{22} - l_{21}^2}, \\
l_{32} &= \frac{1}{l_{22}}\left(a_{32} - l_{21} l_{31}\right), \\
l_{33} &= \sqrt{a_{33} - l_{31}^2 - l_{32}^2}.
\end{align}
Несложно убедиться, что справедливы следующие общие формулы для матрицы произвольной размерности $n \times n$:
\begin{align}
l_{ii} &= \sqrt{a_{ii} - \sum_{j = 1}^{i-1} l_{ij}^2}, \quad i = 1, \dots, n\\
l_{ij} &= \frac{1}{l_{jj}}\left(a_{ij} - \sum_{k = 1}^{j-1} l_{ik} l_{jk}\right), \quad j < i.
\end{align}
Ленточные матрицы
В ряде численных методов результирующая матрица часто является разреженной, то есть обладает преимущественно нулевыми элементами. Одним из важнейших частных случаев разреженных матриц являются ленточные матрицы.
Квадратная матрица $\bm{A} \in \mathbb{R}^{n \times n}$ называется ленточной, если существуют такие $p, q \in \mathbb{Z}$, где $1 < p, q < n$, что $j - i \geq p \Longrightarrow a_{ij} = 0$ и $i - j \geq q \Longrightarrow a_{ij} = 0$. Ширина ленты при этом определяется как $w = p + q - 1$.
В приведенном определении $p$ обозначает число ненулевых диагоналей над главной диагональю (включая ее), в то время как $q$ обозначает число ненулевых диагоналей под главной диагональю (включая ее). Частным случаем ленточной матрицы является трехдиагональная матрица, возникающая, например, при выводе разрещающих уравнений кубического сплайна:
\begin{equation}
\bm{A}
=
\begin{bmatrix}
a_{11} & a_{12} & 0 & \dots & \dots & \dots & 0 \\
a_{21} & a_{22} & a_{23} & \ddots & & & \vdots \\
0 & a_{32} & a_{23} & a_{34} & \ddots & & \vdots \\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\
\vdots & & \ddots & \ddots & \ddots & \ddots & 0 \\
\vdots & & & \ddots & a_{n-1, n-2} & a_{n-1, n-1} & a_{n-1, n} \\
0 & \dots & \dots & \dots & 0 & a_{n, n-1} & a_{nn} \\
\end{bmatrix}.
\end{equation}
Рассмотрим специальный метод решения СЛАУ, имеющих трехдиагональную матрицу –
метод прогонки. Пусть мы имеем следующее матричное уравнение:
\begin{equation}
\begin{bmatrix}
a_{11} & a_{12} & 0 & \dots & \dots & \dots & 0 \\
a_{21} & a_{22} & a_{23} & \ddots & & & \vdots \\
0 & a_{32} & a_{23} & a_{34} & \ddots & & \vdots \\
\vdots & \ddots & \ddots & \ddots & \ddots & \ddots & \vdots \\
\vdots & & \ddots & \ddots & \ddots & \ddots & 0 \\
\vdots & & & \ddots & a_{n-1, n-2} & a_{n-1, n-1} & a_{n-1, n} \\
0 & \dots & \dots & \dots & 0 & a_{n, n-1} & a_{nn} \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
x_3 \\
\vdots \\
\vdots \\
x_{n-1} \\
x_{n} \\
\end{bmatrix}
=
\begin{bmatrix}
b_1 \\
b_2 \\
b_3 \\
\vdots \\
\vdots \\
b_{n-1} \\
b_{n} \\
\end{bmatrix},
\end{equation}
где каждая строка эквивалентна рекуррентному соотношению вида:
\begin{equation}
\label{eq:thomas_original_recur}
a_{i, i-1} x_{i-1} + a_{ii} x_i + a_{i, i+1} x_{i+1} = b_i.
\end{equation}
Очевидно, что методом последовательного исключения можно привести такую СЛАУ к верхней треугольной форме:
\begin{equation}
\begin{bmatrix}
\tilde a_{11} & \tilde a_{12} & 0 & \dots & \dots & 0 \\
0 & \tilde a_{22} & \tilde a_{23} & \ddots & & \vdots \\
\vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\
\vdots & & \ddots & \ddots & \ddots & 0 \\
\vdots & & & \ddots & \tilde a_{n-1, n-1} & \tilde a_{n-1, n} \\
0 & \dots & \dots & \dots & 0 & \tilde a_{nn} \\
\end{bmatrix}
\begin{bmatrix}
x_1 \\
x_2 \\
\vdots \\
\vdots \\
\vdots \\
x_{n-1} \\
x_{n} \\
\end{bmatrix}
=
\begin{bmatrix}
\tilde b_1 \\
\tilde b_2 \\
\vdots \\
\vdots \\
\vdots \\
\tilde b_{n-1} \\
\tilde b_{n} \\
\end{bmatrix}.
\end{equation}
Тогда решение для $x_{i-1}$ может быть выражено через $x_i$:
\begin{equation}
x_{i-1} = \frac{\tilde b_{i-1} - \tilde a_{i-1, i} x_i}{\tilde a_{i-1, i-1}},
\end{equation}
что задает выражение для обратного хода метода Гаусса. Вместо явного вывода неизвестных коэффициентов, мы построим рекуррентное соотношение, позволяющее рекурсивно найти нужные коэффициенты. Для упрощения записи переопределим коэффициенты:
\begin{equation}
\label{eq:thomas_backward_subst}
x_{i-1} = \gamma_i x_i + \beta_i.
\end{equation}
Тогда подстановка в \eqref{eq:thomas_original_recur} дает:
\begin{equation}
\begin{split}
&a_{i, i-1} (\gamma_i x_i + \beta_i) + a_{ii} x_i + a_{i, i+1} x_{i+1} = b_i \\
\Longrightarrow \quad &x_i = \frac{-a_{i, i+1}}{a_{i, i-1} \gamma_i + a_{ii}} x_{i+1} + \frac{b_i - a_{i, i-1} \beta_i}{a_{i, i-1} \gamma_i + a_{ii}}.
\end{split}
\label{eq:thomas_updated_recur}
\end{equation}
Сравнив полученное выражение с \eqref{eq:thomas_backward_subst}, получаем рекуррентные соотношения для определения $\gamma_i$ и $\beta_i$:
\begin{align}
\label{eq:thomas_gamma_recur}
\gamma_{i+1} &= \frac{-a_{i, i+1}}{a_{i, i-1} \gamma_i + a_{ii}}, \\
\label{eq:thomas_beta_recur}
\beta_{i+1} &= \frac{b_i - a_{i, i-1} \beta_i}{a_{i, i-1} \gamma_i + a_{ii}}.
\end{align}
Последовательно вычислив все коэффициенты $\gamma_i$ и $\beta_i$, решение СЛАУ находится с помощью обратного хода метода Гаусса по формуле \eqref{eq:thomas_backward_subst}. Для завершения построения метода прогонки достаточно найти выражения для коэффициентов $\gamma_1, \beta_1$ и неизвестной $x_n$. Рассмотрим уравнение для $x_1$:
\begin{equation}
\begin{split}
&a_{11} x_1 + a_{12} x_2 = b_1 \\
\Longrightarrow \quad &x_1 = - \frac{a_{12}}{a_{11}} x_2 + \frac{b_1}{a_{11}}.
\end{split}
\end{equation}
Тогда сравнение с \eqref{eq:thomas_updated_recur} дает:
\begin{align}
\frac{-a_{12}}{a_{10} \gamma_1 + a_{11}} = -\frac{a_{12}}{a_{11}}, \\
\frac{b_1 - a_{10} \beta_1}{a_{10} \gamma_1 + a_{11}} = \frac{b_{1}}{a_{11}},
\end{align}
из чего следует $\gamma_1 = \beta_1 = 0$ (или, эквивалентно, $a_{10} = 0$). Для определения $x_n$ рассмотрим уравнение, соответствующее последней строке матрицы:
\begin{equation}
\begin{split}
&a_{n, n-1} x_{n-1} + a_{nn} x_n = b_n \\
\Longrightarrow \quad &x_n = - \frac{a_{n, n-1}}{a_{nn}} x_{n-1} + \frac{b_n}{a_{nn}}.
\end{split}
\end{equation}
Подстановка в уравнение \eqref{eq:thomas_backward_subst}, записанное для $i = n$ дает:
\begin{equation}
\begin{split}
&x_{n-1} = \gamma_n x_n + \beta_n \\
\Longrightarrow \quad &x_n = \frac{b_n - a_{n, n-1} \beta_n}{a_{nn} + a_{n, n-1} \gamma_n}.
\end{split}
\end{equation}
Таким образом, используя рекуррентные соотношения \eqref{eq:thomas_backward_subst}, \eqref{eq:thomas_gamma_recur} и \eqref{eq:thomas_beta_recur}, можно найти решение СЛАУ с помощью $O(n)$ операций, что делает метод прогонки предпочтительным перед методом последовательного исключения, требующим $O(n^3)$ операций. Более того, если исходная матрица обладает свойством строгого диагонального преобладания, то метод прогонки является вычислительно устойчивым.