拟牛顿法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)
2. 置H1=In(单位矩阵)，计算出目标函数fx在x1的梯度gk=∇fxk，置k=1
3. 令dk=-Hkgk-1，Hk+1=Hk+p(k)p(k)Tp(k)Tq(k)-Hkq(k)q(k)THkq(k)THkq(k)，p(k)=x(k+1)-xk=λkdk, q(k)=g(k+1)-gk
4. 计算λk=-g1Td(1)d1TAd1，x(k+1)=x(k)+λkd(k)
5. 若gk=0,则停止计算，得x=xk；否则，置k=k+1，返回（3）

 

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
    d[p]=np.dot(-H,g[p])
    print("g="+str(g[p]))
    print("d="+str(d[p]))
def caculateH(p):
    global H,s,f,q
    if p==0:
        s=H
    else:
        f=np.mat(f)
        q=np.mat(q)
        H = s + np.dot(f.T, f) / np.dot(f, q.T) - np.dot(np.dot(np.dot(s, q.T), q), s) / np.dot(np.dot(q, s), q.T)
        s=H
def caculatepq(p):
    global h,w,f,q
    w=np.array(w)
    f = np.dot(w,d[p])
    q = g[p] - g[p-1]
    # print("p=")
    # print(f)
    # print("q=")
    # print(q)
def caculateλ(p):
    global e,d,w
    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)
    w=u[0]
    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,H,v
    x=[]
    getinput()
    g=np.zeros((5000,n))
    H=np.identity(m)
    for i in range(2):
        calulateg(i)
    i = 0
    j = 0
    while judge(i):
        print("第" + str(i + 1) + "次迭代")
        calulateg(i)
        if i!=0:
            caculatepq(i)
        caculateH(i)
        caculated(i)
        caculateλ(i)
        i+=1
    print("最优解：")
    for l in range(m):
        print(f'x{l+1}='+str(format(e[l], '.3f')))
if __name__ == '__main__':
   main()

4.测试样例

 四元：