问题引入¶
实际问题中，经常需要对一个目标函数寻找使输出最大（或最小）的输入。这类问题称为最优化问题（Optimization Problem），其一般形式为：
$$\min_{x \in \mathbb{R}^n} f(x)$$
其中 $f(x)$ 称为目标函数（Objective Function），$x$ 为待求的决策变量。求最大值等价于 $\min -f(x)$。
最优化问题广泛存在于工程、经济、机器学习等领域。本讲将围绕上一讲的直线拟合问题，展示多种求解方法。
准备工作¶
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
import jax
from ipywidgets import interact, IntSlider, FloatSlider, Dropdown
from numba import jit, njit
plt.rcParams['font.family'] = 'Noto Sans CJK SC'
直线拟合的多种方法¶
上一讲做了一些直线拟合（线性回归）。可以用不同的方法来求解这个问题。
探索性数据分析¶
先生成一些含噪声的数据。
In [2]:
# 固定参数生成数据
rng = np.random.default_rng(42)
n = 50
x = np.sort(rng.choice(np.arange(-10, 101), n, replace=False))
y = 5/9 * x - 17.7777777 + 5 * rng.standard_normal(n) # 含噪声的摄氏度数据
x, y
Out[2]:
(array([-5, -4, -3, -2, 0, 4, 5, 6, 11, 17, 18, 20, 25, 26, 27, 30, 31,
38, 39, 42, 44, 46, 47, 50, 53, 55, 57, 58, 61, 62, 63, 64, 65, 68,
69, 71, 74, 75, 76, 80, 81, 85, 86, 87, 88, 92, 93, 94, 97, 99]),
array([-17.16098766, -19.66210457, -17.99884737, -15.73244768,
-25.0635568 , -17.15391156, -17.35186319, -17.63883361,
-13.04237785, -0.8586267 , -12.10693328, -1.82527482,
-12.30323767, -5.00775841, -1.96401237, 1.82000062,
3.00057742, 7.30006959, 2.14526361, 3.24379667,
10.95654615, 6.82125623, 1.95490179, 4.33356401,
7.06940531, 15.26358158, 14.60101765, 17.89687129,
13.97484796, 17.4593652 , 20.35017427, 16.23104516,
20.6172096 , 16.69037037, 18.7402864 , 19.75797727,
17.35413518, 26.32375137, 22.09743282, 26.72913734,
29.62595559, 31.6771004 , 33.32692562, 30.06312821,
28.99461963, 32.93474236, 25.4522168 , 27.20888216,
29.49761313, 32.23598816]))
In [3]:
fig, ax = plt.subplots()
ax.scatter(x, y, c='red', s=15)
ax.set_xlabel('°F')
ax.set_ylabel('°C')
ax.set_title('含噪声的华氏度-摄氏度数据')
plt.show()
交互式数据探索¶
调节下面的滑块，观察数据点数量和噪声强度对拟合结果的影响。
In [4]:
@interact(n_pts=IntSlider(min=3, max=100, step=1, value=50, description='数据点数'),
noise_level=FloatSlider(min=0, max=20, step=0.5, value=5, description='噪声强度'))
def interactive_eda(n_pts, noise_level):
rng_local = np.random.default_rng(42)
xi = np.sort(rng_local.choice(np.arange(-10, 101), n_pts, replace=False))
yi = 5/9 * xi - 17.7777777 + noise_level * rng_local.standard_normal(n_pts)
# 统计学家公式
m_stat = np.cov(xi, yi, ddof=1)[0, 1] / np.var(xi, ddof=1)
b_stat = np.mean(yi) - m_stat * np.mean(xi)
fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(xi, yi, c='red', s=15, label='含噪声数据')
ax.plot(xi, m_stat * xi + b_stat, 'b-', lw=2, label=f'拟合直线: b={b_stat:.2f}, m={m_stat:.4f}')
ax.plot(xi, 5/9 * xi - 17.7777777, 'g--', alpha=0.5, label='理论直线')
ax.set_xlabel('°F')
ax.set_ylabel('°C')
ax.legend(loc='upper left')
ax.set_title(f'n={n_pts}, noise={noise_level:.1f}')
plt.show()
interactive(children=(IntSlider(value=50, description='数据点数', min=3), FloatSlider(value=5.0, description='噪声强度…
最小二乘拟合直线¶
给定数据 $x_i$ 和观测值 $y_i$，最小二乘拟合一条直线意味着找到最佳的 $m$（斜率）和 $b$（截距），使得误差的平方和最小：
$$\min_{m,b} \sum_{i=1}^n \bigl[(b + m x_i) - y_i\bigr]^2$$
直接公式¶
统计学公式¶
用协方差和方差直接计算斜率和截距：
In [5]:
# 统计学公式
m_stat = np.cov(x, y, ddof=1)[0, 1] / np.var(x, ddof=1)
b_stat = np.mean(y) - m_stat * np.mean(x)
print(f'截距 b = {b_stat:.6f}')
print(f'斜率 m = {m_stat:.6f}')
截距 b = -17.416599 斜率 m = 0.531322
线性代数公式¶
更简洁，但需要线性代数知识。好处是可以推广到更复杂的模型。
In [6]:
# 线性代数的最小二乘解
A = np.column_stack([np.ones_like(x), x])
b_la, m_la = np.linalg.lstsq(A, y)[0]
print(f'截距 b = {b_la:.6f}')
print(f'斜率 m = {m_la:.6f}')
截距 b = -17.416599 斜率 m = 0.531322
scipy.optimize.minimize¶
可以使用 SciPy 提供的工具直接求解
$$\min_{b,m} \sum_{i=1}^n [(b + m x_i) - y_i]^2$$
即 $\min_{b,m} \; \mathrm{loss}(b, m)$
In [7]:
# 定义损失函数
def loss(params):
b, m = params
return np.sum((b + m * x - y) ** 2)
# 用默认方法 (BFGS) 优化，初始猜测 [0, 0]
result = optimize.minimize(loss, [0.0, 0.0])
print(f'优化成功: {result.success}')
print(f'截距 b = {result.x[0]:.6f}')
print(f'斜率 m = {result.x[1]:.6f}')
print(f'最小损失 = {result.fun:.6f}')
优化成功: True 截距 b = -17.416598 斜率 m = 0.531322 最小损失 = 689.609046
高阶函数与计算能力¶
优化 $\min_{b,m} \mathrm{loss}(b, m)$ 是一个典型的“函数的函数”——输入是一个函数 loss，输出是最小值的位置或值。
这通常需要大量的计算资源，在计算机的计算能力还不够强大时无法用于求解实际问题。
现代机器学习和更多领域的发展，正是因为计算机的计算能力飞速发展。
使用 SLSQP 方法¶
SLSQP（序列最小二乘规划）是一种支持约束的优化方法。除了目标函数外，还可以指定等式约束和不等式约束：
$$\min_x f(x) \quad \text{s.t.} \quad g_i(x) \geq 0, \; h_j(x) = 0$$
In [8]:
# 使用 SLSQP 方法求解（支持约束的优化方法）
result_slsqp = optimize.minimize(loss, [0.0, 0.0], method='SLSQP')
print(f'无约束: b = {result_slsqp.x[0]:.6f}, m = {result_slsqp.x[1]:.6f}')
# 加入约束：斜率 m 必须等于 5/9（即理论值）
constraints = [{'type': 'eq', 'fun': lambda params: params[1] - 5/9}]
result_eq = optimize.minimize(loss, [0.0, 0.0], method='SLSQP', constraints=constraints)
print(f'等式约束 (m=5/9): b = {result_eq.x[0]:.6f}, m = {result_eq.x[1]:.6f}')
# 加入约束：斜率 m 不能超过 0.5
constraints_ineq = [{'type': 'ineq', 'fun': lambda params: 0.5 - params[1]}]
result_ineq = optimize.minimize(loss, [0.0, 0.0], method='SLSQP', constraints=constraints_ineq)
print(f'不等式约束 (m≤0.5): b = {result_ineq.x[0]:.6f}, m = {result_ineq.x[1]:.6f}')
无约束: b = -17.416598, m = 0.531322 等式约束 (m=5/9): b = -18.611309, m = 0.555556 不等式约束 (m≤0.5): b = -15.872430, m = 0.500000
In [9]:
fig, ax = plt.subplots()
ax.scatter(x, y, c='red', s=15, label='数据')
ax.plot(x, result_slsqp.x[1] * x + result_slsqp.x[0], 'b-', lw=2, alpha=0.7,
label=f'无约束: b={result_slsqp.x[0]:.2f}, m={result_slsqp.x[1]:.4f}')
ax.plot(x, result_eq.x[1] * x + result_eq.x[0], 'm--', lw=2,
label=f'等式约束 m=5/9: b={result_eq.x[0]:.2f}')
ax.plot(x, result_ineq.x[1] * x + result_ineq.x[0], 'c-.', lw=2,
label=f'不等式约束 m≤0.5: b={result_ineq.x[0]:.2f}')
ax.plot(x, 5/9 * x - 17.7777777, 'g:', alpha=0.5, label='理论直线')
ax.set_xlabel('°F')
ax.set_ylabel('°C')
ax.legend(fontsize=8)
ax.set_title('SLSQP 约束优化对比')
plt.show()
梯度¶
上面的优化方法没有显式地提到导数或梯度 (gradients) 信息。梯度对最优化问题的求解非常有用。
对于简单问题，梯度可以手动计算，但对于许多实际问题这是不现实的。
手动计算¶
$$\frac{\partial}{\partial b}\sum_{i=1}^n [(b + m x_i) - y_i]^2 = 2\sum_{i=1}^n [(b + m x_i) - y_i]$$
$$\frac{\partial}{\partial m}\sum_{i=1}^n [(b + m x_i) - y_i]^2 = 2\sum_{i=1}^n x_i[(b + m x_i) - y_i]$$
In [10]:
# 手动计算梯度
def grad_loss(b, m):
"""loss(b, m) 的梯度（手算公式）"""
residual = b + m * x - y
db = 2 * np.sum(residual)
dm = 2 * np.sum(x * residual)
return np.array([db, dm])
grad_loss(0.1, 0.3)
Out[10]:
array([ 611.24203323, 7566.21843351])
有限差分计算¶
In [11]:
# 有限差分近似梯度
eps = 1e-9
f0 = loss([0.1, 0.3])
grad_fd = np.array([
(loss([0.1 + eps, 0.3]) - f0) / eps,
(loss([0.1, 0.3 + eps]) - f0) / eps,
])
grad_fd
Out[11]:
array([ 611.24137574, 7566.21830078])
自动微分¶
对于真实问题，手动求导是不现实的，可以采用自动微分技术。
scipy.optimize.approx_fprime 提供了数值梯度的便捷接口，JAX 提供了更完善的自动微分功能。
注意：自动微分不是有限差分，它在精确度上更优，且不存在选择 ε 的困难。
In [12]:
# 使用 scipy 的数值梯度
grad_auto = optimize.approx_fprime([0.1, 0.3], loss)
# 使用 jax 求自动微分
grad_loss_jax = jax.grad(loss)
grad_jax = grad_loss_jax(np.array([0.1, 0.3]))
print(f'手算梯度: {grad_loss(0.1, 0.3)}')
print(f'有限差分梯度: {grad_fd}')
print(f'SciPy 梯度: {grad_auto}')
print(f'JAX 梯度: {grad_jax}')
手算梯度: [ 611.24203323 7566.21843351] 有限差分梯度: [ 611.24137574 7566.21830078] SciPy 梯度: [ 611.24200439 7566.22094727] JAX 梯度: [ 611.2421 7566.2246]
梯度下降的难点¶
使用梯度下降时，步长的选取是难点，对于复杂函数来说很难找到合适的步长来兼顾所有的数据点。
In [13]:
# 手动梯度下降 + 精确线搜索
b_gd, m_gd = 0.0, 0.0
for _ in range(25):
db, dm = grad_loss(b_gd, m_gd)
# 精确线搜索步长
residual = b_gd + m_gd * x - y
eta = np.sum(residual * (db + dm * x)) / np.sum((db + dm * x) ** 2)
b_gd -= eta * db
m_gd -= eta * dm
print(f'截距 b = {b_gd:.6f}')
print(f'斜率 m = {m_gd:.6f}')
截距 b = -17.416598 斜率 m = 0.531322
随机梯度下降¶
每次随机选取一个数据点（或少量几个），沿着该梯度更新参数，这种方法被称为随机梯度下降（SGD，stochastic gradient descent）。这是机器学习中实际使用的方法。
- 思路：不遍历全体数据，每次仅随机选取单个（或小批量）样本计算梯度并更新。
- 地位：现代机器学习与大规模数据训练中最主流、最有效的策略。
In [14]:
# 随机梯度下降
#@jit # 使用 numba 的 JIT 技术加速函数运行速度
def sgd(n, x, y):
b_sgd, m_sgd = 0.0, 0.0
lr = 0.00002
n_iter = 10_000_000
for t in range(n_iter):
i = np.random.randint(0, n) # 随机选取一个数据点
residual = b_sgd + m_sgd * x[i] - y[i]
db = 2 * residual
dm = 2 * x[i] * residual
b_sgd -= lr * db
m_sgd -= lr * dm
return b_sgd, m_sgd
b_sgd, m_sgd = sgd(n, x, y)
print(f'截距 b = {b_sgd:.6f}')
print(f'斜率 m = {m_sgd:.6f}')
截距 b = -17.410907 斜率 m = 0.547964
Numba：即时编译加速¶
上面的 SGD 函数中使用了 @jit 装饰器。它来自 Numba 库，是一个面向科学计算的 即时（Just-In-Time, JIT）编译器。
工作原理：Numba 在函数第一次被调用时，将其中的 Python 代码编译为机器码。对于包含大量循环和数值运算的函数，编译后的版本可以接近 C 语言的速度。
常用装饰器：
| 装饰器 | 说明 |
|---|---|
@jit |
基础 JIT 编译，兼容性最好 |
@njit |
等价于 @jit(nopython=True)，强制不回退到 Python，速度最快 |
@vectorize |
将标量函数向量化，像 NumPy 的 ufunc 一样使用 |
In [15]:
import time
def python_sum_of_squares(n):
total = 0.0
for i in range(n):
total += i * i
return total
@njit
def numba_sum_of_squares(n):
total = 0.0
for i in range(n):
total += i * i
return total
N = 10_000_000
start = time.perf_counter()
python_sum_of_squares(N)
t_py = time.perf_counter() - start
start = time.perf_counter()
numba_sum_of_squares(N)
t_numba_first = time.perf_counter() - start
start = time.perf_counter()
numba_sum_of_squares(N)
t_numba = time.perf_counter() - start
print(f'纯 Python: {t_py:.4f} s')
print(f'Numba 首次（含编译）: {t_numba_first:.4f} s')
print(f'Numba 编译后: {t_numba:.4f} s')
print(f'加速比: {t_py / t_numba:.1f}x')
纯 Python: 0.5439 s Numba 首次（含编译）: 0.1104 s Numba 编译后: 0.0095 s 加速比: 57.5x
In [16]:
from numba import vectorize, float64
@vectorize([float64(float64, float64)])
def numba_hypot(a, b):
return (a * a + b * b) ** 0.5
a_arr = np.random.randn(1_000_000)
b_arr = np.random.randn(1_000_000)
start = time.perf_counter()
np.hypot(a_arr, b_arr)
t_np = time.perf_counter() - start
start = time.perf_counter()
numba_hypot(a_arr, b_arr)
t_nb_vec = time.perf_counter() - start
print(f'NumPy hypot: {t_np * 1000:.2f} ms')
print(f'Numba @vectorize: {t_nb_vec * 1000:.2f} ms')
NumPy hypot: 2.58 ms Numba @vectorize: 0.46 ms
Numba 使用注意¶
- 首次调用较慢：JIT 编译发生在第一次执行时，后续调用才会快
- 不支持的语法：Numba 的
nopython模式不支持所有 Python 特性（如dict、列表推导式、异常处理等），建议只用于纯数值计算 - NumPy 兼容：大部分 NumPy 函数和数组操作在 Numba 中可以正常使用
- 适用场景：循环密集型、数值运算密集型任务（如 SGD、蒙特卡洛模拟、物理仿真）
交互式 SGD 演示¶
调节学习率和迭代次数，观察随机梯度下降的收敛过程。
In [17]:
@interact(lr=FloatSlider(min=1e-5, max=2e-4, step=1e-5, value=2e-5, readout_format='.5f', description='学习率'),
iterations=IntSlider(min=10000, max=2_000_000, step=50000, value=500000, description='迭代次数'))
def interactive_sgd(lr, iterations):
rng_local = np.random.default_rng(123)
b_s, m_s = 0.0, 0.0
loss_history = []
record_every = max(1, iterations // 200)
for t in range(iterations):
i = rng_local.integers(0, n)
residual = b_s + m_s * x[i] - y[i]
b_s -= lr * 2 * residual
m_s -= lr * 2 * x[i] * residual
if t % record_every == 0:
loss_history.append(loss([b_s, m_s]))
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# 损失收敛曲线
axes[0].plot(loss_history)
axes[0].set_xlabel('迭代步数')
axes[0].set_ylabel('loss')
axes[0].set_ylim(0)
axes[0].set_title('损失函数收敛曲线')
# 拟合结果
axes[1].scatter(x, y, c='red', s=15)
axes[1].plot(x, m_s * x + b_s, 'b-', lw=2, label=f'SGD: b={b_s:.2f}, m={m_s:.4f}')
axes[1].plot(x, 5/9 * x - 17.7777777, 'g--', alpha=0.5, label='理论直线')
axes[1].set_xlabel('°F')
axes[1].set_ylabel('°C')
axes[1].legend(fontsize=8)
axes[1].set_title(f'学习率={lr:.5f}, 迭代={iterations}')
plt.tight_layout()
plt.show()
interactive(children=(FloatSlider(value=2e-05, description='学习率', max=0.0002, min=1e-05, readout_format='.5f',…
3D 损失曲面可视化¶
损失函数 loss(b, m) 是一个二元函数，我们可以用 3D 曲面和等高线来可视化，并观察不同优化算法的路径。
In [19]:
# 损失曲面的 3D 可视化
b_range = np.linspace(-30, 0, 100)
m_range = np.linspace(0, 1, 100)
B, M = np.meshgrid(b_range, m_range)
Z = np.array([[loss([b, m]) for b in b_range] for m in m_range])
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 3D 曲面
ax3d = fig.add_subplot(121, projection='3d')
ax3d.plot_surface(B, M, Z, cmap='viridis', alpha=0.7)
# 标记最优点
ax3d.scatter([b_stat], [m_stat], [loss([b_stat, m_stat])], color='red', s=100, zorder=5)
ax3d.set_xlabel('b (截距)')
ax3d.set_ylabel('m (斜率)')
ax3d.set_zlabel('loss')
ax3d.set_title('损失函数曲面')
# 等高线
ax2d = axes[1]
contour = ax2d.contour(B, M, Z, levels=30, cmap='viridis')
ax2d.scatter([b_stat], [m_stat], color='red', s=100, zorder=5, label='最优解')
ax2d.scatter([0], [0], color='blue', s=50, marker='x', label='初始点')
ax2d.set_xlabel('b (截距)')
ax2d.set_ylabel('m (斜率)')
ax2d.set_title('等高线图')
ax2d.legend()
plt.tight_layout()
plt.show()
本讲小结¶
计算思维：把实际问题建模成最优化问题，使用工具进行求解。
| 方法 | 特点 |
|---|---|
| 统计公式 | 最简单，但仅适用于线性模型 |
线性代数 (lstsq) |
简洁且可推广，但需要线性代数知识 |
scipy.optimize.minimize |
通用优化器，适用于任意非线性问题 |
| SLSQP | 支持约束的优化方法 |
| 梯度下降 | 可手动实现，但步长选择困难 |
| 随机梯度下降 (SGD) | 机器学习的核心算法，每步只用一个数据点 |
最优化是现代计算科学和机器学习的基石。掌握了这些方法，就能理解从简单的直线拟合到神经网络训练背后的统一原理。
使用 numba 可以对 Python 的计算代码进行加速。