深入淺出最優(yōu)化(8) 拉格朗日乘子法

1 拉格朗日乘子法的數(shù)學(xué)背景

當(dāng)使用前面介紹的罰函數(shù)法求解約束問題時，為獲得足夠好的近似解含末，罰參數(shù)需取足夠大的值猜拾，這將導(dǎo)致增廣目標(biāo)函數(shù)的黑森矩陣出現(xiàn)病態(tài)，從而導(dǎo)致數(shù)值計算上的困難佣盒。因此提出拉格朗日乘子法挎袜。

學(xué)高數(shù)的時候我們就學(xué)過等式約束條件下的拉格朗日乘子法。延續(xù)前一節(jié)中對約束最優(yōu)化問題的定義沼撕，則拉格朗日函數(shù)為 $L(x,\mu)=f(x)-\displaystyle\sum_{j\in E}\mu_jh_j(x)$ 宋雏。局部最優(yōu)解的條件是對x求梯度及對所有 $\mu$ 求導(dǎo)均為0。

下面來討論不等式和等式同時約束條件下的拉格朗日乘子法务豺。拉格朗日函數(shù)為 $L(x,\mu)=f(x)-\displaystyle\sum_{i\in I}\lambda_ig_i(x)-\displaystyle\sum_{j\in E}\mu_jh_j(x)$ 磨总。對 $\lambda$ 的情況需要討論：

最優(yōu)解在 $D$ 內(nèi)部時， $g(x^*)>0$ 笼沥，約束條件無效蚪燕，有 $\lambda^*_i=0,i\in I$
最優(yōu)解在 $D$ 邊界上時娶牌， $g(x^*)=0$ ，約束條件有效馆纳。此時對x求梯度設(shè)為0诗良，有 $\nabla f(x^*)=\displaystyle\sum_{i\in I}\lambda_i^*\nabla g_i(x^*)+\sum_{j\in E}\mu^*_j\nabla h_i(x^*)$ , $\nabla f(x^*)$ 指向 $D$ 內(nèi)部（因為邊界處是最優(yōu)解）， $\nabla g_i(x^*)$ 也指向 $D$ 內(nèi)部（因為邊界處為0鲁驶， $D$ 內(nèi)大于0）鉴裹，所以 $\lambda^*_i>0,i\in I$

可見，無論最優(yōu)解在何處钥弯，都有 $\lambda_i^*g_i(x^*)=0,i\in I$ 径荔，這個條件被稱為互補松弛條件。

因此針對不等式和等式同時約束條件下的拉格朗日乘子法脆霎，局部最優(yōu)解的條件可以用如下的KKT條件表示：

$\begin{cases}\nabla_xL(x^*,\lambda^*,\mu^*)=\nabla f(x^*)-\displaystyle\sum_{i\in I}\lambda_i^*\nabla g_i(x^*)-\sum_{j\in E}\mu_j^*\nabla h_j(x^*)=0\\h_j(x^*)=0,j\in E\\g_i(x^*)\geq0,\lambda_i^*\geq0,\lambda_i^*g_i(x^*)=0,i\in I\end{cases}$

滿足KKT條件的點被稱為KKT點总处。

2 拉格朗日乘子法的構(gòu)成

2.1 等式約束問題的乘子法

將外點罰函數(shù)法的思想引入拉格朗日函數(shù)的最優(yōu)化問題。設(shè) $S(x)=\displaystyle\sum_{i\in E}h_i^2(x)$ 睛蛛，構(gòu)造輔助函數(shù) $L_\mu(x,\lambda)=L(x,\lambda)+\frac{1}{2}\mu S(x)=f(x)-\displaystyle\sum_{i\in E}\lambda_ih_i(x)+\frac{1}{2}\mu\sum_{i\in E}h_i^2(x)$

求駐點鹦马，令 $\nabla_xL_\mu(x,\lambda)=\nabla f(x)-\displaystyle\sum_{i\in E}\lambda_i\nabla h_i(x)+\mu\sum_{i\in E}h_i(x)\nabla h_i(x)=0$

根據(jù)KKT條件有 $\nabla f(x^*)-\displaystyle\sum_{i\in E}\lambda_i^*\nabla h_i(x^*)=0$

若 $\nabla_xL_{\mu_k}(x_k,\lambda_k)=\nabla f(x_k)-\displaystyle\sum_{i\in E}\lambda_i^{(k)}\nabla h_i(x_k)+\mu_k\sum_{i\in E}h_i(x_k)\nabla h_i(x_k)=0$

則 $\nabla f(x_k)-\displaystyle\sum_{i\in E}[\lambda_i^{(k)}-\mu_kh_i(x_k)]\nabla h_i(x_k)=0$

為了使得 $\lambda$ 收斂向 $\lambda^*$ ，有 $\lambda_i^{(k+1)}=\lambda_i^{(k)}-\mu_k h_i(x_k),i\in E$
為了使得 $h_j(x)→0,j\in E$ 忆肾，每一步增加罰因子 $\mu_{k+1}=\sigma\mu_k$ （ $\sigma$ 為放大系數(shù)）

