凸优化之有等式约束的优化问题的求解方法

有等式约束的优化问题的求解方法

对数障碍 log barrier \text{log barrier} log barrier

首先, 介绍一下 log barrier \text{log barrier} log barrier 方法:

min ⁡ f 0 ( x ) s.t. f i ( x ) ≤ 0 i = 1 ⋯ m A x = b ⇔ min ⁡ f 0 ( x ) + ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) s.t. A x = b \begin{array}{lll} \min & f_{0}(x) &\\ \text{s.t.} & f_{i}(x) \leq 0 & i=1 \cdots m \\ & Ax = b & \\ \Leftrightarrow & & \\ \min & f_{0}(x)+\sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) \\ \text{s.t.} & A x=b \end{array} mins.t.⇔mins.t.f0(x)fi(x)≤0Ax=bf0(x)+∑i=1m−(1/t)log(−fi(x))Ax=bi=1⋯m

这是一个惩罚的思想. 具体来说就是 f i ( x ) → 0 f_{i}(x) \rightarrow 0 fi(x)→0 时, p ∗ → + ∞ p^{*} \rightarrow+\infty p∗→+∞, 会让目标函数根本无法得到最优解, 就好像有一个“棚栏”一样, 阻止 f i ( x ) f_{i}(x) fi(x) 趋近于 0, 让它只能小于0. 基于 log barrier \text{log barrier} log barrier 方法, 我们之后所有的约束都只讨论等式约束, 不等式约束用 log barrier \text{log barrier} log barrier 去掉. 即所有的凸问题都形如:
min ⁡ f 0 ( x ) s.t. A x = b \begin{array}{lll} \min & f_{0}(x) & \\ \text {s.t. } & Ax=b & \end{array} mins.t. f0(x)Ax=b

由于 − ( 1 / t ) log ⁡ ( − u ) -(1 / t) \log (-u) −(1/t)log(−u) 是 u u u 的单增凸函数, 上式中的目标函数是可微凸函数. 假定恰当的闭性条件成立，则可用 Newton 方法求解该问题.
我们将函数
ϕ ( x ) = − ∑ i = 1 m log ⁡ ( − f i ( x ) ) \phi(x)=-\sum_{i=1}^{m} \log \left(-f_{i}(x)\right) ϕ(x)=−i=1∑mlog(−fi(x))
dom ⁡ ϕ = { x ∈ R n ∣ f i ( x ) < 0 , i = 1 , ⋯ , m } \operatorname{dom} \phi=\left\{x \in \mathbf{R}^{n} \mid f_{i}(x)<0, i=1, \cdots, m\right\} domϕ={x∈Rn∣fi(x)<0,i=1,⋯,m} 称为原问题的对数障碍函数或对数障碍. 其定义域是满足原问题的严格不等式约束的点集. 不管正参数 t t t 取什么值，对于任意 i i i, 当 f i ( x ) → 0 f_{i}(x) \rightarrow 0 fi(x)→0 时，对数障碍函数将趋于无穷大. 并且, 我们可以知道当 t t t 越来越大时, ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) \sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) ∑i=1m−(1/t)log(−fi(x)) 会非常近似于指示函数( indicator function \text{indicator function} indicator function), 即当 f i ( x ) < 0 f_{i}(x)<0 fi(x)<0 时, ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) = 0 \sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) = 0 ∑i=1m−(1/t)log(−fi(x))=0, 当 f i ( x ) = 0 f_{i}(x) = 0 fi(x)=0 时, ∑ i = 1 m − ( 1 / t ) log ⁡ ( − f i ( x ) ) = + ∞ \sum_{i=1}^{m}-(1 / t) \log \left(-f_{i}(x)\right) = +\infty ∑i=1m−(1/t)log(−fi(x))=+∞

有等式约束的优化问题定义

min ⁡ f ( x ) s.t. A x = b (eq0) \begin{array}{ll} \min &f(x) \\ \text{s.t.} &A x = b \end{array} \tag{eq0} mins.t.f(x)Ax=b(eq0)

