Здесь мы опишем методы, использующие специфику при решении задачи \( Ax = b \). В случае, когда \( A \) — симметричная невырожденная матрица, т.е. \( A = A^T \) и \( \det(A) \ne 0 \), существует разложение вида $$ \begin{equation} \tag{4} A = L D L^T, \end{equation} $$ где \( L \) — нижняя унитреугольная матрица, \( D \) — диагональная матрица. В связи с этим работа связанная с получением разложения :eq:$sles-ldl$, составляет половину от того, что требуется для исключения Гаусса. Когда разложение :eq:$sles-ldl$ получено, решение системы \( Ax = b \) может быть найдено посредством решения систем \( Ly = b \) (прямая подстановка), \( Dz = y \) и \( L^T x = z \).
Разложение (4) может быть найдено при помощи исключения Гаусса, вычисляющего \( A = LU \), c последующим определением \( D \) из уравнения \( U = DL^T \). Тем не менее можно использовать интересный альтернативный алгоритм непосредственного вычисления \( L \) и \( D \).
Допустим, что мы знаем первые \( j-1 \) столбцов матрицы \( L \), диагональные элементы \( d_1, d_2, \ldots, d_{j-1} \) матрицы \( D \) для некоторого \( j \), \( 1 \leq j \leq n \). Чтобы получить способ вычисления \( l_{ij} \), \( i = j+1, j+2, \ldots, n \), и \( d_j \) приравняем \( j \)-е столбцы в уравнении \( A = LDL^T \). В частности, $$ \begin{equation} \tag{5} A(1:j,j) = Lv, \end{equation} $$ где $$ v = DL^T e_j = \begin{bmatrix} d_1 l_{j1}\\ \vdots \\ d_{j-1} l_{jj-1}\\ d_j \end{bmatrix}. $$ Следовательно, компоненты \( v_k \), \( k = 1, 2, \ldots, j-1 \) вектора \( v \) могут быть получены простым масштабированием элементов \( j \)-й строки матрицы \( L \). Формула для \( j \)-й компоненты вектора \( v \) получается из \( j \)-го уравнения системы \( L(1:j, 1:j) v = A(1:j, j) \): $$ v_j = a_{jj} - \sum_{k = 1}^{j-1} l_{ik}v_k, $$ Когда мы знаем \( v \), мы вычисляем \( d_j = v_j \). «Нижняя» половина формулы (5) дает уравнение $$ L(j+1:n, 1:j) v(1:j) = A(j+1:n, j), $$ откуда для вычисления \( j \)-го столбца матрицы \( L \) имеем: $$ L(j+1:n, j) = (A(j+1:n,j) - L(j+1:n, 1:j-1)v(1:j-1))/v_j. $$
Реализация.
Для получения \( LDL^T \)-разложения матрицы \( A \) можем написать функцию (сценарий ld.py):
def ld(A):
"""
Для симметричной матрицы A вычисляет нижнюю треугольную
матрицу L и диагональную матрицу D, такие
что A = LDL^T. Элементы a_{ij} замещаются
на l_{ij}, если i > j, и на d_i, если i = j
"""
n = len(A)
LD = np.array(A, float)
for j in range(n):
v = np.zeros(j+1)
v[:j] = LD[j, :j]*LD[range(j), range(j)]
v[j] = LD[j, j] - np.dot(LD[j, :j], v[:j])
LD[j, j] = v[j]
LD[j+1:, j] = (LD[j+1:, j] - np.dot(LD[j+1:, :j], v[:j]))/v[j]
return LD
В этой реализации мы использовали векторизованные вычисления. Разберем некоторые выражения. Строка
v[:j] = LD[j,:j]*LD[range(j),range(j)]
можно заменить следующим циклом:
for i in range(j):
v[i] = LD[j,i]*LD[i,i]
В нашей программе доступ к j диагональным элементам массива A
осуществляется выражением A[range(j),range(j)].
При вычислении v[j] использовалась функция np.dot, которая
вычисляет скалярное произведение векторов.
Отметим также строку
LD[j+1:,j] = (LD[j+1:,j] - np.dot(LD[j+1:,:j],v[:j]))/v[j]
в которой используется срез L[j+1:,j], т.е. элементы с j+1-го до
последнего в j-ом столбце.
Для решения системы $Ax = b$` с использованием \( LDL^T \)-разложения можно написать следующую функцию
def ld_solve(A, b):
"""
Решает систему Ax = b c использованием LDL^T-разложения
"""
LD = ld(A)
b = np.array(b, float)
for i in range(1, len(b)):
b[i] = b[i] - np.dot(LD[i, :i], b[:i])
b[:] = b[:]/LD[range(len(b)), range(len(b))]
for i in range(len(b)-1, -1, -1):
b[i] = (b[i] - np.dot(LD[i+1:, i], b[i+1:]))
return b
Известно, что в случае симметричной положительно определенной матрицы разложение (4) существует и устойчиво. Тем не менее в этом случае можно использовать другое разложение: $$ \begin{equation} \tag{6} A = G G^T \end{equation} $$ известное как разложение Холецкого, а матрицы \( G \) называются треугольниками Холецкого.
Это легко показать, исходя из существования \( LDL^T \) разложения. Так как для симметричной положительно определенной матрицы существует \( A = LDL^T \) и диагональные элементы матрицы \( D \) положительны, то \( G = L\mathrm{diag}(\sqrt{d_{11}}, \sqrt{d_{22}}, \ldots, \sqrt{d_{nn}}) \).
В разделах Метод исключения Гаусса и Методы решения систем с симметричными матрицами процедуры решения систем алгебраических уравнений были связаны с разложением матрицы коэффициентов \( A \). Методы такого типа называются прямыми методами. Противоположностью прямым методам являются итерационные методы. Эти методы порождают последовательность приближенных решений \( \{ x^{(k)} \} \). При оценивании качества итерационных методов в центре внимания вопрос от том, как быстро сходятся итерации \( x^{(k)} \).
Простейшей итерационной схемой, возможно, являются итерации Якоби. Они определяются для матриц с ненулевыми диагональными элементами. Идею метода можно представить, используя запись \( 3 \times 3 \)-системы \( Ax = b \) в следующем виде: $$ \begin{split} x_1 &= (b_1 - a_{12}x_2 - a_{13}x_3) / a_{11}, \\ x_2 &= (b_2 - a_{21}x_1 - a_{23}x_3) / a_{22}, \\ x_3 &= (b_3 - a_{31}x_1 - a_{32}x_2) / a_{33}. \\ \end{split} $$ Предположим, что \( x^{(k)} \) — какое-то приближение к \( x = A^{-1}b \). Чтобы получить новое приближение \( x^{(k+1)} \), естественно взять: $$ \begin{split} x_1^{(k+1)} &= (b_1 - a_{12}x_2^{(k)} - a_{13}x_3^{(k)}) / a_{11}, \\ x_2^{(k+1)} &= (b_2 - a_{21}x_1^{(k)} - a_{23}x_3^{(k)}) / a_{22}, \\ x_3^{(k+1)} &= (b_3 - a_{31}x_1^{(k)} - a_{32}x_2^{(k)}) / a_{33}. \\ \end{split} $$
Эти формулы и определяют итерации Якоби в случае \( n = 3 \). Для произвольных \( n \) мы имеем $$ \begin{equation} \tag{7} x_i^{(k+1)} = \left( b_i - \sum_{j = 1}^{i-1} a_{ij}x_j^{(k)} - \sum_{j = i+1}^{n} a_{ij}x_j^{(k)} \right)/a_{ii}, \quad i = 1, 2, \ldots, n. \end{equation} $$
Заметим, что в итерациях Якоби при вычислении \( x_i^{(k+1)} \) не используется информация, полученная в самый последний момент. Например, при вычислении \( x_2^{(k+1)} \) используется \( x_1^{(k)} \), хотя уже известна компонента \( x_1^{(k+1)} \). Если мы пересмотрим итерации Якоби с тем, чтобы всегда использовать самые последние оценки для \( x_i \), то получим: $$ \begin{equation} \tag{8} x_i^{(k+1)} = \left( b_i - \sum_{j = 1}^{i-1} a_{ij}x_j^{(k+1)} - \sum_{j = i+1}^{n} a_{ij}x_j^{(k)} \right)/a_{ii}, \quad i = 1, 2, \ldots, n. \end{equation} $$ Так определяется то, что называется итерациями Гаусса — Зейделя.
Для итераций Якоби и Гаусса — Зейделя переход от \( x^{(k)} \) к \( x^{(k+1)} \) в сжатой форме описывается в терминах матриц \( L, D \) и \( U \), определяемых следующим образом: $$ \begin{split} L &= \begin{bmatrix} 0 & 0 &\cdots & \cdots & 0 \\ a_{21} & 0 &\cdots & \cdots & 0 \\ a_{31} & a_{32} & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots &\vdots\\ a_{n1} & a_{n2} & \cdots & a_{nn-1} & 0 \end{bmatrix}, \\ D &= \mathrm{diag}(a_{11}, a_{12}, \ldots, a_{nn}), \\ U &= \begin{bmatrix} 0 & a_{12} &a_{13} & \cdots & a_{1n} \\ 0 & 0 & a_{23} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \ddots &\vdots\\ 0 & 0 & \cdots & 0 & a_{n-1n} \\ 0 & 0 & \cdots & 0 & 0 \end{bmatrix}. \end{split} $$ Шаг Якоби имеет вид \( M_J x^{(k+1)} = N_J x^{(k)} + b \), где \( M_J = D \) и \( N_J = -(L+U) \). С другой стороны, шаг Гаусса — Зейделя определяется как \( M_G x^{(k+1)} = N_G x^{(k)} + b \), где \( M_G = (D+L) \) и \( N_G = -U \).
Процедуры Якоби и Гаусса — Зейделя — это типичные представители большого семейства итерационных методов, имеющих вид $$ \begin{equation} \tag{9} M x^{(k+1)} = N x^{(k)} + b, \end{equation} $$ где \( A = M-N \) — расщепление матрицы \( A \). Для практического применения итераций (9) должна «легко» решаться система с матрицей \( M \). Заметим, что для итераций Якоби и Гаусса — Зейделя матрица \( M \) соответственно диагональная и нижняя треугольная.
Сходятся ли итерации (9) к \( x = A^{-1}b \), зависит от собственных значений матрицы \( M^{-1}N \). Определим спектральный радиус произвольной \( n \times n \)-матрицы \( G \) как $$ \rho(G) = \max \{ |\lambda| :\ \lambda \in \lambda(G) \}, $$ тогда если матрица \( M \) невырожденная и \( \rho(M^{-1}N) < 1 \), то итерации \( x^{(k)} \), определенные согласно \( M^{(k+1)} = N x^{(k)} + b \), сходятся к \( x = A^{-1}b \) при любом начальном векторе \( x^{(0)} \).
Метод Гаусса — Зейделя очень привлекателен в силу своей простоты. К несчастью, если спектральный радиус для \( M_G^{-1}N_G \) близок к единице, то метод может оказаться непозволительно медленным из-за того, что ошибки стремятся к нулю как \( \rho(M_G^{-1}N_G)^k \). Чтобы исправить это, возьмем \( \omega \in \mathbb{R} \) и рассмотрим следующую модификацию шага Гаусса — Зейделя: $$ \begin{equation} \tag{10} x_i^{(k+1)} = \omega \left( b_i - \sum_{j=1}^{i-1} a_{ij}x_{j}^{(k+1)} - \sum_{j={i+1}}^{n} a_{ij}x_{j}^{(k)} \right)/a_{ii} + (1 - \omega) x_i^{(k)}. \end{equation} $$ Так определяется метод последовательной верхней релаксации (SOR — Successive Over Relaxation). В матричных обозначениях шаг SOR выглядит как $$ M_\omega x^{(k+1)} = N_\omega x^{(k)} + \omega b, $$ где \( M_\omega = D + \omega L \) и \( N_\omega = (1-\omega) D - \omega U \). Для небольшого числа специфических задач значения реалксационного параметра \( \omega \), минимизируещего \( \rho(M_\omega^{-1}N_\omega) \), является известным. В более сложных задачах, однако, для того чтобы определить подходящее \( \omega \), может возникнуть необходимость в выполнении весьма трудного анализа собственных значений.
Трудность, связанная с SOR и такого же типа методами, заключается в том, что они зависят от параметров, правильный выбор которых иногда бывает затруднителен. Например, для того чтобы чебышевское ускорение было успешным, нам нужны хорошие оценки для наибольшего и наименьшего собственных значений соответствующей итерационной матрицы \( М^{-1}N \). Если эта матрица не устроена по-особому, то получение их в аналитическом виде, скорее всего, невозможно, а вычисление дорого.
Вывод метода связан с минимизацией функционала: $$ \varphi(x) = \frac{1}{2} x^TAx - x^Tb, $$ где \( b \in \mathbb{R}^n \) и матрица \( A \) предполагается положительно определенной и симметричной. Минимальное значение \( \varphi \) равно \( -b^TA^{-1}b/2 \) и достигается при \( x = A^{-1}b \). Таким образом, минимизация \( \varphi \) и решение системы \( Ax = b \) --- эквивалентные задачи.
Одной из самых простых стратегий минимизации функционала \( \varphi \) является метод наискорейшего спуска. В текущей точке \( x_c \) функция \( \varphi \) убывает наиболее быстро в направлении антиградиента \( \nabla\varphi(x_c) = b - A x_c \). Мы называем \( r_c = b - Ax_c \) невязкой вектора \( x_c \). Если невязка ненулевая, то \( \varphi(x_c + \alpha r_c) < \varphi(x_c) \) для некоторого положительного \( \alpha \) (будем называть этот параметр поправкой). В методе наискорейшего спуска (с точной минимизацией на прямой) мы берем поправку $$ \alpha = \frac{r_c^Tr_c}{r_c^T A r_c}, $$ дающую минимум для \( \varphi(x_c + \alpha r_c) \). Итерационный процесс запишется следующим образом $$ \begin{split} \alpha_k &= \frac{r_{k-1}^Tr_{k-1}}{r_{k-1}^T A r_{k-1}},\\ x_k &= x_{k-1} + \alpha_kr_{k-1},\\ r_k &= b - Ax_k, \quad k = 1, 2, \ldots \end{split} $$ при начальных векторах \( x_0 = 0 \), \( r_0 = b \).
К несчастью, скорость сходимости может быть недопустимо медленной, если число обусловленности \( \kappa(A) = \lambda_1(A)/\lambda_2(A) \) большое. В этом случае линии уровня для \( \varphi \) являются сильно вытянутыми гиперэллипсоидами, а минимизация соответствует поиску самой нижней точки на относительно плоском дне крутого оврага. При наискорейшем спуске мы вынуждены переходить с одной стороны оврага на другую вместо того, чтобы спуститься к его дну. Направления градиента, возникающие при итерациях, являются слишком близкими; это и замедляет продвижение к точке минимума.
Чтобы избежать ловушек при наискорейшем спуске, мы рассмотрим последовательную минимизацию \( \varphi \) вдоль какого-либо множества направлений \( \{p_1, p_2, \ldots\} \), которые не обязаны соответствовать невязкам \( \{r_0, r_1, \ldots\} \). Легко показать, что минимум \( \varphi(x_{k-1} + \alpha p_k) \) по \( \alpha \) дает $$ \alpha_k = \frac{p_k^Tr_{k-1}}{p_k^TAp_k}. $$
Для того, чтобы обеспечить уменьшение функционала \( \varphi \), мы должны потребовать, чтобы \( p_k \) не был ортогонален к \( r_{k-1} \). Проблема состоит в том, как выбирать эти векторы, чтобы гарантировать глобальную сходимость и в то же время обойти ловушки наискорейшего спуска.
Как было сказано выше направления спуска \( p_k \) нужно выбирать так, чтобы они не были ортогональны к невязкам \( r_{k-1} \), т.е. \( p_kr_{k-1} \ne 0 \). Кроме того, метод сопряженных градиентов основан на том, что требуется, чтобы направление \( p_k \) было \( A \)-сопряженным по отношению к \( p_1, p_2, \ldots, p_{k-1} \), т.е. \( p_m^T Ap_k = 0 \) для \( m = 1, 2, \ldots, k-1 \).
Поскольку наша цель — осуществить быстрое сокращение величины невязок, естественно выбирать в качестве \( p_k \) вектор, который ближе всего к \( r_{k-1} \) среди векторов, \( A \)-сопряженных с \( p_1, p_2, \ldots, p_{k-1} \).
Для получения таких направлений спуска и нахождения приближенного решения используется метод сопряженных градиентов. Ниже представлен код функции, реализующий данный алгоритм (файл cg.py)
def cg(A, b, tol, it_max):
it = 0
x = 0
r = np.copy(b)
r_prev = np.copy(b)
rho = np.dot(r, r)
p = np.copy(r)
while (np.sqrt(rho) > tol*np.sqrt(np.dot(b, b)) and it < it_max):
it += 1
if it == 1:
p[:] = r[:]
else:
beta = np.dot(r, r)/np.dot(r_prev, r_prev)
p = r + beta*p
w = np.dot(A, p)
alpha = np.dot(r, r)/np.dot(p, w)
x = x + alpha*p
r_prev[:] = r[:]
r = r - alpha*w
rho = np.dot(r, r)
return x, it