有限内存BFGS以及非精确牛顿法

最新推荐文章于 2023-04-18 21:42:36 发布
最新推荐文章于 2023-04-18 21:42:36 发布
阅读量2.5k 收藏 13
点赞数 3
分类专栏： 优化方法
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/forrestguang/article/details/111055075
版权

有限内存BFGS

BFGS算法的每一步迭代皆形如：
x k + 1 = x k − α k H k g k x_{k+1}=x_{k} - \alpha_{k}H_{k}g_{k} xk+1​=xk​−αk​Hk​gk​
H k + 1 = V k T H k V k + ρ k s k s k T , ρ k = 1 y k T s k , V k = I − ρ k y k s k T , s k = x k + 1 − x k , y k = g k + 1 − g k . H_{k+1} = V_{k}^{T}H_{k}V_{k} + \rho_{k}s_{k}s_{k}^{T},\\ \rho_k = \frac{1}{y_{k}^{T}s_{k}},V_{k}=I - \rho_{k}y_{k}s_{k}^{T},\\ s_{k} = x_{k+1}-x_{k},y_{k}=g_{k+1}-g_{k}. Hk+1​=VkT​Hk​Vk​+ρk​sk​skT​,ρk​=ykT​sk​1​,Vk​=I−ρk​yk​skT​,sk​=xk+1​−xk​,yk​=gk+1​−gk​.
在这里插入图片描述
在这里插入图片描述
具体介绍参考链接添加链接描述

非精确牛顿法

对于非精确牛顿法，直接考虑下降方向的选取，以往的牛顿方法直接求解下面这个方程：

G k d = − g k G_{k}d=-g_{k} Gk​d=−gk​
具体的考虑是，通过某种方法得到一个近似解 d k d_k dk​，使得 ∥ G k d k + g k ∥ ≤ η k ∥ g k ∥ \parallel G_{k}d_{k}+g_{k} \parallel \leq \eta_{k} \parallel g_{k} \parallel ∥Gk​dk​+gk​∥≤ηk​∥gk​∥，这里我们使用G-S迭代得到 d k d_k dk​， η \eta η的选取针对不同的函数不太一样，具体参考下文 η \eta η的选取参考下面两种方式。
η 1 \eta_1 η1​表示 η k = m i n { 0.5 , ∥ g k ∥ } \eta_k = min\{0.5,\parallel g_k \parallel \} ηk​=min{0.5,∥gk​∥}， η 2 \eta_2 η2​表示 η k = m a x { γ ∥ g k ∥ / ∥ g k − 1 ∥ , γ η k − 1 } \eta_k = max\{\gamma \parallel g_k \parallel / \parallel g_{k-1} \parallel, \gamma \eta_{k-1}\} ηk​=max{γ∥gk​∥/∥gk−1​∥,γηk−1​}。这个 γ = 0.2 \gamma = 0.2 γ=0.2，且G-S迭代最多进行100次。

测试函数

