Проекционные методы решения СЛАУ — класс итерационных методов, в которых решается задача проектирования неизвестного вектора на некоторое пространство оптимально относительно другого некоторого пространства.
Рассмотрим СЛАУ A x = b , {\displaystyle Ax=b,} где A {\displaystyle A} - квадратная матрица размерности n . {\displaystyle n.} Пусть K {\displaystyle K} и L {\displaystyle L} - два m {\displaystyle m} -мерных подпространства пространства R n . {\displaystyle R^{n}.} Необходимо найти такой вектор x ~ ∈ K {\displaystyle {\tilde {x}}\in K} , чтобы b − A x ~ ⊥ L , {\displaystyle b-A{\tilde {x}}\perp L,} т.е. выполнялось условие:
∀ l ∈ L : ( A x ~ , l ) = ( b , l ) , {\displaystyle \forall l\in L:(A{\tilde {x}},l)=(b,l),}
называемое условием Петрова-Галёркина.
Если известно начальное приближение x 0 {\displaystyle x_{0}} , то тогда решение должно проектироваться на аффинное пространство x 0 + K . {\displaystyle x_{0}+K.} Представим x ~ = x 0 + δ {\displaystyle {\tilde {x}}=x_{0}+\delta } и обозначим невязку начального приближения как r 0 = b − A x 0 . {\displaystyle r_{0}=b-Ax_{0}.}
b − A x ~ = b − A ( x 0 + δ ) = b − A x 0 − A δ = ( b − A x 0 ) − A δ = r 0 − A δ . {\displaystyle b-A{\tilde {x}}=b-A(x_{0}+\delta )=b-Ax_{0}-A\delta =(b-Ax_{0})-A\delta =r_{0}-A\delta .}
Тогда постановку задачи можно сформулировать следующим образом: Необходимо найти такое δ ∈ K , {\displaystyle \delta \in K,} чтобы r 0 − A δ ⊥ L , {\displaystyle r_{0}-A\delta \perp L,} т.е. выполнялось условие:
x ~ = x 0 + δ , δ ∈ K {\displaystyle {\tilde {x}}=x_{0}+\delta ,\ \delta \in K}
∀ l ∈ L : ( r 0 − A δ , l ) = 0 {\displaystyle \forall l\in L:(r_{0}-A\delta ,l)=0}
Введём матричные базисы в пространствах K {\displaystyle K} и L : {\displaystyle L:}
V = [ v 1 , v 2 , . . . , v m ] {\displaystyle V=[v_{1},v_{2},...,v_{m}]} - матрица размера n × m {\displaystyle n\times m} составленная из базисных векторов-столбцов пространства K . {\displaystyle K.} W = [ w 1 , w 2 , . . . , w m ] {\displaystyle W=[w_{1},w_{2},...,w_{m}]} - матрица размера n × m {\displaystyle n\times m} составленная из базисных векторов-столбцов пространства L . {\displaystyle L.}
Тогда δ = V y {\displaystyle \delta \ =Vy} и вектор-решение x ~ {\displaystyle {\tilde {x}}} может быть записан:
x ~ = x 0 + V y , {\displaystyle {\tilde {x}}=x_{0}+Vy,}
где y ∈ R m {\displaystyle y\in R^{m}} - вектор коэффициентов.
Тогда выражение ( r 0 − A δ , l ) = 0 {\displaystyle (r_{0}-A\delta \ ,l)=0} может быть переписано в виде:
W T ( r 0 − A δ ) = 0 ¯ , {\displaystyle W^{T}(r_{0}-A\delta \ )={\bar {0}},}
откуда W T A V y = W T r 0 {\displaystyle W^{T}AVy=W^{T}r_{0}} и
y = ( W T A V ) − 1 W T r 0 , {\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0},}
Таким образом решение должно уточняться в соответствии с формулой:
x 1 = x 0 + V ( W T A V ) − 1 W T r 0 , {\displaystyle x_{1}=x_{0}+V(W^{T}AV)^{-1}W^{T}r_{0},}
Общий вид любого метода проекционного класса:
Делать, пока не найдено решение.
Выбираем пару подпространств K {\displaystyle K} и L . {\displaystyle L.} Построение для K {\displaystyle K} и L {\displaystyle L} базисов V = [ v 1 , v 2 , . . . , v m ] {\displaystyle V=[v_{1},v_{2},...,v_{m}]} и W = [ w 1 , w 2 , . . . , w m ] . {\displaystyle W=[w_{1},w_{2},...,w_{m}].} r 0 = b − A x 0 . {\displaystyle r_{0}=b-Ax_{0}.} y = ( W T A V ) − 1 W T r 0 . {\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0}.} x 1 = x 0 + V y . {\displaystyle x_{1}=x_{0}+Vy.} Выбор пространств K {\displaystyle K} и L {\displaystyle L} и способ построения для них базисов полностью определяет вычислительную схему метода.
В случае когда пространства K {\displaystyle K} и L {\displaystyle L} одномерны, их матричные базисы являются векторами: V = [ v ] {\displaystyle V=[v]} и W = [ w ] {\displaystyle W=[w]} и выражение x ~ = x 0 + V y , {\displaystyle {\tilde {x}}=x_{0}+Vy,} можно переписать как
x k + 1 = x k + γ k v k , {\displaystyle x_{k+1}=x_{k}+\gamma _{k}\ v_{k},}
где γ k {\displaystyle \gamma _{k}\ } - неизвестный коэффициент, который легко находится из условия ортогональности r k − A ( γ k v k ) ⊥ w k : {\displaystyle r_{k}-A(\gamma _{k}\ v_{k})\perp w_{k}:}
( r k − γ k A v k , w k ) = ( r k , w k ) − γ k ( A v k , w k ) = 0 ¯ , {\displaystyle (r_{k}-\gamma _{k}\ Av_{k},w_{k})=(r_{k},w_{k})-\gamma _{k}\ (Av_{k},w_{k})={\bar {0}},}
откуда γ k = ( r k , w k ) ( A v k , w k ) . {\displaystyle \gamma _{k}\ ={\frac {(r_{k},w_{k})}{(Av_{k},w_{k})}}.}
Методы с выбором одномерных подпространств K {\displaystyle K} и L {\displaystyle L} :
В практических задачах методы использующие одномерные пространства K {\displaystyle K} и L {\displaystyle L} обладают достаточно медленной сходимостью.
Методы Крыловского типа (или методы подпространства Крылова ) - это методы для которых в качестве подпространства K {\displaystyle K} выбирается подпространство Крылова:
K m ( r 0 , A ) = s p a n { r 0 , A r 0 , A 2 r 0 , . . . , A m − 1 r 0 } , {\displaystyle {\mathcal {K}}_{m}(r_{0},A)=span\{r_{0},Ar_{0},A^{2}r_{0},...,A^{m-1}r_{0}\},}
где r 0 = b − A x 0 {\displaystyle r_{0}=b-Ax_{0}} - невязка начального приближения. Различные версии методов подпространства Крылова обуславливаются выбором подпространства L . {\displaystyle L.}
С точки зрения теории аппроксимации , приближения x ~ , {\displaystyle {\tilde {x}},} полученные в методах подпространства Крылова имеют форму
A − 1 b ≈ x ~ = x 0 + q m − 1 ( A ) r 0 , {\displaystyle A^{-1}b\approx {\tilde {x}}=x_{0}+q_{m-1}(A)r_{0},}
где q m − 1 {\displaystyle q_{m-1}} - полином степени m − 1. {\displaystyle m-1.} Если положить x 0 = 0 , {\displaystyle x_{0}=0,} , то
A − 1 b ≈ q m − 1 ( A ) b . {\displaystyle A^{-1}b\approx q_{m-1}(A)b.}
Другими словами, A − 1 b {\displaystyle A^{-1}b} аппроксимируется q m − 1 ( A ) b . {\displaystyle q_{m-1}(A)b.}
Хотя выбор подпространства L {\displaystyle L} и не оказывает влияния на тип полиномиальной аппроксимации, он оказывает существенное влияние на эффективность метода. На сегодняшний день известны 2 способа выбора подпространства L , {\displaystyle L,} дающие наиболее эффективные результаты:
L = K {\displaystyle L=K} и L = A K {\displaystyle L=AK} L = K m ( r ~ 0 , A T ) {\displaystyle L={\mathcal {K}}_{m}({\tilde {r}}_{0},A^{T})} L = K {\displaystyle L=K} и L = A K {\displaystyle L=AK} В силу положительной определённости матрицы A {\displaystyle A} функционал Φ 1 ( x ) {\displaystyle \Phi _{1}(x)} достигает своего минимума при x = x ~ {\displaystyle x={\tilde {x}}} и является строго выпуклым. При этом
Φ 1 ( x ) = ( A ( x − x ~ ) , x − x ~ ) = ( A x , x ) − ( A x ~ , x ) − ( A x , x ~ ) + ( A x ~ , 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}})=} = ( A x , x ) − ( A x , x ~ ) − ( b , x ) + ( b , x ~ ) = x T A x − x ~ T A x − b T x + b T 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}}.} В силу симметричности матрицы A {\displaystyle A} справедливо x ~ T A x = b T A ( A − 1 ) x = b T x , {\displaystyle {\tilde {x}}^{T}Ax=b^{T}A(A^{-1})x=b^{T}x,} и функционал равен
Φ 1 ( x ) = x T A x − 2 b T x + b T x ~ . {\displaystyle \Phi _{1}(x)=x^{T}Ax-2b^{T}x+b^{T}{\tilde {x}}.} По условию теоремы K = L , {\displaystyle K=L,} следовательно V = W . {\displaystyle V=W.} Функционал Φ 1 ( x ) {\displaystyle \Phi _{1}(x)} является строго выпуклым. Таким образом сформулированная в условии задача минимизации сводится к нахождению
y = arg min y Φ 1 ( x 0 + V y ) . {\displaystyle y=\arg \min _{y}\Phi _{1}(x_{0}+Vy).} Рассмотрим эту задачу. В силу выпуклости достаточно найти стационарную точку функционала Ψ ( y ) = Φ 1 ( x 0 + V y ) , {\displaystyle \Psi (y)=\Phi _{1}(x_{0}+Vy),} т.е. решить систему r Ψ ( y ) = 0. {\displaystyle {\mathcal {r}}\Psi (y)=0.}
Ψ ( y ) = ( x 0 + V y ) T A ( x 0 + V y ) − 2 b T ( x 0 + V y ) + b T x ~ = {\displaystyle \Psi (y)=(x_{0}+Vy)^{T}A(x_{0}+Vy)-2b^{T}(x_{0}+Vy)+b^{T}{\tilde {x}}=} = ( x 0 T A − b T ) x 0 − b T x 0 + 2 ( x 0 T A − b T ) V y + y T ( V T A V ) y + b T 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}}=} = y T ( V T A V ) y − r 0 T x 0 − b T x 0 − 2 r 0 T V y + b T 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}}.} Градиент этого функционала равен r Ψ 1 ( y ) = 2 V T A V y − 2 V T r 0 . {\displaystyle {\mathcal {r}}\Psi _{1}(y)=2V^{T}AVy-2V^{T}r_{0}.} Приравнивая его к нулю, получим
y = ( V T A V ) − 1 V T r 0 , {\displaystyle y=(V^{T}AV)^{-1}V^{T}r_{0},} что в точности совпадает с выражением y = ( W T A V ) − 1 W T r 0 , {\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0},} если положить в нём V = W . {\displaystyle V=W.}
Подставив в формулу y = ( W T A V ) − 1 W T r 0 {\displaystyle y=(W^{T}AV)^{-1}W^{T}r_{0}} соотношение для базисов W = A V , {\displaystyle W=AV,} получим:
y = ( V T A T A V ) − 1 V T A T r 0 . {\displaystyle y=(V^{T}A^{T}AV)^{-1}V^{T}A^{T}r_{0}.} Это означает что рассматриваемая ситуация эквивалентна выбору L = K {\displaystyle L=K} для симметризованной системы A T A x = A T b . {\displaystyle A^{T}Ax=A^{T}b.}
Учитывая соотношение
∥ x − x ~ ∥ A T A 2 = ( A T A ( x − x ~ ) , ( x − x ~ ) ) 2 = ( A ( x − x ~ ) , A ( x − x ~ ) ) 2 = ( r x , r x ) 2 =∥ r x ∥ 2 2 {\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}} и применяя к такой системе предыдущую теорему получим сформулированное в условии утверждение.
L = K m ( r ~ 0 , A T ) {\displaystyle L={\mathcal {K}}_{m}({\tilde {r}}_{0},A^{T})} Для построения каждого нового вектора v k {\displaystyle v_{k}} алгоритм ортогонализации Арнольди требует нахождения ( k − 1 ) {\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. Прямые методы Итерационные методы Общее