由此得到等式約束問題乘子法的步驟：

選定初始點 $x_0\in R^n$ 荸频，初始乘子估計 $\lambda_1$ ，初始罰因子 $\mu_1>0$ 难菌，常數(shù) $\sigma>1$ 试溯， $\beta\in(0,1)$ ，精度 $\epsilon>0$ 郊酒，置 $k=1$
構(gòu)造增廣目標(biāo)函數(shù) $L_{\mu_k}(x,\lambda_k)=L(x,\lambda_k)+\frac{1}{2}\mu_kS(x)$
以 $x_{k-1}$ 為初始點求解無約束問題 $minL_{\mu_{k}}(x,\lambda_k)$ 遇绞，其解為 $x_{k}$
若 $S(x_k)^{\frac{1}{2}}\leq\epsilon$ ，則得解 $x_k$ 燎窘，停止迭代
若 $\frac{S(x_k)^{\frac{1}{2}}}{S(x_{k-1})^{\frac{1}{2}}}\leq\beta$ 成立摹闽，則令 $\mu_{k+1}=\mu_k$ ，否則 $\mu_{k+1}=\sigma\mu_k$
令 $\lambda_i^{(k+1)}=\lambda_i^{(k)}-\mu_kh_i(x_k),i\in E$ 褐健，令 $k=k+1$ 付鹿，轉(zhuǎn)步2

2.2 一般約束問題的乘子法

引入松弛變量z將不等式問題化為等價的等式約束問題： $g_i(x)\geq 0,i\in I\Rightarrow g_i(x)-z_i^2=0,i\in I$

構(gòu)造增廣拉格朗日函數(shù)： $\overline{L_\mu}(x,z,\lambda)=f(x)-\displaystyle\sum_{i\in I}\lambda_i[g_i(x)-z_i^2]+\frac{\mu}{2}\sum_{i\in I}[g_i(x)-z_i^2]^2$

配方后可得： $\overline{L_\mu}(x,z,\lambda)=f(x)-\displaystyle\sum_{i\in I}\{\frac{\mu}{2}[z_i^2-\frac{1}{\mu}(\mu g_i(x)-\lambda_i)]^2-\frac{\lambda^2_i}{2\mu}\}$

將該式對 $z_i$ 求偏導(dǎo)，令偏導(dǎo)為0蚜迅，有 $2z_i\{\lambda_i-\mu[g_i(x)-z_i^2]\}=0$

因此 $z_i^2=\frac{1}{\mu}max\{0,\mu g_i(x)-\lambda_i\}$ 舵匾，代入消去z得 $L_\mu(x,\lambda)=f(x)+\frac{1}{2\mu}\displaystyle\sum_{i\in I}(max^2\{0,\lambda_i-\mu g_i(x)\}-\lambda_i^2)$

則 $\nabla_xL(x,\lambda)=\nabla f(x)-\displaystyle\sum_{i\in I}(max\{0,\lambda_i-\mu g_i(x)\})\nabla g_i(x)$

根據(jù)KKT條件有 $\nabla f(x^*)-\displaystyle\sum_{i\in E}\lambda_i^*\nabla g_i(x^*)=0$

為了使得 $\lambda$ 收斂向 $\lambda^*$ ，根據(jù)KKT條件谁不，有 $\lambda_i^{(k+1)}=max\{\lambda_i^{(k)}-\mu g_i(x_k),0\},i\in E$

可以得到一般約束條件下的增廣目標(biāo)函數(shù)和 $\lambda$ 的遞推公式：

$L_\mu(x,\lambda)=f(x)+\frac{1}{2\mu}\displaystyle\sum_{i\in I}(max^2\{0,\lambda_i-\mu g_i(x)\}-\lambda_i^2)-\displaystyle\sum_{j\in E}\lambda_jh_j(x)+\frac{1}{2}\mu\sum_{j\in E}h_j^2(x)$

$\begin{cases}\lambda_i^{(k+1)}=max\{\lambda_i^{(k)}-\mu g_i(x_k),0\},i\in I\\\lambda_j^{(k+1)}=\lambda_j^{(k)}-\mu h_j(x_k),j\in E\end{cases}$

由此得到一般約束問題乘子法的步驟：

選定初始點 $x_0\in R^n$ 坐梯，初始乘子估計 $\lambda_1$ ，初始罰因子 $\mu_1>0$ 刹帕，常數(shù) $\sigma>1$ 吵血， $\beta\in(0,1)$ 谎替，精度 $\epsilon>0$ ，置 $k=1$
構(gòu)造增廣目標(biāo)函數(shù)
以 $x_{k-1}$ 為初始點求解無約束問題 $minL_{\mu_{k}}(x,\lambda_k)$ 蹋辅，其解為 $x_{k}$
若 $(\displaystyle\sum_{j\in E}h_j^2(x_k))^{\frac{1}{2}}+(\displaystyle\sum_{i\in I}min^2\{g_i(x_k),0\})^{\frac{1}{2}}\leq\epsilon$ 钱贯，則得解 $x_k$ ，停止迭代
若 $\frac{(\displaystyle\sum_{j\in E}h_j^2(x_k))^{\frac{1}{2}}+(\displaystyle\sum_{i\in I}min^2\{g_i(x_k),0\})^{\frac{1}{2}}}{(\displaystyle\sum_{j\in E}h_j^2(x_{k-1}))^{\frac{1}{2}}+(\displaystyle\sum_{i\in I}min^2\{g_i(x_{k-1}),0\})^{\frac{1}{2}}}\leq\beta$ 成立侦另，則令 $\mu_{k+1}=\mu_k$ 秩命，否則 $\mu_{k+1}=\sigma\mu_k$
計算 $\lambda^{(k+1)}$ ，令 $k=k+1$ 褒傅，轉(zhuǎn)步2