这里注意我们取 X = ( x 0 , … , x n − 1 ) X = (x_0,\ldots,x_{n-1}) X=(x0​,…,xn−1​)，如果 X ∈ R n X \in R^n X∈Rn。\
f = ∑ i = 1 m f i 2 ( x 0 , x 2 , … , x n − 1 ) f = \sum_{i=1}^{m}f_{i}^{2}(x_0,x_2,\ldots,x_{n-1}) f=∑i=1m​fi2​(x0​,x2​,…,xn−1​)
∙ \bullet ∙ Penalty I function:
α = 1 0 − 5 , \alpha = 10^{-5}, α=10−5,
f i ( X ) = α 1 2 ( x i − 1 ) , 0 ≤ i ≤ n − 1 , f n ( X ) = ( ∑ i = 0 n − 1 x i 2 ) − 0.25. x i n i t = ( 1 , 2 , … , n ) . f_{i}(X) = \alpha^{\frac{1}{2}}(x_i - 1),0\leq i \leq n - 1, f_{n}(X) = (\sum_{i=0}^{n-1} x_{i}^2) - 0.25. \\ x_{init} = (1,2,\ldots,n). fi​(X)=α21​(xi​−1),0≤i≤n−1,fn​(X)=(∑i=0n−1​xi2​)−0.25.xinit​=(1,2,…,n).
∙ \bullet ∙ Trigonometric function:
f i ( X ) = n − ∑ i = 0 n − 1 c o s ( x i ) + ( i + 1 ) ( 1 − c o s ( x i ) ) − s i n ( x i ) , 0 ≤ i ≤ n − 1 , x i n i t = ( 1 n , … , 1 n ) . f_{i}(X) = n - \sum_{i = 0}^{n - 1} cos(x_i) + (i+1)(1 - cos(x_i)) - sin(x_i) ,0\leq i \leq n - 1, x_{init} = (\frac{1}{n},\ldots,\frac{1}{n}). fi​(X)=n−∑i=0n−1​cos(xi​)+(i+1)(1−cos(xi​))−sin(xi​),0≤i≤n−1,xinit​=(n1​,…,n1​).
∙ \bullet ∙ Extended Freudenstein and Roth function:
f i ( X ) = ( x i + x i + 1 ( 5 x i + 1 − x i + 1 2 − 2 ) − 13 ) 2 + ( x i + x i + 1 ( x i + 1 2 + x i + 1 − 14 ) − 29 ) 2 , 0 ≤ i ≤ n − 2 , x i n i t = ( − 2 , − 2 , … , − 2 ) . f_{i}(X) = (x_i + x_{i+1}(5x_{i+1} - x_{i+1}^2 - 2) - 13)^2 + (x_i + x_{i+1}(x_{i+1}^2 + x_{i+1} - 14) - 29)^2 ,\\ 0\leq i \leq n - 2,\\ x_{init} = (-2,-2,\ldots,-2). fi​(X)=(xi​+xi+1​(5xi+1​−xi+12​−2)−13)2+(xi​+xi+1​(xi+12​+xi+1​−14)−29)2,0≤i≤n−2,xinit​=(−2,−2,…,−2).
∙ \bullet ∙ Extended Rosenbrock function:
f 2 i = 10 ∗ ( x 2 i + 1 − x 2 i 2 ) , f 2 i + 1 = 1 − x 2 i . x i n i t = ( − 1.2 , 1 , … , − 1.2 , 1 ) . f_{2i} = 10*(x_{2i+1}-x_{2i}^2),\\ f_{2i+1} = 1 - x_{2i}.\\ \\ x_{init} = (-1.2,1,\ldots,-1.2,1). f2i​=10∗(x2i+1​−x2i2​),f2i+1​=1−x2i​.xinit​=(−1.2,1,…,−1.2,1).

有限内存BFGS代码

首先这里需要重点注意定义函数的方式，绝对要避免for循环，采取张量运算可以提高运行速度。
先拿penalty函数举例子：

import numpy as np
import torch
import time
def FF(X):
    n = X.shape[1]
    al = 1e-5
    return al*((X - 1)**2).sum() + ((X**2).sum() - 0.25)**2
def GG(X):
    n = X.shape[1]
    al = 1e-5
    return 2*al*(X - 1) + 4*X*((X**2).sum() - 0.25)
def HH(X):
    n = X.shape[1]
    al = 1e-5
    return 8*X.T@X + 4*((X**2).sum() - 0.25)*np.eye(n,n)
def gold(FF,GG,d,X):
    rho = 0.3
    
    am = 0;aM = 1;al = 0.5*(am + aM)
    for i in range(100):
        rig1 = FF(X) + rho*al*GG(X)@d.T
        rig2 = FF(X) + (1 - rho)*al*GG(X)@d.T
        lef = FF(X + al*d)
        if lef <= rig1:
            if lef >= rig2:
                break
            else:
                am = al
                al = am + 0.6*(aM - am)
        else:
            aM = al
            al = am + 0.6*(aM - am)
    al_k = al
    #print(i)
    return al_k,2*(i + 1)
def BFGS(FF,GG,HH,X,m,N):
    tic = time.time()
    eps = 1e-5
    n = X.shape[1]
    fun_iter = 0
    X_new = np.zeros([1,n])
    S = np.zeros([m,n])
    Y = np.zeros([m,n])
    alpha = np.zeros(m)
    for i in range(N):
        if i <= m:
            if i == 0:
                al_k,ite = gold(FF,GG,-GG(X),X)
                fun_iter += ite
                X_new = X - GG(X)*al_k
            else:
                s = X_new - X
                y = GG(X_new) - GG(X)
                S[i - 1:i,:] = s
                Y[i - 1:i,:] = y
                #d = np.linalg.solve(HH(X_new),-GG(X_new).T)
                d = -GG(X).T
                al_k,ite = gold(FF,GG,d.T,X_new)
                #print(np.linalg.norm(X))
                X = X_new.copy()
                X_new = X_new - al_k*GG(X_new)
                fun_iter += ite
                #print(i,al_k,np.linalg.norm(X_new - X))
                #print(i,al_k,np.linalg.norm(GG(X_new),np.linalg.norm(GG(X_new)))
                #print(i,al_k,np.linalg.norm(GG(X_new)))
        else:
            s = X_new - X
            y = GG(X_new) - GG(X)
            S[0:m - 1,:] = S[1:m,:]
            Y[0:m - 1,:] = Y[1:m,:]
            S[m - 1:m,:] = s
            Y[m - 1:m,:] = y
            tmp = (s*y).sum()/(y*y).sum()
            q = GG(X_new)
            for j in range(m - 1,- 1,-1):
                alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                q = q - alpha[j]*Y[j:j + 1,:]
            r = tmp*q
            for j in range(0,m,1):
                beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                r = r + S[j:j + 1,:]*(alpha[j] - beta)
            al_k,ite = gold(FF,GG,-r,X_new)
            X = X_new.copy()
            X_new = X_new + (-r)
        left = np.linalg.norm(GG(X))
        rig = max(1,np.linalg.norm(X))*eps
        if left < rig:
            break
    ela = time.time() - tic
    print('iteration:%d,grad_2:%.3e,fun_iter = %d,time = %.2f'%(i + 1,(GG(X)**2).sum(),fun_iter,ela))
    return X_new
