拟牛顿法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(单位矩阵),计算出目标函数fxx1的梯度gk=fxk,置k=1
  3. dk=-Hkgk-1Hk+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)d1TAd1x(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.测试样例

 四元:

 

  • 0
    点赞
  • 3
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

别管我啦就是说

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值