Проекционные методы решения СЛАУ

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

Постановка задачи

Рассмотрим СЛАУ $~ Ax=b,$ где $~A$ - квадратная матрица размерности $~n.$ Пусть $~K$ и $~L$ - два $~m$-мерных подпространства пространства $~R^n.$ Необходимо найти такой вектор $\tilde{x} \in K$, чтобы $b-A\tilde{x} \perp L,$ т.е. выполнялось условие:

$\forall l\in L: (A\tilde{x} ,l)=(b,l),$

называемое условием Петрова-Галёркина.

Если известно начальное приближение $~x_0$, то тогда решение должно проектироваться на аффинное пространство $~ x_0 +K.$ Представим $\tilde{x} =x_0+ \delta$ и обозначим невязку начального приближения как $~r_0 = b- Ax_0.$

$~b-A \tilde{x} =b-A( x_0 + \delta ) =b- Ax_0 - A \delta =(b- Ax_0)- A \delta = r_0 -A \delta .$

Тогда постановку задачи можно сформулировать следующим образом: Необходимо найти такое $\delta \in K,$ чтобы $r_0 -A \delta \perp L ,$ т.е. выполнялось условие:

$~\tilde{x}=x_0+ \delta,\ \delta \in K$

$\forall l\in L:(r_0-A \delta ,l)=0$

Общий подход к построению проекционных методов

Введём матричные базисы в пространствах $~K$ и $~L:$

$~V=[ v_1 , v_2 , ... , v_m ]$ - матрица размера $~ n \times m$ составленная из базисных векторов-столбцов пространства $~K.$ $~W=[ w_1 , w_2 , ... , w_m ]$ - матрица размера $~ n \times m$ составленная из базисных векторов-столбцов пространства $~L.$

Тогда $\delta\ = Vy$ и вектор-решение $\tilde{x}$ может быть записан:

$\tilde{x} = x_0 + Vy,$

где $~y\in R^m$ - вектор коэффициентов.

Тогда выражение $(r_0-A \delta\ ,l)=0$ может быть переписано в виде:

$W^T (r_0-A \delta\ )= \bar{0},$

откуда $~ W^T A Vy = W^T r_0$ и

$~ y = (W^T A V)^{-1} W^T r_0,$

Таким образом решение должно уточняться в соответствии с формулой:

$~ x_1 = x_0 + V(W^T A V)^{-1} W^T r_0,$

Общий вид любого метода проекционного класса:

Делать, пока не найдено решение.

1. Выбираем пару подпространств $~K$ и $~L.$
2. Построение для $~K$ и $~L$ базисов $~V=[ v_1 , v_2 , ... , v_m ]$ и $~W=[ w_1 , w_2 , ... , w_m ].$
3. $~r_0 = b - Ax_0.$
4. $~y = (W^T A V)^{-1} W^T r_0.$
5. $~x_1 = x_0 + Vy.$

Выбор пространств $~K$ и $~L$ и способ построения для них базисов полностью определяет вычислительную схему метода.

Выбор подпространств K и L

Случай одномерных подпространств K и L

В случае когда пространства $~K$ и $~L$ одномерны, их матричные базисы являются векторами: $~V=[v]$ и $~W=[w]$ и выражение $\tilde{x} = x_0 + Vy,$ можно переписать как

$x_{k+1} = x_k + \gamma_k\ v_k,$

где $\gamma_k\$ - неизвестный коэффициент, который легко находится из условия ортогональности $r_k - A( \gamma_k\ v_k ) \perp w_k :$

$(r_k - \gamma_k\ Av_k,w_k) =(r_k , w_k) - \gamma_k\ (Av_k,w_k)= \bar{0},$

откуда $\gamma_k\ = \frac {(r_k , w_k)}{(Av_k,w_k)}.$

Методы с выбором одномерных подпространств $~K$ и $~L$ :

• Метод наискорейшего спуска
• Метод наискорейшего уменьшения невязки
• Метод Гаусса-Зейделя
• Методы ABS-класса

В практических задачах методы использующие одномерные пространства $~K$ и $~L$ обладают достаточно медленной сходимостью.

Методы Крыловского типа

Методы Крыловского типа (или методы подпространства Крылова) - это методы для которых в качестве подпространства $~K$ выбирается подпространство Крылова:

$~ \mathcal{K}_m (r_0,A) = span \mathcal {f} r_0, Ar_0, A^2 r_0, ..., A^{m-1} r_0 \mathcal {g} ,$

где $~r_0=b-Ax_0$ - невязка начального приближения. Различные версии методов подпространства Крылова обуславливаются выбором подпространства $~L.$

С точки зрения теории аппроксимации, приближения $\tilde{x},$ полученные в методах подпространства Крылова имеют форму

$~ A^{-1}b \approx \tilde{x}=x_0+q_{m-1}(A)r_0,$

где $~q_{m-1}$ - полином степени $~m-1.$ Если положить $~x_0=0,$, то