n = 1000
m = 15
X = np.zeros([1,n])
for i in range(n):
    X[:,i] = i + 1
X_new = BFGS(FF,GG,HH,X,m,100)
#print(X_new)
print(FF(X_new))

上面的FF,GG,HH表示函数，梯度和Hessen矩阵，定义采取张量，使用goldstein搜索步长（wolf原则本人代码太差，发现不好调节参数），BFGS就是算法。
Extended Freudenstein and Roth function

import torch
import numpy as np
import time

def FF(X):
    n = X.shape[1]
    x = X[:,0:n - 1];y = X[:,1:n]
    return (((x + y*((5 - y)*y - 2) - 13)**2 + (x + y*((y + 1)*y - 14) - 29)**2)**2).sum()
def GG(X):
    n = X.shape[1]
    x = X[:,0:n - 1];y = X[:,1:n]
    vec = np.zeros([1,n])
    lef = (x + y*((5 - y)*y - 2) - 13);rig = (x + y*((y + 1)*y - 14) - 29)
    tmp = (x + y*((5 - y)*y - 2) - 13)**2 + (x + y*((y + 1)*y - 14) - 29)**2
    vec[:,0:n - 1] = 4*tmp*(lef + rig)
    vec[:,1:n] = 4*tmp*(lef*(-3*y**2 + 10*y - 2) + rig*(3*y**2 + 2*y - 14))
    return vec
def HH(X):
    n = X.shape[1]
    x = X[:,0:n - 1];y = X[:,1:n]
    mat = np.zeros([n,n])
    lef = (x + y*((5 - y)*y - 2) - 13);rig = (x + y*((y + 1)*y - 14) - 29)
    tmp = (x + y*((5 - y)*y - 2) - 13)**2 + (x + y*((y + 1)*y - 14) - 29)**2
    tmp_x = 2*(lef + rig);tmp_y = 2*(lef*(-3*y**2 + 10*y - 2) + rig*(3*y**2 + 2*y - 14))
    m_xx = 2*(tmp_x)**2 + 4*tmp
    m_yy = 2*(tmp_y)**2 + 4*tmp*((-3*y**2 + 10*y - 2)**2 + (3*y**2 + 2*y - 14)**2 + lef*(-6*y - 10) + rig*(6*y + 2))
    m_xy = 2*tmp_x*tmp_y + 2*tmp*(-3*y**2 + 10*y - 2 + 3*y**2 + 2*y - 14)
    mat[0:n - 1,0:n - 1] = np.diag(m_xx[0,:])
    mat[1:n,1:n] += np.diag(m_yy[0,:])
    mat += np.diag(m_xy[0,:],1) + np.diag(m_xy[0,:],-1)
    return mat
def gold(FF,GG,d,X):
    rho = 0.3
    
    am = 0;aM = 1;al = 0.5*(am + aM)
    for i in range(100):
        rig1 = FF(X) + rho*al*GG(X)@d.T
        rig2 = FF(X) + (1 - rho)*al*GG(X)@d.T
        lef = FF(X + al*d)
        if lef <= rig1:
            if lef >= rig2:
                break
            else:
                am = al
                al = am + 0.6*(aM - am)
        else:
            aM = al
            al = am + 0.6*(aM - am)
    al_k = al
    #print(i)
    return al_k,2*(i + 1)
