《机器学习数学基础》补充资料：拉格朗日乘子法

瑞士数学家欧拉（Leonhard Euler，1707-1783）的大名，如雷贯耳——欧拉，是按德文发音翻译。欧拉不仅是公认的十八世纪最伟大的数学家，还是目前史上最多产的数学家。所著的书籍及论文多达 886 部（篇），涉及非常多领域。平均每年发表约 800 页的高水平学术论文，而且和贝多芬类似，他失明以后发表论文的速度更快了！

如果你生活在十八世纪，就可以把你遇到的数学难题寄到德国柏林科学院，交给非常厉害的欧拉来解决，他几乎什么难题都能解出来。但当时有一个条件极值问题，欧拉长期没能解开。

条件极值问题是这样的，比如说我们想求 $z=xy$ 的极值，却有条件限制：$x^2+xy-y^3=2$。也就是说，我们要从那些满足方程式 $x^2+xy-y^3=2$ 的点中，选取点代入 $z=xy$，看看哪些点代入后会有极大值或极小值。

当时有一位初生牛犊不怕虎的青年，约瑟夫-路易・拉格朗日（Joseph-Louis Lagrange）。他 17 岁以前对数学都没有兴趣，有一天无意中读了英国埃德蒙・哈雷（Edmund Halley，就是命名一颗彗星的那个哈雷）的文章后，开始对数学产生兴趣，并开始学习数学。而在他 19 岁时，便寄了一封信给欧拉，告诉欧拉条件极值问题其实可以轻松解决。欧拉看完他的信后，抱着头说：“天啊，这么简单我怎么没想到！” 随后，拉格朗日声名大噪，当上了皇家炮兵学校的教授，一位 19 岁的教授。日后他也做出了许多非常杰出的研究，拿破仑称赞他是 “数学上高耸的金字塔”。

这个故事告诉我们，即使本来对数学没有兴趣也没关系，我们还是有可能轻松学好数学的。为了帮助你成为现代的拉格朗日，我们先来学习如何求解条件极值，这个方法称为拉格朗日乘子法（Lagrange multiplier method）。

我们先将目标函数设为 $z=f(x,y)=xy$，要求它的极值。至于条件限制，先设为 $g(x,y)=x^2+y^2=1$​。也就是说，我们现在要从单位圆上取点，来找出哪些点会使 $z=xy$ 取得极值。

我们先在坐标平面上，画出 $x^2+y^2=1$。然后在同一张图上，画出目标函数 $z=xy$ 的等高线图。在下图中，画出了 $z=0.1,0.2,0.4,0.6,0.8$​ 时的等高线。

现在要从 $x^2+y^2=1$ 上取点代入 $z=xy$ ，那么取的点，必然是 $x^2+y^2=1$ 和某条等高线的交点。在 $z=0.1,0.2,0.4$ 这几条等高线，都还与条件限制 $x^2+y^2=1$ 有交点。但 $z=0.6,0.8$ 这两条，就超出范围了，我将它们标为红色，表示那不是我们关注的。

将某一条等高线，越往某个方向移动， $z$ 值就会呈现某种趋势，可能是越来越大或者越来越小。在这个例子中，等高线越往圆外， $z$ 值就越大。我们现在要找极值，所以就拼命地将等高线往圆外拉，一直拉到不能再拉为止。什么叫做 “不能再拉” 呢？就是仍然与 $x^2+y^2=1$ 有交点，但再拉就会超出范围了。什么样的状态，会是仍有交点，但再拉就会超出范围呢？就是相切的时候！

所以，我们要找的那些点，就是条件限制 $x^2+y^2=1$ 与目标函数 $z=xy$​ 的等高线相切的那些点。至于，什么叫做 “相切” 呢？相切即是，它们在某个点的法向量是平行的！而该点称为切点。于是，我们现在只要分别求条件限制 $g(x,y)=x^2+y^2=1$ 与各等高线 $f(x,y)=xy=c_k$ 的法向量，并解出使它们平行的条件即可。

如果对一个隐函数取梯度后，再代入点 $P$ ，就会得到它在点 $P$ 处的法向量。而我们的条件限制 $g(x,y)=x^2+y^2=1$ 与各等高线 $f(x,y)=xy=c_k$ ，确实都是隐函数的形式。至于两个向量平行，就是其中一个向量为另一个向量的常数倍。例如 $(3,2)$ 与 $(-6,-4)$ 平行，可写成 $(-6,-4)=-2(3,2)$ 。

