数值计算之 拉格朗日乘子法初探

前言

LM算法的置信域方法中，通过在优化函数后添加一个优化量约束来提升迭代过程中的稳定性。

拉格朗日乘子法(Lagrange Multiplier Method)是非线性优化中的常用方法，置信域约束实际上就是拉格朗日乘子法的应用。

本篇将学习记录拉格朗日乘子法的原理和使用。

注意：之所以是初探，是因为我还没有总结过凸优化的基础内容，现在不会深入研究拉格朗日乘子法需要满足的限定条件。

等式约束优化

对于优化函数$f:{\mathbf{R}}^n\to \mathbf{R}$，求$f$在某个可行域$\mathbf{x}\in {\mathbf{R}}^{\mathbf{n}}$下的极值点。如果可行域通过等式$g(\mathbf{x})=0$表达，则记为：
$$\begin{gathered} \underset{\mathbf{x}}{\min }f(\mathbf{x}) \\ s.t.g(\mathbf{x})=0 \end{gathered}$$
假设$f(\mathbf{x}),g(\mathbf{x})$都连续可导，则定义拉格朗日函数，获得极值条件：
$$\begin{gathered} L(\mathbf{x},\lambda )=f(\mathbf{x})+\lambda g(\mathbf{x}) \\ \underset{\mathbf{x},\lambda }{\min }L(\mathbf{x},\lambda )=\underset{\mathbf{x},\lambda }{\min }(f(\mathbf{x})+\lambda g(\mathbf{x})) \\ \to \{\begin{array}{l}{\nabla }_{\mathbf{x}}L(\mathbf{x},\lambda )={\nabla }_{\mathbf{x}}f(\mathbf{x})+\lambda {\nabla }_{\mathbf{x}}g(\mathbf{x})=0\\ {\nabla }_{\lambda }L(\mathbf{x},\lambda )=g(\mathbf{x})=0\end{array} \end{gathered}$$
于是，等式约束优化就变成了无约束优化问题。

不过，拉格朗日函数的构造有一些突兀，这里给出一个几何解释。

等式约束的几何解释

下图是求极大值的例子，蓝色虚线表示优化函数$f(x,y)$，$f(x,y)=d_1,d_2,\dots$表示该函数的一个个等高线，$d_1>d_2>d_3$，蓝色小箭头表示该处的梯度；红色线表示约束方程$g(x,y)=c$，红色小箭头表示该处的梯度。需要注意，该图像是从三维空间投影到xy平面上的结果。

等式约束条件下，$x,y$的值只能在红线上取。优化函数的极值必然会在蓝线与红线相切的切点获得，这是因为：蓝色梯度是优化函数增长最快的方向，如果极值点是切点附近的点，那么向切点运动，函数值仍然会增加，这就不符合极大值点的定义了。求极小值的情况也同理。

给出切点处是极小值点的数学形式，即：
$$\begin{gathered} {\nabla }_{x,y}f(x,y)=-\lambda {\nabla }_{x,y}g(x,y) \\ g(x,y)=0 \end{gathered}$$
第一个式子表示两个函数$f,g$处的梯度平行（说$f,g$相切并不严谨，因为图示只是在xy平面上的函数投影），第二个式子表示约束条件决定的可行域。

不等式约束优化

优化函数不变，约束条件变成了不等式$g(\mathbf{x})\le 0$，则可行域从曲线变成了一个曲面，如下图所示，灰色部分代表可行域。

此时需要考虑两种情况。

第一种情况如图左，如果极值点落在可行域边界内，即$g(\mathbf{x})<0$，则约束条件并没有影响优化函数本来的极值点位置${\mathbf{x}}^{\ast }$，也就是说拉格朗日函数$L(\mathbf{x},\lambda )$与$f(\mathbf{x})$的极值是等价的，$\lambda =0$，得到极值条件：
$$\begin{gathered} {\nabla }_{\mathbf{x}}L(\mathbf{x},\lambda )=0 \\ \lambda =0 \\ g(\mathbf{x})<0 \end{gathered}$$

第二种情况如图右，如果极值点落在可行域边界上，即$g(\mathbf{x})=0$，则与等式约束优化的极值条件相同：
$$\begin{gathered} {\nabla }_{\mathbf{x}}L(\mathbf{x},\lambda )=0 \\ \lambda \ge 0 \\ g(\mathbf{x})=0 \end{gathered}$$
其中，$\lambda \ge 0$是因为，在求极小值时，等式约束条件是${\nabla }_{\mathbf{x}}f(\mathbf{x})=-\lambda \nabla g_{\mathbf{x}}(\mathbf{x})$，其中${\nabla }_{\mathbf{x}}f(\mathbf{x})$指向可行域内（如果指向可行域外，那么在可行域内移动时，优化函数值还会继续下降，不符合极小值定义），${\nabla }_{\mathbf{x}}g(\mathbf{x})$指向可行域外（梯度方向指向$g(\mathbf{x})>0$）。这就限定了$\lambda \ge 0$。$\lambda =0$意味着可行域边界通过了${\nabla }_{\mathbf{x}}f(\mathbf{x})=0$的位置。

将两种情况综合考虑，获得不等式约束优化的KKT(Karush-Kuhn-Tucker)条件：
$$\begin{gathered} {\nabla }_{\mathbf{x}}L(\mathbf{x},\lambda )={\nabla }_{\mathbf{x}}f(\mathbf{x})+\lambda {\nabla }_{\mathbf{x}}g(\mathbf{x})=0 \\ g(\mathbf{x})\le 0 \\ \lambda \ge 0 \\ \lambda g(\mathbf{x})=0 \end{gathered}$$

多约束优化

上面的等式与不等式约束可以共同推广到多个约束的优化情况：
$$\begin{gathered} L(\mathbf{x},\lambda )=f(\mathbf{x})+\sum _{i=1}^m{\lambda }_i{g}_i(\mathbf{x})+\sum _{j=1}^n{\mu }_j{h}_j(\mathbf{x}) \\ s.t.g_i(\mathbf{x})\le 0,h_j(\mathbf{x})=0 \end{gathered}$$

其KKT条件为：
$$\begin{gathered} {\nabla }_{\mathbf{x}}L(\mathbf{x},\lambda )={\nabla }_{\mathbf{x}}f(\mathbf{x})+\sum _{i=1}^m{\lambda }_i\nabla g_i(\mathbf{x})+\sum _{j=1}^n{\mu }_j\nabla h_j(\mathbf{x})=0 \\ {g}_i(\mathbf{x})\le 0 \\ {h}_j(\mathbf{x})=0 \\ {\lambda }_i\ge 0 \\ {\lambda }_i{g}_i(\mathbf{x})=0 \end{gathered}$$

KKT条件

上面给出的KKT条件，实际上是约束优化问题有解时满足的必要条件。也就是说，非极值也可能满足KKT条件。因此，只有满足某些正则条件时，KKT条件才能正确使用，比如最常用的线性无关条件LICQ：

对于极值点落在边界上的不等式约束条件在极值点处的梯度与等式约束在极值点处的梯度线性无关时。

由于这一部分涉及很多凸优化内容，因此等我学习到凸优化时，再来深入探究KKT条件。

代码示例

这里给出一个基于python sympy的等式约束拉格朗日乘子法示例。由于sympy不能单独完成不等式约束，需要添加其它模块，以后会补充不等式约束代码。

import numpy as np
import math
import scipy
import sympy


def solve_lagarange(F, G=0, H=0):
    x, y, l, u = sympy.symbols('x, y, l, u')
    L = F + l * G + u * H
    px, py, pl, pu = sympy.derive_by_array(L, [x, y, l, u])
    res = sympy.solve([px, py, pl, pu], x, y, l, u)
    print(res)


def eq_limit_solve():
    a, b = sympy.symbols('x, y')
    optim_function = non_linear_func_2(a, b)
    eq_limit_function = limite_func_1(a, b)
    solve_lagarange(optim_function, 0, eq_limit_function)


def non_linear_func(x, y):
    fxy = 0.5 * (x ** 2 + y ** 2)
    return fxy


def non_linear_func_2(x, y):
    fxy = x*x + 2*y*y + 2*x*y + 3*x - y - 2
    return fxy


def non_linear_func_3(x, y):
    fxy = 0.5 * (x ** 2 - y ** 2)
    return fxy


def non_linear_func_4(x, y):
    fxy = x**4 + 2*y**4 + 3*x**2*y**2 + 4*x*y**2 + x*y + x + 2*y + 0.5
    return fxy


def limite_func_1(x, y):
    fxy = 2 * x + y - 5
    return fxy


def limite_func_2(x, y):
    fxy = x**2 + y**2 - 9
    return fxy


if __name__ == '__main__':
    a, b = sympy.symbols('x, y')
    eq_limit_solve()