はじめに
線形方程式 $A\boldsymbol{x}=\boldsymbol{b}$ は工学の様々な分野で登場するため、これをコンピュータで計算する手法は盛んに研究されてきました。そういった手法の一種として頻繁に用いられるのが共役勾配法(CG法)です。この記事では、CG法を理解する前段階として、変分原理、反復法、最急降下法の解説を行います。
CG法の解説は次の記事で行っています。
目次
変分原理
数学的なざっくり解説
反復法
最急降下法
参考資料
変分原理
次のような形で表される、ベクトル $\boldsymbol{x}$ に関する線形方程式を考えます。
A\boldsymbol{x}=\boldsymbol{b} \tag{1}
このような線形方程式は工学の様々な場面で登場しますが、これを数値的に解くにはどうしたらいいでしょう?まず、次のように変形します。
A\boldsymbol{x}-\boldsymbol{b}=0 \tag{2}
右辺がゼロになってくれた方がいろいろやりやすいからです。次に、関数 $f(\boldsymbol{x})$ を次のように定義します。
f(\boldsymbol{x}) \equiv \frac{1}{2}(\boldsymbol{x}, A\boldsymbol{x})-(\boldsymbol{b}, \boldsymbol{x}) \tag{3}
なお、 $(\boldsymbol{a}, \boldsymbol{b})$ は $\boldsymbol{a}$ と$\boldsymbol{b}$ の内積を表します。すると、
f'(\boldsymbol{x}) = A\boldsymbol{x}-\boldsymbol{b} \tag{4}
となり、元の線形方程式(1)は結局のところ
f'(\boldsymbol{x}) = 0 \tag{5}
へと変形されます。
さて、微分値がゼロとなるということはその関数は局所的な最小値、もしくは最大値となっているということですね。ここで、数学的説明をすっ飛ばします(後でちゃんと解説します)が、工学的分野では大抵の場合、
f'(\boldsymbol{x}) = 0 \Leftrightarrow f(\boldsymbol{x})が最小 \tag{6}
が成り立ちます。従って、与えられた線形方程式(1)の解を探すには、関数 $f(\boldsymbol{x})$ を最小とする $\boldsymbol{x}$ を求めればよいわけです。
このように、ある方程式の解が、ある汎関数を最小とする点として特徴付けされる場合、そのことを「変分原理」と呼びます。
数学的なざっくり解説
ここでは、先ほどすっ飛ばした、式(6)の解説をざっくりと行います。少し数学チックになってしまうので、特に気にならない方はすっ飛ばして頂いて構いません。
さて、関数 $f(\boldsymbol{x})$ の定義式(3)を見ると、これは変数ベクトル $\boldsymbol{x}$ の「2次形式」となっていることがお分かりでしょうか。次のように変形した方が分かりやすいかもしれません。
f(\boldsymbol{x}) \equiv \boldsymbol{x} ^\top A \boldsymbol{x} - (\boldsymbol{b}, \boldsymbol{x}) \tag{7}
2次形式が始めましての方は、こちらのサイトなどを参考にしてください。
2次形式で表される関数は、スカラーの2次関数と同様に微分値がゼロとなる点は最大値か最小値のどちらかでしかないことが知られています。
また、行列 $A$ が「正定値行列」と呼ばれる行列だった場合、微分値がゼロとなる点は最大値ではなく、最小値となります。工学的に扱う行列はほとんどが正定値行列ですので、式(6)が成り立つわけです。
あまり深くは触れられませんでしたが、より深く学習したい方への指針は示せたかと思います。
反復法
$f(\boldsymbol{x})$ を最小とする $\boldsymbol{x}$ を反復法で求めるにはどうしたらいいでしょう?
第 $k$ ステップの $\boldsymbol{x}_k$ に対し、第 $k+1$ ステップの $\boldsymbol{x}_{k+1}$ を次式で定めることにします。
\boldsymbol{x}_{k+1} = \boldsymbol{x}_k + \alpha_k\boldsymbol{p}_k \tag{8}
つまり、修正の方向を定めるのが $\boldsymbol{p}_k$、その方向への修正量を定めるのが $\alpha_k$ というわけです。そこで、以下では$\boldsymbol{p}_k$ を「修正方向ベクトル」、$\alpha_k$ を「修正量」と名付けます。
では、$\alpha_k$ と $p_k$ はどのように定めればよいのでしょうか?
これから式変形が続きます。頑張ってついてきてください。まず、式(8)を用いると、$f(\boldsymbol{x}_{k+1})$ が次のように求まります。
\begin{align}
f(\boldsymbol{x}_{k+1}) &= \frac{1}{2}(\boldsymbol{x}_k + \alpha_k\boldsymbol{p}_k, A(\boldsymbol{x}_k + \alpha_k\boldsymbol{p}_k))
- (\boldsymbol{b}, \boldsymbol{x}_k + \alpha_k\boldsymbol{p}_k) \\
&= \frac{1}{2}(\boldsymbol{x}_k, A\boldsymbol{x}_k)
- (\boldsymbol{b}, \boldsymbol{x}_k)
+ \frac{1}{2}\alpha_k(\boldsymbol{x}_k, A\boldsymbol{p}_k)
+ \frac{1}{2}\alpha_k(\boldsymbol{p}_k, A\boldsymbol{x}_k)
+ \frac{1}{2}\alpha_k^2(A\boldsymbol{p}_k, \boldsymbol{p}_k)
- \alpha_k(\boldsymbol{b}, \boldsymbol{p}_k) \tag{9}
\end{align}
ここで、最初の2項は $f(\boldsymbol{x}_k)$ そのものです。また、$A$ が対称行列であれば $(\boldsymbol{x}_k, A\boldsymbol{p}_k)=(\boldsymbol{p}_k, A\boldsymbol{x}_k)$ となります。「ざっくり解説」で工学的なほとんどの場合では $A$ が「正定値行列」であるとお話ししましたが、正定値行列はすべて対称行列です。従って、この変形も大抵の場合で問題ないことになります。すると、$f(\boldsymbol{x}_{k+1})$ は次のように変形できます。
f(\boldsymbol{x}_{k+1})
= \frac{(A\boldsymbol{p}_k, \boldsymbol{p}_k)}{2}\alpha_k^2
+ (A\boldsymbol{x}_k - \boldsymbol{b}, \boldsymbol{p}_k)\alpha_k
+ f(\boldsymbol{x}_k) \tag{10}
ここで、残差 $\boldsymbol{r}_k$ を次のように定義します。
\boldsymbol{r}_k \equiv 0 - f'(\boldsymbol{x}_k) = \boldsymbol{b} - A\boldsymbol{x}_k \tag{11}
これを用いれば、$f(\boldsymbol{x}_{k+1})$ が次式へと変形されます。
f(\boldsymbol{x}_{k+1})
= \frac{(A\boldsymbol{p}_k, \boldsymbol{p}_k)}{2}\alpha_k^2
- (\boldsymbol{r}_k, \boldsymbol{p}_k)\alpha_k
+ f(\boldsymbol{x}_k) \tag{12}
式変形は以上で終了です。お疲れ様でした。
さて、式(12)を $\alpha_k$ に関する降べきの順に整理したことにお気づきでしょうか。実は、これには意図があります。今、私たちの目的は $f(\boldsymbol{x}_{k+1})$ を最小にすることです。この時、$\alpha_k$ とその他の変数は独立していますから、$\alpha_k$ 以外の変数を固定し、$\alpha_k$ だけを動かせばいいのです。つまり、
\begin{align}
\frac{\partial f(\boldsymbol{x}_{k+1})}{\partial \alpha_k} &= 0 \\
(A\boldsymbol{p}_k, \boldsymbol{p}_k)\alpha_k - (\boldsymbol{r}_k, \boldsymbol{p}_k) &= 0 \\
\alpha_k &= \frac{(\boldsymbol{r}_k, \boldsymbol{p}_k)}{(A\boldsymbol{p}_k, \boldsymbol{p}_k)} \tag{13}
\end{align}
となります。
さて、修正量 $\alpha_k$ は式(13)のように求めればよいです。では、修正方向ベクトル $\boldsymbol{p}_k$ はどのように求めたらいいのでしょうか?
修正方向ベクトル $\boldsymbol{p}_k$ の定め方には主に2つの流儀があり、それぞれ「最急降下法」と「共役勾配法(CG法)」と呼ばれています。まずは、最急降下法について解説します。
最急降下法
$f'(\boldsymbol{x})$ は $f(\boldsymbol{x})$ の勾配、すなわち $f(\boldsymbol{x})$ が最も急速に増加する方向を示しています。裏を返せば、$-f'(\boldsymbol{x})$ は $f(\boldsymbol{x})$ が最も急速に減少する方向を示しています。これを修正方向ベクトル $\boldsymbol{p}$ として採用するのが「最急降下法」です。すなわち、
\begin{align}
\boldsymbol{p}_k &= -f'(\boldsymbol{x}_k) \\
&= \boldsymbol{b} - A\boldsymbol{x}_k \\
&=\boldsymbol{r}_k \tag{14}
\end{align}
とするわけですね。
しかし、最急降下法は実はそれほど速くないのです。なぜなら、各反復での修正方向として「最急降下」方向を取るという戦略は、過去の反復の際の修正方向データを使っていないからです。
$f(\boldsymbol{x})$ は2次形式ですが、これはスカラーの2次関数とよく似た性質を持ちます。従って、ある方向に最小値を探索したら、その方向への探索はもう必要ないのです。
最急降下法の問題を改善し、修正方向を選ぶ際に過去の修正方向と重ならないようにするのが共役勾配法(CG法)です。
CG法については次の記事で解説していますので、ぜひご覧ください。
参考資料
戸川隼人(1977)『シリーズ新しい応用の数学 17 共役勾配法』教育出版株式会社
この記事の大部分はこの資料から勉強しました。丸々一冊がCG法に関する内容という非常に贅沢な本です。展開も付いていきやすいです。本腰を入れて勉強したい方にはぜひおすすめしたい良書です。