def BFGS(FF,GG,HH,X,m,N):
    tic = time.time()
    eps = 1e-5
    n = X.shape[1]
    fun_iter = 0
    X_new = np.zeros([1,n])
    S = np.zeros([m,n])
    Y = np.zeros([m,n])
    alpha = np.zeros(m)
    for i in range(N):
        if i <= m:
            if i == 0:
                al_k,ite = gold(FF,GG,-GG(X),X)
                fun_iter += ite
                X_new = X - GG(X)*al_k
            else:
                s = X_new - X
                y = GG(X_new) - GG(X)
                S[i - 1:i,:] = s
                Y[i - 1:i,:] = y
                tmp = (s*y).sum()/(y*y).sum()
                q = GG(X_new)
                for j in range(i - 1,-1,-1):
                    alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                    q = q - alpha[j]*Y[j:j + 1,:]
                r = tmp*q
                for j in range(0,i,1):
                    beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                    r = r + (alpha[j] - beta)*S[j:j + 1,:]
                al_k,ite = gold(FF,GG,-r,X_new)
                fun_iter += ite
                X = X_new.copy()
                X_new = X_new + al_k*( - r)
        else:
            s = X_new - X
            y = GG(X_new) - GG(X)
            S[0:m - 1,:] = S[1:m,:]
            Y[0:m - 1,:] = Y[1:m,:]
            S[m - 1:m,:] = s
            Y[m - 1:m,:] = y
            tmp = (s*y).sum()/(y*y).sum()
            q = GG(X_new)
            for j in range(m - 1,- 1,-1):
                alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                q = q - alpha[j]*Y[j:j + 1,:]
            r = tmp*q
            for j in range(0,m,1):
                beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                r = r + S[j:j + 1,:]*(alpha[j] - beta)
            al_k,ite = gold(FF,GG,-r,X_new)
            X = X_new.copy()
            X_new = X_new + (-r)
        left = np.linalg.norm(GG(X))
        rig = max(1,np.linalg.norm(X))*eps
        if left < rig:
            break
    ela = time.time() - tic
    print('iteration:%d,grad_2:%.3e,fun_iter = %d,time = %.2f'%(i + 1,(GG(X)**2).sum(),fun_iter,ela))
    return X_new
n = 1000

m = 5
X = -2*np.ones([1,n])

X_new = BFGS(FF,GG,HH,X,m,200)
#print(X_new)
print(FF(X_new))

Trigonometric function

import torch
import numpy as np
import time
def FF(X):
    n = X.shape[1]
    return ((n - np.cos(X).sum() + np.arange(1,n + 1)*(1 - np.cos(X)) - np.sin(X))**2).sum()
def GG(X):
    n = X.shape[1]
    #vec = np.zeros_like(X)
    tmp = n - np.cos(X).sum() + np.arange(1,n + 1)*(1 - np.cos(X)) - np.sin(X)
    return 2*np.sin(X)*tmp.sum() + 2*tmp*(np.arange(1,n + 1)*np.sin(X) - np.cos(X))
def HH(X):
    n = X.shape[1]
    mat = np.zeros([n,n])
    tmp = n - np.cos(X).sum() + np.arange(1,n + 1)*(1 - np.cos(X)) - np.sin(X)
    dia = 2*np.cos(X)*tmp.sum() + 2*np.sin(X)*((np.arange(1,n + 1) + n)*np.sin(X) - np.cos(X)) +\
    2*((np.arange(1,n + 1) + 1)*np.sin(X) - np.cos(X))*(np.arange(1,n + 1)*np.sin(X) - np.cos(X)) + \
    + 2*tmp*(np.sin(X) + np.arange(1,n + 1)*np.cos(X))
    mat = 2*np.sin(X).T*((np.arange(1,n + 1) + n)*np.sin(X) - np.cos(X)) + \
    (np.arange(1,n + 1)*np.sin(X) - np.cos(X)).T*2*np.sin(X)
    mat.flat[::(n + 1)] = dia
    return mat
def gold(FF,GG,d,X):
    rho = 0.3
    
    am = 0;aM = 1;al = 0.5*(am + aM)
    for i in range(100):
        rig1 = FF(X) + rho*al*GG(X)@d.T
        rig2 = FF(X) + (1 - rho)*al*GG(X)@d.T
        lef = FF(X + al*d)
        if lef <= rig1:
            if lef >= rig2:
                break
            else:
                am = al
                al = am + 0.6*(aM - am)
        else:
            aM = al
            al = am + 0.6*(aM - am)
    al_k = al
    #print(i)
    return al_k,2*(i + 1)
