共役勾配法

出典: フリー百科事典『ウィキペディア(Wikipedia)』
移動: 案内検索
線型方程式の二次形式を最小化するための、最適なステップサイズによる最急降下法(緑)の収束と共役勾配法(赤)の収束の比較。共役勾配法は、厳密にはn次の係数行列に対して高々nステップで収束する(ここではn=2)。

共役勾配法(きょうやくこうばいほう、: conjugate gradient method、CG法とも呼ばれる)は対称正定値行列を係数とする連立一次方程式を解くためのアルゴリズムである。反復法として利用され、コレスキー分解のような直接法では大きすぎて取り扱えない、大規模な疎行列を解くために利用される。そのような問題は偏微分方程式などを数値的に解く際に常に現れる。

共役勾配法は、エネルギー最小化などの最適化問題を解くために用いることもできる。

双共役勾配法は、共役勾配法の非対称問題への拡張である。また、非線形問題を解くために、さまざまな非線形共役勾配法が提案されている。

詳説[編集]

対称正定値行列Aを係数とするn元連立一次方程式

Ax = b

の解をx*とする。

直接法としての共役勾配法[編集]

非零ベクトルuv

 \mathbf{u}^{\mathrm{T}} \mathbf{A} \mathbf{v} = \mathbf{0}

を満たすとき、uvAに関して共役であるという。Aは対称正定値なので、左辺から内積

 \langle \mathbf{u},\mathbf{v} \rangle_\mathbf{A} := \langle \mathbf{A}^{\mathrm{T}} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{A} \mathbf{u}, \mathbf{v}\rangle = \langle \mathbf{u}, \mathbf{A}\mathbf{v} \rangle = \mathbf{u}^{\mathrm{T}} \mathbf{A} \mathbf{v}

を定義することができる。この内積に関して2つのベクトルが直交するなら、それらのベクトルは互いに共役である。 この関係は対称で、uvに対して共役なら、vuに対して共役である(この場合の「共役」は複素共役と無関係であることに注意)。

{pk}をn個の互いに共役なベクトル列とする。pk基底Rnを構成するので、Ax = bの解x*をこの基底で展開すると、

 \mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{p}_i

と書ける。ただし係数は

 \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i  \mathbf{A} \mathbf{p}_i = \mathbf{b}.
 \mathbf{p}_k^{\mathrm{T}} \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i\mathbf{p}_k^{\mathrm{T}} \mathbf{A} \mathbf{p}_i= \mathbf{p}_k^{\mathrm{T}} \mathbf{b}.
 \alpha_k = \frac{\mathbf{p}_k^{\mathrm{T}} \mathbf{b}}{\mathbf{p}_k^{\mathrm{T}} \mathbf{A} \mathbf{p}_k} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\langle \mathbf{p}_k,  \mathbf{p}_k\rangle_\mathbf{A}} = \frac{\langle \mathbf{p}_k, \mathbf{b}\rangle}{\,\,\,\|\mathbf{p}_k\|_\mathbf{A}^2}.

で与えられる。

この結果は、上で定義した内積を考えるのが最も分かりやすいと思われる。

以上から、Ax = bを解くための方法が得られる。すなわち、まずn個の共役な方向を見つけ、それから係数αkを計算すればよい。

反復法としての共役勾配法[編集]

共役なベクトル列pkを注意深く選ぶことにより、一部のベクトルからx*の良い近似を得られる可能性がある。そこで、共役勾配法を反復法として利用することを考える。こうすることで、nが非常に大きく、直接法では解くのに時間がかかりすぎるような問題にも適用することができる。

x*の初期値をx0 = 0 とする。x*二次形式

 f(\mathbf{x}) = \frac12 \mathbf{x}^{\mathrm{T}} \mathbf{A}\mathbf{x} - \mathbf{b}^{\mathrm{T}} \mathbf{x} , \quad \mathbf{x}\in\mathbf{R}^n.

を最小化する一意な解であることに注意し、最初の基底ベクトルp1x = x0でのfの勾配Ax0b=−bとなるように取る。このとき、基底の他のベクトルは勾配に共役である。そこで、この方法を共役勾配法と呼ぶ。

rkkステップ目での残差

 \mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_k

とする。rkx = xkでのfの負の勾配であることに注意されたい。最急降下法rkの方向に進む解法である。pkは互いに共役でなければならないので、rkに最も近い方向を共役性を満たすように取る。これは

 \mathbf{p}_{k+1} = \mathbf{r}_k - \frac{\mathbf{p}_k^{\mathrm{T}} \mathbf{A} \mathbf{r}_k}{\mathbf{p}_k^{\mathrm{T}}\mathbf{A} \mathbf{p}_k} \mathbf{p}_k

のように表すことができる(記事冒頭の図を参照)。

アルゴリズム[編集]

以上の方法を簡素化することにより、Aが実対称正定値である場合にAx = bを解くための以下のアルゴリズムを得る。初期ベクトルx0は近似解もしくは0とする。

r_0 := b - A x_0 \,
p_0 := r_0 \,
k := 0 \,

repeat

\alpha_k := \frac{r_k^\top r_k}{p_k^\top A p_k}  \,
x_{k+1} := x_k + \alpha_k p_k \,
r_{k+1} := r_k - \alpha_k A p_k \,
if rk+1 is "sufficiently small" then exit loop end if
\beta_k := \frac{r_{k+1}^\top r_{k+1}}{r_k^\top r_k}  \,
p_{k+1} := r_{k+1} + \beta_k p_k \,
k := k + 1 \,

end repeat
The result is x_{k+1} \,

Octaveでの共役勾配法の記述例[編集]

Gnu Octaveで書くとは以下のようになる。

 function [x] = conjgrad(A,b,x0)

    r = b - A*x0;
    w = -r;
    z = A*w;
    a = (r'*w)/(w'*z);
    x = x0 + a*w;
    B = 0;

    for i = 1:size(A)(1);
       r = r - a*z;
       if( norm(r) < 1e-10 )
            break;
       endif
       B = (r'*z)/(w'*z);
       w = -r + B*w;
       z = A*w;
       a = (r'*w)/(w'*z);
       x = x + a*w;
    end
 end

前処理[編集]

前処理行列とは、Aと同値なP-1A (PT)-1条件数Aより小さく、Ax=bよりP-1A (PT)-1 x′ =P-1b′の方が容易に解けるような正定値行列 P.PTを指す。前処理行列の生成には、ヤコビ法ガウス・ザイデル法、対称SOR法などが用いられる。

最も単純な前処理行列は、Aの対角要素のみからなる対角行列である。これはヤコビ前処理または対角スケーリングとして知られている。対角行列は逆行列の計算が容易かつメモリも消費しない点で、入門用として優れた方法である。より洗練された方法では、κ(A)の減少による収束の高速化とP-1の計算に要する時間とのトレードオフを考えることになる。

正規方程式に対する共役勾配法[編集]

ATAは任意の行列に対して対称(半)正定値となるので、係数行列をATA、右辺をATbとする正規方程式を解くことにより、共役勾配法を任意のn×m行列に対して適用することができる(CGNR法)。

ATAx = ATb

反復法としては、ATAを明示的に保持する必要がなく、行列ベクトル積、転置行列ベクトル積を計算すればよいので、Aが疎行列である場合にはCGNR法は特に有効である。ただし、条件数κ(ATA)がκ(A2)に等しいことから収束は遅くなる傾向があり、前処理行列を使用するCGLS、LSQRなどの解法が提案されている。LSQRはAが悪条件である場合に最も数値的に安定な解法である。

関連項目[編集]

参考文献[編集]

外部リンク[編集]