瑞士数学家欧拉(Leonhard Euler,1707-1783)的大名,如雷贯耳——欧拉,是按德文发音翻译。欧拉不仅是公认的十八世纪最伟大的数学家,还是目前史上最多产的数学家。所著的书籍及论文多达 886 部(篇),涉及非常多领域。平均每年发表约 800 页的高水平学术论文,而且和贝多芬类似,他失明以后发表论文的速度更快了!
如果你生活在十八世纪,就可以把你遇到的数学难题寄到德国柏林科学院,交给非常厉害的欧拉来解决,他几乎什么难题都能解出来。但当时有一个条件极值问题,欧拉长期没能解开。
条件极值问题是这样的,比如说我们想求 z=xyz = xyz=xy 的极值,却有条件限制:x2+xy−y3=2x^2 + xy - y^3 = 2x2+xy−y3=2。也就是说,我们要从那些满足方程式 x2+xy−y3=2x^2 + xy - y^3 = 2x2+xy−y3=2 的点中,选取点代入 z=xyz = xyz=xy,看看哪些点代入后会有极大值或极小值。
当时有一位初生牛犊不怕虎的青年,约瑟夫-路易・拉格朗日(Joseph-Louis Lagrange)。他 17 岁以前对数学都没有兴趣,有一天无意中读了英国埃德蒙・哈雷(Edmund Halley,就是命名一颗彗星的那个哈雷)的文章后,开始对数学产生兴趣,并开始学习数学。而在他 19 岁时,便寄了一封信给欧拉,告诉欧拉条件极值问题其实可以轻松解决。欧拉看完他的信后,抱着头说:“天啊,这么简单我怎么没想到!” 随后,拉格朗日声名大噪,当上了皇家炮兵学校的教授,一位 19 岁的教授。日后他也做出了许多非常杰出的研究,拿破仑称赞他是 “数学上高耸的金字塔”。
这个故事告诉我们,即使本来对数学没有兴趣也没关系,我们还是有可能轻松学好数学的。为了帮助你成为现代的拉格朗日,我们先来学习如何求解条件极值,这个方法称为拉格朗日乘子法(Lagrange multiplier method)。
我们先将目标函数设为 z=f(x,y)=xyz = f(x, y) = xyz=f(x,y)=xy,要求它的极值。至于条件限制,先设为 g(x,y)=x2+y2=1g(x, y) = x^2 + y^2 = 1g(x,y)=x2+y2=1。也就是说,我们现在要从单位圆上取点,来找出哪些点会使 z=xyz = xyz=xy 取得极值。
我们先在坐标平面上,画出 x2+y2=1x^2 + y^2 = 1x2+y2=1。然后在同一张图上,画出目标函数 z=xyz = xyz=xy 的等高线图。在下图中,画出了 z=0.1,0.2,0.4,0.6,0.8z = 0.1, 0.2, 0.4, 0.6, 0.8z=0.1,0.2,0.4,0.6,0.8 时的等高线。
现在要从 x2+y2=1x^2 + y^2 = 1x2+y2=1 上取点代入 z=xyz = xyz=xy ,那么取的点,必然是 x2+y2=1x^2 + y^2 = 1x2+y2=1 和某条等高线的交点。在 z=0.1,0.2,0.4z = 0.1, 0.2, 0.4z=0.1,0.2,0.4 这几条等高线,都还与条件限制 x2+y2=1x^2 + y^2 = 1x2+y2=1 有交点。但 z=0.6,0.8z = 0.6, 0.8z=0.6,0.8 这两条,就超出范围了,我将它们标为红色,表示那不是我们关注的。
将某一条等高线,越往某个方向移动, zzz 值就会呈现某种趋势,可能是越来越大或者越来越小。在这个例子中,等高线越往圆外, zzz 值就越大。我们现在要找极值,所以就拼命地将等高线往圆外拉,一直拉到不能再拉为止。什么叫做 “不能再拉” 呢?就是仍然与 x2+y2=1x^2 + y^2 = 1x2+y2=1 有交点,但再拉就会超出范围了。什么样的状态,会是仍有交点,但再拉就会超出范围呢?就是相切的时候!
所以,我们要找的那些点,就是条件限制 x2+y2=1x^2 + y^2 = 1x2+y2=1 与目标函数 z=xyz = xyz=xy 的等高线相切的那些点。至于,什么叫做 “相切” 呢?相切即是,它们在某个点的法向量是平行的!而该点称为切点。于是,我们现在只要分别求条件限制 g(x,y)=x2+y2=1g(x, y) = x^2 + y^2 = 1g(x,y)=x2+y2=1 与各等高线 f(x,y)=xy=ckf(x, y) = xy = c_kf(x,y)=xy=ck 的法向量,并解出使它们平行的条件即可。
如果对一个隐函数取梯度后,再代入点 PPP ,就会得到它在点 PPP 处的法向量。而我们的条件限制 g(x,y)=x2+y2=1g(x, y) = x^2 + y^2 = 1g(x,y)=x2+y2=1 与各等高线 f(x,y)=xy=ckf(x, y) = xy = c_kf(x,y)=xy=ck ,确实都是隐函数的形式。至于两个向量平行,就是其中一个向量为另一个向量的常数倍。例如 (3,2)(3, 2)(3,2) 与 (−6,−4)(-6, -4)(−6,−4) 平行,可写成 (−6,−4)=−2(3,2)(-6, -4) = -2(3, 2)(−6,−4)=−2(3,2) 。
所以总结以上,我们需要求解:
∇f(x,y)=λ∇g(x,y)
\nabla f(x, y) = \lambda \nabla g(x, y)
∇f(x,y)=λ∇g(x,y)
也就是
{fx(x,y)=λgx(x,y)fy(x,y)=λgy(x,y)\begin{cases} f_x(x, y) = \lambda g_x(x, y) \\ f_y(x, y) = \lambda g_y(x, y) \end{cases}
{fx(x,y)=λgx(x,y)fy(x,y)=λgy(x,y)
解出符合条件的 (x,y)(x, y)(x,y) 。至于其中的 λ\lambdaλ ,称为拉格朗日乘子(Lagrange multiplier)。
为什么要用 λ\lambdaλ 这个符号,可能是因为拉格朗日(Lagrange)的首字母 LLL 对应的希腊字母是 λ\lambdaλ 。
另外也别忘了,我们所找出的点,一定要满足条件限制 g(x,y)g(x, y)g(x,y) 。所以准确地说,应该是:
Lagrange 乘子法(Lagrange multiplier method)
在条件限制 g(x,y)=kg(x,y)=kg(x,y)=k 的条件下,求目标函数 f(x,y)f(x,y)f(x,y) 的极值。先求出所有的 (x,y,λ)(x, y, \lambda)(x,y,λ) 满足:
{fx(x,y)=λgx(x,y)fy(x,y)=λgy(x,y)g(x,y)=k \begin{cases} f_x(x,y)=\lambda g_x(x,y) \\f_y(x, y)=\lambda g_y(x,y) \\g(x,y)=k \end{cases} ⎩⎨⎧fx(x,y)=λgx(x,y)fy(x,y)=λgy(x,y)g(x,y)=k
然后解出这些 (x,y)(x,y)(x,y),代入到目标函数 f(x,y)f(x,y)f(x,y) 中,得到最大的即为极大值,最小的即为极小值。
例
求 f(x,y)=x2+y2f(x,y) = x^2 + y^2f(x,y)=x2+y2 的极值,条件为 x2−2x+y2−4y=0x^2 - 2x + y^2 - 4y = 0x2−2x+y2−4y=0。
解:
设约束条件为 g(x,y)=x2−2x+y2−4y=0g(x,y) = x^2 - 2x + y^2 - 4y = 0g(x,y)=x2−2x+y2−4y=0,根据拉格朗日乘数法,列联立方程组:
{2x=λ(2x−2)(梯度分量相等,fx=λgx)2y=λ(2y−4)(梯度分量相等,fy=λgy)x2−2x+y2−4y=0(满足约束条件)
\begin{cases} 2x = \lambda (2x - 2) \quad \text{(梯度分量相等,$f_x = \lambda g_x$)} \\ 2y = \lambda (2y - 4) \quad \text{(梯度分量相等,$f_y = \lambda g_y$)} \\ x^2 - 2x + y^2 - 4y = 0 \quad \text{(满足约束条件)} \end{cases}
⎩⎨⎧2x=λ(2x−2)(梯度分量相等,fx=λgx)2y=λ(2y−4)(梯度分量相等,fy=λgy)x2−2x+y2−4y=0(满足约束条件)
步骤 1:消去拉格朗日乘子 λ\boldsymbol{\lambda}λ
从前两个方程解 λ\lambdaλ(若分母非零):
-
若 2x−2≠02x - 2 \neq 02x−2=0(即 x≠1x \neq 1x=1),则 λ=2x2x−2=xx−1\lambda = \dfrac{2x}{2x - 2} = \dfrac{x}{x - 1}λ=2x−22x=x−1x;
-
若 2y−4≠02y - 4 \neq 02y−4=0(即 y≠2y \neq 2y=2),则 λ=2y2y−4=yy−2\lambda = \dfrac{2y}{2y - 4} = \dfrac{y}{y - 2}λ=2y−42y=y−2y。
令两式相等:
xx−1=yy−2
\dfrac{x}{x - 1} = \dfrac{y}{y - 2}
x−1x=y−2y
交叉相乘化简:
x(y−2)=y(x−1) ⟹ xy−2x=xy−y ⟹ y=2x
x(y - 2) = y(x - 1) \implies xy - 2x = xy - y \implies y = 2x
x(y−2)=y(x−1)⟹xy−2x=xy−y⟹y=2x
步骤 2:代入约束条件求解 (x,y)(x, y)(x,y)
将 y=2xy = 2xy=2x 代入约束条件 x2−2x+y2−4y=0x^2 - 2x + y^2 - 4y = 0x2−2x+y2−4y=0:
x2−2x+(2x)2−4⋅(2x)=0x2−2x+4x2−8x=05x2−10x=05x(x−2)=0
\begin{align*} 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{align*}
x2−2x+(2x)2−4⋅(2x)x2−2x+4x2−8x5x2−10x5x(x−2)=0=0=0=0
解得:
-
x=0x = 0x=0,对应 y=2⋅0=0y = 2 \cdot 0 = 0y=2⋅0=0;
-
x=2x = 2x=2,对应 y=2⋅2=4y = 2 \cdot 2 = 4y=2⋅2=4。
步骤 3:验证 λ=0\boldsymbol{\lambda = 0}λ=0 的特殊情况
若 λ=0\lambda = 0λ=0,代入前两个方程:
-
2x=0 ⟹ x=02x = 0 \implies x = 02x=0⟹x=0;
-
2y=0 ⟹ y=02y = 0 \implies y = 02y=0⟹y=0。
验证约束条件:02−2⋅0+02−4⋅0=00^2 - 2 \cdot 0 + 0^2 - 4 \cdot 0 = 002−2⋅0+02−4⋅0=0,满足条件,故 (0,0)(0, 0)(0,0) 也是解。
步骤 4:计算函数值并判断极值
将两点代入 f(x,y)=x2+y2f(x,y) = x^2 + y^2f(x,y)=x2+y2:
-
f(0,0)=02+02=0f(0, 0) = 0^2 + 0^2 = 0f(0,0)=02+02=0;
-
f(2,4)=22+42=20f(2, 4) = 2^2 + 4^2 = 20f(2,4)=22+42=20。
几何意义补充(辅助理解)
约束条件配方为:(x−1)2+(y−2)2=5(x - 1)^2 + (y - 2)^2 = 5(x−1)2+(y−2)2=5
这是以 (1,2)(1, 2)(1,2) 为圆心、5\sqrt{5}5 为半径的圆。
-
f(x,y)f(x,y)f(x,y) 表示点 (x,y)(x,y)(x,y) 到原点的距离的平方。
-
原点到圆心 (1,2)(1, 2)(1,2) 的距离为 12+22=5\sqrt{1^2 + 2^2} = \sqrt{5}12+22=5(等于圆的半径),说明 原点在圆上(对应点 (0,0)(0,0)(0,0),距离最小为 0);
-
点 (2,4)(2, 4)(2,4) 是圆上离原点最远的点(距离平方为 20)。
因此,f(0,0)=0f(0, 0) = 0f(0,0)=0 是极小值,f(2,4)=20f(2, 4) = 20f(2,4)=20 是极大值。
拉格朗日乘子法的流程并不困难,就算你不明白它背后的原理,只是把这个方法死记下来也能使用。难就难在解联立方程的时候,技巧千变万化,没有固定的方法。
程序实现
以下是使用Python实现拉格朗日乘子法的代码。
方法1:符号解法(使用SymPy)
import sympy as spdef 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 = ffor 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**2constraints = [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 minimizedef lagrange_numeric():"""数值解法实现拉格朗日乘子法"""# 目标函数def f(X):x, y = Xreturn 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 npdef 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 = 0grad_eq = gradf(x) for i, grad_con in enumerate(grad_constraints):grad_eq += lambdas[i] * grad_con(x)# 约束方程: gi(x) = 0con_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.xreturn 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]}")
《机器学习数学基础》一书在各大平台有售,请认准作者、出版社