pytorch如何求解非线性优化 python非线性约束优化

目录

• 迭代法求零点

• 基本思想
• 具体做法
• 几何含义
• 重要定理

• 迭代法求解无约束优化问题

• 1. 最速下降 (SD) 法 (负梯度方法)

• 梯度和 Hesse 矩阵
• SD 法
• 一维精确线搜索
• Python 实现

• 2. Newton 法

无约束优化问题就是没有任何的约束限制的优化问题, 如求最小值 $\min \ f(x)$, 其中 $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$. 求解无约束优化问题的迭代算法有最速下降 (SD) 法和 Newton 法等.

迭代法求零点

基本思想

$f (x) = 0 \xLeftrightarrow{\text{等价变换}} x=\varphi(x) \ (\text { 迭代函数 })$

$f(x) \text { 的零点 } x^{*} \Longleftrightarrow \varphi(x) \text { 的不动点 } x^{*}$

不动点迭代: $x_{k+1}=\varphi\left(x_{k}\right)$

具体做法

从一个给定的初值 $x_{0}$ 出发, 计算 $x_{1}=\varphi\left(x_{0}\right), x_{2}=\varphi\left(x_{1}\right), \ldots,$ 若 $\left\{x_{k}\right\}_{k=0}^{\infty}$ 收玫, 即存在 $x^{*}$ 使得 $\lim_{k \rightarrow \infty} x_{k}=x^{*}$, 则由 $\varphi$ 的连续性和 $\lim_{k \rightarrow \infty} x_{k+1}=\lim_{k \rightarrow \infty} \varphi\left(x_{k}\right)$ 可得 $x^{*}=\varphi\left(x^{*}\right)$, 即 $x^{*}$ 是 $\varphi$ 的不 动点, 也就是 $f(x)$

几何含义

求曲线 $y=\varphi(x)$ 与直线 $y=x$

重要定理

定理 (不动点原理, 压缩映像定理) 设 $\varphi(x)$ 在 $[a, b]$ 上连续, 且一阶导数连续, 若:
(1) $a \leq \varphi(x) \leq b$ 对 一切 $x \in[a, b]$ 都成立; (封闭性)
(2) 存在 $0 \leq L<1$, 使得 $\left|\varphi^{\prime}(x)\right| \leq L$ 对 $\forall x \in[a, b]$ 成立, (压缩性)
则函数 $f(x)=x-\varphi(x)$ 在 $[a, b]$ 中有唯一的零点 $x^{*}$, 即 $\varphi(x)$ 存在唯一的的不动点 $x^{*}$, s.t. $\varphi\left(x^{*}\right)=x^{*}$.

证明思路 : 存在性: $f(a)=a-\varphi(a) \leq 0, \ f(b)=b-\varphi(b) \geq 0$ $\Longrightarrow$ $f(x)$ 在 $[a, b]$ 上有零点.
唯一性: 用反证法, 假设存在 $x^{*}, y^{*} \in[a, b]$

