共役勾配法

Aが対称行列のときにAx=bを解く、反復法の最も基本的なやつです。
Aの状態数が悪い場合は、不完全コレスキー分解等を用いて前処理を行う方法があります。

from numpy import *
import copy

# 共役勾配法によりAx=bを解く(Aは正値対称行列)
def ConjugateGradientMethod(amat, b, eps, nloop):
   x = matrix(zeros((amat.shape[1],1)))
   r = b - amat*x
   p = copy.copy(r)
   for k in range(nloop):
      if linalg.norm(r) < eps:
           break
       alpha = (r.T*r)/(p.T*amat*p)
       x += alpha[0,0]*p
       rold = copy.copy(r)
       r -= alpha[0,0]*amat*p
       beta = (r.T*r)/(rold.T*rold)
       p = r + beta[0,0]*p
   return x

A = matrix([[1.0, 0.5],
             [0.5, 1.0]])
b = matrix([[0.5],[0.0]])
x = ConjugateGradientMethod(A, b, 1.0e-5, 2)
print x
最終更新:2011年04月30日 02:41