共役勾配法

ファイル:Conjugate gradient illustration.svg

線型方程式の二次形式を最小化するための、最適なステップサイズによる最急降下法（緑）の収束と共役勾配法（赤）の収束の比較。共役勾配法は、厳密にはn次の係数行列に対して高々nステップで収束する（ここではn=2）。

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

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

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

詳説

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

Ax = b

の解をx_*とする。

直接法としての共役勾配法

非零ベクトルu、vが

[math] \mathbf{u}^{\mathrm{T}} \mathbf{A} \mathbf{v} = \mathbf{0}[/math]

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

[math] \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}[/math]

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

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

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

と書ける。ただし係数は

[math] \mathbf{A}\mathbf{x}_* = \sum^{n}_{i=1} \alpha_i \mathbf{A} \mathbf{p}_i = \mathbf{b}.[/math]
[math] \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}.[/math]
[math] \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}. [/math]

で与えられる。

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

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

反復法としての共役勾配法

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

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

[math] 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. [/math]

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

r_kをkステップ目での残差

[math] \mathbf{r}_k = \mathbf{b} - \mathbf{Ax}_k [/math]

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

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

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

アルゴリズム

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

[math]r_0 = b - A x_0[/math]
[math]p_0 = r_0[/math]

for (k = 0; ; k++) 
    [math]\alpha_k = \frac{r_k^{\rm T} p_k}{p_k^{\rm T} A p_k}[/math]
    [math]x_{k+1} = x_k + \alpha_k p_k[/math]
    [math]r_{k+1} = r_k - \alpha_k A p_k[/math]

    if [math]r_{k + 1}[/math] が十分に小さい then
        break

    [math]\beta_k = \frac{r_{k+1}^{\rm T} r_{k+1}}{r_k^{\rm T} r_k}[/math]
    [math]p_{k+1} = r_{k+1} + \beta_k p_k[/math]
結果は [math]x_{k+1}[/math]

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

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

正規方程式に対する共役勾配法

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

A^TAx = A^Tb

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

関連項目

• 双共役勾配法 (BICG)
• 非線形共役勾配法

参考文献

• Hestenes, Magnus R.; Stiefel, Eduard (1952-12). “Methods of Conjugate Gradients for Solving Linear Systems”. Journal of Research of the National Bureau of Standards 49 (6). http://nvl.nist.gov/pub/nistpubs/jres/049/6/V49.N06.A08.pdf. 
• Atkinson, Kendell A. (1988). “Section 8.9”, An introduction to numerical analysis, 2nd, John Wiley and Sons. ISBN 0-471-50023-2. 
• Avriel, Mordecai (2003). Nonlinear Programming: Analysis and Methods.. Dover Publishing.. ISBN 0-486-43227-0. 
• Golub, Gene H. (1996). “Chapter 10”, Matrix computations, 2nd, Johns Hopkins University Press.. ISBN 0-8018-5414-8. 

外部リンク

• Conjugate Gradient Method by Nadir Soualem.
• Preconditioned conjugate gradient method by Nadir Soualem.
• An Introduction to the Conjugate Gradient Method Without the Agonizing Pain by Jonathan Richard Shewchuk.
• Iterative methods for sparse linear systems by Yousef Saad
• LSQR: Sparse Equations and Least Squares by Christopher Paige and Michael Saunders.