使用 Tensorflow 中的有限差分计算四阶导数的准确性答案

【问题标题】：Accuracy in Calculating Fourth Derivative using Finite Differences in Tensorflow使用 Tensorflow 中的有限差分计算四阶导数的准确性
【发布时间】：2021-11-06 12:46:30
【问题描述】：

我正在编写一个小代码来使用张量流中的有限差分法计算四阶导数。如下：

def action(y,x):
    #spacing between points.
    h = (x[-1] - x[0]) / (int(x.shape[0]) - 1)
    
    #fourth derivative 
    dy4 = (y[4:] - 4*y[3:-1] + 6*y[2:-2] - 4*y[1:-3] + y[:-4])/(h*h*h*h)

    return dy4

x = tf.linspace(0.0, 30, 1000)
y = tf.tanh(x)
dy4 = action(y,x)

sess = tf.compat.v1.Session()
plt.plot(sess.run(dy4))

结果如下图：

但是，如果我使用基本相同的代码但只使用 numpy，结果会更清晰：

def fourth_deriv(y, x):
    h = (x[-1] - x[0]) / (int(x.shape[0]) - 1)
    dy = (y[4:] - 4*y[3:-1] + 6*y[2:-2] - 4*y[1:-3] + y[:-4])/(h*h*h*h)
    return dy

x = np.linspace(0.0, 30, 1000)
test = fourth_deriv(np.tanh(x), x)
plt.plot(test)

这给出了：

这里有什么问题？起初我在想点之间的间隔可能太小而无法给出准确的计算，但显然，如果 numpy 可以很好地处理它，情况并非如此。

【问题讨论】：

标签： python numpy tensorflow

【解决方案1】：

这个问题与浮点类型的选择有关。

tf.linspace 自动选择 tf.float32 作为其类型，而
np.linspace 创建一个 float64 数组，该数组具有更高的精度。

进行如下修改：

start = tf.constant(0.0, dtype = tf.float64)
end = tf.constant(30.0, dtype = tf.float64)
x = tf.linspace(start, end, 1000)

导致出现平滑图：

进一步值得注意的是，Tensorflow 确实包含自动微分，这对于机器学习训练至关重要，因此经过充分测试 - 您可以使用梯度磁带来访问它并评估四阶导数，而不会使用有限的数值微分不精确区别：

with tf.compat.v1.Session() as sess2:
  x = tf.Variable(tf.linspace(0, 30, 1000))
  sess2.run(tf.compat.v1.initialize_all_variables())
  with tf.GradientTape() as t4:
    with tf.GradientTape() as t3:
      with tf.GradientTape() as t2:
        with tf.GradientTape() as t1:
          y = tf.tanh(x)

        der1 = t1.gradient(y, x)
      der2 = t2.gradient(der1, x)
    der3 = t3.gradient(der2, x)
  der4 = t4.gradient(der3, x)
  print(der4)

  plt.plot(sess2.run(der4))

这种方法的精度远远优于使用有限差分方法所能达到的精度。下面的代码比较了auto diff的准确率和有限差分法的准确率：

x = np.linspace(0.0, 30, 1000)
sech = 1/np.cosh(x)
theoretical = 16*np.tanh(x) * np.power(sech, 4) - 8*np.power(np.tanh(x), 3)*np.power(sech,2)

finite_diff_err = theoretical[2:-2] - from_finite_diff
autodiff_err = theoretical[2:-2] - from_autodiff[2:-2]

print('Max err with autodiff: %s' % np.max(np.abs(autodiff_err)))
print('Max err with finite difference: %s' % np.max(np.abs(finite_diff_err)))

line, = plt.plot(np.log10(np.abs(autodiff_err)))
line.set_label('Autodiff log error')
line2, = plt.plot(np.log10(np.abs(finite_diff_err)))
line2.set_label('Finite difference log error')
plt.legend()

并产生以下输出：

Max err with autodiff: 3.1086244689504383e-15
Max err with a finite difference: 0.007830900165363808

和下面的图（两条线在 X 轴上大约 600 处重叠）：

【讨论】：

感谢您的精彩回答！也非常感谢有关 tensorflow 梯度的注释。
其实还可以再问一个问题。是否有在 tensorflow 中进行积分的功能？就像你输入一个 x 值和 y 值的数组并返回积分一样？
@Tabin 我不相信有，但我可能弄错了。我的理解是，符号和半符号（类似于 autodiff）机制对于微分来说是简单的并且总是可以解决的，但是非常神秘并且并不总是可以以封闭形式来解决以进行集成。（xkcd.com/2117 是对此的相关喜剧）。 math.stackexchange.com/a/109907/31402 似乎证实了这种算法并不广为人知。
@Tabin "functionality to do an integral" 没有。原因是因为函数的任意组合可以使用链式法则进行微分。所以在 tensorflow 中，如果每个函数都知道它的导数，那么任何函数组合都可以轻松地数值计算它的导数。不存在用于整合功能组合的类似物。

【解决方案2】：

您的问题是 Tensorflow 中的默认精度是 32 位，但 numpy 中的默认精度是 64 位。一个简单的解决方案是按如下方式替换您的 linspace：

x = tf.linspace(tf.constant(0.0, dtype='float64'), 
                tf.constant(30, dtype='float64'), 1000)

【讨论】：