所以总结以上，我们需要求解：
$$\nabla f(x,y)=\lambda \nabla g(x,y)$$
也就是
$$\{\begin{array}{l}{f}_x(x,y)=\lambda g_x(x,y)\\ f_y(x,y)=\lambda g_y(x,y)\end{array}$$
解出符合条件的 $(x,y)$ 。至于其中的 $\lambda$ ，称为拉格朗日乘子（Lagrange multiplier）。

为什么要用 $\lambda$ 这个符号，可能是因为拉格朗日（Lagrange）的首字母 $L$ 对应的希腊字母是 $\lambda$ 。

另外也别忘了，我们所找出的点，一定要满足条件限制 $g(x,y)$ 。所以准确地说，应该是：

Lagrange 乘子法（Lagrange multiplier method）

在条件限制 $g(x,y)=k$ 的条件下，求目标函数 $f(x,y)$ 的极值。先求出所有的 $(x,y,\lambda )$ 满足：
$$\{\begin{array}{l}{f}_x(x,y)=\lambda g_x(x,y)\\ f_y(x,y)=\lambda g_y(x,y)\\ g(x,y)=k\end{array}$$
然后解出这些 $(x,y)$，代入到目标函数 $f(x,y)$ 中，得到最大的即为极大值，最小的即为极小值。

例

求 $f(x,y)=x^2+y^2$ 的极值，条件为 $x^2-2x+y^2-4y=0$。

解：

设约束条件为 $g(x,y)=x^2-2x+y^2-4y=0$，根据拉格朗日乘数法，列联立方程组：
$$\{\begin{array}{l}2x=\lambda (2x-2)\text{（梯度分量相等，fx=λgx）}\\ 2y=\lambda (2y-4)\text{（梯度分量相等，fy=λgy）}\\ x^2-2x+y^2-4y=0\text{（满足约束条件）}\end{array}$$
步骤 1：消去拉格朗日乘子 $\mathbf{\lambda }$

从前两个方程解 $\lambda$​（若分母非零）：

• 若 $2x-2\ne 0$（即 $x\ne 1$），则 $\lambda =\frac{2x}{2x-2}=\frac{x}{x-1}$；

• 若 $2y-4\ne 0$（即 $y\ne 2$），则 $\lambda =\frac{2y}{2y-4}=\frac{y}{y-2}$。

令两式相等：
$$\frac{x}{x-1}=\frac{y}{y-2}$$
交叉相乘化简：
$$x(y-2)=y(x-1)⟹xy-2x=xy-y⟹y=2x$$
步骤 2：代入约束条件求解 $(x,y)$

将 $y=2x$ 代入约束条件 $x^2-2x+y^2-4y=0$：
$$\begin{array}{rl}{x}^2-2x+(2x{)}^2-4\cdot (2x)& =0\\ x^2-2x+4x^2-8x& =0\\ 5x^2-10x& =0\\ 5x(x-2)& =0\end{array}$$
解得：

• $x=0$，对应 $y=2\cdot 0=0$；

• $x=2$，对应 $y=2\cdot 2=4$。

步骤 3：验证 $\mathbf{\lambda }=0$ 的特殊情况

若 $\lambda =0$，代入前两个方程：

• $2x=0⟹x=0$；

• $2y=0⟹y=0$。

验证约束条件：$0^2-2\cdot 0+0^2-4\cdot 0=0$，满足条件，故 $(0,0)$ 也是解。

步骤 4：计算函数值并判断极值

将两点代入 $f(x,y)=x^2+y^2$：

• $f(0,0)=0^2+0^2=0$；

• $f(2,4)=2^2+4^2=20$。

几何意义补充（辅助理解）

约束条件配方为：$(x-1{)}^2+(y-2{)}^2=5$

这是以 $(1,2)$ 为圆心、$\sqrt{5}$ 为半径的圆。

• $f(x,y)$ 表示点 $(x,y)$ 到原点的距离的平方。

• 原点到圆心 $(1,2)$ 的距离为 $\sqrt{1^2+2^2}=\sqrt{5}$（等于圆的半径），说明 原点在圆上（对应点 $(0,0)$，距离最小为 0）；

• 点 $(2,4)$ 是圆上离原点最远的点（距离平方为 20）。

因此，$f(0,0)=0$ 是极小值，$f(2,4)=20$​ 是极大值。

拉格朗日乘子法的流程并不困难，就算你不明白它背后的原理，只是把这个方法死记下来也能使用。难就难在解联立方程的时候，技巧千变万化，没有固定的方法。

程序实现

以下是使用Python实现拉格朗日乘子法的代码。

方法1：符号解法（使用SymPy）

import sympy as sp

def lagrange_multiplier(f, constraints, variables):
    """
    符号解法实现拉格朗日乘子法
    :param f: 目标函数
    :param constraints: 等式约束列表 [g1, g2, ...]
    :param variables: 变量列表 [x1, x2, ...]
    :return: 极值点及乘子组成的字典
    """
    # 创建拉格朗日乘子符号
    lambdas = sp.symbols(f'λ1:{len(constraints)+1}')
    
    # 构造拉格朗日函数
    L = f
    for i, con in enumerate(constraints):
        L += lambdas[i] * con
    
    # 生成方程组
    equations = []
    for var in variables:
        equations.append(sp.diff(L, var))
    equations.extend(constraints)
    
    # 求解方程组
    all_vars = list(variables) + list(lambdas)
    solutions = sp.solve(equations, all_vars, dict=True)
    return solutions

# 示例：求 f(x,y) = x² + 2y² 在约束 x + y - 1 = 0 下的极值
if __name__ == "__main__":
    # 定义变量
    x, y = sp.symbols('x y')
    
    # 定义目标函数和约束
    f = x**2 + 2*y**2
    constraints = [x + y - 1]  # 等式约束列表
    
    # 调用拉格朗日乘子法
    solutions = lagrange_multiplier(f, constraints, [x, y])
    
    # 打印结果
    print("符号解法结果：")
    for i, sol in enumerate(solutions):
        print(f"解{i+1}:")
        print(f"  x = {sol[x]}, y = {sol[y]}, λ = {sol[sp.Symbol('λ1')]}")
        print(f"  目标函数值: {f.subs({x: sol[x], y: sol[y]})}\n")

输出示例：

符号解法结果：
解1:
  x = 2/3, y = 1/3, λ = -4/3
  目标函数值: 2/3

方法2：数值解法（使用SciPy）

import numpy as np
from scipy.optimize import minimize

def lagrange_numeric():
    """
    数值解法实现拉格朗日乘子法
    """
    # 目标函数
    def f(X):
        x, y = X
        return x**2 + 2*y**2
    
    # 约束条件（字典形式）
    cons = (
        {'type': 'eq', 'fun': lambda X: X[0] + X[1] - 1}
    )
    
    # 初始猜测点
    initial_guess = np.array([0.0, 0.0])
    
    # 使用SLSQP算法求解
    result = minimize(f, initial_guess, constraints=cons, method='SLSQP')
    
    return result

# 示例求解
if __name__ == "__main__":
    res = lagrange_numeric()
    
    print("\n数值解法结果：")
    print(f"最优解: x = {res.x[0]:.6f}, y = {res.x[1]:.6f}")
    print(f"目标函数值: {res.fun:.6f}")
    print(f"是否成功: {res.success}")
    print(f"迭代次数: {res.nit}")

输出示例：

数值解法结果：
最优解: x = 0.666667, y = 0.333333
目标函数值: 0.666667
是否成功: True
迭代次数: 2

方法3：KKT系统求解法（多约束通用）

from scipy.optimize import fsolve
import numpy as np

def lagrange_kkt(f, gradf, constraints, grad_constraints, initial_guess):
    """
    求解KKT系统的数值方法（支持多约束）
    :param f: 目标函数（未使用，仅为接口一致）
    :param gradf: 目标函数梯度
    :param constraints: 约束函数列表 [g1, g2, ...]
    :param grad_constraints: 约束梯度列表 [∇g1, ∇g2, ...]
    :param initial_guess: 初始猜测 [x1, x2, ..., λ1, λ2, ...]
    :return: 解向量
    """
    n_vars = len(initial_guess) - len(constraints)
    
    def equations(vars):
        # 拆分变量和乘子
        x = vars[:n_vars]
        lambdas = vars[n_vars:]
        
        # 梯度方程: ∇f + Σλi∇gi = 0
        grad_eq = gradf(x) 
        for i, grad_con in enumerate(grad_constraints):
            grad_eq += lambdas[i] * grad_con(x)
        
        # 约束方程: gi(x) = 0
        con_eq = [con(x) for con in constraints]
        
        return np.concatenate([grad_eq, con_eq])
    
    return fsolve(equations, initial_guess)

# 示例：求解多约束问题
if __name__ == "__main__":
    # 问题：min x² + y² + z²
    # 约束: 
    #   x + y - 1 = 0
    #   y + z - 2 = 0
    
    # 定义梯度函数
    gradf = lambda X: np.array([2*X[0], 2*X[1], 2*X[2]])
    
    # 定义约束及其梯度
    constraints = [
        lambda X: X[0] + X[1] - 1,
        lambda X: X[1] + X[2] - 2
    ]
    grad_constraints = [
        lambda X: np.array([1, 1, 0]),
        lambda X: np.array([0, 1, 1])
    ]
    
    # 初始猜测 [x, y, z, λ1, λ2]
    initial_guess = np.array([0.0, 0.0, 0.0, 0.0, 0.0])
    
    # 求解KKT系统
    solution = lagrange_kkt(None, gradf, constraints, grad_constraints, initial_guess)
    
    print("\nKKT系统解法结果：")
    print(f"x = {solution[0]:.6f}, y = {solution[1]:.6f}, z = {solution[2]:.6f}")
    print(f"λ1 = {solution[3]:.6f}, λ2 = {solution[4]:.6f}")
    print(f"目标函数值: {sum(x**2 for x in solution[:3]):.6f}")

输出示例：

KKT系统解法结果：
x = 0.000000, y = 1.000000, z = 1.000000
λ1 = 0.000000, λ2 = -2.000000
目标函数值: 2.000000

方法4：深度学习方法

import numpy as np
import tensorflow as tf
from tensorflow import keras
from scipy.optimize import minimize

# 定义神经网络架构
def create_optimizer_model(input_dim, output_dim):
    """创建用于预测拉格朗日乘子法解的神经网络模型
    
    参数:
        input_dim: 输入维度 (目标函数参数+约束参数)
        output_dim: 输出维度 (解向量维度)
    """
    model = keras.Sequential([
        keras.layers.Input(shape=(input_dim,)),
        keras.layers.Dense(64, activation='relu'),
        keras.layers.Dense(64, activation='relu'),
        keras.layers.Dense(output_dim)  # 输出解向量
    ])
    return model

# 生成训练数据
def generate_training_data(num_samples, problem_dim):
    """生成训练数据: 随机问题参数和对应的数值解
    
    参数:
        num_samples: 样本数量
        problem_dim: 问题维度 (变量数+约束数)
    """
    # 随机生成目标函数参数和约束参数
    problem_params = np.random.rand(num_samples, problem_dim)
    
    # 存储数值解
    solutions = np.zeros((num_samples, problem_dim))
    
    # 为每个样本计算数值解
    for i in range(num_samples):
        # 使用scipy计算数值解 (替代原solve_lagrange)
        res = minimize(
            fun=lambda x: np.sum((x - problem_params[i, :problem_dim//2])**2),  # 示例目标函数
            x0=np.zeros(problem_dim//2),  # 初始猜测
            constraints={'type': 'eq', 'fun': lambda x: np.sum(x) - problem_params[i, -1]},
            method='SLSQP'
        )
        solutions[i, :len(res.x)] = res.x
    
    return problem_params, solutions

# 自定义损失函数
def custom_loss(y_true, y_pred):
    """均方误差损失函数"""
    return tf.reduce_mean(tf.square(y_true - y_pred))

# 主程序
if __name__ == "__main__":
    # 参数设置
    problem_dim = 6  # 问题维度 (例如: 3个变量 + 3个约束参数)
    output_dim = 3   # 输出维度 (解向量维度, 等于变量数)
    num_samples = 1000  # 训练样本数量
    
    # 创建神经网络模型
    optimizer_model = create_optimizer_model(problem_dim, output_dim)
    optimizer_model.compile(optimizer='adam', loss=custom_loss)
    
    # 生成训练数据
    problem_params, solutions = generate_training_data(num_samples, problem_dim)
    
    # 训练神经网络
    print("开始训练神经网络...")
    history = optimizer_model.fit(
        problem_params, 
        solutions, 
        epochs=50, 
        batch_size=32,
        validation_split=0.2
    )
    
    print("\n训练完成!")
    print(f"最终训练损失: {history.history['loss'][-1]:.4f}")
    print(f"最终验证损失: {history.history['val_loss'][-1]:.4f}")
    
    # 测试神经网络预测
    test_params = np.random.rand(1, problem_dim)
    prediction = optimizer_model.predict(test_params)
    print(f"\n测试输入: {test_params[0]}")
    print(f"预测解: {prediction[0]}")

《机器学习数学基础》一书在各大平台有售，请认准作者、出版社