C++ 中的快速自动微分
4.67/5 (10投票s)
一篇关于快速自动微分 (FAD) 的文章，这是计算函数导数的另一种方法。
引言
快速自动微分 (FAD) 是除众所周知的符号微分和有限差分方法之外，计算函数导数的另一种方法。虽然不如这两种方法流行，但 FAD 可以很好地补充它们。它专门用于区分仅仅由以您喜欢的任何编程语言编写的代码片段定义的函数。
在这里，我想概述这种方法的优点，并描述我自己在 C++ 中实现的 FAD 算法。
FAD 可以帮助在您的程序中嵌入微分功能。您实际上不需要了解它的工作原理就能有效地使用它。（尽管如此，您可能仍然应该了解微分是什么以及如何使用函数的导数。）当您需要对迭代或递归定义的函数进行微分，或者当相应的代码片段包含分支等时，此方法真正大放异彩。
这里 CodeProject 上有一篇关于如何实现符号微分的出色文章。它还包含了一些记号约定和有用函数的微分公式，所以我不会在这里重复。
这里 dy/dx 表示 y 关于 x 的导数（如果 y 还依赖于除 x 之外的参数，则为偏导数）。
FAD 基础
要理解自动微分的工作原理，您只需要知道两件事。幸运的是，它们很容易理解。（尽管我再说一遍，它们并不是完全有必要才能使用 FAD。）第一个是所谓的微分的“链式法则”。
假设我们有一个复合函数 y(x) = f(g(x))，并且 g(x) 和 f(t) 分别在 x == x0 和 g(x0) 处可微。那么 y 在 x == x0 处也可微，其导数为
dy/dx = df/dt * dg/dx,
其中 dg/dx 是为 x == x0 计算的，df/dt 是为 t == g(x0) 计算的。如果反过来 g() 函数是复合的（例如 g(x) = h(k(x))），我们可以重复应用此规则，扩展“导数链”，直到我们得到非复合的函数（“基本函数”）
dy/dx = df/dt * dh/du * dk/dx,
其中 x == x0, t == g(x0), u == k(x0)。
请注意，f()、h() 和 k() 函数是“基本函数”（如 sin、cos、exp 等），因此我们已经知道了它们导数的精确公式。
此规则可以推广到多于一个参数的函数。在这种情况下，公式稍微复杂一些，但原理保持不变。
现在回想一下，许多对浮点值进行操作的计算过程无非就是对某个复合函数进行求值。对于给定的输入数据集，相应的程序执行一系列“基本”运算，如 +, -, *, /, sin, cos, tan, exp, log 等。结论是简单明了的。如果我们能够像处理复合函数一样处理这个序列，那么我们就能计算出它的导数，无论它有多长。快速自动微分正是这样做的。（有关更多详细信息，请参阅[1] 综述。近年来，有关此主题的出版物很多。FAD 背后的思想实际上并不新。）
一些 FAD 算法需要某种方式存储此序列以供进一步处理，而其他算法则不需要。用于保存此序列的最流行的数据结构之一是“计算图”（与上面描述的链式法则一起理解 FAD 的第二个重要概念）。在计算图中 (CG)，节点对应于基本运算，而弧线将这些运算与其参数连接起来。请参见下文的图，该图显示了函数 z(x) = f(g(h(x), k(x))) 的示例 CG。（"in" 表示任何输入操作。）
如果代码片段包含分支、循环、递归等，则不同输入的 CG 可能不同。这就是为什么 CG 不是控制流图或抽象语法树 (AST) 或类似的东西。
自动微分算法在函数求值的同一地点计算函数的导数（与符号方法不同）。一旦记录了 CG，就可以遍历它，从参数到感兴趣的函数（在上例中从 x 到 z，即“前向”），或反向遍历（“后向”）来计算导数。因此，FAD 算法有两个主要类别，即“前向”和“后向”。
“前向” FAD 在一次 CG 遍历中计算每个函数相对于一个给定参数的导数。（尽管某些前向算法不需要实际存储 CG。）
“后向” FAD 在一次 CG 遍历中计算一个给定函数 的整个梯度（相对于每个独立变量的导数）。这在处理具有许多参数的少数函数时可能很有益。例如，在矩阵运算等中可能会出现此类函数。
“前向”类似 FAD 的算法也可以用来估计计算过程中累积的函数值不准确性。正如您可能知道的，浮点值通常以有限的相对精度存储。例如，以双精度 IEEE 兼容浮点数存储的 1.0 和 (1.0 + 1.0e-16) 通常被认为是相同的。这些“舍入误差会累积并导致计算值存在不确定性。舍入误差估计 (REE) 算法估计这些不确定性（使用类似于区间算术的方法）。
C++ 实现
在 C++ 中实现 FAD 算法的一种常见方法是为特殊的“活动变量”类使用重载的运算符和函数。（当然，还有其他方法。其中一些方法严重依赖于解析，这超出了本文的范围。您可以在 [2] 中找到有关这些方法的信息。）
考虑以下代码片段
double x0, y0;
double x, y, f;
// get x0 and y0 from somewhere (ask the user, for example)
//...
x = x0; // initialize the arguments (x and y)
y = y0;
f = sin(y * cos(x));
为了使用 FAD 算法，您需要将感兴趣的变量（x、y 和 f）的类型从 double 更改为 FAD 库提供的 CADouble（“活动双精度”）。您还需要标记活动代码段的开始和结束（定义要微分的函数以及其参数初始化的代码片段）。然后您将得到以下结果
double x0, y0;
CADouble x, y, f; // instead of "double x, y, f;"
CActiveSection as;
// mark the beginning of the 'active section'
as.begin(); // start recording operations in a computational graph (CG)
x = x0; // initialize the arguments (x and y)
y = y0;
f = sin(y * cos(x));
as.end(); // stop recording, the active section ends
这些代码片段看起来非常相似，不是吗？在这两种情况下，都会计算 f(x, y)，但在第二种情况下，还会记录基本运算的相应序列（CG）。稍后可以使用 CG 来计算 f() 的导数，如下所示
double out_x, out_y; // here the derivative values will be stored // compute df/dx partial derivative using forward mode FAD as.forward(x, f, out_x); // compute df/dy partial derivative using forward mode FAD as.forward(y, f, out_y);
为 CADouble 实例重载了有用的运算符和函数，因此这些操作不仅执行其通常的含义，而且还会将自身记录在 CG 中。
就是这样。只需稍微更改您的代码，即可嵌入微分功能。当然，天下没有免费的午餐。使用 FAD 可能会使受影响的代码段的执行速度减慢 2 到 400 倍（甚至更多，具体取决于您使用的 FAD 算法、您的硬件等）。遍历大型 CG 会涉及大量内存访问操作，这些操作当然可能比简单地执行一系列数学运算慢。
请注意，不需要实际存储 CG 的方法（所谓的“快速前向”和“快速 REE”算法）通常速度更快，所需内存也更少。
使用代码（“Reed 库”）
我编写的“Reed 库”提供以下内容
- 活动变量和活动段的类；
- 基于 CG 的 FAD 方法：前向和后向（均用于计算一个或多个双精度参数的函数的一阶导数）；
- 基于 CG 的舍入误差估计 (REE) 方法；
- “快速前向”微分方法和“快速 REE”算法的实现，这些算法不需要 CG；
- 一些实用工具（CG 缓存，如果可能，在不重建 CG 的情况下更新以适应新的参数值，等等）。
（“Reed” 是“Rounding Error Estimation and Differentiation”的缩写。事实证明，在我的工作中使用 FAD 时，REE 功能是强制性的，而 FAD 似乎是可选的。这就是为什么 REE 在这个名字中排在第一位。FAD 在我身边也多次很有帮助。）
您可以在提供的 Doxygen 生成的文档中找到类和方法的详细参考（请参阅 /reed/doc/reed library v01.chm）。我在这里只提及要点。
库中的大多数有用功能都位于 namespace reed 中。
目前支持的基本运算包括
- 算术运算（
+, -, *, /），以及复合运算，如+=, -=, *=, /=； fmin(x, y)和fmax(x, y)- 两个值（x和y）的最小值和最大值；- 条件赋值（如果
x为正，则condAssign(x, y) == y，否则为0）； abs(x)-x的绝对值；sin(x), cos(x), tan(x);exp(x), log(x);sqrt(x).
上述运算中有一些不是处处可微的（例如，abs(x) 在 x == 0 处没有导数）。在这些点求值时，仍会计算某些值，系统不会崩溃，但当然，这个值是不可信的。（也许我会在未来版本中添加一个标志或类似的东西来指示这种情况。不幸的是，我现在没有时间。）
关系运算（==, >, < 等）也为 CADouble 和 CADoubleF<...> 重载。在基于 CG 的方法中，它们用于检测 CG 结构是否需要根据参数的新值进行更改。在“快速”方法中，它们仅返回相应浮点值比较的结果。
可以通过 setValue() 和 getValue() 方法访问活动变量的值。
reed::CADouble 类的实例可以在活动段内和活动段外创建。但是（出于效率原因）您只能在活动段内为它们赋值并使用重载的运算符。未能这样做将在调试模式下导致断言，并在发布模式下导致未定义行为。
为了存储 CG，我选择创建一个专门的固定大小对象内存分配器，它比使用默认的 new 和 delete 快得多。请参阅 reed/mem_pool.h。您可以在例如 [3] 中找到有关构建此类分配器的信息。
所有基于 CG 的 FAD 算法都实现为 reed::CActiveSection 类的成员。其中包括
int forward(CADouble arg, CADouble func, double& out);
使用前向模式 FAD 计算函数 func 相对于参数 arg 的导数，并将结果返回到 out。返回值目前不重要。
int reverse(CADouble arg, CADouble func, double& out);
使用后向模式执行与上述相同的操作。
int ree_lower(std::vector<CADouble>& arg, std::vector<double>& err);
计算误差的下界估计。默认情况下，参数应以机器精度存储（即，对于 IEEE 兼容的规范化双精度值，相对精度约为 1.1e-16）。您可以通过设置 arg 中变量的初始绝对误差值来覆盖此设置。这些值包含在 err 向量中。ree_lower() 返回后，您可以通过 CADouble::getAcc() 方法检索任何变量的累积误差值，例如 double e = f.getAcc()。
还有其他重载的此方法版本可能很有用。有关详细信息，请参阅本文档和本文提供的示例。
要使用这些方法，您应该在代码中 #include reed/adouble.h 和 reed/asection.h，并将您的项目链接到 reed/lib/reed.lib 库（或调试版本的 reed_d.lib）。请参阅提供的示例（“simple”、“iter”和“gradient”）。
“快速”前向和 REE 算法实现在 reed/adouble_f.h 中。您无需任何活动段或 CG 即可使用它们。也不需要 reed.lib。只需在项目中 #include reed/adouble_f.h，将 util.cpp 添加到其中，即可完成所有操作。“iter_f”示例显示了如何使用这些“快速”算法。
特殊情况
如前所述，您可以使用 FAD 算法（包括基于 CG 和“快速”版本）对迭代定义的函数、带分支的函数等进行微分，前提是函数在感兴趣的点处可微。考虑以下示例（有关完整实现，请参阅示例 zip 存档中的“iter_f”项目）。假设 y = y(x) 定义如下
double x0 = ...; // set the argument value or ask the user for it size_t m = ...; // Number of iterations to perform. Let the user define it too. // Note: the value of m can be unknown at the compile-time. FAD will handle this. double x, y; x = x0; //argument y = x; for (size_t i = 0; i < m; i++) { if (y < 0.5) { y = 4.0 * y * (1.0 - y); } else if (y > 0.5) { y = 1.0 - 0.75 * y * y; } else { y = -1.0; std::cout << "1/2 has been hit!\n"; // shouldn't get here break; } }
需要在 x == x0 处计算 dy/dx？没问题。生成的代码片段使用“快速”前向方法，但您也可以使用基于 CG 的方法。更改类型如上所述，设置初始导数值（dx/dx == 1.0），您将得到
double x0 = ...; // set the argument value or ask the user for it size_t m = ...; // Number of iterations to perform. Let the user define it too. CADoubleFF x, y; // for REE use CADoubleRee instead. // Active sections are not necessary for "fast forward" and "fast REE" // algorithms because these algorithms do not require CG to be stored. x = x0; //argument x.setAcc(1.0); // Accumulated derivative is 1.0 in x (dx/dx). By default it is zero. y = x; for (size_t i = 0; i < m; i++) { if (y < 0.5) { y = 4.0 * y * (1.0 - y); } else if (y > 0.5) { y = 1.0 - 0.75 * y * y; } else { y = -1.0; std::cout << "1/2 has been hit!\n"; // shouldn't get here break; } }
dy/dx 的值将与 x == x0 处的 y(x) 一起计算。然后，您可以通过 y.getAcc() 和 y.getValue() 分别检索这些值。很简单，不是吗？
示例说明
“FADcpp_examples” Visual Studio 解决方案包含 5 个项目。它仅在 MS Visual C++ 7.1 下进行了测试。（由于我使用的一些模板技巧，它可能无法在 6.0 或更早版本中运行。尽管我尚未在 6.0 下进行测试。）
“FADcpp_examples/reed”目录包含“Reed 库”本身（代码和文档）。您可以使用此解决方案的 “reed_lib” 项目重新构建该库。
其他 4 个项目是使用该库的示例
- “simple” - 一个简单的示例。计算
f = sin(y * cos(x))的导数，x 和 y 由用户提供。使用前向和后向基于 CG 的算法。也演示了 REE。还展示了如何重用现有的 CG 来计算另一个点的导数（如果可能）。 - “iter”。微分迭代定义的函数（见上文）。用户提供参数和迭代次数。
- “iter_f”。与“iter”相同，只是使用了“快速”前向和 REE 算法而不是基于 CG 的算法。
- “gradient”。使用后向 FAD 算法计算二次函数
f(x) = (x * Ax)的梯度，其中A是一个 3x3 矩阵，x = (x0, x1, x2)。在这种情况下，仅使用一次 CG 遍历。（如果使用前向算法，则需要 3 次 CG 遍历）。
结论
本文概述了所谓的快速自动微分 (FAD)。当您需要在程序中嵌入微分功能和/或处理具有分支、循环和递归等的函数时，FAD 会很有帮助。这样，自动微分就可以补充符号微分。后者主要在函数“封闭表达式”的公式已知时使用。如果不是，FAD 通常可以提供帮助。希望我实现的快速自动微分和舍入误差估计算法能对您有所帮助。
试试吧！欢迎提出任何反馈。
参考文献和推荐阅读
- M.Iri. 自动微分和舍入误差估计的历史。自动微分算法：理论、实现和应用，A.Griewank 和 G.F.Corliss，编，SIAM，费城，PA，1991。
- C.Bischof, M.Buecker. 计算计算机程序的导数。
- D.Bulka, D.Mayhew. 高效 C++ 性能编程技术。Addison Wesley，1999。
历史
- v.0.1 本文和示例的早期版本可用。