Gaussian elimination, compact scheme for tridiagonal matrices and its modifications

Main author: А.В.Фролов

1 Compact scheme of Gaussian elimination for tri-diagonal matrices and its modifications

The triangular decomposition of a tri-diagonal matrix is based on the formulas of the same decomposition for general matrices . Tri-diagonal matrices that occur in applied problems are usually diagonally dominant. Consequently, there is no need to use permutations for ensuring the numerical stability of the decomposition. It can be shown, by applying, for instance, the Cauchy–Binet formula to the minors of the product of two matrices :

[math] A \begin{pmatrix} \alpha_{1} & \alpha_{2} & \cdots & \alpha_{k} \\ \beta_{1} & \beta_{2} & \cdots & \beta_{k} \\ \end{pmatrix} = \sum_{\begin{smallmatrix} \gamma_{1} \lt \gamma_{2} \lt \cdots \lt \gamma_{k} \end{smallmatrix}} L \begin{pmatrix} \alpha_{1} & \alpha_{2} & \cdots & \alpha_{k} \\ \gamma_{1} & \gamma_{2} & \cdots & \gamma_{k} \\ \end{pmatrix} U\begin{pmatrix} \gamma_{1} & \gamma_{2} & \cdots & \gamma_{k} \\ \beta_{1} & \beta_{2} & \cdots & \beta_{k} \\ \end{pmatrix} [/math]

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

[math] A = \begin{bmatrix} a_{11} & a_{12} & 0 & \cdots & \cdots & 0 \\ a_{21} & a_{22} & a_{23}& \cdots & \cdots & 0 \\ 0 & a_{32} & a_{33} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & a_{n-1 n-2} & a_{n-1 n-1} & a_{n-1 n} \\ 0 & \cdots & \cdots & 0 & a_{n n-1} & a_{n n} \\ \end{bmatrix} [/math]

[math]LU[/math]-разложение будет содержать две двухдиагональные матрицы:

[math] L = \begin{bmatrix} l_{11} & 0 & 0 & \cdots & \cdots & 0 \\ l_{21} & l_{22} & 0 & \cdots & \cdots & 0 \\ 0 & l_{32} & l_{33} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & l_{n-1 n-2} & l_{n-1 n-1} & 0 \\ 0 & \cdots & \cdots & 0 & l_{n n-1} & l_{n n} \\ \end{bmatrix}, U = \begin{bmatrix} u_{11} & u_{12} & 0 & \cdots & \cdots & 0 \\ 0 & u_{22} & u_{23}& \cdots & \cdots & 0 \\ 0 & 0 & u_{33} & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & u_{n-1 n-1} & u_{n-1 n} \\ 0 & \cdots & \cdots & 0 & 0 & u_{n n} \\ \end{bmatrix} [/math]

Поэтому формулы компактной схемы метода Гаусса и её модификаций существенно упрощаются.

1.1 [math]LU[/math]-разложение

1.1.1 Компактная схема метода Гаусса для трёхдиагональной матрицы

Компактная схема метода Гаусса для трёхдиагональной матрицы вычисляет такое [math]LU[/math]-разложение на две двухдиагональные матрицы, в котором у нижней двухдиагональной матрицы [math]L[/math] на диагонали стоят только единицы:

[math] L = \begin{bmatrix} 1 & 0 & 0 & \cdots & \cdots & 0 \\ l_{21} & 1 & 0 & \cdots & \cdots & 0 \\ 0 & l_{32} & 1 & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & l_{n-1 n-2} & 1 & 0 \\ 0 & \cdots & \cdots & 0 & l_{n n-1} & 1 \\ \end{bmatrix} [/math]

Формулы метода следующие:

[math]u_{11} = a_{11}[/math]

[math]u_{i i+1} = a_{i i+1}[/math], [math]\quad i \in [1, n-1][/math]

[math]l_{i+1 i} = a_{i+1 i} / u_{i i} [/math], [math]\quad i \in [1, n-1] [/math]

[math]u_{ii} = a_{ii} - l_{i i-1} u_{i-1 i}[/math], [math]\quad i \in [2, n] [/math]

Подробно её свойства описаны на соответствующей странице. В неизменённом виде не подлежит распараллеливанию, алгоритм чисто последовательный.Существует и блочная версия разложения - для блочно-трёхдиагональных матриц. Компактная схема метода Гаусса легко подвергается модификации, в которой не в матрице [math]L[/math], а в матрице [math]U[/math] на диагонали только единицы. Это используется, например, в качестве части метода прогонки.

1.1.2 Алгоритм Стоуна для [math]LU[/math]-разложения трёхдиагональной матрицы

Алгоритм сдваивания Стоуна для LU-разложения трёхдиагональной матрицы является частью алгоритма Стоуна для решения трёхдиагональных СЛАУ и является первым параллельным алгоритмом [math]LU[/math]-разложения трёхдиагональной матрицы. Математическая его основа - линейная рекуррентная формула для главных миноров раскладываемой матрицы. После записи этой рекуррентной формулы в матричном виде оказывается, что можно воспользоваться ассоциативностью операции перемножения матриц. Стоун использовал для распараллеливания схему сдваивания. Однако область устойчивости этого метода гораздо уже, чем у компактной схемы метода Гаусса, несмотря на то, что искомые матрицы - те же самые. Блочная схема метода для блочно-трёхдиагональных матриц не разрабатывалась, хотя теоретически и может существовать, из-за неустойчивости схемы.

