共役勾配法

対称正定値行列を係数とする連立一次方程式を解くためのアルゴリズム

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

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

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

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

詳説

編集

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

Ax = b

の解をx*とする。

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

編集

非零ベクトルuv

 

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

 

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

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

 

と書ける。ただし係数は

 
 
 

で与えられる。

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

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

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

編集

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

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

 

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

rkkステップ目での残差

 

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

 

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

アルゴリズム

編集

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

 
 

for (k = 0; ; k++) 
     
     
     

    if   が十分に小さい then
        break

     
     
結果は  

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を指す[4]。前処理行列の生成には、ヤコビ法ガウス・ザイデル法、対称SOR法などが用いられる[15][16]

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

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

編集

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

ATAx = ATb

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

関連項目

編集

脚注

編集

出典

編集
  1. ^ a b c 山本哲朗『数値解析入門』(増訂版)サイエンス社〈サイエンスライブラリ 現代数学への入門 14〉、2003年6月。ISBN 4-7819-1038-6 
  2. ^ a b c d e 森正武. 数値解析 第2版. 共立出版.
  3. ^ a b c d e 数値線形代数の数理とHPC, 櫻井鉄也, 松尾宇泰, 片桐孝洋編(シリーズ応用数理 / 日本応用数理学会監修, 第6巻)共立出版, 2018.8
  4. ^ a b c d e f g 皆本晃弥. (2005). UNIX & Informatioin Science-5 C 言語による数値計算入門.
  5. ^ 田端正久; 偏微分方程式の数値解析, 2010. 岩波書店.
  6. ^ 登坂宣好, & 大西和榮. (2003). 偏微分方程式の数値シミュレーション. 東京大学出版会.
  7. ^ Zworski, M. (2002). Numerical linear algebra and solvability of partial differential equations. Communications in mathematical physics, 229(2), 293-307.
  8. ^ Gill, P. E., Murray, W., & Wright, M. H. (1991). Numerical linear algebra and optimization (Vol. 1, p. 74). Redwood City, CA: Addison-Wesley.
  9. ^ Gilbert, J. C., & Nocedal, J. (1992). Global convergence properties of conjugate gradient methods for optimization. SIAM Journal on optimization, 2(1), 21-42.
  10. ^ Steihaug, T. (1983). The conjugate gradient method and trust regions in large scale optimization. SIAM Journal on Numerical Analysis, 20(3), 626-637.
  11. ^ Black, Noel and Moore, Shirley. "Biconjugate Gradient Method." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/BiconjugateGradientMethod.html
  12. ^ Dai, Y. H. (2010). Nonlinear conjugate gradient methods. Wiley Encyclopedia of Operations Research and Management Science.
  13. ^ Hager, W. W., & Zhang, H. (2006). A survey of nonlinear conjugate gradient methods. Pacific journal of Optimization, 2(1), 35-58.
  14. ^ Dai, Y., Han, J., Liu, G., Sun, D., Yin, H., & Yuan, Y. X. (2000). Convergence properties of nonlinear conjugate gradient methods. SIAM Journal on Optimization, 10(2), 345-358.
  15. ^ Eisenstat, S. C. (1981). Efficient implementation of a class of preconditioned conjugate gradient methods. SIAM Journal on Scientific and Statistical Computing, 2(1), 1-4.
  16. ^ Kaasschieter, E. F. (1988). Preconditioned conjugate gradients for solving singular systems. Journal of Computational and Applied Mathematics, 24(1-2), 265-275.
  17. ^ Black, Noel and Moore, Shirley. "Conjugate Gradient Method on the Normal Equations." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. mathworld.wolfram.com/ConjugateGradientMethodontheNormalEquations.html
  18. ^ Bjorck, A. (1996). Numerical methods for least squares problems (Vol. 51). SIAM.
  19. ^ Paige, C. and Saunders, M. "LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares." ACM Trans. Math. Soft. 8, 43-71, 1982.
  20. ^ Paige, C. C., & Saunders, M. A. (1982). Algorithm 583: LSQR: Sparse linear equations and least squares problems. ACM Transactions on Mathematical Software (TOMS), 8(2), 195-209.

参考文献

編集

外部リンク

編集