共轭梯度法python实现

目录

1.简介

2.程序设计思路

3.程序代码

4.测试样例

---

1.简介

通过采用共轭梯度类似的解题步骤，用python编程实现，做到考虑问题minfx=12xTAx+bTx+c(其中 A为对称矩阵，c为常数)，输入初始区间A,b,c,以及初始点，输出最优解。

2.程序设计思路

共轭梯度法的步骤：

考虑问题minfx=12xTAx+bTx+c(其中 A为对称矩阵，c为常数)

1. 任意给定一个初始点x(1),置k=1
2. 计算出目标函数fx在这点的梯度gk=∇fxk，若gk=0,则停止计算，得x=xk；否则，进行下一步
3. 令dk=-gk+βk-1dk-1，当k=1时，βk-1=0，当k>1时，βk=gk+12gk2
4. 计算λk=-g1Td(1)d1TAd1，x(k+1)=x(k)+λkd(k)
5. 若k=n，则停止计算，得x=xk+1；否则，置k=k+1，返回步骤（2）

3.程序代码

'''
共轭梯度法
'''
import numpy as np
from fractions import Fraction as f
def getinput():
    global matrix,b,c,m,n,e
    string=input('''求解的函数 min 1/2XT*A*X+bXT+C
    请输入对称正定矩阵A,以;作为每一行的分隔：
    ''')
    a = [list(map(eval,row.split())) for row in string.split(';')]
    matrix = np.mat(a)
    m, n = matrix.shape
    string2 = input('''请输入B:''')
    b=list(map(eval,string2.split(' ')))
    string3=input('''请输入C:''')
    c=int(string3)
    string4=input('''取初始点，例如：（5，4），则输入“5 4”''')
    e=list(map(eval,string4.split(' ')))
    e=np.array(e)
    x=[]
    print('\n输入的待求解函数为')
    for i in range(m):
        for j in range(n):
            if matrix[i,j]!=0:
                x.append(f'{matrix[i,j]}*x_{i+1}*x_{j+1}')
    for i in range(len(b)):
        if b[i]!=0:
            x.append(f'{b[i]}*x_{i+1}')
    if c!=0:
        x.append(f'{c}')
    print('min ' + ' + '.join(x))
    # print('\n初始点为：['+str(e[0])+','+str(e[1])+']')
def calulateg(p):
    global g
    for i in range(m):
        for j in range(n):
            g[p][i] = 0
    for i in range(m):
        for j in range(n):
            # print(i)
            g[p][i] += 2 * matrix[i, j] * e[j]
    g[p][i] += b[i]
def caculated(p):
    global d
    if p==0:
        d = -g
    else:
        d[p]=-g[p]+np.dot(h[0],d[p-1])
    print("d=="+str(d[p]))
def caculateβ(p):
    global h
    #h=np.dot(np.dot(d[p],matrix),g[p+1])/np.dot(np.dot(d[p],matrix),d[p])
    h=np.dot(g[p],g[p])/np.dot(g[p-1],g[p-1])
    h = np.array(h)
    h = h.reshape(-1)
    print("β="+str(h[0]))
def caculateλ(p):
    global e,d
    u=np.dot(-1,np.dot(g[p],d[p]))/np.dot(np.dot(d[p],np.dot(2,matrix)),d[p])
    u=np.array(u)
    u=u.reshape(-1)
    e=np.float64(e)
    d=np.float64(d)
    print("λ=" + str(u[0]))
    for i in range(len(e)):
        e[i]+=u[0]*d[p][i]
def judge(p):
    if np.dot(g[p],g[p])==0:
        return 0
    else:
        return 1
def main():
    global g
    x=[]
    getinput()
    g=np.zeros((5000,n))

    for i in range(2):
        calulateg(i)
    i = 0
    j = 0
    while judge(i):
        print("第" + str(i + 1) + "次迭代")
        calulateg(i)
        if i!=0:
            caculateβ(i)
        caculated(i)
        caculateλ(i)
        i+=1
    print("最优解：")
    for l in range(m):
        print(f'x{l+1}='+str(format(e[l], '.2f')))
if __name__ == '__main__':
   main()

4.测试样例

四元：

五元：