Подробно свойства алгоритма Стоуна описаны на соответствующей странице.

1.1.3 Последовательно-параллельный алгоритм для [math]LU[/math]-разложения трёхдиагональной матрицы

Последовательно-параллельный алгоритм для LU-разложения трёхдиагональной матрицы разработан для распараллеливания нахождения того же [math]LU[/math]-разложения трёхдиагональной матрицы, что получается из компактной схемы метода Гаусса. Как и метод Стоуна, использует ассоциативность умножения матриц, но использует также нормировку в последовательных ветвях вычислений, что даёт схеме большую область устойчивости, чем у схемы Стоуна. Для блочно-трёхдиагональных матриц существует блочная версия метода, для которой, однако, необходима невырожденность блоков не только на главной, но и на одной из побочных диагоналей матрицы.

1.2 [math]LDU[/math]-разложение

В случае, если существует [math]LU[/math]-разложение на две двухдиагональные матрицы, в котором нижняя двухдиагональная матрица [math]L[/math] имеет на главной диагонали только единицы, легко можно найти и [math]LDU[/math]-разложение, в котором и верхняя двухдиагональная матрица также имеет на главной диагонали только единицы:

[math] U = \begin{bmatrix} 1 & u_{12} & 0 & \cdots & \cdots & 0 \\ 0 & 1 & u_{23}& \cdots & \cdots & 0 \\ 0 & 0 & 1 & \cdots & \cdots & 0 \\ \vdots & \vdots & \ddots & \ddots & \ddots & 0 \\ 0 & \cdots & \cdots & 0 & 1 & u_{n-1 n} \\ 0 & \cdots & \cdots & 0 & 0 & 1 \\ \end{bmatrix} [/math]

Формулы видоизменённой компактной схемы метода Гаусса для нахождения этого разложения выглядят следующим образом:

[math]d_{11} = a_{11} [/math]

[math]u_{i i+1} = a_{i i+1} / d_{i i}, \quad i \in [1, n-1] [/math]

[math]l_{i+1 i} = a_{i+1 i} / d_{i i} , \quad i \in [1, n-1] [/math]

[math]d_{ii} = a_{ii} - l_{i i-1} a_{i-1 i}, \quad i \in [2, n] [/math]

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

Получаемое разложение легче использовать для решения трёхдиагональных СЛАУ, поскольку в обоих подстановках деление не будет перемежаться с операциями умножения-сложения/вычитания.

2 Компактная схема треугольного разложения для трёхдиагональной эрмитовой матрицы и её модификации

Эрмитовость трёхдиагональной матрицы позволяет вычислять такое её разложение, которое можно хранить в меньшей памяти с учетом симметрии разложения.

2.1 [math]LL^*[/math]-разложение

Для вычисления [math]LL^*[/math]-разложения матрицы логично использовать метод Холецкого, формулы которого после коррекции будут адаптированы под трёхдиагональность исходной матрицы и двухдиагональность результата:

[math]l_{11} = \sqrt{a_{11}}, [/math]

[math]l_{i+1 i} = \frac{a_{i+1 i}}{l_{ii}}, \quad i \in [1, n-1], [/math]

[math]l_{ii} = \sqrt{a_{ii} - |l_{i i-1}|^2}, \quad i \in [2, n]. [/math]

Аналогичные изменения можно провести в методах Стоуна и последовательно-параллельного [math]LU[/math]-разложения. Однако на практике использование схемы с квадратными корнями для трёхдиагональных матриц не нужно, поскольку есть более быстрое [math]LDL^*[/math]-разложение.

2.2 [math]LDL^*[/math]-разложение

[math]LDL^*[/math]-разложение для трёхдиагональной эрмитовой матрицы - самое экономное из последовательных её разложений. Формулы разложения имеют вид:

[math]d_{11} = a_{11} [/math]

[math]l_{i+1 i} = a_{i+1 i} / d_{i i} , \quad i \in [1, n-1] [/math]

[math]d_{ii} = a_{ii} - l_{i i-1} a_{i i-1}^*, \quad i \in [2, n] [/math]

Аналогично, можно видоизменить и формулы метода Стоуна и последовательно-параллельного метода разложения. Как и для [math]LDU[/math]-разложения, при решении СЛАУ с трёхдиагональной матрицей с помощью [math]LDL^*[/math]-разложения легче решаются полученные двухдиагональные СЛАУ.

3 Разложения для трёхдиагональных матриц специального вида

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

Разложения для трёхдиагональных матриц специального вида можно найти на специальной отдельной странице.

4 Использование разложений

Обычно разложения трёхдиагональных матриц в произведение двухдиагональных (и, возможно, ещё одной диагональной) используют для решения СЛАУ с трёхдиагональной матрицей, полученной из какой-либо физической модели. Поэтому зачастую полезность того или иного метода разложения следует сравнивать не только с другими методами разложения трёхдиагональных матриц, но и с различными методами решения СЛАУ с такими матрицами, в том числе и не использующими такие разложения - например, такими, как метод циклической редукции или метод окаймления.