共役勾配法

GPUによる科学計算の例として、たびたび出てくる共役勾配法(Conjugate Gradient Method)を、復習してみた。

http://daisy.math.sci.ehime-u.ac.jp/users/tsuchiya/math/linear/cg/

に、ソースコードがあったので、これを参考にアルゴリズムを整理すれば、Ax=bのxを求めようとした場合、

1. x=0として、残差ベクトルr=b-Ax、修正ベクトルp=b、を初期値とする
2. alpha=pr/pAp
3. x=x+alpha*p, r=r-alpha*Ap
4. もし、rが一定値以下なら、終了
5. beta=-rAp/pAp
6. p=r+beta*p
7. 2に戻る

という感じ。演算自体は単純なので、すぐに実装できたが、なぜに修正ベクトルの更新が、p=r+beta*p で良いのかが、まだ理解しきれていないので、どうも気持ち悪い。