PyTorch与NumPy二阶导数计算的数值差异问题排查

对比PyTorch与NumPy的二阶导数数值计算：采用完全相同的手动差分逻辑，NumPy计算结果符合预期；但PyTorch计算时二阶导数出现数值爆炸，一阶导数存在大量局部振荡。而使用PyTorch的autograd功能计算则无此问题，以下为具体代码、结果及原因分析。

NumPy实现代码与结果 #

import numpy as np
import matplotlib.pyplot as plt

dt  = 1e-6
t   = np.arange(0,10,dt)

u_SZ   = np.cos(2*t)

# ---------------------
# --- NUMPY ---
# ---------------------
t_SZ_c    = t[1:-1]

u_SZ_t    = u_SZ[2:]
u_SZ_c    = u_SZ[1:-1]
u_SZ_b    = u_SZ[:-2]
u_t_SZ_t  = (u_SZ_t - u_SZ_c) / dt
u_t_SZ_b  = (u_SZ_c - u_SZ_b) / dt
u_t_SZ_c  = (u_t_SZ_t + u_t_SZ_b) / 2
u_tt_SZ_c = (u_t_SZ_t - u_t_SZ_b) / dt

ax = plt.subplot(131)
ax.plot(t_SZ_c, u_SZ_c)
ax = plt.subplot(132)
ax.plot(t_SZ_c, u_t_SZ_c)
ax = plt.subplot(133)
ax.plot(t_SZ_c, u_tt_SZ_c)
plt.show(block = True)

PyTorch实现代码与结果 #

dt = 1e-6
t_SZ  = torch.arange(0,10,dt).view(-1,1).requires_grad_(True)
u_SZ = torch.cos(2*t_SZ); 
t_SZ_c    = t_SZ[1:-1]

u_SZ_t    = u_SZ[2:]
u_SZ_c    = u_SZ[1:-1]
u_SZ_b    = u_SZ[:-2]
u_t_SZ_t  = (u_SZ_t - u_SZ_c) / dt
u_t_SZ_b  = (u_SZ_c - u_SZ_b) / dt
u_t_SZ_c  = (u_t_SZ_t + u_t_SZ_b) / 2
u_tt_SZ_c = (u_t_SZ_t - u_t_SZ_b) / dt

t_SZ_c = (t_SZ_c.detach().numpy())
u_SZ_c = (u_SZ_c.detach().numpy())
u_t_SZ_c = (u_t_SZ_c.detach().numpy())
u_tt_SZ_c =(u_tt_SZ_c.detach().numpy())


ax = plt.subplot(131)
ax.plot(t_SZ_c, u_SZ_c)
ax = plt.subplot(132)
ax.plot(t_SZ_c, u_t_SZ_c)
ax = plt.subplot(133)
ax.plot(t_SZ_c, u_tt_SZ_c)
plt.show(block = True)

原因解析 #

1. 浮点精度差异：
   NumPy的np.arange默认生成**float64（双精度）数组，而PyTorch的torch.arange在CPU环境下默认生成float32（单精度）**张量。手动二阶差分对数值精度极其敏感：两次差分操作会将原始浮点数的舍入误差放大1/dt²倍（此处dt=1e-6，误差放大1e12倍）。float32的有效位数仅约7位，远不足以抵消这种误差放大，最终导致二阶导数数值爆炸、一阶导数出现振荡。

2. 计算图追踪的间接影响：
   虽然最终调用detach()脱离计算图，但设置requires_grad_(True)后，PyTorch在计算过程中会为梯度追踪保留中间张量的精度状态，进一步加剧了单精度下的误差累积。而NumPy无计算图追踪，直接进行双精度数值计算，误差可控。

3. autograd的优势：
   PyTorch的autograd采用解析导数（反向传播时直接计算导数的解析表达式），完全避免了数值差分带来的精度损耗，因此结果符合预期。

验证方案 #

将PyTorch张量改为float64精度，即可得到与NumPy一致的结果：

dt = 1e-6
# 指定dtype为float64
t_SZ  = torch.arange(0,10,dt, dtype=torch.float64).view(-1,1).requires_grad_(True)
u_SZ = torch.cos(2*t_SZ)
# 后续差分逻辑不变

内容的提问来源于stack exchange，提问作者Nena