计算它的 KKT \text{KKT} KKT 条件:
x ∗ 是最优解 , v ∗ 是拉格朗日乘子 { A x ∗ = b primal feasibility ∇ f ( x ∗ ) + A T v ∗ = 0 stationarity x^{*} 是最优解, v^{*} 是拉格朗日乘子 \\ \left\{ \begin{array}{l} A x^{*}=b \quad \text{primal feasibility} \\ \nabla f\left(x^{*}\right)+A^{T} v^{*}=0 \quad \text{stationarity} \end{array} \right. x∗是最优解,v∗是拉格朗日乘子{Ax∗=bprimal feasibility∇f(x∗)+ATv∗=0stationarity

解约束优化问题就不得不用到 KKT \text{KKT} KKT 条件了，解 KKT \text{KKT} KKT 条件(如果可解的话), 那么相当于一次性把原问题和对偶问题的最优解都解出来了. 目前有一个从 KKT \text{KKT} KKT 条件出发得出的算法, 就叫做拉格朗日乘子法]

当目标函数 f ( x ) f(x) f(x) 是二次函数时, KKT \text{KKT} KKT 条件是线性方程组

当目标函数 f ( x ) f(x) f(x) 是二次函数时, KKT \text{KKT} KKT 条件是线性方程组.例如
min 1 2 x T P x + q T x + r s.t. A x = b \begin{array}{ll} \text{min} &\frac{1}{2} x^{T} P x+q^{T} x+r \\ \text{s.t.} &A x = b \end{array} mins.t.21xTPx+qTx+rAx=b

其中 P ∈ S + n , A P \in S_{+}^{n}, A P∈S+n,A 一般是不可逆的(这样 x x x 才会是有一个解集的), 否则可逆矩阵的方程组有唯一解就谈不上要优化了.
计算这个问题的 KKT \text{KKT} KKT 条件为:
{ A x ∗ = b P x ∗ + q + A T v ∗ = 0 ⇔ [ P A T A 0 ] ⋅ [ x ∗ v ∗ ] = [ − q b ] \left\{\begin{array}{ll}A x^{*} =b \\ P x^{*} +q+A^{T} v^{*}=0\end{array}\right. \Leftrightarrow \left[\begin{array}{cc}P & A^{T} \\ A & 0\end{array}\right] \cdot\left[\begin{array}{l}x^{*} \\ v^{*}\end{array}\right]=\left[\begin{array}{c}-q \\ b\end{array}\right] {Ax∗=bPx∗+q+ATv∗=0⇔[PAAT0]⋅[x∗v∗]=[−qb]
我们会发现，这时候 KKT \text{KKT} KKT 条件也是一个线性的方程组, 主要原因就是 ∇ f ( x ∗ ) \nabla f\left(x^{*}\right) ∇f(x∗) 是线性的, 这个时候就可以用Gauss-Seidel等解线性方程组的方法求解最优值了. 这样的情况不需要什么下降算法, 本质就是求解线性方程组. 当然, 你也可以转化为求无约束优化的对偶问题.

当目标函数 f ( x ) f(x) f(x) 是高次函数时, KKT \text{KKT} KKT 条件不是线性的方程组

当 ∇ f ( x ∗ ) \nabla f\left(x^{*}\right) ∇f(x∗) 不是一个线性函数时,则有
{ A x ∗ = b ∇ f ( x ∗ ) + A T v ∗ = 0 \left\{\begin{array}{l}A x^{*}=b \\ \nabla f\left(x^{*}\right)+A^{T} v^{*}=0\end{array}\right. {Ax∗=b∇f(x∗)+ATv∗=0

如果 ∇ f ( x ∗ ) \nabla f\left(x^{*}\right) ∇f(x∗) 不是一个关于 x ∗ x^{*} x∗ 的线性函数时, 自然会想到要找一个办法把它变成线性函数. 例如使用近似的思想,如果函数二次可微的话, 可以对函数进行二阶展开.
假设现在有一个可行点 x k x^{k} xk, x k x^{k} xk 可以使上式 ( e q 0 ) (eq0) (eq0) 成立,注意, x k x^{k} xk 只是可行解,而不一定是最优解.
现在我们构建一个问题 ( e q 1 ) (eq1) (eq1) ( e q 1 ) (eq1) (eq1): 在 x k x^{k} xk 的邻域附近寻找一个新的点 x k + d x^{k}+d xk+d,使得 f ( x k + d ) f(x^{k}+d) f(xk+d) 尽可能地小. 问题 ( e q 1 ) (eq1) (eq1) ( e q 1 ) (eq1) (eq1) 一定有解, 最坏情况下, d = 0 d=0 d=0, ( e q 1 ) (eq1) (eq1) 的解与 ( e q 0 ) (eq0) (eq0) 的解相同. 令 x k + 1 = x k + d x^{k+1} = x^{k} + d xk+1=xk+d, 然后我们再将 x k + 1 x^{k+1} xk+1 当做 x k x^{k} xk, 再在 x k + 1 x^{k+1} xk+1 邻域内寻找问题 ( e q 1 ) (eq1) (eq1) ( e q 1 ) (eq1) (eq1) 的解, 如此往复, 则我们就可以得到原问题 ( e q 0 ) (eq0) (eq0) 的解(或者近似解).

min ⁡ d f ( x k + d ) s.t. A ( x k + d ) = b ⇔ min ⁡ d f ( x k + d ) s.t. A d = 0 因为 A x k = b , A ( x k + d ) = b 所以 A d = 0 (eq1) \begin{array}{ll} \underset{d}{\min} & f\left(x^{k}+d\right) \\ \text{s.t.} & A\left(x^{k}+d\right)=b \end{array} \Leftrightarrow \begin{array}{ll} \underset{d}{\min} & f\left(x^{k}+d\right) \\ \text {s.t.} &A d=0 \\ \end{array} \\ 因为A x^{k}=b, A\left(x^{k}+d\right)=b 所以A d=0 \\ \tag{eq1} dmins.t.f(xk+d)A(xk+d)=b⇔dmins.t.f(xk+d)Ad=0因为Axk=b,A(xk+d)=b所以Ad=0(eq1)
则我们现在把原问题(求解全局最优解)等价为新的优化问题(在 x k x^{k} xk的邻域附近寻找一个新的点 x k + d x^{k}+d xk+d,使得 f ( x k + d ) f(x^{k}+d) f(xk+d) 尽可能地小).
这是一个关于 d d d 的新的优化问题. 为了让原问题 KKT \text{KKT} KKT 线性, 先对问题 ( e q 1 ) (eq1) (eq1)目标函数做二阶泰勒展开
f ( x k + d ) = f ( x k ) + ∇ f ( x k ) T d + 1 2 d T ( ∇ 2 f ( x k ) ) d f\left(x^{k}+d\right)=f\left(x^{k}\right)+\nabla f\left(x^{k}\right)^{T} d+\frac{1}{2} d^{T}\left(\nabla^{2} f\left(x^{k}\right)\right) d f(xk+d)=f(xk)+∇f(xk)Td+21dT(∇2f(xk))d
则问题 ( e q 1 ) (eq1) (eq1)又被近似等价为问题 ( e q 2 ) (eq2) (eq2)
min ⁡ d f ( x k ) + ∇ f ( x k ) T d + 1 2 d T ∇ 2 f ( x k ) T d s.t. A d = 0 (eq2) \begin{array}{ll} \underset{d}{\min} & f\left(x^{k}\right)+\nabla f\left(x^{k}\right)^{T} d+\frac{1}{2} d^{T} \nabla^{2} f\left(x^{k}\right)^{T} d \\ \text {s.t.} & A d=0 \end{array} \\ \tag{eq2} dmins.t.f(xk)+∇f(xk)Td+21dT∇2f(xk)TdAd=0(eq2)
求解这个问题 ( e q 2 ) (eq2) (eq2)的 KKT \text{KKT} KKT 条件为:
[ ∇ 2 f ( x k ) A T A 0 ] ⋅ [ d ∗ v ∗ ] = [ − ∇ f ( x k ) 0 ] \left[\begin{array}{cc}\nabla^{2} f\left(x^{k}\right) & A^{T} \\ A & 0\end{array}\right] \cdot\left[\begin{array}{l}d^{*} \\ v^{*}\end{array}\right]=\left[\begin{array}{c}-\nabla f\left(x^{k}\right) \\ 0\end{array}\right] [∇2f(xk)AAT0]⋅[d∗v∗]=[−∇f(xk)0]
便又回到解线性方程组了. 但是这里要注意的是, 我们这样近似解出来的是 d k d^{k} dk, 也就是用上一次的迭代点沿着 d d d 方向展开了得到的 d k d^{k} dk, 这一个过程我们进行了二阶展开这样的近似操作, 并不是求解出了原问题的最优解但是我们可以发现, 解出来的 d k d^{k} dk 满足 A d k = 0 A d^{k}=0 Adk=0, 也就是说 x k + d k x^{k}+d^{k} xk+dk 起码是可行点,满足 A ( x k + d k ) = b A\left(x^{k}+d^{k}\right)=b A(xk+dk)=b .
到这里我们就会思考，可不可以用这个方向去下降呢? 答案是最好先别. 因为 x k + d k x^{k}+d^{k} xk+dk 只是可行点, 不一定相对于 x k x^{k} xk 是下降点. 所以我们可以给 d k d^{k} dk 乘上一个步长, 对步长做一次线搜索，再来当作下降方向来用. 由此我们得到带等式约束的牛顿法:
{ [ ∇ 2 f ( x k ) A T A 0 ] ⋅ [ d ∗ v ∗ ] = [ − ∇ f ( x k ) 0 ] d k = d ∗ α k = arg ⁡ min ⁡ α ≥ 0 f ( x k + α d k ) x k + 1 = x k + α k d k , x k + 1 始终满足 A x k + 1 = b \left\{ \begin{array}{l} \left[\begin{array}{cc}\nabla^{2} f\left(x^{k}\right) & A^{T} \\ A & 0\end{array}\right] \cdot\left[\begin{array}{l}d^{*} \\ v^{*}\end{array}\right]=\left[\begin{array}{c}-\nabla f\left(x^{k}\right) \\ 0\end{array}\right] \\ d^{k}=d^{*} \\ \alpha^{k}=\arg \underset{\alpha \geq 0}{\min} f \left(x^{k}+\alpha d^{k}\right) \\ x^{k+1}=x^{k}+\alpha^{k} d^{k}, \; x^{k+1}始终满足A x^{k+1} = b \\ \end{array} \right. ⎩⎪⎪⎪⎪⎨⎪⎪⎪⎪⎧[∇2f(xk)AAT0]⋅[d∗v∗]=[−∇f(xk)0]dk=d∗αk=argα≥0minf(xk+αdk)xk+1=xk+αkdk,xk+1始终满足Axk+1=b
所以，对于带约束优化问题, KKT \text{KKT} KKT 条件还是蛮重要的

拉格朗日乘子法Method of Lagrangian Multiplier, 拉格朗日法Lagrangian Method

拉格朗日乘子法/拉格朗日法是通过某种方法将原变量和对偶变量看成一个大的变量集合体, 然后确定一个优化方向, 同时优化原变量和对偶变量

{ x k + 1 = x k − α k ( ∇ f ( x k ) + A T v k ) v k + 1 = v k + α k ( A x k − b ) ⇔ [ x k + 1 v k + 1 ] = [ x k v k ] + α k [ − ( ∇ f ( x k ) + A T v k ) A x k − b ] \left\{\begin{array}{l}x^{k + 1}=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \\ v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k}-b\right)\end{array}\right. \\ \Leftrightarrow \\ \left [\begin{array}{c} x^{k + 1} \\ v^{k+1} \end{array}\right ] = \left [\begin{array}{c} x^{k} \\ v^{k} \end{array}\right]+ \alpha^{k} \left [\begin{array}{c} -\left( \nabla f (x^{k})+A^{T} v^{k} \right) \\ A x^{k}-b \end{array}\right] {xk+1=xk−αk(∇f(xk)+ATvk)vk+1=vk+αk(Axk−b)⇔[xk+1vk+1]=[xkvk]+αk[−(∇f(xk)+ATvk)Axk−b]

第二个式子乘子的迭代部分, x k , v k x^{k},v^{k} xk,vk 肯定是可以换成上一步已经计算好的 x k + 1 , v k + 1 x^{k+1},v^{k+1} xk+1,vk+1. 很明显,拉格朗日寻找下一个更新点的方法如上式,找到一个方向 [ − ( ∇ f ( x k ) + A T v k ) A x k − b ] \left [\begin{array}{c} -\left( \nabla f (x^{k})+A^{T} v^{k} \right) \\ A x^{k}-b \end{array}\right] [−(∇f(xk)+ATvk)Axk−b], 再寻找合适的步长 α k \alpha^{k} αk, 然后进行更新.

从鞍点角度理解拉格朗日法

因为
L ( x , v ) = f ( x ) + v T ( A x − b ) L(x, v)=f(x)+v^{T}(A x-b) L(x,v)=f(x)+vT(Ax−b)
所以鞍点
{ ( x ∗ , v ∗ ) = arg ⁡ max ⁡ v min ⁡ x L ( x , v ) ( x ∗ , v ∗ ) = arg ⁡ min ⁡ x max ⁡ v L ( x , v ) \begin{aligned} \left \{ \begin{array}{ll} (x^{*}, v^{*}) &=\arg \underset{v}{\max} \; \underset{x}{\min} \; L(x, v) \\ (x^{*}, v^{*}) &=\arg \underset{x}{\min} \; \underset{v}{\max} \; L(x, v) \end{array} \right. \end{aligned} {(x∗,v∗)(x∗,v∗)=argvmaxxminL(x,v)=argxminvmaxL(x,v)
我们知道，如果强对偶性成立，那么鞍点就是最优点. 此时
x ∗ = arg ⁡ min ⁡ x L ( x , v ∗ ) ( 已知 v ∗ 求 x ∗ ) v ∗ = arg ⁡ max ⁡ v L ( x ∗ , v ) ( 已知 x ∗ 求 v ∗ ) 上述求 x ∗ , v ∗ 的过程是无约束优化问题求极值的方法 , 可以使用梯度下降或者梯度上升的算法 . x^{*}=\arg \underset{x}{\min} L\left(x, v^{*}\right)\left(\right. 已知 v^{*} 求 \left.x^{*}\right) \\ v^{*}=\arg \underset{v}{\max} L\left(x^{*}, v\right)\left(\right.已知 x^{*}求 \left.v^{*}\right) \\ 上述求 x^{*}, v^{*}的过程是无约束优化问题求极值的方法,可以使用梯度下降或者梯度上升的算法. x∗=argxminL(x,v∗)(已知v∗求x∗)v∗=argvmaxL(x∗,v)(已知x∗求v∗)上述求x∗,v∗的过程是无约束优化问题求极值的方法,可以使用梯度下降或者梯度上升的算法.

如果已知某一个分量的最优点, 求另一个分量的最优点本质上就是求拉格朗日函数对于单个变量的最优解. 接下来分类讨论下:

当已知 v ∗ v^{*} v∗ 求 x ∗ x^{*} x∗ 时:
用梯度下降法迭代求解 L ( x , v ∗ ) L\left(x, v^{*}\right) L(x,v∗) . 因为 − ∇ L ( x , v ∗ ) = − ∇ f ( x ) − A T v ∗ -\nabla L\left(x, v^{*}\right)=-\nabla f(x)-A^{T} v^{*} −∇L(x,v∗)=−∇f(x)−ATv∗, 所以我们的迭代算法为 x k + 1 = x k + α k ( − ∇ f ( x k ) − A T v ∗ ) x^{k+1}=x^{k}+\alpha^{k}\left(-\nabla f\left(x^{k}\right)-A^{T} v^{*}\right) xk+1=xk+αk(−∇f(xk)−ATv∗) 但是我们不可能真的已知 v ∗ v^{*} v∗ 是什么, 所以我们用 v k v^{k} vk 近似. 这便是拉格朗日法的关于 x x x 的迭代部分： x k + 1 = x k − α k ( ∇ f ( x k ) + A T v k ) x^{k+1}=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) xk+1=xk−αk(∇f(xk)+ATvk)
当已知 x ∗ x^{*} x∗ 求 v ∗ v^{*} v∗ 时:
因为我们求极大, 所以上面的迭代方向就是梯度方向,同样地 x ∗ x^{*} x∗ 用 x k x^{k} xk 替代. v k + 1 = v k + α k ( A x k − b ) v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k}-b\right) vk+1=vk+αk(Axk−b)
拉格朗日法实际上是在求 L ( x , v ) L(x, v) L(x,v) 的鞍点, 来得出原问题的最优值.

从罚函数的角度理解拉格朗日法(最好理解的方式)

我们把带等式约束的原问题的 KKT \text{KKT} KKT 条件:
{ A x ∗ = b ∇ f ( x ∗ ) + A T v ∗ = 0 \left\{\begin{array}{l}A x^{*}=b \\ \nabla f\left(x^{*}\right)+A^{T} v^{*}=0\end{array}\right. \\ {Ax∗=b∇f(x∗)+ATv∗=0
改写成一个罚函数的样子, 去当作一个新的优化问题来解:
min ⁡ P ( x , v ) = 1 2 ∥ A x − b ∥ 2 2 + 1 2 ∥ ∇ f ( x ) + A T v ∥ 2 2 ( 一般为非凸问题 ) (eq3) \min P(x, v)=\frac{1}{2}\|A x-b\|_{2}^{2}+\frac{1}{2}\left\|\nabla f(x)+A^{T} v\right\|_{2}^{2} \quad (一般为非凸问题) \\ \tag{eq3} minP(x,v)=21∥Ax−b∥22+21∥∥∇f(x)+ATv∥∥22(一般为非凸问题)(eq3)
特点: 这个新的无约束优化问题的解就是满足 KKT \text{KKT} KKT 条件的点. 某种意义上说，相当于把原问题转化成了无约束优化问题, 我们可以去求这个新问题的目标函数的负梯度方向，然后梯度下降求最优解:

− ∇ P ( x k , v k ) = [ − ∇ x P ( x k , v k ) , − ∇ v P ( x k , v k ) ] = − [ A T ( A x k − b ) + ∇ 2 f ( x k ) ( ∇ f ( x k ) + A T v k ) A ( ∇ f ( x k ) + A T v k ) ] \begin{array}{ll} &-\nabla P\left(x^{k}, v^{k}\right) \\ &=\left[\begin{array}{l}-\nabla_{x} P\left(x^{k}, v^{k}\right), -\nabla_{v} P\left(x^{k}, v^{k}\right)\end{array}\right] \\ &=-\left[ \begin{array}{c} A^{T}\left(A x^{k}-b\right)+\nabla^{2} f\left(x^{k}\right)\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \\ A\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \end{array} \right] \end{array}\\ −∇P(xk,vk)=[−∇xP(xk,vk),−∇vP(xk,vk)]=−[AT(Axk−b)+∇2f(xk)(∇f(xk)+ATvk)A(∇f(xk)+ATvk)]

拉格朗日方法中迭代方向的正确性证明

我们选择 − ∇ P ( x k , v k ) -\nabla P\left(x^{k}, v^{k}\right) −∇P(xk,vk) 为梯度方向,如果拉格朗日法所确定的迭代方向 d k d^{k} dk与
− ∇ P ( x k , v k ) -\nabla P\left(x^{k}, v^{k}\right) −∇P(xk,vk) 夹角小于 π 2 \frac{\pi}{2} 2π,那么我们认为迭代方向 d k d^{k} dk 能够使函数值下降.
d k = [ − ( ∇ f ( x k ) + A T v k ) A x k − b ] d^{k}=\left [ \begin{array}{c} -\left( \nabla f (x^{k})+A^{T} v^{k} \right) \\ A x^{k}-b \end{array} \right] dk=[−(∇f(xk)+ATvk)Axk−b]

( d k ) T ( − ∇ P ( x k , v k ) ) = ( ∇ f ( x k ) + A T v k ) T A T ( A x k − b ) + ( ∇ f ( x k ) + A T v k ) T ∇ 2 f ( x k ) ( ∇ f ( x k ) + A T v k ) − ( A T x k − b ) T A ( ∇ f ( x k ) + A T v k ) = ( ∇ f ( x k ) + A T v k ) T ∇ 2 f ( x k ) ( ∇ f ( x k ) + A T v k ) \begin{aligned} &(d^{k})^{T} \left(-\nabla P(x^{k}, v^{k})\right)\\ &=\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right)^{T} A^{T}\left(A x^{k}-b\right) \\ &\qquad + \left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right)^{T} \nabla^{2} f\left(x^{k}\right) \left(\nabla f\left(x^{k}\right)+A^{T} v^{k} \right) \\ &\qquad - (A^{T}x^{k}-b)^{T}A(\nabla f\left(x^{k}\right)+A^{T} v^{k})\\ &=\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right)^{T} \nabla^{2} f\left(x^{k}\right)\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) \end{aligned} (dk)T(−∇P(xk,vk))=(∇f(xk)+ATvk)TAT(Axk−b)+(∇f(xk)+ATvk)T∇2f(xk)(∇f(xk)+ATvk)−(ATxk−b)TA(∇f(xk)+ATvk)=(∇f(xk)+ATvk)T∇2f(xk)(∇f(xk)+ATvk)

上边这个式子(二次型), 当 ∇ 2 f ( x k ) ≻ 0 \nabla^{2} f\left(x^{k}\right) \succ 0 ∇2f(xk)≻0 且 ( ∇ f ( x k ) + A T v k ) \left(\nabla f\left(x^{k}\right)+A^{T} v^{k}\right) (∇f(xk)+ATvk) 不全为 0 时, 会大于等于0 . 即: { ∇ f ( x k ) + A T v k ≠ 0 ∇ 2 f ( x k ) ≻ 0 \left\{\begin{array}{l}\nabla f\left(x^{k}\right)+A^{T} v^{k} \neq 0 \\ \nabla^{2} f\left(x^{k}\right) \succ 0\end{array}\right. {∇f(xk)+ATvk=0∇2f(xk)≻0 时, d k d^{k} dk 是 P ( u k , v k ) P\left(u^{k}, v^{k}\right) P(uk,vk) 的一个下降方向(和负梯度方向同向). 而上面这个条件是成立的. 因为原问题凸, 所以Hessian矩阵半正定. 而当随着迭代 x k → x ∗ x^{k} \rightarrow x^{*} xk→x∗ 的时候, ∇ 2 f ( x ∗ ) = 0 \nabla^{2} f\left(x^{*}\right)=0 ∇2f(x∗)=0, 那么大部分情况下 ∇ 2 f ( x k ) \nabla^{2} f\left(x^{k}\right) ∇2f(xk) 都是正定的. 另外, 只要 v k , x k v^{k}, x^{k} vk,xk 还没有到达最优解, ∇ f ( x k ) + A T v k ≠ 0 \nabla f\left(x^{k}\right)+A^{T} v^{k} \neq 0 ∇f(xk)+ATvk=0 就一直满足, 当 v ∗ → v ∗ , x k → x ∗ v^{*} \rightarrow v^{*}, x^{k} \rightarrow x^{*} v∗→v∗,xk→x∗ 时, ∇ f ( x k ) + A T v k = 0 \nabla f\left(x^{k}\right)+A^{T} v^{k}=0 ∇f(xk)+ATvk=0, 此时 x k + 1 = x k , v k + 1 = v k x^{k+1}=x^{k}, v^{k+1}=v^{k} xk+1=xk,vk+1=vk,算法也就停止更新了. 对拉格朗日法进行理论分析很有意义, 但是这个算法实际运行起来速度很慢, 所以为了提高一下效率, 诞生了增广拉格朗日法.

增广拉格朗日法

增广拉格朗日法是交替优化原变量和对偶变量,并以对偶变量优化为主
增广拉格朗日法相比与拉格朗日法, 会构建一个新的函数, 称为增广拉格朗日函数. 增广拉格朗日函数是由拉格朗日函数加上一个罚项得到的, 记为:
L c ( x , v ) = f ( x ) + v T ( A x − b ) + c 2 ∥ A x − b ∥ 2 2 L_{c}(x, v)=f(x)+v^{T}(A x-b)+\frac{c}{2}\|A x-b\|_{2}^{2} Lc(x,v)=f(x)+vT(Ax−b)+2c∥Ax−b∥22
它实际上是下面这个优化问题 ( e q 4 ) (eq4) (eq4) 的拉格朗日函数:
min ⁡ f ( x ) + c 2 ∥ A x − b ∥ 2 2 s.t. A x = b (eq4) \begin{array}{ll} \min & f(x)+\frac{c}{2}\|A x-b\|_{2}^{2} \\ \text{s.t.} & A x=b \end{array} \tag{eq4} mins.t.f(x)+2c∥Ax−b∥22Ax=b(eq4)
也就是把原问题的等式约束变为一个惩罚项再放到目标函数里面去, 得到的优化问题的拉格朗日函数, 被称为原问题拉格朗日函数的增广拉格朗日函数.
接下来很容易可以验证, 上面 ( e q 4 ) (eq4) (eq4) 这个问题和开头的标准问题 ( e q 0 ) (eq0) (eq0) 的最优解相同, 两个问题是等价的.
对于问题 ( e q 0 ) (eq0) (eq0), 其拉格朗日函数为： L ( x , v ) = f ( x ) + v T ( A x − b ) L(x, v)=f(x)+v^{T}(A x-b) L(x,v)=f(x)+vT(Ax−b), 假设拉格朗日函数的最优解为 ( x ∗ , v ∗ ) \left(x^{*}, v^{*}\right) (x∗,v∗), 则 L ( x ∗ , v ∗ ) = 0 L\left(x^{*}, v^{*}\right)=0 L(x∗,v∗)=0, 即: ∇ x { f ( x ∗ ) + v ∗ T ( A x ∗ − b ) } = 0 \nabla_{x}\left\{f\left(x^{*}\right)+v^{* T}\left(A x^{*}-b\right)\right\}=0 ∇x{f(x∗)+v∗T(Ax∗−b)}=0
对于问题 ( e q 4 ) (eq4) (eq4), 其拉格朗日函数为: L c ( x , v ) = f ( x ) + v T ( A x − b ) + c 2 ∥ A x − b ∥ 2 2 L_{c}(x, v)=f(x)+v^{T}(A x-b)+\frac{c}{2}\|A x-b\|_{2}^{2} Lc(x,v)=f(x)+vT(Ax−b)+2c∥Ax−b∥22
假设拉格朗日函数的最优解为 ( x 1 ∗ , v 1 ∗ ) \left(x_{1}^{*}, v_{1}^{*}\right) (x1∗,v1∗), 则 L c ( x 1 ∗ , v 1 ∗ ) = 0 L_{c}\left(x_{1}^{*}, v_{1}^{*}\right)=0 Lc(x1∗,v1∗)=0, 即:
∇ x { f ( x 1 ∗ ) + v 1 ∗ T ( A x 1 ∗ − b ) } + c A T ( A x 1 ∗ − b ) = 0 \nabla_{x}\left\{f\left(x_{1}^{*}\right)+v_{1}^{* T}\left(A x_{1}^{*}-b\right)\right\}+c A^{T}\left(A x_{1}^{*}-b\right)=0 ∇x{f(x1∗)+v1∗T(Ax1∗−b)}+cAT(Ax1∗−b)=0
又因为 A x 1 ∗ − b = 0 A x_{1}^{*}-b=0 Ax1∗−b=0, 所以实际上上式为 ∇ x { f ( x 1 ∗ ) + v 1 ∗ T ( A x 1 ∗ − b ) } = 0 \nabla_{x}\left\{f\left(x_{1}^{*}\right)+v_{1}^{* T}\left(A x_{1}^{*}-b\right)\right\}=0 ∇x{f(x1∗)+v1∗T(Ax1∗−b)}=0
这和 ( e q 0 ) (eq0) (eq0) 的结论是一样的. 所以说 ( x ∗ , v ∗ ) = ( x 1 ∗ , v 1 ∗ ) \left(x^{*}, v^{*}\right)=\left(x_{1}^{*}, v_{1}^{*}\right) (x∗,v∗)=(x1∗,v1∗), 即 ( e q 0 ) ⇔ ( e q 4 ) (eq0) \Leftrightarrow(eq4) (eq0)⇔(eq4) .
既然两个问题是等价的，也就是它们各自的拉格朗日函数的最优值是等价的. 那么, 自然可以想到用 ∇ L c ( x , v ) \nabla L_{c}(x, v) ∇Lc(x,v) 去替换掉 ∇ L ( x , v ) \nabla L(x, v) ∇L(x,v) .
使用拉格朗日法来解决问题 ( e q 4 ) (eq4) (eq4)(问题 ( e q 0 ) (eq0) (eq0)的增广拉格朗日函数就是问题 ( e q 4 ) (eq4) (eq4)的拉格朗日函数)
{ x k + 1 = x k − α k ∇ x L c ( x k , v k ) = x k − α k ( ∇ f ( x k ) + A T v k + c A T ( A x k − b ) ) v k + 1 = v k + α k ∇ v L c ( x k , v k ) = v k + α k ( A x k − b ) \left\{\begin{array}{l}x^{k+1}=x^{k}-\alpha^{k} \nabla_{x} L_{c}\left(x^{k}, v^{k}\right)=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}+c A^{T}\left(A x^{k}-b\right)\right) \\ v^{k+1}=v^{k}+\alpha^{k} \nabla_{v} L_{c}\left(x^{k}, v^{k}\right)=v^{k}+\alpha^{k}\left(A x^{k}-b\right)\end{array}\right. {xk+1=xk−αk∇xLc(xk,vk)=xk−αk(∇f(xk)+ATvk+cAT(Axk−b))vk+1=vk+αk∇vLc(xk,vk)=vk+αk(Axk−b)
但是上面这组迭代式子, 有个东西没有利用起来, 就是第一个式子会更新好的 x k + 1 x^{k+1} xk+1 . 在第二式子计算 v k + 1 v^{k+1} vk+1 时用的还是 x k x^{k} xk 的信息, 这是把它替换成新的 x k + 1 x^{k+1} xk+1, 那么整个算法变为:
{ x k + 1 = x k − α k ( ∇ f ( x k ) + A T v k + c A T ( A x k − b ) ) v k + 1 = v k + α k ( A x k + 1 − b ) \left\{\begin{array}{l}x^{k+1}=x^{k}-\alpha^{k}\left(\nabla f\left(x^{k}\right)+A^{T} v^{k}+c A^{T}\left(A x^{k}-b\right)\right) \\ v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k+1}-b\right)\end{array}\right. {xk+1=xk−αk(∇f(xk)+ATvk+cAT(Axk−b))vk+1=vk+αk(Axk+1−b)
如果增广拉格朗日函数 L c ( x , v ) L_{c}(x, v) Lc(x,v) 性质比较理想或者比较容易被优化的话, 不一定非要用梯度下降法, 比如二阶可微也可以用牛顿法, 等等. 所以统一改成用 arg ⁡ min ⁡ \arg \min argmin 的形式来表达上面这个算法为:
{ x k + 1 = arg ⁡ min ⁡ x L c ( x , v k ) v k + 1 = v k + α k ( A x k + 1 − b ) \left\{\begin{array}{l}x^{k+1}=\arg \underset{x}{\min} L_{c}\left(x, v^{k}\right) \\ v^{k+1}=v^{k}+\alpha^{k}\left(A x^{k+1}-b\right)\end{array}\right. {xk+1=argxminLc(x,vk)vk+1=vk+αk(Axk+1−b)
我们还可以把第二个式子的步长替换为罚参数, 得到最终的增广拉格朗日法：
{ x k + 1 = arg ⁡ min ⁡ x L c ( x , v k ) v k + 1 = v k + c ( A x k + 1 − b ) \left\{\begin{array}{l}x^{k+1}=\arg \underset{x}{\min} L_{c}\left(x, v^{k}\right) \\ v^{k+1}=v^{k}+c\left(A x^{k+1}-b\right)\end{array}\right. {xk+1=argxminLc(x,vk)vk+1=vk+c(Axk+1−b)

因此最终的增广拉格朗日法 Augmented Lagrangian Method \text{Augmented Lagrangian Method} Augmented Lagrangian Method为
x k + 1 = arg ⁡ min ⁡ x L c ( x , v k ) v k + 1 = v k + c ( A x k + 1 − b ) \begin{array}{l} x^{k+1}=\arg \underset{x}{\min} L_{c}\left(x, v^{k}\right) \\ v^{k+1}=v^{k}+c\left(A x^{k+1}-b\right) \end{array} xk+1=argxminLc(x,vk)vk+1=vk+c(Axk+1−b)
步长 α k \alpha^{k} αk 都是要自己找的, 比如精确搜索或非精确搜索.

注：这里 c c c 要自己选定, 但是选定 c c c 肯定比去线搜索 α \alpha α 简单多了. 特别的，如果选择 c = 1 c=1 c=1, 基本上可以保证增广拉格朗日算法收敘. 而且, 从收敘性分析可以知道, c c c 也可以是递增序列，不用担心不收敘.

算法性质:

(1). 若 v = v ∗ v=v^{*} v=v∗,则 ∀ c > 0 , x ∗ = arg ⁡ min ⁡ x L c ( x , v ∗ ) \forall c>0, x^{*}=\arg \underset{x}{\min} L_{c}\left(x, v^{*}\right) ∀c>0,x∗=argxminLc(x,v∗), 只要 v → v ∗ , c v \rightarrow v^{*}, c v→v∗,c 的值只要大于0就可以得到最优解 x ∗ x^{*} x∗
(2). 若 c → + ∞ c \rightarrow+\infty c→+∞,则 ∀ v , x ∗ = arg ⁡ min ⁡ x L c ( x , v ) \forall v, \quad x^{*}=\arg \underset{x}{\min} L_{c}(x, v) ∀v,x∗=argxminLc(x,v), 步长 c c c 可以选择递增序列 { c k } \left\{c^{k}\right\} {ck}, 此时一定收敛.假设 c = + ∞ c=+\infty c=+∞

参考资料

[1] Convex Optimization – Boyd and \nuandenberghe
[2] 中科大-凸优化
[2] 知乎-落日翻车-凸优化总结