python实现共轭梯度法

共轭梯度法（Conjugate Gradient）是介于最速下降法与牛顿法之间的一个方法，它仅需利用一阶导数信息，但克服了最速下降法收敛慢的缺点，又避免了牛顿法需要存储和计算Hesse矩阵并求逆的缺点，共轭梯度法不仅是解决大型线性方程组最有用的方法之一，也是解大型非线性最优化最有效的算法之一。

   这里初始化A为对角线元素全为2且横纵坐标差值绝对值为1时为-1的矩阵，b为全1矩阵，求解x

  算法实现过程：

import numpy as np
A = np.zeros((100, 100))
for i in range(100): #generate A
    for j in range(100):
        if (i == j):
            A[i, j] = 2
        if (abs(i - j) == 1):
            A[i, j] = A[j, i] = -1
b = np.ones((100, 1))  #generate b
print("Conjugate Gradient x:")
x=np.zeros((100,1)) # 初始值x0

r=b-np.dot(A,x)
p=r  #p0=r0
#while np.linalg.norm(np.dot(A, x) - b) / np.linalg.norm(b) >= 10 ** -6:
for i in range(100):
    r1=r
    a=np.dot(r.T,r)/np.dot(p.T,np.dot(A,p))
    x = x + a * p    #x(k+1)=x(k)+a(k)*p(k)
    r=b-np.dot(A,x)  #r(k+1)=b-A*x(k+1)
    q = np.linalg.norm(np.dot(A, x) - b) / np.linalg.norm(b)
    if q<10**-6:
        break
    else:
        beta=np.linalg.norm(r)**2/np.linalg.norm(r1)**2
        p=r+beta*p  #p(k+1)=r(k+1)+beta(k)*p(k)

print(x)
print("done Conjugate Gradient!")