$\left\{\begin{array}{l}x^{*}=\varphi\left(x^{*}\right) \\ y^{*}=\varphi\left(y^{*}\right)\end{array}\right.$

$\Longrightarrow$ $\left|\boldsymbol{x}^{*}-\boldsymbol{y}^{*}\right|=\left|\varphi\left(\boldsymbol{x}^{*}\right)-\varphi\left(\boldsymbol{y}^{*}\right)\right|=\left|\varphi^{\prime}(\xi)\right| \cdot\left|\boldsymbol{x}^{*}-\boldsymbol{y}^{*}\right| \leq {L}\left|\boldsymbol{x} *-\boldsymbol{y}^{*}\right|$, 推出矛盾.

迭代法求解无约束优化问题

无约束优化问题就是没有任何的约束限制的优化问题, 如求最小值 $\min \ f(x)$, 其中 $f: \mathbb{R}^{n} \rightarrow \mathbb{R}$.

1. 最速下降 (SD) 法 (负梯度方法)

梯度和 Hesse 矩阵

$\nabla f(x)=\left[\begin{array}{llll} \frac{\partial f(x)}{\partial x_{1}} & \frac{\partial f(x)}{\partial x_{2}} & \cdots & \frac{\partial f(x)}{\partial x_{n}} \end{array}\right]^{\mathrm{T}}=g(x)$

$\nabla^{2} f(x)=\left[\begin{matrix} \frac{\partial^{2} f(x)}{\partial x_{1}^{2}} & \frac{\partial^{2} f(x)}{\partial x_{1} \partial x_{2}} & \cdots & \frac{\partial^{2} f(x)}{\partial x_{1} \partial x_{n}} \\ \frac{\partial^{2} f(x)}{\partial x_{2} \partial x_{1}} & \frac{\partial^{2} f(x)}{\partial x_{2}^{2}} & \cdots & \frac{\partial^{2} f(x)}{\partial x_{2} \partial x_{n}} \\ \vdots & \vdots & & \vdots \\ \frac{\partial^{2} f(x)}{\partial x_{n} \partial x_{1}} & \frac{\partial^{2} f(x)}{\partial x_{n} \partial x_{2}} & \cdots & \frac{\partial^{2} f(x)}{\partial x_{n}^{2}} \end{matrix}\right]=G(x)$

SD 法

步 1: 给出 $x_{0} \in \mathbb{R}^{n}, \varepsilon>0, k:=0$;
步 2: 若满足终止条件 $\left\|\nabla f\left(x_{k}\right)\right\| \leq \varepsilon$, 则迭代停止;
步 3: 计算 $d_{k}=-g_{k}$ ;
步 4: 一维精确线搜索求 $\alpha_{k}$;
步 5: $x_{k+1}:=x_{k}+\alpha_{k} d_{k}, k:=k+1$, 转步 2.

其中一维精确线搜索确定最优步长的方法如下

一维精确线搜索

$\min _{\alpha>0} f\left(x_{k}-\alpha g_{k}\right)$

其中
$\begin{aligned} f\left(x_{k}-\alpha g_{k}\right) &=\frac{1}{2}\left(x_{k}-\alpha g_{k}\right)^{\mathrm{T}} G\left(x_{k}-\alpha g_{k}\right)+b^{\mathrm{T}}\left(x_{k}-\alpha g_{k}\right) \\ &=\frac{1}{2} g_{k}^{\mathrm{T}} G g_{k} \alpha^{2}-g_{k}^{\mathrm{T}} g_{k} \alpha+f\left(x_{k}\right) \end{aligned}$

得最速下降方法的步长
$\alpha_{k}=\frac{g_{k}^{\mathrm{T}} g_{k}}{g_{k}^{\mathrm{T}} G g_{k}}$

Python 实现

#最速下降法 (负梯度方法)
import numpy as np #导入numpy模块
import matplotlib.pyplot as plt #导入matplotlib.pyplot模块
from sympy import*
from scipy.interpolate import griddata
import numpy.linalg as LA

#目标函数
def nf(x1,x2):
    y=1/2*(21*x1*x1+4*x2*x1+4*x1*x2+15*x2*x2)+2*x1+3*x2+10
    return y

#目标函数求梯度
def ndfx(x1,x2):
    fx=np.array([diff(nf(x1,x2),x1),diff(nf(x1,x2),x2)])
    #分别就目标函数对x1，x2一阶偏导, array数组
    # sympy.diff(func,x,n) 求导
    # func是要求导的函数; x是要对其求导的变量; n是求导阶数, 缺省为1(可选)
    return fx

#初始化参数
n=100   #设置最大迭代次数
x=np.zeros((n,2)) #迭代点xk
x[0,:]=[-30,100]
xk=[]
xk.append(x[0,:])
d=np.zeros((n,2)) #下降方向dk
df=np.zeros((n,2)) #梯度, 建立array数组的形式 

#--------------------------------------

def main():
    x1,x2=var('x1 x2') # sympy.var 将字符串变成变量
    eps=1e-5 #精度要求
    k=1
    
    #设置初始迭代点: 
    A=np.arange(4*n).reshape(n,2,2)
    fx=ndfx(x1,x2)#梯度

    alpha=0.1
    for i in range(n-2):
        #二阶导数矩阵
        #计算一阶导数值
        f1x=fx[0].subs({x1:x[i][0],x2:x[i][1]})#截取xk
        f2x=fx[1].subs({x1:x[i][0],x2:x[i][1]})

        df[i,:]=[f1x,f2x]
        b=-df[i,:].T#负梯度方向
        #subs函数只能应用于sympy类型的符号函数，并非array数组类型
        A[i,:,:]=[[21,4],[4,15]]#SD方法  
        alpha=(b*b)/(21*b[0]*b[0]+4*b[1]*b[0]+4*b[0]*b[1]+15*b[1]*b[1])#步长公式
        d[i,:]=alpha*b#负梯度方向
        x[i+1,:]=x[i,:]+d[i,:]
        xk.append(x[i+1,:])#存储迭代点
        if np.abs(x[i+1][0]-x[i][0])<eps and np.abs(x[i+1][1]-x[i][1])<eps :
             break   
        i=i+1
        print(nf(x[i][0],x[i][1]))
    print('-----------------------------------')
    print('最后迭代 The result is:',i)
    print('最优解:x*=',x[i,:])
    print('最优目标值:y=%.3f'%(   nf(x[i][0],x[i][1])  )  )   

#--------------------------------------

main()


plt.rc('font',size=16);plt.rc('font',family='SimHei');plt.rc('text',usetex=True)

xh0=np.arange(-100,100.1,0.1)#区域横坐标
yh0=np.arange(-100,100.1,0.1)#区域纵坐标   
xh,yh=np.meshgrid(xh0,yh0)
contr=plt.contour(xh,yh,nf(xh,yh),20)
plt.clabel(contr)#等高线
plt.xlabel('$x_1$');plt.ylabel('$x_2$',rotation=90)
xk=np.reshape(xk,(1,-1))
xk=np.reshape(xk,(int(xk.size/2),2))
xtk=np.array(xk.T[0,:])
ytk=np.array(xk.T[1,:])
#print(xtk)
#print(ytk)
plt.plot(xtk,ytk,c='r',marker='H')
plt.show()

2. Newton 法

步 1: 给出 $x_{0} \in \mathbb{R}^{n}, \varepsilon>0, k:=0$;
步 2: 若满足终止准则 $\left\|\nabla f\left(x_{k}\right)\right\| \leq \varepsilon$, 则输出有关信息, 停止迭代;
步 3: 由 $G_{k} d=-g_{k}$ $\Rightarrow$ $d_{k}=-G_{k}^{-1} g_{k}$;
步 4: $x_{k+1}:=x_{k}+d_{k}, k:=k+1$, 转步 2.

# Newton 法
import numpy as np #导入numpy模块
import matplotlib.pyplot as plt #导入matplotlib.pyplot模块
from sympy import*
from scipy.interpolate import griddata
import numpy.linalg as LA

#目标函数
def nf(x1,x2):
    y=3*x1**2+3*x2**2-x1**2*x2
    #y=1/2*(21*x1*x1+4*x2*x1+4*x1*x2+15*x2*x2)+2*x1+3*x2+10
    return y

#目标函数求梯度 (一阶偏导)
def ndfx(x1,x2):
    fx=np.zeros(2)#二维向量
    fx=[diff(nf(x1,x2),x1),diff(nf(x1,x2),x2)]#分别就目标函数对x1，x2一阶偏导，列表类型
    fx=np.array(fx)#转换为array数组类型
    return fx

#目标函数求 Hesse 矩阵 (二阶偏导)
def ndfxx(x1,x2):
    fx=ndfx(x1,x2)#梯度
    fxx=[[],[]]
    fxx[0]=[diff(fx[0],x1),diff(fx[0],x2)]
    fxx[1]=[diff(fx[1],x1),diff(fx[1],x2)]
    fxx=np.array(fxx) #转换为数组类型
    return fxx

#初始化参数
n=100   #设置最大迭代次数
x=np.zeros((n,2))#迭代点xk
x[0,:]=[-2,4]#设置初始迭代点
d=np.zeros((n,2))#下降方向dk
df=np.zeros((n,2))#梯度，建立array数组的形式
A=np.arange(4*n).reshape(n,2,2)#建立array数组的形式
xk=[]
xk.append(x[0,:])

#--------------------------------------

def main():
    x1,x2=var('x1 x2') #转化为sympy变量
    eps=1e-6 #精度要求
    fx=ndfx(x1,x2)#梯度
    z1=ndfxx(x1,x2)#hesse矩阵
    for i in range(n-2):
        #计算一阶导数值
        f1x=fx[0].subs({x1:x[i][0],x2:x[i][1]})#梯度在xk点值
        f2x=fx[1].subs({x1:x[i][0],x2:x[i][1]})
        df[i,:]=[f1x,f2x]
        #subs函数只能应用于sympy类型的符号函数, 并非array数组类型
        #为二阶导数矩阵每个元素循环赋值
        z2=np.arange(4).reshape(2,2)#建立array数组的形式
        for j in range(2):
            for k in range(2):
                z2[j][k]=z1[j][k].subs({x1:x[i][0],x2:x[i][1]})#Hesse矩阵不同分量取值
        #A[i,:,:]=z1 #Gk Hesse矩阵在xk点处值
        #----------------------------------
        #d[i,:]=-np.dot(LA.inv(A[i,:,:]),df[i,:]) #G*d=-g
        d[i,:]=-np.dot(LA.inv(z2[:,:]),df[i,:]) #G*d=-g 改成z2形式之后Newton有二阶速度
        #基本newton
        alpha=1.0 #全步长=1              
        x[i+1,:]=x[i,:]+alpha*d[i,:]
        xk.append(x[i+1,:])
        if LA.norm(df[i,:])<eps :
            break
        i=i+1
    print('-----------------------------------')
    print('最后迭代 The result is:',i)
    print('最优解:x*=',x[i,:])
    print('最优目标值:y=%.3f'%(   nf(x[i][0],x[i][1])  )  )      
        
#--------------------------------------

if __name__ == '__main__':
    main()

plt.rc('font',size=16);plt.rc('font',family='SimHei');plt.rc('text',usetex=True)
xh0=np.arange(-6,6.01,0.01)#区域横坐标
yh0=np.arange(-6,6.01,0.01)#区域纵坐标     
xh,yh=np.meshgrid(xh0,yh0)
contr=plt.contour(xh,yh,nf(xh,yh),20)
plt.clabel(contr)#等高线
plt.xlabel('$x1$')#x,y坐标轴
plt.ylabel('$x2$',rotation=90)
xk=np.reshape(xk,(1,-1))
xk=np.reshape(xk,(int(xk.size/2),2))
xtk=np.array(xk.T[0,:])
ytk=np.array(xk.T[1,:])
#print(xtk)
#print(ytk)
plt.plot(xtk,ytk,c='r',marker='H')
plt.show()