Проекционные методы решения СЛАУ — класс итерационных методов, в которых решается задача проектирования неизвестного вектора на некоторое пространство оптимально относительно другого некоторого пространства.
Рассмотрим СЛАУ
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.
Прямые методы Итерационные методы Общее