def BFGS(FF,GG,HH,X,m,N):
    tic = time.time()
    eps = 1e-5
    n = X.shape[1]
    fun_iter = 0
    X_new = np.zeros([1,n])
    S = np.zeros([m,n])
    Y = np.zeros([m,n])
    alpha = np.zeros(m)
    for i in range(N):
        if i <= m:
            if i == 0:
                al_k,ite = gold(FF,GG,-GG(X),X)
                fun_iter += ite
                X_new = X - GG(X)*al_k
            else:
                s = X_new - X
                y = GG(X_new) - GG(X)
                S[i - 1:i,:] = s
                Y[i - 1:i,:] = y
                tmp = (s*y).sum()/(y*y).sum()
                q = GG(X_new)
                for j in range(i - 1,-1,-1):
                    alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                    q = q - alpha[j]*Y[j:j + 1,:]
                r = tmp*q
                for j in range(0,i,1):
                    beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                    r = r + (alpha[j] - beta)*S[j:j + 1,:]
                al_k,ite = gold(FF,GG,-r,X_new)
                fun_iter += ite
                X = X_new.copy()
                X_new = X_new + al_k*( - r)
        else:
            s = X_new - X
            y = GG(X_new) - GG(X)
            S[0:m - 1,:] = S[1:m,:]
            Y[0:m - 1,:] = Y[1:m,:]
            S[m - 1:m,:] = s
            Y[m - 1:m,:] = y
            tmp = (s*y).sum()/(y*y).sum()
            q = GG(X_new)
            for j in range(m - 1,- 1,-1):
                alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                q = q - alpha[j]*Y[j:j + 1,:]
            r = tmp*q
            for j in range(0,m,1):
                beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                r = r + S[j:j + 1,:]*(alpha[j] - beta)
            al_k,ite = gold(FF,GG,-r,X_new)
            X = X_new.copy()
            X_new = X_new + (-r)
        left = np.linalg.norm(GG(X))
        rig = max(1,np.linalg.norm(X))*eps
        if left < rig:
            break
    ela = time.time() - tic
    print('iteration:%d,grad_2:%.3e,fun_iter = %d,time = %.2f'%(i + 1,(GG(X)**2).sum(),fun_iter,ela))
    return X_new
n = 1000
m = 15
X = np.ones([1,n])/n 
X_new = BFGS(FF,GG,HH,X,m,200)
#print(X_new)
print(FF(X_new))

Extended Rosenbrock function

import torch
import numpy as np
import time



def FF(X):
    n = X.shape[1]
    X_even = X[:,0:n:2];X_odd = X[:,1:n:2]
    return 100*((X_odd - X_even**2)**2).sum() + ((1 - X_even)**2).sum()
def GG(X):
    n = X.shape[1]
    vec = np.zeros_like(X)
    X_even = X[:,0:n:2];X_odd = X[:,1:n:2]
    vec[:,0:n:2] = -400*X_even*(X_odd - X_even**2) - 2*(1 - X_even)
    vec[:,1:n:2] = 200*(X_odd - X_even**2)
    return vec
def HH(X):
    n = X.shape[1];rank = int(n/2)
    X_even = X[:,0:n:2];X_odd = X[:,1:n:2]
    mat = np.zeros([n,n])
    diag = np.zeros(n)
    diag[0:n:2] = -400*(X_odd - X_even**2) + 800*X_even**2 + 2
    diag[1:n:2] = 200*np.ones([1,rank])
    xie = np.zeros(n - 1)
    xie[0:n - 1:2] = -400*X_even
    mat = np.diag(diag) + np.diag(xie,1) + np.diag(xie,-1)
    return mat
def gold(FF,GG,d,X):
    rho = 0.3
    
    am = 0;aM = 1;al = 0.5*(am + aM)
    for i in range(100):
        rig1 = FF(X) + rho*al*GG(X)@d.T
        rig2 = FF(X) + (1 - rho)*al*GG(X)@d.T
        lef = FF(X + al*d)
        if lef <= rig1:
            if lef >= rig2:
                break
            else:
                am = al
                al = am + 0.6*(aM - am)
        else:
            aM = al
            al = am + 0.6*(aM - am)
    al_k = al
    #print(i)
    return al_k,2*(i + 1)
def BFGS(FF,GG,HH,X,m,N):
    tic = time.time()
    eps = 1e-5
    n = X.shape[1]
    fun_iter = 0
    X_new = np.zeros([1,n])
    S = np.zeros([m,n])
    Y = np.zeros([m,n])
    alpha = np.zeros(m)
    for i in range(N):
        if i <= m:
            if i == 0:
                al_k,ite = gold(FF,GG,-GG(X),X)
                fun_iter += ite
                X_new = X - GG(X)*al_k
            else:
                s = X_new - X
                y = GG(X_new) - GG(X)
                S[i - 1:i,:] = s
                Y[i - 1:i,:] = y
                tmp = (s*y).sum()/(y*y).sum()
                q = GG(X_new)
                for j in range(i - 1,-1,-1):
                    alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                    q = q - alpha[j]*Y[j:j + 1,:]
                r = tmp*q
                for j in range(0,i,1):
                    beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                    r = r + (alpha[j] - beta)*S[j:j + 1,:]
                al_k,ite = gold(FF,GG,-r,X_new)
                fun_iter += ite
                X = X_new.copy()
                X_new = X_new + al_k*( - r)
        else:
            s = X_new - X
            y = GG(X_new) - GG(X)
            S[0:m - 1,:] = S[1:m,:]
            Y[0:m - 1,:] = Y[1:m,:]
            S[m - 1:m,:] = s
            Y[m - 1:m,:] = y
            tmp = (s*y).sum()/(y*y).sum()
            q = GG(X_new)
            for j in range(m - 1,- 1,-1):
                alpha[j] = (S[j:j + 1,:]*q).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                q = q - alpha[j]*Y[j:j + 1,:]
            r = tmp*q
            for j in range(0,m,1):
                beta = (Y[j:j + 1,:]*r).sum()/(S[j:j + 1,:]*Y[j:j + 1,:]).sum()
                r = r + S[j:j + 1,:]*(alpha[j] - beta)
            al_k,ite = gold(FF,GG,-r,X_new)
            X = X_new.copy()
            X_new = X_new + (-r)
        left = np.linalg.norm(GG(X))
        rig = max(1,np.linalg.norm(X))*eps
        if left < rig:
            break
    ela = time.time() - tic
    print('iteration:%d,grad_2:%.3e,fun_iter = %d,time = %.2f'%(i + 1,(GG(X)**2).sum(),fun_iter,ela))
    return X_new
