从0开始的机器学习（笔记系列）——导数代码系列
PS
这个系列纯粹笔者学习留下的笔记，寻思下发出来接触一些新东西
内容大纲：
- 数值求导方法：前向差分、中心差分、复步法（complex-step）
1. 导数的极限定义与直观理解
函数 ( f(x) ) 在点 ( x ) 处的导数定义为极限：
直观上，导数表示函数在该点的瞬时变化率，几何意义上对应于函数曲线在该点的切线斜率。
在数值计算中，我们无法让 ( h ) 真正趋近于 0，因此通常使用有限差分来近似导数。不同的差分方案具有不同的截断误差阶：
- 前向差分（Forward difference）
- 中心差分（Central difference）：
其截断误差为 ( O(h^2) )，在实践中通常比前向差分更精确。 - 复步法（Complex-step method）（非常推荐）：
如果函数 ( f ) 可以在复平面上解析延拓，则有该方法几乎不存在差分消去误差，在极小步长下依然保持良好的数值稳定性。
2. 基本求导法则（常用公式）
给定常数 ，以及可导函数 ，常见的一阶求导法则如下：
常数规则：
幂函数法则：
和差法则：
常数倍法则：
乘法（积）法则：
商法则：
链式法则（复合函数）：
若 ，则
下面我们将在代码中，通过数值差分与 NumPy 计算，验证并演示这些求导法则在实际计算中的效果。
# 测试以下能不能用
print("Hello Differ!")
# 导入必要库并定义数值求导的常用方法（NumPy 实现）
import numpy as np
import matplotlib.pyplot as plt
# 前向差分
defforward_diff(f, x, h=1e-6):
return (f(x + h) - f(x)) / h
# 中心差分
defcentral_diff(f, x, h=1e-6):
return (f(x + h) - f(x - h)) / (2*h)
# 复步法 (complex-step)
defcomplex_step_diff(f, x, h=1e-20):
# 要求 f 能接受复数输入并返回复数输出（NumPy 函数通常可以）
return np.imag(f(x + 1j*h)) / h
# 二阶中心差分 (近似二阶导)
defsecond_central_diff(f, x, h=1e-5):
return (f(x + h) - 2*f(x) + f(x - h)) / (h*h)
# 向量化（便于在数组上批量计算）
v_forward_diff = np.vectorize(lambda xi, f=None: forward_diff(lambda t: f(t), xi), otypes=[float])
# 注意：上面的简易 vectorize 示例我们将在具体使用时用不同方式传 f
print('helper functions ready')
3. 常见函数的导数示例与数值检验
我们用一些典型函数来对比解析导数与数值导数：
下面用 NumPy 的数值方法（中心差分与复步法）进行检验。
# 定义几个函数并比较解析导数与数值导数
deff1(x): return x**3
deff1p_analytic(x): return3*x**2
deff2(x): return np.sin(x)
deff2p_analytic(x): return np.cos(x)
deff3(x): return np.exp(2*x)
deff3p_analytic(x): return2*np.exp(2*x)
xs = np.linspace(-2.0, 2.0, 9) # 一小段点用于展示
h = 1e-6
for fun, anafun, name in [(f1, f1p_analytic, 'x^3'), (f2, f2p_analytic, 'sin(x)'), (f3, f3p_analytic, 'exp(2x)')]:
print('Function:', name)
for x in xs:
num_central = central_diff(fun, x, h=h)
num_complex = complex_step_diff(fun, x, h=1e-20)
analytic = anafun(x)
print(f' x={x:6.3f} analytic={analytic:10.6f} central={num_central:10.6f} complex={num_complex:10.6f}')
print()
4. 绘制函数与导数我们用可视化来直观比较原函数与其导数（用中心差分和复步法）。
# 绘制示例：f(x) = sin(x) + 0.5 x^2
deffunc(x): return np.sin(x) + 0.5 * x**2
deffunc_p_analytic(x): return np.cos(x) + x
x = np.linspace(-6, 6, 400)
dx = x[1] - x[0]
# 数值导数（中心差分）在向量上用 np.gradient 更方便（等距网格）
num_grad = np.gradient(func(x), dx) # np.gradient 在等距网格上近似中心差分
# 复步法逐点计算（更高精度）
num_complex = np.array([complex_step_diff(func, xi, h=1e-20) for xi in x])
plt.figure(figsize=(8,4))
plt.plot(x, func(x), label='f(x)') # no explicit color
plt.plot(x, num_grad, label='num derivative (np.gradient)')
plt.plot(x, num_complex, label='num derivative (complex-step)', linestyle='--')
plt.plot(x, func_p_analytic(x), label='analytic derivative', linestyle=':')
plt.legend()
plt.title('Function and its derivative')
plt.xlabel('x')
plt.grid(True)
plt.show()
5. 用导数求极值（数值流程）一维情形常用流程：
- 找到一阶导数 (f'(x)=0) 的候选点（数值上可用根搜索或查找导数符号变化）
- 用二阶导数 (f''(x)) 检验：若 (f''(x)>0) 为局部极小，若 (f''(x)<0) 为局部极大。
下面用数值方法演示：先用导数的符号变化（zero-crossings）找到候选点，再用二阶中心差分判断。
# 极值寻找（用符号变化 + 二阶检验）的简单实现
deffind_critical_points(f, x_grid, h=1e-5):
# 使用 np.gradient 得到一阶导数的网格值
dx = x_grid[1] - x_grid[0]
fprime = np.gradient(f(x_grid), dx)
# 找到符号变化点（简单方法）
sign = np.sign(fprime)
zero_crossings = np.where(np.roll(sign, -1) * sign <= 0)[0] # 索引
candidates = []
for idx in zero_crossings:
# 拟合候选点在 x_grid[idx]
xc = x_grid[idx]
f2 = second_central_diff(f, xc, h=h)
candidates.append((xc, fprime[idx], f2))
return candidates
# 测试函数：f(x) = x^3 - 3x + 1 （有多个极值）
deftest_fun(x): return x**3 - 3*x + 1
xg = np.linspace(-3, 3, 10001)
crit = find_critical_points(test_fun, xg, h=1e-6)
crit[:10] # 显示几个候选点
上面得到的候选点中，我们再看二阶导数符号来判断极值类型。
for xc, fp, f2 in crit:
typ = 'min'if f2 > 0else ('max'if f2 < 0else'inflection/uncertain')
print(f"x≈{xc:8.5f}, f'(x)≈{fp:8.3e}, f\"(x)≈{f2:8.3e} -> {typ}")
6. 二阶导数与多元情形的 Hessian
对一维函数，二阶导数 描述函数在该点附近的曲率信息；
对多元函数，所有二阶偏导数组成 Hessian 矩阵，用于刻画局部二次结构：
在数值计算中，可以通过对一阶偏导数再次使用中心差分的方法，逐项逼近 Hessian 矩阵中的各个二阶偏导数。
# 数值计算二维函数的 Hessian（有限差分）
defnumeric_gradient_2d(f, x, y, h=1e-6):
# 返回 (df/dx, df/dy) 在点 (x,y)
fx = (f(x+h, y) - f(x-h, y)) / (2*h)
fy = (f(x, y+h) - f(x, y-h)) / (2*h)
return fx, fy
defnumeric_hessian_2d(f, x, y, h=1e-4):
f_xx = (f(x+h, y) - 2*f(x, y) + f(x-h, y)) / (h*h)
f_yy = (f(x, y+h) - 2*f(x, y) + f(x, y-h)) / (h*h)
f_xy = (f(x+h, y+h) - f(x+h, y-h) - f(x-h, y+h) + f(x-h, y-h)) / (4*h*h)
return np.array([[f_xx, f_xy],
[f_xy, f_yy]])
# 测试二维函数
deff2d(x, y): return np.sin(x) * np.exp(-y**2) + x*y
pt = (0.5, -0.3)
grad = numeric_gradient_2d(f2d, *pt)
hess = numeric_hessian_2d(f2d, *pt)
print('point', pt)
print('numeric gradient:', grad)
print('numeric hessian:\n', hess)
# Hessian 的特征值用于判定鞍点/最小/最大
eigvals = np.linalg.eigvals(hess)
print('hessian eigenvalues:', eigvals)
7. 偏导数、方向导数与梯度场可视化
对于多元函数 ：
偏导数：
梯度：
方向导数：
在单位向量 方向上的导数为
下面将在一个二维网格上计算并可视化梯度向量场，并在指定点给出方向导数的数值示例。
# 网格上的梯度场与方向导数示例
X = np.linspace(-2, 2, 41)
Y = np.linspace(-2, 2, 41)
XX, YY = np.meshgrid(X, Y)
ZZ = f2d(XX, YY)
# 使用 np.gradient 计算标量场的偏导数（注意网格间距）
dx = X[1] - X[0]
dy = Y[1] - Y[0]
dZ_dX, dZ_dY = np.gradient(ZZ, dy, dx) # 注意参数顺序：先沿 axis=0 的 spacing，再沿 axis=1 的 spacing
plt.figure(figsize=(6,5))
plt.contourf(XX, YY, ZZ, levels=30)
plt.quiver(XX, YY, dZ_dX, dZ_dY) # 梯度向量场
plt.title('Scalar field and its gradient vectors')
plt.xlabel('x'); plt.ylabel('y')
plt.show()
# 在某点计算方向导数
pt = (0.8, 0.2)
grad_pt = numeric_gradient_2d(f2d, *pt)
u = np.array([1.0, 1.0])
u = u / np.linalg.norm(u) # 单位向量
directional = np.dot(grad_pt, u)
print('point', pt, 'gradient', grad_pt, 'unit direction', u, 'directional derivative', directional)
8. NumPy 在数值求导方面的经典 API 与建议
常用 NumPy 函数 / 技巧：
np.gradient：在等距格点上使用中心差分近似一阶导数，支持多维数组，非常适合标量场梯度的离散近似。np.diff：计算相邻元素之差，常用于构造差分格式，但需要注意结果长度会变短以及边界的处理方式。np.pad：在差分计算中进行边界填充（如 mirror、edge 等策略），有助于统一数组尺寸。np.vectorize、np.apply_along_axis：便于在数组上应用标量函数（注意 vectorize 只是语法糖，并不会带来性能加速）。- 复步法（complex-step）：当函数能够接受复数输入时，这是数值求导中非常稳定且精度极高的方法（具体实现见上文）。
np.linalg.eigvals / np.linalg.eig：用于分析 Hessian 矩阵的特征值，从而判断驻点的极值类型。
外部 / 更专业的库（补充说明）：
SymPySciPy：包括 scipy.optimize 中的优化与根查找工具，以及部分数值微分接口（注意其中一些函数可能已被弃用）。numdifftools：专门的数值求导工具包，支持高阶导数、Jacobian、Hessian 的数值计算及误差估计。autograd / JAX：基于自动微分（AD）的求导工具，在机器学习与科学计算中非常常见，适合需要高效、稳定梯度计算的场景。
