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

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

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

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

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

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

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

${\displaystyle b-A{\tilde {x}}=b-A(x_{0}+\delta )=b-Ax_{0}-A\delta =(b-Ax_{0})-A\delta =r_{0}-A\delta .}$

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

${\displaystyle {\tilde {x}}=x_{0}+\delta ,\ \delta \in K}$

${\displaystyle \forall l\in L:(r_{0}-A\delta ,l)=0}$

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

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

${\displaystyle V=[v_{1},v_{2},...,v_{m}]}$ - матрица размера ${\displaystyle n\times m}$ составленная из базисных векторов-столбцов пространства ${\displaystyle K.}$ ${\displaystyle W=[w_{1},w_{2},...,w_{m}]}$ - матрица размера ${\displaystyle n\times m}$ составленная из базисных векторов-столбцов пространства ${\displaystyle L.}$

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

${\displaystyle {\tilde {x}}=x_{0}+Vy,}$

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

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

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

откуда ${\displaystyle W^{T}AVy=W^{T}r_{0}}$ и

${\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0},}$

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

${\displaystyle x_{1}=x_{0}+V(W^{T}AV)^{-1}W^{T}r_{0},}$

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

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

1. Выбираем пару подпространств ${\displaystyle K}$ и ${\displaystyle L.}$
2. Построение для ${\displaystyle K}$ и ${\displaystyle L}$ базисов ${\displaystyle V=[v_{1},v_{2},...,v_{m}]}$ и ${\displaystyle W=[w_{1},w_{2},...,w_{m}].}$
3. ${\displaystyle r_{0}=b-Ax_{0}.}$
4. ${\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0}.}$
5. ${\displaystyle x_{1}=x_{0}+Vy.}$

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

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

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

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

${\displaystyle x_{k+1}=x_{k}+\gamma _{k}\ v_{k},}$

где ${\displaystyle \gamma _{k}\ }$ - неизвестный коэффициент, который легко находится из условия ортогональности ${\displaystyle r_{k}-A(\gamma _{k}\ v_{k})\perp w_{k}:}$

${\displaystyle (r_{k}-\gamma _{k}\ Av_{k},w_{k})=(r_{k},w_{k})-\gamma _{k}\ (Av_{k},w_{k})={\bar {0}},}$

откуда ${\displaystyle \gamma _{k}\ ={\frac {(r_{k},w_{k})}{(Av_{k},w_{k})}}.}$

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

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

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

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

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

${\displaystyle {\mathcal {K}}_{m}(r_{0},A)=span\{r_{0},Ar_{0},A^{2}r_{0},...,A^{m-1}r_{0}\},}$

где ${\displaystyle r_{0}=b-Ax_{0}}$ - невязка начального приближения. Различные версии методов подпространства Крылова обуславливаются выбором подпространства ${\displaystyle L.}$

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

${\displaystyle A^{-1}b\approx {\tilde {x}}=x_{0}+q_{m-1}(A)r_{0},}$

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

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

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

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

• ${\displaystyle L=K}$ и ${\displaystyle L=AK}$
• ${\displaystyle L={\mathcal {K}}_{m}({\tilde {r}}_{0},A^{T})}$

${\displaystyle L=K}$ и ${\displaystyle L=AK}$

Теорема

Если матрица А симметрична и положительно определена, то задача проектирования СЛАУ ${\displaystyle Ax=b}$ на любое подпространство ${\displaystyle K}$ ортогонально к подпространству ${\displaystyle L=K}$ эквивалентна задаче минимизации функционала

${\displaystyle \Phi _{1}(x)=\parallel x-{\tilde {x}}\parallel _{A}^{2},}$

где ${\displaystyle \parallel x\parallel _{A}^{2}=(Ax,x).}$

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

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

${\displaystyle \Phi _{1}(x)=(A(x-{\tilde {x}}),x-{\tilde {x}})=(Ax,x)-(A{\tilde {x}},x)-(Ax,{\tilde {x}})+(A{\tilde {x}},{\tilde {x}})=}$ ${\displaystyle =(Ax,x)-(Ax,{\tilde {x}})-(b,x)+(b,{\tilde {x}})=x^{T}Ax-{\tilde {x}}^{T}Ax-b^{T}x+b^{T}{\tilde {x}}.}$

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

${\displaystyle \Phi _{1}(x)=x^{T}Ax-2b^{T}x+b^{T}{\tilde {x}}.}$

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

${\displaystyle y=\arg \min _{y}\Phi _{1}(x_{0}+Vy).}$

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

${\displaystyle \Psi (y)=(x_{0}+Vy)^{T}A(x_{0}+Vy)-2b^{T}(x_{0}+Vy)+b^{T}{\tilde {x}}=}$ ${\displaystyle =(x_{0}^{T}A-b^{T})x_{0}-b^{T}x_{0}+2(x_{0}^{T}A-b^{T})Vy+y^{T}(V^{T}AV)y+b^{T}{\tilde {x}}=}$ ${\displaystyle =y^{T}(V^{T}AV)y-r_{0}^{T}x_{0}-b^{T}x_{0}-2r_{0}^{T}Vy+b^{T}{\tilde {x}}.}$

Градиент этого функционала равен ${\displaystyle {\mathcal {r}}\Psi _{1}(y)=2V^{T}AVy-2V^{T}r_{0}.}$ Приравнивая его к нулю, получим

${\displaystyle y=(V^{T}AV)^{-1}V^{T}r_{0},}$

что в точности совпадает с выражением ${\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0},}$ если положить в нём ${\displaystyle V=W.}$

Теорема

Если матрица А невырождена, то задача проектирования СЛАУ ${\displaystyle Ax=b}$ на любое подпространство ${\displaystyle K}$ ортогонально к подпространству ${\displaystyle L=AK}$ эквивалентна задаче минимизации функционала

${\displaystyle \Phi _{2}(x)=\parallel r_{x}\parallel _{2}^{2}.}$

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

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

${\displaystyle y=(V^{T}A^{T}AV)^{-1}V^{T}A^{T}r_{0}.}$

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

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

${\displaystyle \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}}$

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

${\displaystyle L={\mathcal {K}}_{m}({\tilde {r}}_{0},A^{T})}$

Для построения каждого нового вектора ${\displaystyle v_{k}}$ алгоритм ортогонализации Арнольди требует нахождения ${\displaystyle (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.