n = 1000
rank = int(n/2)
m = 5
X = np.zeros([1,n])
for i in range(rank):
    X[:,2*i] = -1.2
    X[:,2*i + 1] = 1
X_new = BFGS(FF,GG,HH,X,m,200)
#print(X_new)
print(FF(X_new))

关于上面的代码，重点是函数的定义

非精确牛顿法代码

上面有一个三角函数的Hessen矩阵不正定，G-S迭代不收敛，那里需要添加一个对角矩阵修改：

import torch
import numpy as np
import time
def FF(X):
    n = X.shape[1]
    return ((n - np.cos(X).sum() + np.arange(1,n + 1)*(1 - np.cos(X)) - np.sin(X))**2).sum()
def GG(X):
    n = X.shape[1]
    #vec = np.zeros_like(X)
    tmp = n - np.cos(X).sum() + np.arange(1,n + 1)*(1 - np.cos(X)) - np.sin(X)
    return 2*np.sin(X)*tmp.sum() + 2*tmp*(np.arange(1,n + 1)*np.sin(X) - np.cos(X))
def HH(X):
    n = X.shape[1]
    mat = np.zeros([n,n])
    tmp = n - np.cos(X).sum() + np.arange(1,n + 1)*(1 - np.cos(X)) - np.sin(X)
    dia = 2*np.cos(X)*tmp.sum() + 2*np.sin(X)*((np.arange(1,n + 1) + n)*np.sin(X) - np.cos(X)) +\
    2*((np.arange(1,n + 1) + 1)*np.sin(X) - np.cos(X))*(np.arange(1,n + 1)*np.sin(X) - np.cos(X)) + \
    + 2*tmp*(np.sin(X) + np.arange(1,n + 1)*np.cos(X))
    mat = 2*np.sin(X).T*((np.arange(1,n + 1) + n)*np.sin(X) - np.cos(X)) + \
    (np.arange(1,n + 1)*np.sin(X) - np.cos(X)).T*2*np.sin(X)
    mat.flat[::(n + 1)] = dia
    return mat
def GS(HH,GG,X,eta,N):
    R = HH(X);b = -GG(X).T
    n = X.shape[1]
    nu = max(abs(np.diag(R))) + 1
    A = R + nu*np.eye(n,n)
    b_norm = (b**2).mean()
    x0 = b
    x = x0.copy()
    for k in range(N):
        x[0,0] = (b[0,0] - (A[0:1,1:n]@x0[1:n,0]).item())/A[0,0]
        for i in range(1,n):
            x[i,0] = (b[i,0] - (A[i:i + 1,0:i]@x[0:i]).item() - (A[i: i + 1,i + 1:n]@x0[i + 1:n]).item())/A[i,i]
        x0 = x
        r_norm = ((R@x - b)**2).mean()
        if r_norm < eta*b_norm:
            break
    #print(k)
    return x
def INewton(FF,GG,HH,X,N):
    tic = time.time()
    eps = 1e-5
    fun_iter = 0
    gama = 0.2
    eta = np.zeros(2)
    X_old = X.copy()
    for i in range(N):
        if i == 0:
            eta[0] = 1e-2
            d = GS(HH,GG,X_old,eta[0],100)
            fun_iter += 1
            X = X_old + d.T
            eta[1] = gama*np.linalg.norm(X)/np.linalg.norm(X_old)
            X_old = X
        else:
            if gama*eta[0] > 0.1:
                eta[1] = max(eta[1],gama*eta[0])
                d = GS(HH,GG,X_old,eta[1],100)
                fun_iter += 1
                X = X_old + d.T
                eta[0] = eta[1]
                eta[1] = gama*np.linalg.norm(X)/np.linalg.norm(X_old)
                X_old = X
            else:
                d = GS(HH,GG,X_old,eta[1],100)
                fun_iter += 1
                X = X_old + d.T
                eta[0] = eta[1]
                eta[1] = gama*np.linalg.norm(X)/np.linalg.norm(X_old)
                X_old = X
            
        if (i + 1)%50 == 0:
            print('the iteration:%d'%(i + 1))
        tmp_lf = (GG(X)**2).sum()
        tmp_ri = max(1,(X**2).sum())*eps
        if tmp_lf < tmp_ri:
            break
    ela = time.time() - tic
    print('iteration:%d,grad_2:%.3e,fun_iter = %d,time = %.2f'%(i + 1,(GG(X)**2).sum(),fun_iter,ela))
    print('the end:%d'%(i))
    return X

