記事の背景
機械学習(数理最適化)と制御との関係の理解に向け勉強中です．その一環として，学生時代から触れている最適制御問題を対象とし，その中でも特にモデル予測制御について調べています．
本記事では，モデル予測制御の実装方法の一つで，高速な求解を特徴とするC/GMRES法を理解する前段として，GMRES法について説明したいと思います．
#後日，モデル予測制御に関する記事も投稿予定です．
はじめに
以下のような連立1次方程式を解く(つまり，$x$を求める)ことを考える．
$$Ax = b,\ A \in {\mathbb R}^{n\times n},\ x,b\in {\mathbb R}^n \tag{1}$$
この連立1次方程式を解く方法は，直接法(e.g. LU分解法，Gaussの消去法)と反復法(e.g.ニュートン法，GMRES法)とに大別される．各手法の主な特徴を以下に記す．
- 直接法
- 安定な計算(初期解に依存しない)
- 記憶容量・計算時間を多く要する
- 反復法
- 初期解に依存した解
- 省メモリ
以下では，"大規模連立1次方程式に対する一般化残差法について"を参考にして，反復法の一般化最小残差法(GMRES法)について勉強した内容をまとめる．
GMRES法
まず初めに，GMRES法のアルゴリズムの概観を述べ，続いてその概観を踏まえてアルゴリズムの要所について説明する．
##目次
1.アルゴリズム：概観
2.アルゴリズム：詳細
#1. アルゴリズム:概観
GMRES法とは，方程式(1)の解$x$に対する$k$ステップ目における近似解$x_k$での方程式(1)の残差，すなわち，$r_k = b-Ax_k$を用いて近似解を構成する方法である．GMRES法のアルゴリズムの疑似コードを以下に記す．
※"大規模連立1次方程式に対する一般化残差法について"より引用
※Markdown上で,texのAlgorithmicのようなものの使い方が分からず，引用させて頂きました.python等でプログラムを作成したら，転載しようと思います．
GMRES法の要点：
1. Arnoldi法による$m$($n\geqq m$)次元部分空間(Krylov部分空間)の生成
(=係数行列$A$を利用したベクトル生成 + Gram-Schmidtの正規直交化)
2. 生成された部分空間上で最適な解の探索
GMRES法のメリット：
元々の方程式(1)の求解は$n$次元空間から解を探す問題であったが，GMRES法では$ m $次元部分空間に探索範囲を限ることで計算量を減らすことができる． ($ m $の適切な与え方が分かりませんでした... 残差$r_k$がある閾値以内になる次元で反復計算を打切るのかな...)
※2021/9/17追記：妥当な次元数$m$を与える指標としては，Hankelノルムにより近似精度を測って定めることが多いそうです．また，それに関係して，制御工学のモデル低次元化という分野の一手法としてKrylov部分空間法が適用されているそうです．余力があれば，その記事も書いてみようと思います．
#2. アルゴリズム:詳細
基本的には，GMRES法は上記のステップを経て近似解を取得する方法であるが，恐らく聞きなじみの無い単語がいくつかあるので，本セクションでは，①Krylov部分空間法および②Arnoldi法($\simeq$ Gram-Schmidtの正規直交化法)の説明をして，最後に③GMRES法のアルゴリズムをまとめようと思う．
①Krylov部分空間法
Krylov部分空間とは，行列$A$のべき乗と初期残差ベクトル$r_0$の積によって張られるベクトル空間
$$\mathcal{K}_{\small k} (A;r_0) = \mbox{span} \{r_0,Ar_0,\cdots,A^{k-1}r_0\}\tag{2}$$と定義される．このベクトル空間に含まれるベクトル$z_k \in \mathcal{K}{\small k} (A;r_0)$を用いて近似解を得る解法のことをKrylov部分空間法と呼ぶ．では，定義は分かったが，この考え方がどのように生まれたのかについて，以下で説明する．
$k$ステップ目の方程式(1)の近似解$x_k$を，前ステップでの近似解$x_{k-1}$と残差$r_{k-1}$を用いて以下で定義する．
$$ x_k = x_{k-1}+r_{k-1} \tag{3}$$$$ r_k = b - Ax_k \tag{4}$$(4)を(3)に代入して$x_k$を書き下してみると，$$\left\{ \begin{align*} x_1&=x_0+r_0&\\ x_2&=x_0+2r_0-Ar_0\\ x_3&=x_0+3\underline{r_0}-3\underline{Ar_0}+\underline{A^2r_0} \\ &\vdots
\end{align*} \right. \tag{5}$$
という関係式を得る．ここで，$x_3$の右辺を見ると，初期ベクトル$ x_0 $と下線部の線形結合で解が表現されていることが分かる．
下線部のベクトルが線形独立であれば基底として捉えることができ，近似解$ x_k $は$x_k = x_0 +z_k$ (ここで，$z_k$ は下線部の線形結合，すなわちKrylov部分空間$\mathcal{K}{\small k} (A;r_0)$の元)として表現できる．このようにして，Krylov部分空間法の考え方に至る．
(注意：(5)式に従って解を構築していっても解に収束するわけではない．各基底ベクトルに対して適切な重みづけを行う必要がある．実際，$x_k = x_{k-1}+r_{k-1} = x_{k-1}+b-Ax_{k-1}= (I-A)x_{k-1}+b,$ (ただし，$I$は$n$次単位行列)であり，この数列が収束するには，$(I-A)$の固有値の絶対値がすべて1未満である必要があるし，収束先が方程式の解となるわけではない．(と思う...))
②Arnoldi法
係数行列$A$のべき乗と$r_0$の積により生成される基底が線形独立となる保証はなく，その線形和は数値計算上不安定になってしまう．そこで，Gram-Schmidtの正規直交化法などを用いてKrylov部分空間の正規直交基底を生成する必要がある．
ここで，Arnoldi法の疑似コードを再掲する．
この疑似コードから，Arnoldi法は，係数行列$A$を利用したベクトル生成+Gram-Schmidtの正規直交化をしているだけであることが分かる．ちなみに，Gram-Schmidtの正規直交化法では，あるベクトル列$\{v_1,v_2,\dots,v_k \}$が与えられた時に，それらのベクトルから正規直交基底を作成する．一方で，Arnoldi法では係数行列$A$を利用することで，逐次，ベクトル列を作成していく(Arnoldi法の疑似コードの3ステップ目)．
#長くなってしまったので，Gram-Schmidtの正規直交化法についての解説は別の記事にしたいと思います．
では，この手順にどんな価値があるのかについて考える．Arnoldi法のアルゴリズムを行列表現にしてみる．Arnoldi法のアルゴリズムで新しい基底ベクトルを生成する流れをまとめると，$$ h_{j+1,j} v_{j+1} = A v_j - \sum^{j}_{i=1} h_{i,j} v_i \tag{6}$$となる．左辺を右辺の$\sum$にまとめると，
$$ Av_j = \sum^{j+1}_{i=1} h_{i,j} v_i \tag{7}$$と表現できるので，Arnoldi法によって生成された正規直交基底を集めた$n \times m$行列を$V_m(= (v_1,v_2,\cdots, v_m)\in\mathbb{R}^{n\times m})$と定義すると，式(7)から，
$$ \begin{eqnarray*} Av_j = \sum^{j+1}_{i=1} h_{i,j} v_i &=& (v_1,v_2,\cdots,v_{j+1})\left(h_{1,j},\cdots,h_{j+1,j}\right)^T \\
&=& (v_1,v_2,\cdots,v_{j+1},\cdots,v_{m+1})\left(h_{1,j},\cdots,h_{j+1,j},\cdots,h_{m+1,j}\right)^T \hspace{0.3cm} (\because i>j+1 の時，h_{i,j} = 0 ) \\
&=& V_{m+1}\left(h_{1,j},\cdots,h_{j+1,j},\cdots,h_{m+1,j}\right)^T \end{eqnarray*} \tag{8}$$
と書くことができる．したがって，
$$ AV_m = V_{m+1}\bar{H}_{m} \tag{9}$$
$$ \mbox{ただし，} \bar{H}_{m}= \begin{bmatrix}
h_{1,1} & h_{1,2} & \cdots & h_{1,m-2}& h_{1,m-1} & h_{1,m}\
h_{2,1}& h_{2,2} & \cdots & h_{2,m-2}& h_{2,m-1} & h_{2,m}\
0 & h_{3,2} & \cdots & h_{3,m-2}& h_{3,m-1} & h_{3,m}\
& \ddots & \cdots & \vdots& \vdots & \vdots\
0 & 0 & \cdots & 0& h_{m,m-1} & h_{m,m}\
0 & 0 & \cdots & 0 & 0 & h_{m+1,m}\\
\end{bmatrix} \in \mathbb{R}^{m+1 \times m}\tag{10}$$
なお，式(10)の一番下の行を除いた$m\times m$行列をヘッセンベルグ行列$H_m$と呼ぶ．
③GMRES法
Arnoldi法によって生成されたKrylov部分空間の内で最も良い解を探索したい．そこでまず，評価関数を定める．
$$\| r_m \|_2 = \min_{x_m \in \mathcal{K}{\small k}(A;r_0)} \| b - Ax_m \|_2 \tag{11}$$
を最小にする$x_m$をKrylov部分空間から求める．ここで，係数ベクトル$y=(y_1,y_2,\cdots,y_m)^T \in \mathbb{R}^m$を用いて，Krylov部分空間内のベクトル$x_m$を表現すると，$$x_m = x_0 + \sum_{i=1}^mv_iy_i = x_0+V_my\tag{12}$$となるので，式(11)の最適化問題は，係数ベクトル$y$についての最適化問題に置き換わる．以下のように式変形をすると，
$$ \begin{eqnarray*} b - Ax_m &=& b-A(x_0+V_my) = \underbrace{(b-Ax_0)}_{r_0}+AV_my \\ &=& \|r_0 \|_2v_1 -AV_my = \|r_0 \|_2V_{m+1}e_1 -V_{m+1}\bar{H}_my = V_{m+1}(\|r_0 \|_2e_1 -\bar{H}_my) \end{eqnarray*} \tag{13}$$
となるので，元々の最適化問題(式(11))に代入すると，以下を得る．
$$\begin{eqnarray*} \| r_m \|_2 &=& \underbrace{\min_{x_m \in \mathcal{K}{\small k}(A;r_0)} \| b - Ax_m \|_2}_{(i)} &=& \underbrace{\min_y \left\| V_{m+1}(\|r_0 \|_2e_1 -\bar{H}_my)\right\|_2}_{(ii)} \end{eqnarray*} \tag{11'}$$
式(11')の($i$)ではm次元部分空間の中から最も残差の小さいベクトル$x_m$を持ってくる問題だったが，($ii$)まで変形すると，係数ベクトル$y$を求める問題に落ちていることが分かる．($V_{m+1}$は$y$に依存しないため，$\| \|r_0 \|_2e_1 -\bar{H}_my \|_2$の最小化問題を解けばよい．) そして，ここで重要なのが，$\bar{H}_m$が上三角行列に近い形になっていることから計算が容易となる点である．(後退代入，つまり一番下の行から順次，方程式を解くことで$y_i$の値を求めていき，1,2行目については連立方程式を解く形でできるかと思います．)
GMRES法のアルゴリズムをまとめると以下のようになる．(再掲)
1. Arnoldi法による$m$($n\geqq m$)次元部分空間(Krylov部分空間)の生成
(=係数行列$A$を利用したベクトル生成 + Gram-Schmidtの正規直交化)
2. 生成された部分空間上で最適な解の探索
(=$m$次元部分空間での逆行列計算で係数を最適化)
#Givens回転という方法で，解くべき問題をさらに簡単にする方法が参考文献には書いてあるので，参考文献をご確認ください．
まとめ
連立一次方程式の反復解法であるGMRES法について説明してみました．
執筆に関して，動画や図を加えて，もう少し読みやすい文章にしないといけないなと反省しています...
GMRES法の内容に関しては，残差を用いて定義された数列を基底ベクトルの線形結合と見なすという，私には新鮮な視点だったのでとても面白いと感じました．また，計算を自分の手で行ってみればアルゴリズムの理解もそんなに難しいものでは無いかなと思いましたので，皆様もぜひ，式変形から行ってみることをお勧めいたします．
# 文献を見てみると式が多いため気が引けますが，手を動かし始めたら意外とすんなり頭に入ると思います．