Linear systems

Система линейных алгебраических уравнений (СЛАУ) из \(m\) уравнений с \(n\) неизвестными имеет вид

\[\begin{split} \begin{align*} a_{11}x_1 &+ a_{12}x_2+ \ldots +a_{1n}x_n = b_1 \\ a_{21}x_1 &+ a_{22}x_2 + \ldots +a_{2n}x_n = b_2 \\ \ldots &\ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots \\ a_{m1}x_1 &+ a_{m2}x_2 +\ldots + a_{mn}x_n = b_m \\ \end{align*}, \end{split}\]

где коэффициенты \(a_{ij}\) и числа в правой части \(b_i\) считаются известными, а найти надо неизвестные переменные \(x_j\). Записывать систему линейных уравнений в таком громоздком виде не очень удобно, вместо этого чаще используют матричную запись. Коэффициенты \(a_{ij}\) естественным образом образуют матрицу \(\boldsymbol A\in\mathbb R^{m\times n}\), правая часть — вектор \(\boldsymbol b \in \mathbb R^m\), а искомые переменные — вектор \(\boldsymbol x \in \mathbb R^n\). И тогда СЛАУ переписывается в чрезвычайно простом и коротком виде

\[ \boldsymbol{Ax} = \boldsymbol b. \]

In NumPy a linear system of equations with \(m=n\) can be solved by np.linalg.solve:

import numpy as np
from scipy.linalg import pascal
A = pascal(4, kind='lower')
A

array([[1, 0, 0, 0],
       [1, 1, 0, 0],
       [1, 2, 1, 0],
       [1, 3, 3, 1]], dtype=uint64)

b = np.array([1, -1, 1, -1])
np.linalg.solve(A, b)

array([ 1., -2.,  4., -8.])

Gaussian elimintaion

Если матрица СЛАУ диагональная, \(\boldsymbol A = \mathrm{diag}\{a_1, \ldots, a_n\}\), причём все \(a_k\ne 0\), то решить систему не составляет труда: \(x_k = \frac{b_k}{a_k}\). В этом случае решение единственно. Если же \(a_k = 0\) при некотором \(k\), то наличие решения зависит от значения \(b_k\): при \(b_k\ne 0\) система несовместна (решений нет), а при \(b_k = 0\) решений бесконечно много, так как любое значение \(x_k\) удовлетворяет системе уравнений.

Метод Гаусса заключается в преобразовании квадратной матрицы сначала к верхнему треугольному виду, а затем и к диагональному, с помощью элементарных шагов, не меняющих решение СЛАУ.

Transform to upper triangular form

Идея проста: если \(a_{11} \ne 0\), занулим все остальные коэффициенты первого столбца матрицы \(\boldsymbol A\) с помощью основных элементарных преобразований

\[ a_{ij} := a_{ij} - \frac{a_{1j}}{a_{11}}a_{i1}, \quad 2 \leqslant i \leqslant n, \quad 1\leqslant j \leqslant n. \]

В векторном виде их можно записать как

\[ \boldsymbol a_{i}^\top := \boldsymbol a_{i}^\top - \frac{a_{i1}}{a_{11}} \boldsymbol a_1^\top, \quad 2 \leqslant i \leqslant n. \tag{1} \]

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

Единственная проблема данного алгоритма в том, что некоторые диагональные элементы \(a_{ii}\) могут оказаться нулевыми. Если это произошло, то поищем ненулевой элемент \(a_{ij}\) ниже в том же столбце. В случае успеха обменяем \(i\)-ю и \(j\)-ю строки. Такой обмен эквивалентен умножению слева на матрицу перестановки \(\boldsymbol P_{ij}\), и он не изменяет решение системы.

Если же \(a_{ij} = 0\) при всех \(j \geqslant i\), то сделать с \(i\)–м столбцом уже ничего нельзя, так как матрица оказалась вырожденной. Подробнее о вырожденных матрицах см. ниже, пока ограничимся СЛАУ с невырожденными матрицами. Такие матрицы методом Гаусса приводятся к верхнетреугольному виду с ненулевыми элементами на главной диагонали.

Оценим алгоритмическую сложность такого преобразования. Элементарное преобразование (1) требует \(O(n)\) арфметических операций, при занулении элементов \(i\)-го столбца надо произвести \(n-i\) таких преобразований. Перестановка строк требует не более \(O(n)\) операций. Поэтому суммарная сложность равна

\[ O\bigg(\sum\limits_{i=1}^n n(n-i) \bigg) + O(n) = O(n^3). \]

Backward pass

В результате прямого хода метода Гаусса невырожденная квадратная матрица \(\boldsymbol A\) превращается в верхнюю треугольную матрицу

\[\begin{split} \boldsymbol U = \begin{pmatrix} u_{11} & u_{12} & u_{13} & \ldots & u_{1,n-1} & u_{1n}\\ 0 & u_{22} & u_{23} & \ldots & u_{2,n-1} & u_{2n}\\ 0 & 0 & u_{33} & \ldots & u_{3, n-1} & u_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & \ldots & u_{n-1, n-1} & u_{n-1,n} \\ 0 & 0 & 0 & \ldots & 0 & u_{nn} \\ \end{pmatrix} \end{split}\]

с ненулевыми элементами на главной диагонали. Система линейных уравнений \(\boldsymbol{Ax} = \boldsymbol b\) при этом переводится в эквивалентную ей систему \(\boldsymbol{Ux} = \boldsymbol c\). Для решения последней системы используют обратный ход метода Гаусса, который состоит в последовательном исключении неизвестных, начиная с последней строки:

\[ u_{nn}x_n = c_n \Rightarrow x_n = \frac{c_n}{u_{nn}}, \]

\[ u_{ii} x_i + \sum\limits_{j=i+1}^n u_{ij}x_j = c_i \Rightarrow x_i = \frac 1{u_{ii}}\Big(c_i - \sum\limits_{j=i+1}^n u_{ij}x_j\Big), \quad i=n-1, \ldots, 1. \]

Вычисление переменной \(x_i\) требует \(O(n + 1 -i)\) арифметических операций, поэтому общая сложность обратного хода метода Гаусса составляет

\[ O\bigg(\sum\limits_{i=1}^n (n+1-i) \bigg) = O(n^2). \]

Таким образом, суммарная сложность прямого и обратного хода метода Гаусса равна \(O(n^3) + O(n^2) = O(n^3)\).