def Inewton(FF,GG,HH,X,N):
    tic = time.time()
    eps = 1e-5
    fun_iter = 0
    X_old = X.copy()
    for i in range(N):
        eta = min(1e-4,(GG(X_old)**2).sum())
        
        d = GS(HH,GG,X_old,eta,1000)
        fun_iter += 1
        X = X_old + d.T
        #print(np.linalg.norm(X - X_old))
        X_old = X
        if (i + 1)%50 == 0:
            print('the iteration:%d'%(i + 1))
        tmp_lf = (GG(X)**2).sum()
        tmp_ri = max(1,(X**2).sum())*eps
        if tmp_lf < tmp_ri:
            break
    ela = time.time() - tic
    print('iteration:%d,grad_2:%.3e,fun_iter = %d,time = %.2f'%(i + 1,(GG(X)**2).sum(),fun_iter,ela))
    print('the end:%d'%(i))
    return X
n = 1000
X = np.ones([1,n])/n
X_new = Inewton(FF,GG,HH,X,400)#这里只需替换为X_new = INewton(FF,GG,HH,X,500)就得到另一种修正
#print(X_new)
print(FF(X_new))

重点
对于不同的函数，G-S迭代跳出循环的条件和最大迭代次数需要调整，不然得到的不是下降方向。

