从0开始的机器学习（笔记系列）——导数代码系列

PS

这个系列纯粹笔者学习留下的笔记，寻思下发出来接触一些新东西

内容大纲：

1. 导数的极限定义与直观理解
2. 基本求导法则（含公式）
3. 常见函数的导数示例
4. 数值求导方法：前向差分、中心差分、复步法（complex-step）
5. 用导数求极值：一阶条件与二阶检验 + 数值求根
6. 二阶导数与 Hessian（多元情形）
7. 偏导数、方向导数、梯度场与可视化
8. NumPy 相关经典 API 与实用建议

1. 导数的极限定义与直观理解

函数 ( f(x) ) 在点 ( x ) 处的导数定义为极限：

$$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}$$

直观上，导数表示函数在该点的瞬时变化率，几何意义上对应于函数曲线在该点的切线斜率。

在数值计算中，我们无法让 ( h ) 真正趋近于 0，因此通常使用有限差分来近似导数。不同的差分方案具有不同的截断误差阶：

• 前向差分（Forward difference）：
  $$f^{\prime }(x)\approx \frac{f(x+h)-f(x)}{h},$$
  其截断误差为 ( O(h) )。

• 中心差分（Central difference）：
  $$f^{\prime }(x)\approx \frac{f(x+h)-f(x-h)}{2h},$$
  其截断误差为 ( O(h^2) )，在实践中通常比前向差分更精确。

• 复步法（Complex-step method）（非常推荐）：
  如果函数 ( f ) 可以在复平面上解析延拓，则有
  $$f^{\prime }(x)\approx \frac{\mathrm{Im}\left(f(x+ih)\right)}{h}.$$
  该方法几乎不存在差分消去误差，在极小步长下依然保持良好的数值稳定性。

2. 基本求导法则（常用公式）

给定常数 $c$，以及可导函数 $u(x),v(x)$，常见的一阶求导法则如下：

• 常数规则：
  $(c{)}^{\prime }=0$

• 幂函数法则：
  $(x^n{)}^{\prime }=n{x}^{n-1}$

• 和差法则：
  $(u+v{)}^{\prime }=u^{\prime }+v^{\prime }$

• 常数倍法则：
  $(cu{)}^{\prime }=c{u}^{\prime }$

• 乘法（积）法则：
  $(uv{)}^{\prime }=u^{\prime }v+u{v}^{\prime }$

• 商法则：
  ${\left(\frac{u}{v}\right)}^{\prime }=\frac{u^{\prime }v-u{v}^{\prime }}{v^2}$

• 链式法则（复合函数）：
  若 $y=u(v(x))$，则
  $(u\circ v{)}^{\prime }=(u^{\prime }\circ v)\cdot v^{\prime }$

下面我们将在代码中，通过数值差分与 NumPy 计算，验证并演示这些求导法则在实际计算中的效果。

# 测试以下能不能用
print("Hello Differ!")

# 导入必要库并定义数值求导的常用方法（NumPy 实现）
import numpy as np
import matplotlib.pyplot as plt

# 前向差分
def forward_diff(f, x, h=1e-6):
    return (f(x + h) - f(x)) / h

# 中心差分
def central_diff(f, x, h=1e-6):
    return (f(x + h) - f(x - h)) / (2*h)

# 复步法 (complex-step)
def complex_step_diff(f, x, h=1e-20):
    # 要求 f 能接受复数输入并返回复数输出（NumPy 函数通常可以）
    return np.imag(f(x + 1j*h)) / h

# 二阶中心差分 (近似二阶导)
def second_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. 常见函数的导数示例与数值检验

我们用一些典型函数来对比解析导数与数值导数：

• $f(x)=x^3$ ：$f^{\prime }(x)=3x^2$
• $g(x)=\sin (x)$ ：$g^{\prime }(x)=\cos (x)$
• $h(x)=e^{2x}$ ：$h^{\prime }(x)=2e^{2x}$

下面用 NumPy 的数值方法（中心差分与复步法）进行检验。

# 定义几个函数并比较解析导数与数值导数
def f1(x): return x**3
def f1p_analytic(x): return 3*x**2

def f2(x): return np.sin(x)
def f2p_analytic(x): return np.cos(x)

def f3(x): return np.exp(2*x)
def f3p_analytic(x): return 2*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
def func(x): return np.sin(x) + 0.5 * x**2
def func_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. 用导数求极值（数值流程）一维情形常用流程：

1. 找到一阶导数 (f’(x)=0) 的候选点（数值上可用根搜索或查找导数符号变化）
2. 用二阶导数 (f’‘(x)) 检验：若 (f’‘(x)>0) 为局部极小，若 (f’'(x)<0) 为局部极大。

下面用数值方法演示：先用导数的符号变化（zero-crossings）找到候选点，再用二阶中心差分判断。

# 极值寻找（用符号变化 + 二阶检验）的简单实现
def find_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 （有多个极值）
def test_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 > 0 else ('max' if f2 < 0 else 'inflection/uncertain')
    print(f"x≈{xc:8.5f}, f'(x)≈{fp:8.3e}, f\"(x)≈{f2:8.3e} -> {typ}")

6. 二阶导数与多元情形的 Hessian

对一维函数，二阶导数 $f^{\prime \prime }(x)$ 描述函数在该点附近的曲率信息；
对多元函数，所有二阶偏导数组成 Hessian 矩阵，用于刻画局部二次结构：

$$H_{ij}=\frac{{\partial }^{2}f}{\partial x_i\partial x_j}$$

在数值计算中，可以通过对一阶偏导数再次使用中心差分的方法，逐项逼近 Hessian 矩阵中的各个二阶偏导数。

# 数值计算二维函数的 Hessian（有限差分）
def numeric_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

def numeric_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]])

# 测试二维函数
def f2d(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. 偏导数、方向导数与梯度场可视化

对于多元函数 $f(x,y)$：

• 偏导数：
  $\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}$

• 梯度：
  $\nabla f=\left(\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\right)$

• 方向导数：
  在单位向量 $\mathbf{u}$ 方向上的导数为
  $D_{\mathbf{u}}f=\nabla f\cdot \mathbf{u}$

下面将在一个二维网格上计算并可视化梯度向量场，并在指定点给出方向导数的数值示例。

# 网格上的梯度场与方向导数示例
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 矩阵的特征值，从而判断驻点的极值类型。

外部 / 更专业的库（补充说明）：

• SymPy：用于符号微分，适合解析导数计算与公式推导。
• SciPy：包括 scipy.optimize 中的优化与根查找工具，以及部分数值微分接口（注意其中一些函数可能已被弃用）。
• numdifftools：专门的数值求导工具包，支持高阶导数、Jacobian、Hessian 的数值计算及误差估计。
• autograd / JAX：基于自动微分（AD）的求导工具，在机器学习与科学计算中非常常见，适合需要高效、稳定梯度计算的场景。