共轭梯度法（CG） python实现

共轭梯度法是求解对称正定线性方程组的首选方法

下面给出代码

import numpy as np

def CG(x,A,b,M,epsilon):
    '''
    :param x0: 初值
    :param A: Ax=b的系数矩阵
    :param b: 常数向量
    :param M: 对大迭代次数
    :param epsilon: 误差最大限度 
    :return: x的数值解
    '''
    r = b - np.dot(A,x)
    beta = np.linalg.norm(r,2)
    if beta < epsilon:
        print(x)
    else:
        for m in range(M):
            rho = np.dot(np.transpose(r),r)
            if m == 0:
                p = r
            elif m > 0:
                mju = rho / rho0
                p = r + np.dot(mju,p)
            q = np.dot(A,p)
            Xi = rho / np.dot(np.transpose(p),q)
            x = x + np.dot(Xi,p)
            r = r - np.dot(Xi,q)
            relres = np.linalg.norm(r,2) / beta
            if relres < epsilon:
                break
            else:
                rho0 = rho
        if relres < epsilon:
            return x
        else:
            return 0


from sklearn.datasets import make_spd_matrix
A =make_spd_matrix(10, random_state=42)#随机生成对阵正定矩阵
m = 1000
x = np.zeros(10, dtype = float)
b = np.array([2., 3., 3.,3., 3.,3., 3., 3., 3., 2.])
f=CG(x,A,b,m,0.0000001)
print(f)     

下面验证代码的正确性：

由随机取得的A第一行与得到的x的第一列做乘积求和

A1 =[ 0.89102384  0.72750936 -0.43020245  0.01204145 -0.90973822  0.2562303
  -0.08113697  0.14219761  0.23773209  1.28879948]

x1 =[40.28357866 15.254577   84.6691862  78.080485   53.79086054 37.54488589
 52.00030268 47.05693282 46.05108637 12.71619   ]

得到的结果为1.999999516893963可以看做b的第一个值2

故代码是正确的

算法来源于 矩阵计算讲义 潘建瑜  算法7.7