1 拉格朗日乘子法的数学背景
当使用前面介绍的罚函数法求解约束问题时,为获得足够好的近似解,罚参数需取足够大的值,这将导致增广目标函数的黑森矩阵出现病态,从而导致数值计算上的困难。因此提出拉格朗日乘子法。
学高数的时候我们就学过等式约束条件下的拉格朗日乘子法。延续前一节中对约束最优化问题的定义,则拉格朗日函数为 L ( x , μ ) = f ( x ) − ∑ j ∈ E μ j h j ( x ) L(x,\mu)=f(x)-\displaystyle\sum_{j\in E}\mu_jh_j(x) L(x,μ)=f(x)−j∈E∑μjhj(x)。局部最优解的条件是对x求梯度及对所有 μ \mu μ求导均为0。
下面来讨论不等式和等式同时约束条件下的拉格朗日乘子法。拉格朗日函数为 L ( x , μ ) = f ( x ) − ∑ i ∈ I λ i g i ( x ) − ∑ j ∈ E μ j h j ( x ) L(x,\mu)=f(x)-\displaystyle\sum_{i\in I}\lambda_ig_i(x)-\displaystyle\sum_{j\in E}\mu_jh_j(x) L(x,μ)=f(x)−i∈I∑λigi(x)−j∈E∑μjhj(x)。对 λ \lambda λ的情况需要讨论:
- 最优解在 D D D内部时, g ( x ∗ ) > 0 g(x^*)>0 g(x∗)>0,约束条件无效,有 λ i ∗ = 0 , i ∈ I \lambda^*_i=0,i\in I λi∗=0,i∈I
- 最优解在 D D D边界上时, g ( x ∗ ) = 0 g(x^*)=0 g(x∗)=0,约束条件有效。此时对x求梯度设为0,有 ∇ f ( x ∗ ) = ∑ i ∈ I λ i ∗ ∇ g i ( x ∗ ) + ∑ j ∈ E μ j ∗ ∇ h i ( x ∗ ) \nabla f(x^*)=\displaystyle\sum_{i\in I}\lambda_i^*\nabla g_i(x^*)+\sum_{j\in E}\mu^*_j\nabla h_i(x^*) ∇f(x∗)=i∈I∑λi∗∇gi(x∗)+j∈E∑μj∗∇hi(x∗), ∇ f ( x ∗ ) \nabla f(x^*) ∇f(x∗)指向 D D D内部(因为边界处是最优解), ∇ g i ( x ∗ ) \nabla g_i(x^*) ∇gi(x∗)也指向 D D D内部(因为边界处为0, D D D内大于0),所以 λ i ∗ > 0 , i ∈ I \lambda^*_i>0,i\in I λi∗>0,i∈I
可见,无论最优解在何处,都有 λ i ∗ g i ( x ∗ ) = 0 , i ∈ I \lambda_i^*g_i(x^*)=0,i\in I λi∗gi(x∗)=0,i∈I,这个条件被称为互补松弛条件。
因此针对不等式和等式同时约束条件下的拉格朗日乘子法,局部最优解的条件可以用如下的KKT条件表示:
{ ∇ x L ( x ∗ , λ ∗ , μ ∗ ) = ∇ f ( x ∗ ) − ∑ i ∈ I λ i ∗ ∇ g i ( x ∗ ) − ∑ j ∈ E μ j ∗ ∇ h j ( x ∗ ) = 0 h j ( x ∗ ) = 0 , j ∈ E g i ( x ∗ ) ≥ 0 , λ i ∗ ≥ 0 , λ i ∗ g i ( x ∗ ) = 0 , i ∈ I \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} ⎩ ⎨ ⎧∇xL(x∗,λ∗,μ∗)=∇f(x∗)−i∈I∑λi∗∇gi(x∗)−j∈E∑μj∗∇hj(x∗)=0hj(x∗)=0,j∈Egi(x∗)≥0,λi∗≥0,λi∗gi(x∗)=0,i∈I
满足KKT条件的点被称为KKT点。
2 拉格朗日乘子法的构成
2.1 等式约束问题的乘子法
将外点罚函数法的思想引入拉格朗日函数的最优化问题。设 S ( x ) = ∑ i ∈ E h i 2 ( x ) S(x)=\displaystyle\sum_{i\in E}h_i^2(x) S(x)=i∈E∑hi2(x),构造辅助函数 L μ ( x , λ ) = L ( x , λ ) + 1 2 μ S ( x ) = f ( x ) − ∑ i ∈ E λ i h i ( x ) + 1 2 μ ∑ i ∈ E h i 2 ( x ) 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) Lμ(x,λ)=L(x,λ)+21μS(x)=f(x)−i∈E∑λihi(x)+21μi∈E∑hi2(x)
求驻点,令 ∇ x L μ ( x , λ ) = ∇ f ( x ) − ∑ i ∈ E λ i ∇ h i ( x ) + μ ∑ i ∈ E h i ( x ) ∇ h i ( x ) = 0 \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 ∇xLμ(x,λ)=∇f(x)−i∈E∑λi∇hi(x)+μi∈E∑hi(x)∇hi(x)=0
根据KKT条件有 ∇ f ( x ∗ ) − ∑ i ∈ E λ i ∗ ∇ h i ( x ∗ ) = 0 \nabla f(x^*)-\displaystyle\sum_{i\in E}\lambda_i^*\nabla h_i(x^*)=0 ∇f(x∗)−i∈E∑λi∗∇hi(x∗)=0
若 ∇ x L μ k ( x k , λ k ) = ∇ f ( x k ) − ∑ i ∈ E λ i ( k ) ∇ h i ( x k ) + μ k ∑ i ∈ E h i ( x k ) ∇ h i ( x k ) = 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 ∇xLμk(xk,λk)=∇f(xk)−i∈E∑λi(k)∇hi(xk)+μki∈E∑hi(xk)∇hi(xk)=0
则 ∇ f ( x k ) − ∑ i ∈ E [ λ i ( k ) − μ k h i ( x k ) ] ∇ 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 ∇f(xk)−i∈E∑[λi(k)−μkhi(xk)]∇hi(xk)=0
- 为了使得 λ \lambda λ收敛向 λ ∗ \lambda^* λ∗,有 λ i ( k + 1 ) = λ i ( k ) − μ k h i ( x k ) , i ∈ E \lambda_i^{(k+1)}=\lambda_i^{(k)}-\mu_k h_i(x_k),i\in E λi(k+1)=λi(k)−μkhi(xk),i∈E
- 为了使得 h j ( x ) → 0 , j ∈ E h_j(x)→0,j\in E hj(x)→0,j∈E,每一步增加罚因子 μ k + 1 = σ μ k \mu_{k+1}=\sigma\mu_k μk+1=σμk( σ \sigma σ为放大系数)
由此得到等式约束问题乘子法的步骤:
- 选定初始点 x 0 ∈ R n x_0\in R^n x0∈Rn,初始乘子估计 λ 1 \lambda_1 λ1,初始罚因子 μ 1 > 0 \mu_1>0 μ1>0,常数 σ > 1 \sigma>1 σ>1, β ∈ ( 0 , 1 ) \beta\in(0,1) β∈(0,1),精度 ϵ > 0 \epsilon>0 ϵ>0,置 k = 1 k=1 k=1
- 构造增广目标函数 L μ k ( x , λ k ) = L ( x , λ k ) + 1 2 μ k S ( x ) L_{\mu_k}(x,\lambda_k)=L(x,\lambda_k)+\frac{1}{2}\mu_kS(x) Lμk(x,λk)=L(x,λk)+21μkS(x)
- 以 x k − 1 x_{k-1} xk−1为初始点求解无约束问题 m i n L μ k ( x , λ k ) minL_{\mu_{k}}(x,\lambda_k) minLμk(x,λk),其解为 x k x_{k} xk
- 若 S ( x k ) 1 2 ≤ ϵ S(x_k)^{\frac{1}{2}}\leq\epsilon S(xk)21≤ϵ,则得解 x k x_k xk,停止迭代
- 若 S ( x k ) 1 2 S ( x k − 1 ) 1 2 ≤ β \frac{S(x_k)^{\frac{1}{2}}}{S(x_{k-1})^{\frac{1}{2}}}\leq\beta S(xk−1)21S(xk)21≤β成立,则令 μ k + 1 = μ k \mu_{k+1}=\mu_k μk+1=μk,否则 μ k + 1 = σ μ k \mu_{k+1}=\sigma\mu_k μk+1=σμk
- 令 λ i ( k + 1 ) = λ i ( k ) − μ k h i ( x k ) , i ∈ E \lambda_i^{(k+1)}=\lambda_i^{(k)}-\mu_kh_i(x_k),i\in E λi(k+1)=λi(k)−μkhi(xk),i∈E,令 k = k + 1 k=k+1 k=k+1,转步2
2.2 一般约束问题的乘子法
引入松弛变量z将不等式问题化为等价的等式约束问题: g i ( x ) ≥ 0 , i ∈ I ⇒ g i ( x ) − z i 2 = 0 , i ∈ I g_i(x)\geq 0,i\in I\Rightarrow g_i(x)-z_i^2=0,i\in I gi(x)≥0,i∈I⇒gi(x)−zi2=0,i∈I
构造增广拉格朗日函数: L μ ‾ ( x , z , λ ) = f ( x ) − ∑ i ∈ I λ i [ g i ( x ) − z i 2 ] + μ 2 ∑ i ∈ I [ g i ( x ) − z i 2 ] 2 \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 Lμ(x,z,λ)=f(x)−i∈I∑λi[gi(x)−zi2]+2μi∈I∑[gi(x)−zi2]2
配方后可得: L μ ‾ ( x , z , λ ) = f ( x ) − ∑ i ∈ I { μ 2 [ z i 2 − 1 μ ( μ g i ( x ) − λ i ) ] 2 − λ 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}\} Lμ(x,z,λ)=f(x)−i∈I∑{2μ[zi2−μ1(μgi(x)−λi)]2−2μλi2}
将该式对 z i z_i zi求偏导,令偏导为0,有 2 z i { λ i − μ [ g i ( x ) − z i 2 ] } = 0 2z_i\{\lambda_i-\mu[g_i(x)-z_i^2]\}=0 2zi{λi−μ[gi(x)−zi2]}=0
因此 z i 2 = 1 μ m a x { 0 , μ g i ( x ) − λ i } z_i^2=\frac{1}{\mu}max\{0,\mu g_i(x)-\lambda_i\} zi2=μ1max{0,μgi(x)−λi},代入消去z得 L μ ( x , λ ) = f ( x ) + 1 2 μ ∑ i ∈ I ( m a x 2 { 0 , λ i − μ g i ( x ) } − λ i 2 ) 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) Lμ(x,λ)=f(x)+2μ1i∈I∑(max2{0,λi−μgi(x)}−λi2)
则 ∇ x L ( x , λ ) = ∇ f ( x ) − ∑ i ∈ I ( m a x { 0 , λ i − μ g i ( x ) } ) ∇ g i ( x ) \nabla_xL(x,\lambda)=\nabla f(x)-\displaystyle\sum_{i\in I}(max\{0,\lambda_i-\mu g_i(x)\})\nabla g_i(x) ∇xL(x,λ)=∇f(x)−i∈I∑(max{0,λi−μgi(x)})∇gi(x)
根据KKT条件有 ∇ f ( x ∗ ) − ∑ i ∈ E λ i ∗ ∇ g i ( x ∗ ) = 0 \nabla f(x^*)-\displaystyle\sum_{i\in E}\lambda_i^*\nabla g_i(x^*)=0 ∇f(x∗)−i∈E∑λi∗∇gi(x∗)=0
为了使得 λ \lambda λ收敛向 λ ∗ \lambda^* λ∗,根据KKT条件,有 λ i ( k + 1 ) = m a x { λ i ( k ) − μ g i ( x k ) , 0 } , i ∈ E \lambda_i^{(k+1)}=max\{\lambda_i^{(k)}-\mu g_i(x_k),0\},i\in E λi(k+1)=max{λi(k)−μgi(xk),0},i∈E
可以得到一般约束条件下的增广目标函数和 λ \lambda λ的递推公式:
L μ ( x , λ ) = f ( x ) + 1 2 μ ∑ i ∈ I ( m a x 2 { 0 , λ i − μ g i ( x ) } − λ i 2 ) − ∑ j ∈ E λ j h j ( x ) + 1 2 μ ∑ j ∈ E h j 2 ( x ) 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) Lμ(x,λ)=f(x)+2μ1i∈I∑(max2{0,λi−μgi(x)}−λi2)−j∈E∑λjhj(x)+21μj∈E∑hj2(x)
{ λ i ( k + 1 ) = m a x { λ i ( k ) − μ g i ( x k ) , 0 } , i ∈ I λ j ( k + 1 ) = λ j ( k ) − μ h j ( x k ) , j ∈ E \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} {λi(k+1)=max{λi(k)−μgi(xk),0},i∈Iλj(k+1)=λj(k)−μhj(xk),j∈E
由此得到一般约束问题乘子法的步骤:
- 选定初始点 x 0 ∈ R n x_0\in R^n x0∈Rn,初始乘子估计 λ 1 \lambda_1 λ1,初始罚因子 μ 1 > 0 \mu_1>0 μ1>0,常数 σ > 1 \sigma>1 σ>1, β ∈ ( 0 , 1 ) \beta\in(0,1) β∈(0,1),精度 ϵ > 0 \epsilon>0 ϵ>0,置 k = 1 k=1 k=1
- 构造增广目标函数
- 以 x k − 1 x_{k-1} xk−1为初始点求解无约束问题 m i n L μ k ( x , λ k ) minL_{\mu_{k}}(x,\lambda_k) minLμk(x,λk),其解为 x k x_{k} xk
- 若 ( ∑ j ∈ E h j 2 ( x k ) ) 1 2 + ( ∑ i ∈ I m i n 2 { g i ( x k ) , 0 } ) 1 2 ≤ ϵ (\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 (j∈E∑hj2(xk))21+(i∈I∑min2{gi(xk),0})21≤ϵ,则得解 x k x_k xk,停止迭代
- 若 ( ∑ j ∈ E h j 2 ( x k ) ) 1 2 + ( ∑ i ∈ I m i n 2 { g i ( x k ) , 0 } ) 1 2 ( ∑ j ∈ E h j 2 ( x k − 1 ) ) 1 2 + ( ∑ i ∈ I m i n 2 { g i ( x k − 1 ) , 0 } ) 1 2 ≤ β \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 (j∈E∑hj2(xk−1))21+(i∈I∑min2{gi(xk−1),0})21(j∈E∑hj2(xk))21+(i∈I∑min2{gi(xk),0})21≤β成立,则令 μ k + 1 = μ k \mu_{k+1}=\mu_k μk+1=μk,否则 μ k + 1 = σ μ k \mu_{k+1}=\sigma\mu_k μk+1=σμk
- 计算 λ ( k + 1 ) \lambda^{(k+1)} λ(k+1),令 k = k + 1 k=k+1 k=k+1,转步2
3 实战测试
对于上节深入浅出最优化(7) 罚函数法中提出的约束最优化问题, x 1 , x 2 , x 3 x_1,x_2,x_3 x1,x2,x3的初值均在 [ 0 , 4 ] [0,4] [0,4]的范围内随机生成,总共生成100组起点。统计迭代成功(在1000步内得到最优解且单次步长搜索迭代次数不超过1000次)的样本的平均迭代步数、平均迭代时间和得到的最优解及开销函数最小值。
迭代步数 | 迭代时间 | 最优解 | 函数最小值 |
---|---|---|---|
17 | 24.3729s | x 1 = 1.1051 x 2 = 1.1969 x 3 = 1.5352 x_1=1.1051~x_2=1.1969~x_3=1.5352 x1=1.1051 x2=1.1969 x3=1.5352 | 0.03257 0.03257 0.03257 |
代码实现
本博客所有代码在https://github.com/HarmoniaLeo/optimization-in-a-nutshell开源,如果帮助到你,请点个star,谢谢这对我真的很重要!
你可以在上面的GitHub链接或本文集的第一篇文章深入浅出最优化(1) 最优化问题概念与基本知识中找到Function.py和lagb.py
使用共轭梯度PRP+法的拉格朗日乘子法
import numpy as np
from Function import Function #定义法求导工具
from lagb import * #线性代数工具库
from scipy import linalg
n=3 #x的长度
mu=2
lj=np.ones(1) #λj初值,长度等于等式限制条件个数
li=np.ones(2) #λi初值,长度等于不等式限制条件个数
def func(x): #目标函数,x是一个包含所有参数的列表
return (x[0]-1)**2+(x[0]-x[1])**2+(x[1]-x[2])**2
def hj(x): #构造数组h,第j位是第j+1个等式限制条件计算的值,x是一个包含所有参数的列表
return np.array([x[0]*(1+x[1]**2)+x[2]**4-4-3*np.sqrt(2)])
def gi(x): #构造数组g,第i位是第i+1个不等式限制条件计算的值,x是一个包含所有参数的列表
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)