共役勾配法

提供: testwiki
2025年3月10日 (月) 15:23時点におけるimported>Oyyo37による版 (校正)
(差分) ← 古い版 | 最新版 (差分) | 新しい版 → (差分)
ナビゲーションに移動 検索に移動
線型方程式の二次形式を最小化するための、最適なステップサイズによる最急降下法（緑）の収束と共役勾配法（赤）の収束の比較。共役勾配法は、厳密にはn次の係数行列に対して高々nステップで収束する（ここではn=2）。

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

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

テンプレート:仮リンクは、共役勾配法の非対称問題への拡張である[11]。

また、非線形問題を解くために、さまざまな非線形共役勾配法が提案されている[12][13][14][15]。

詳説

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

テンプレート:Indent

の解をx*とする。

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

非零ベクトルu、vが テンプレート:Indent を満たすとき、u、vはAに関して共役であるという[2][3][4]。Aは対称正定値なので、左辺から内積 テンプレート:Indent を定義することができる。この内積に関して2つのベクトルが直交するなら、それらのベクトルは互いに共役である。 この関係は対称で、uがvに対して共役なら、vもuに対して共役である（この場合の「共役」は複素共役と無関係であることに注意）。

{pk}をn個の互いに共役なベクトル列とする。pkは基底Rnを構成するので、Ax = bの解x*をこの基底で展開すると、 テンプレート:Indent と書ける。ただし係数は テンプレート:Indent で与えられる。

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

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

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

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

x*の初期値をx0 = 0 とする。x*が二次形式 テンプレート:Indent を最小化する一意な解であることに注意し、最初の基底ベクトルp1をx = x0でのfの勾配Ax0−b=−bとなるように取る。このとき、基底の他のベクトルは勾配に共役である。そこで、この方法を共役勾配法と呼ぶ[2][3][4]。

rkをkステップ目での残差 テンプレート:Indent とする。rkはx = xkでのfの負の勾配であることに注意されたい。最急降下法はrkの方向に進む解法である。pkは互いに共役でなければならないので、rkに最も近い方向を共役性を満たすように取る。これは テンプレート:Indent のように表すことができる（記事冒頭の図を参照）。

アルゴリズム

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

r0=b−Ax0
p0=r0

for (k = 0; ; k++) 
    αk=rkTpkpkTApk
    xk+1=xk+αkpk
    rk+1=rk−αkApk

    if rk+1 が十分に小さい then
        break

    βk=rk+1Trk+1rkTrk
    pk+1=rk+1+βkpk
結果は xk+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を指す[4]。前処理行列の生成には、ヤコビ法、ガウス・ザイデル法、対称SOR法などが用いられる[16][17]。

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

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

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

テンプレート:Indent

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

脚注

テンプレート:脚注ヘルプ

出典

テンプレート:Reflist

参考文献

テンプレート:Refbegin

テンプレート:Refend

関連項目

外部リンク

テンプレート:Linear algebra テンプレート:最適化アルゴリズム テンプレート:Authority control

  1. ↑ 1.0 1.1 1.2 テンプレート:Cite book
  2. ↑ 2.0 2.1 2.2 2.3 2.4 森正武. 数値解析 第2版. 共立出版.
  3. ↑ 3.0 3.1 3.2 3.3 3.4 数値線形代数の数理とHPC, 櫻井鉄也, 松尾宇泰, 片桐孝洋編（シリーズ応用数理 / 日本応用数理学会監修, 第6巻）共立出版, 2018.8
  4. ↑ 4.0 4.1 4.2 4.3 4.4 4.5 4.6 皆本晃弥. (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. テンプレート:URL
  12. ↑ Loyce Adams, and J. L. Nazareth (Eds.): Linear and Nonlinear Conjugate Gradients-Related Methods, SIAM, ISBN 0-89871-376-5 (1996).
  13. ↑ Dai, Y. H. (2010). Nonlinear conjugate gradient methods. Wiley Encyclopedia of Operations Research and Management Science.
  14. ↑ Hager, W. W., & Zhang, H. (2006). A survey of nonlinear conjugate gradient methods. Pacific journal of Optimization, 2(1), 35-58.
  15. ↑ 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.
  16. ↑ 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.
  17. ↑ Kaasschieter, E. F. (1988). Preconditioned conjugate gradients for solving singular systems. Journal of Computational and Applied Mathematics, 24(1-2), 265-275.
  18. ↑ Black, Noel and Moore, Shirley. "Conjugate Gradient Method on the Normal Equations." From MathWorld--A Wolfram Web Resource, created by Eric W. Weisstein. テンプレート:URL
  19. ↑ Bjorck, A. (1996). Numerical methods for least squares problems (Vol. 51). SIAM.
  20. ↑ Paige, C. and Saunders, M. "LSQR: An Algorithm for Sparse Linear Equations and Sparse Least Squares." ACM Trans. Math. Soft. 8, 43-71, 1982.
  21. ↑ 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.
「https://ja.wiki.beta.math.wmflabs.org/w/index.php?title=共役勾配法&oldid=5127」から取得