限制内存Broyden-Fletcher-Goldfarb-Shanno算法（Limited-memory BFGS，L-BFGS）
qq_25426809的博客
05-08 608
限制内存Broyden-Fletcher-Goldfarb-Shanno算法（Limited-memory BFGS，L-BFGS）是一种用于求解大型无约束优化问题的迭代方法。它是BFGS（Broyden-Fletcher-Goldfarb-Shanno）算法的内存优化版本，适用于处理大规模问题。L-BFGS通过近似海森矩阵的逆来更新解的估计值，但与BFGS不同，它仅存储有限数量的过去梯度和更新向量，从而减少存储和计算开销。
matlab二阶锥规划,二阶锥规划的光滑非精确牛顿法 Smoothing Inexact Newton Method for the Second Order Cone Programming...
weixin_34910865的博客
03-20 1015
Advances in Applied MathematicsVol.04No.03(2015), Article ID:15898,6pages10.12677/AAM.2015.43033Smoothing InexactNewton Method for the Second Order Cone ProgrammingLi Dong1, Siqi Xu2, Jingen Yang11Co...
LBFGS机器学习优化算法
06-23
用C实现的liblbfgs的代码，有样例程序
【转】解决函数最优化问题算法之 L-BFGS算法
Toby的博客
07-22 1113
原文： 一文读懂L-BFGS算法 原文链接：https://blog.csdn.net/weixin_39445556/article/details/84502260 L-BFGS算法：有限内存中进行BFGS算法,L是limited memory的意思. BFGS是四位数学家名字的首字母缩写，Broyden,Fletcher,Goldfarb,Shanno.所以算法的名字就叫做BFGS算法 牛顿法求根问题 牛顿法求函数的驻点 牛顿法求驻点的本质 BFGS算法 L-BFGS算法 牛顿法求根问题
L-BFGS算法（有限内存BFGS算法）
kaka19880812的专栏
08-02 2810
有限内存BFGS求解泊松方程
forrestguang的博客
12-14 475
构建优化函数 利用差分法得到AX=b，求解F(x)=0.5∗xTAx−bTxF(x)=0.5*x^{T}Ax-b^{T}xF(x)=0.5∗xTAx−bTx 首先明确里面的A是对称正定矩阵，然后我们开始构造A,b。 import torch import numpy as np import time def UU(X, order,prob):#X表示(x,t) if prob==1: temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])
内存有限的问题
另一个空间
02-13 565
虽然《苹果往事》这本由外行翻译的书有很多的遗憾，我还是觉得很值得一看，当然，直接去Folklore网站看更好。 从中能窥得很多苹果的趣事，开发Macintosh的趣事。我最喜欢《We're Not Hackers!》一篇，虽然它的Overall Rating在我打了五星后还是只有3.93。 我感兴趣的是，他们总是在处理内存有限的问题。 现在很少开发人员关注这个，无论磁盘空间还是内存空间，都早
BFGS算法（拟牛顿法）.docx
02-26
拟牛顿法是一种在数值最优化领域广泛应用的迭代方法，主要用来寻找函数的局部极小值。...然而，对于大规模问题，存储和更新Hessian近似矩阵可能成为瓶颈，这时可以考虑使用更节省存储的L-BFGS（有限内存BFGS）算法。
L-BFGS算法/Broyden族/BFGS算法/阻尼牛顿法的Python实现代码
sinat_28371057的博客
03-21 5717
下面定义了三个Python语言编写的函数：函数表达式fun，梯度向量gfun，和海森矩阵hess。这三个表达式在后面各个算法的实现中会用到。 # 函数表达式fun fun = lambda x:100*(x[0]**2-x[1])**2 + (x[0]-1)**2 # 梯度向量 gfun gfun = lambda x:np.array([400*x[0]*(x[0]**2-x[1])+2*(x[0]-1), -200*(x[0]**2-x[1])]) # 海森矩阵 hess hess = lamb
基于LR的优化方法：梯度下降法，随机梯度下降法，牛顿法，LBFGS，BFGS.zip
最新发布
08-25
本资料包主要探讨了五种常见的优化方法，分别是梯度下降法、随机梯度下降法、牛顿法以及两种有限内存的拟牛顿法——LBFGS和BFGS。这些方法在逻辑回归（Logistic Regression, LR）和其他复杂模型如神经网络中广泛使用...
论文研究-求解非线性方程组的改进不精确雅可比牛顿法.pdf
09-07
在分析不精确雅可比牛顿法的基础上，进一步研究了不精确雅可比矩阵在精确解附近奇异的求解方法。利用雅可比矩阵与函数自身，在不增加新的计算量前提下，得到改进的求解非线性方程组的不精确雅可比牛顿算法。数值结果表明，改进后算法与原不精确雅可比牛顿法具有相同的计算效率，而且在使用上更为方便，有效。
从牛顿法到L-BFGS的算法演变
Eason的博客
09-10 4360
前言(本文主要学习自该博主的文章：http://blog.csdn.net/itplus，以下是本人的笔记，主要记录了结论部分，省略了推导的部分。对具体推导过程有兴趣的同学请访问原博主的博客~)   拟牛顿法是求解非线性优化问题最有效的方法之一，其中DFP方法，BFGS方法以及L-BFGS方法都是重要的拟牛顿法。我们现在考虑如下无约束的极小化问题： >minxf(x),其中x=(x1,x2
BFGS算法：高效全局优化与拟牛顿法详解
BFGS算法在处理非线性最小化问题时展现出了优势，尤其是在梯度下降法和牛顿法之间找到了平衡。虽然它在计算效率上不如梯度下降法那样简单，但相比牛顿法，其内存消耗和计算负担更小，适合处理中等规模的数据。然而，...
L-BFGS-B(Limited-memory Broyden–Fletcher–Goldfarb–Shanno )算法理解_附代码
没有最好，只有更好
04-27 8924
L-BFGS-B(Limited-memory Broyden–Fletcher–Goldfarb–Shanno )算法理解
各种牛顿法、拟牛顿法
热门推荐
nameofcsdn的博客
05-23 1万+
...
一文读懂L-BFGS算法
thedha3ege的博客
06-14 967
本章我们来学习L-BFGS算法.L-BFGS是机器学习中解决函数最优化问题比较常用的手段,本文主要包括以下六部分: 1-L-BFGS算法简介 2-牛顿法求根问题 3-牛顿法求函数的驻点 4-牛顿法求驻点的本质 5-BFGS算法 ...
大规模无约束优化
Fullofstars
04-18 418
大规模无约束优化方法
数值优化（Numerical Optimization）学习系列-大规模无约束最优化（Large-Scale Unconstrained Optimization）
fangqingan_java的专栏
12-27 5102
概述 当最优化问题参数个数增加，求解问题所需要的时间和空间复杂度会增加。计算时间和空间是一个权衡，只需要存储一阶梯度时，时间复杂度可能为超线性；如果利用Hessian矩阵可以达到二次收敛，但是需要o(n2)o(n^2)的空间复杂度。 另外对于拟牛顿算法所得到的Hessian矩阵式稠密的，不能利用到稀疏矩阵的一些性质。针对以上问题本小节给出求解大规模无约束最优化问题的一些思路，主要内容
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

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

打赏作者

谨慎付费（看不懂试读博客不要订阅）

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

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

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

打赏作者

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

抵扣说明：

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

余额充值