3 實戰(zhàn)測試

對于上節(jié)深入淺出最優(yōu)化(7) 罰函數(shù)法中提出的約束最優(yōu)化問題硫麻， $x_1,x_2,x_3$ 的初值均在 $[0,4]$ 的范圍內(nèi)隨機生成，總共生成100組起點樊卓。統(tǒng)計迭代成功（在1000步內(nèi)得到最優(yōu)解且單次步長搜索迭代次數(shù)不超過1000次）的樣本的平均迭代步數(shù)、平均迭代時間和得到的最優(yōu)解及開銷函數(shù)最小值杠河。

迭代步數(shù)	迭代時間	最優(yōu)解	函數(shù)最小值
17	24.3729s	$x_1=1.1051~x_2=1.1969~x_3=1.5352$	$0.03257$

代碼實現(xiàn)

本博客所有代碼在https://github.com/HarmoniaLeo/optimization-in-a-nutshell開源碌尔，如果幫助到你，請點個star券敌，謝謝這對我真的很重要唾戚！

你可以在上面的GitHub鏈接或本文集的第一篇文章深入淺出最優(yōu)化(1) 最優(yōu)化問題概念與基本知識中找到Function.py和lagb.py

使用共軛梯度PRP+法的拉格朗日乘子法

import numpy as np
from Function import Function   #定義法求導(dǎo)工具
from lagb import *  #線性代數(shù)工具庫
from scipy import linalg

n=3 #x的長度
mu=2
lj=np.ones(1)   #λj初值，長度等于等式限制條件個數(shù)
li=np.ones(2)   #λi初值待诅，長度等于不等式限制條件個數(shù)

def func(x):    #目標(biāo)函數(shù)叹坦，x是一個包含所有參數(shù)的列表
    return (x[0]-1)**2+(x[0]-x[1])**2+(x[1]-x[2])**2

def hj(x):  #構(gòu)造數(shù)組h，第j位是第j+1個等式限制條件計算的值卑雁，x是一個包含所有參數(shù)的列表
    return np.array([x[0]*(1+x[1]**2)+x[2]**4-4-3*np.sqrt(2)])

def gi(x):  #構(gòu)造數(shù)組g募书，第i位是第i+1個不等式限制條件計算的值，x是一個包含所有參數(shù)的列表
    return np.array([x[0]+10,-x[0]+10])

def myFunc(x):
    return  func(x)+\
        1/(2*mu)*np.sum(np.power(np.where(li-mu*gi(x)>0,li-mu*gi(x),0),2)-np.power(li,2))-\
            np.sum(lj*hj(x))+0.5*mu*np.sum(np.power(hj(x),2))

def cdt(x):
    return np.sqrt(np.sum(np.power(hj(x),2)))+np.sqrt(np.sum(np.power(np.where(gi(x)>0,0,gi(x)),2)))

sigma2=1.5  #放大因子
e2=0.001
beta=0.5
x=np.array([2.0,2.0,2.0])   #初值點
k1=0
while cdt(x)>=e2:
    e=0.001
    beta1=1
    sigma=0.4
    rho=0.55
    tar=Function(myFunc)
    k=0
    d=-tar.grad(x)
    x1=x
    while tar.norm(x)>e:
        a=1
        if not (tar.value(x+a*d)<=tar.value(x)+rho*a*dot(turn(tar.grad(x)),d) and \
            np.abs(dot(turn(tar.grad(x+a*d)),d))>=sigma*dot(turn(tar.grad(x)),d)):
            a=beta1
            while tar.value(x+a*d)>tar.value(x)+rho*a*dot(turn(tar.grad(x)),d):
                a*=rho
            while np.abs(dot(turn(tar.grad(x+a*d)),d))<sigma*dot(turn(tar.grad(x)),d):
                a1=a/rho
                da=a1-a
                while tar.value(x+(a+da)*d)>tar.value(x)+rho*(a+da)*dot(turn(tar.grad(x)),d):
                    da*=rho
                a+=da
        lx=x
        x=x+a*d
        beta=np.max((dot(turn(tar.grad(x)),tar.grad(x)-tar.grad(lx))/(tar.norm(lx)**2),0))  #PRP+
        d=-tar.grad(x)+beta*d
        k+=1
        print(k1,k)
    if cdt(x)/cdt(x1)>beta:
        mu*=sigma2
    k1+=1
print(x)