$~ A^{-1}b \approx q_{m-1}(A)b.$

Другими словами, $~A^{-1}b$ аппроксимируется $~q_{m-1}(A)b.$

Хотя выбор подпространства $~L$ и не оказывает влияния на тип полиномиальной аппроксимации, он оказывает существенное влияние на эффективность метода. На сегодняшний день известны 2 способа выбора подпространства $~L,$ дающие наиболее эффективные результаты:

• $~L=K$ и $~L=AK$
• $~L=\mathcal{K}_m (\tilde{r}_0,A^T)$

$~L=K$ и $~L=AK$

Теорема. Если матрица А симметрична и положительно определена, то задача проектирования СЛАУ $~Ax=b$ на любое подпространство $~K$ ортогонально к подпространству $~L=K$ эквивалентна задаче минимизации функционала $\Phi_1(x) = \parallel x- \tilde{x} \parallel_A^2,$ где $\parallel x \parallel_A^2 = (Ax,x).$

Доказательство  

В силу положительной определённости матрицы $~A$ функционал $~ \Phi_1(x)$ достигает своего минимума при $x= \tilde{x}$ и является строго выпуклым. При этом

$\Phi_1(x)=(A(x- \tilde{x}),x- \tilde{x})= (Ax,x)-(A \tilde{x},x)-(Ax, \tilde{x})+(A \tilde{x},\tilde{x})=$ $=(Ax,x)-(Ax,\tilde{x})-(b,x)+(b,\tilde{x})= x^TAx-\tilde{x}^TAx-b^Tx+b^T\tilde{x}.$

В силу симметричности матрицы $~A$ справедливо $~\tilde{x}^TAx=b^TA(A^{-1})x=b^Tx,$ и функционал равен

$\Phi_1(x)=x^TAx-2b^Tx+b^T\tilde{x}.$

По условию теоремы $~K=L,$ следовательно $~V=W.$ Функционал $~ \Phi_1(x)$ является строго выпуклым. Таким образом сформулированная в условии задача минимизации сводится к нахождению

$~y=\arg\min_y \Phi_1(x_0+Vy).$

Рассмотрим эту задачу. В силу выпуклости достаточно найти стационарную точку функционала $~ \Psi(y)=\Phi_1(x_0+Vy),$ т.е. решить систему $~\mathcal {r} \Psi(y)=0.$

$~\Psi(y)=(x_0+Vy)^TA(x_0+Vy)-2b^T(x_0+Vy)+b^T\tilde{x}=$ $~=(x_0^TA-b^T)x_0-b^Tx_0+2(x_0^TA-b^T)Vy+y^T(V^TAV)y+b^T\tilde{x}=$ $~=y^T(V^TAV^T)y-r_0^Tx_0-b^Tx_0-2r_0^TVy+b^T\tilde{x}.$

Градиент этого функционала равен $\mathcal {r} \Psi_1(x)=2V^TAVy-2V^Tr_0.$ Приравнивая его к нулю, получим

$~y=(V^TAV)^{-1}V^Tr_0,$

что в точности совпадает с выражением $~y=(W^TAV)^{-1}W^Tr_0,$ если положить в нём $~V=W.$

Теорема. Если матрица А невырождена, то задача проектирования СЛАУ $~Ax=b$ на любое подпространство $~K$ ортогонально к подпространству $~L=AK$ эквивалентна задаче минимизации функционала $\Phi_2(x) = \parallel r_x \parallel_2^2.$

Доказательство  

Подставив в формулу $~y=(W^T AV)^{-1} W^T r_0$ соотношение для базисов $~W=AV,$ получим:

$~y=( V^T A^T AV)^{-1} V^T A^T r_0.$

Это означает что рассматриваемая ситуация эквивалентна выбору $~L=K$ для симметризованной системы $~A^T Ax=A^T b.$

Учитывая соотношение

$\parallel x- \tilde{x} \parallel_{A^T A}^2 = (A^T A (x- \tilde{x}),(x- \tilde{x}))_2=(A(x- \tilde{x}),A(x- \tilde{x} ))_2=(r_x,r_x)_2= \parallel r_x \parallel_2^2$

и применяя к такой системе предыдущую теорему получим сформулированное в условии утверждение.

$~L=\mathcal{K}_m (\tilde{r}_0,A^T)$

Для построения каждого нового вектора $~v_k$ алгоритм ортогонализации Арнольди требует нахождения $~ (k-1)$ скалярных произведений и столько же операций линейного комбинирования.

Литература

• Saad Y. Iterative methods for sparse linear systems. — 2nd edition. — SIAM Society for Industrial & Applied Mathematics, 2003. — С. 477. — ISBN 0898715342
• Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. — Новосибирск: НГТУ, 2000. — С. 70.
• Голуб Дж., Ван Лоун Ч. Матричные вычисления. — Москва: Мир, 1999.
• Ильин В.П. Методы неполной факторизации для решения линейных систем. — Москва: Физматлит, 1995.

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