C/C++中的衍生品？

Question

我有一些表达式，比如x^2+y^2，我想用它们进行一些数学计算。其中一件事是要对这些表达式进行偏导数。

例如，如果 f(x,y) = x^2 + y^2 ，那么对于 x 的偏导数将会是 2x，对于 y 的偏导数则是 2y 。

我使用了一个有限差分法编写的小型函数，但是遇到了很多浮点精度问题。例如，我得到的结果是1.99234而不是2。是否有任何支持符号微分的库？还有其他建议吗？

- idontwanttortfm

---

自己编写可能不是一个好主意。我修过一门计算机代数的研究生课程，但事实上还没有完全实现所有功能。简单的微分可以通过应用数学规则，然后尝试进行替换和评估来完成...在C/C++中很难做到。 - Calyth

请参考 https://dev59.com/inRB5IYBdhLWcg3wbGtB， 或在问题中打上 derivative 标签。 - Daniel Daranas

Answer 1

我已经在几种不同的语言中实现了这样的库，但遗憾的是没有C。如果你只处理多项式（求和、乘积、变量、常数和指数），那么它就很容易。三角函数也不算太难。任何更复杂的东西，你可能最好花时间去掌握别人的库。

如果您决定自己开发，我有一些建议，可以简化您的生活：

• 使用不可变数据结构（纯函数式数据结构）表示表达式。

• 使用Hans Boehm 的垃圾收集器来管理内存。

• 为了表示线性求和，使用有限映射（例如二叉搜索树）将每个变量映射到其系数。

如果您愿意将Lua嵌入到C代码中，并在那里进行计算，我已经将我的Lua代码放在http://www.cs.tufts.edu/~nr/drop/lua上。其中一个较好的功能是它可以接受符号表达式，对其进行微分，并将结果编译成Lua。当然，你会发现完全没有文档 :-(

- Norman Ramsey

Answer 2

如果您正在进行数值微分（“在x = x0处评估f(x)的导数”），并且您预先知道方程式（即非用户输入），那么我建议使用FADBAD++。它是一个C ++模板库，用于使用自动微分解决数值导数。 它非常快速和准确。

- eduffy

---

谢谢您推荐这个库。我已经仔细检查并为OP发布了答案。它看起来是一个非常好的库。再次感谢。 - CroCo

Answer 3

在数值微分中，正确性（即最小化误差）的实现可能非常棘手。要开始，请查看Numerical Recipies关于数值导数的部分。

对于免费的符号数学软件包，您应该看看GiNaC。您也可以查看SymPy, 一个自包含的、纯Python符号数学软件包。您会发现SymPy更容易探索，因为您可以在Python命令行中交互地使用它。

在商业端，Mathematica和Maple都有C API。您需要安装/许可程序的版本才能使用库，但两者都可以轻松地进行所需的符号微分。

- Barry Wark

Answer 4

您可以通过两种简单的方法提高数值微分的准确性：

1. 使用更小的delta。您似乎使用了约为1e-2的值。从1e-8开始，测试是否使得更小的值有益或有害。显然，您不能太接近机器精度——双精度约为1e-16。

2. 使用中心差分而不是前向（或后向）差分。 即df_dx=(f(x+delta)-f(x-delta))/(2.0*delta) 由于抵消更高截断项的原因，中心差分估计的误差是前向差分的delta阶。请参见http://en.wikipedia.org/wiki/Finite_difference。

- Ian

Answer 5

如果这确实是您想要使用的函数类型，编写一个类库将非常容易。从单个项开始，包括系数和指数。创建一个多项式，它由一组项组成。

如果您为感兴趣的数学方法定义一个接口（例如，add/sub/mul/div/differentiate/integrate），那么您正在考虑GoF组合模式。Term和Polynomial都将实现该接口。Polynomial将简单地迭代其集合中的每个Term。

- duffymo

Answer 6

利用已有的软件包要比编写自己的软件包更容易，但是如果你决心要编写自己的软件包，并且准备花一些时间学习 C++ 的一些黑暗角落，那么你可以使用Boost中的Boost.Proto来设计你自己的库。

基本上，Boost.Proto 允许你将任何有效的C++表达式（如x * x + y * y）转换为表达式模板——使用嵌套的结构体表示该表达式的语法树，然后通过在稍后调用 proto::eval() 对该解析树执行任意计算。默认情况下，proto::eval() 用于按照直接运行的方式评估该树，但没有理由不能修改每个函数或运算符的行为，以取代符号导数。

虽然这可能是一个极其复杂的解决方案，但它仍然比尝试使用C ++模板元编程技术来滚动自己的表达式模板要容易得多。

- j_random_hacker

---

距离我第一次听说boost-proto已经过去了4年，但我仍然不清楚boost-proto提供了什么价值。 - user678269

@user678269：有两个原因：（1）表现力。您可以使用它轻松地创建AST数据结构，例如通过在代码中逐字输入C++表达式“x * x + y * y”来表示一个AST；如果您已经编写了用于操作AST的代码（例如用于评估、编译、漂亮打印或（对于数学表达式）求导/积分/展开/简化），则可以应用此代码。... - j_random_hacker

(2) 效率。如果你只想评估表达式，通常的 C++ 计算会不必要地创建和销毁许多临时变量；通过构建并立即评估 AST，你通常可以做得更好，AST 的行为有点像一个定制的构造函数，只执行一次内存分配（有关更多详细信息，请参见第一个谷歌搜索结果）。 - j_random_hacker

Answer 7

很抱歉在6年后提出这个问题。然而，我正在寻找一个适合我的项目的库，我看到@eduffy建议使用FADBAD++。我已经阅读了文档并回来回答你的问题。我觉得我的答案会有所帮助，因此，以下代码适用于你的情况。

#include <iostream>
#include "fadiff.h"

using namespace fadbad;

F<double> func(const F<double>& x, const F<double>& y)
{
    return x*x + y*y;
}

int main()
{
    F<double> x,y,f;     // Declare variables x,y,f
    x=1;                 // Initialize variable x
    x.diff(0,2);         // Differentiate with respect to x (index 0 of 2)
    y=1;                 // Initialize variable y
    y.diff(1,2);         // Differentiate with respect to y (index 1 of 2)
    f=func(x,y);         // Evaluate function and derivatives

    double fval=f.x();   // Value of function
    double dfdx=f.d(0);  // Value of df/dx (index 0 of 2)
    double dfdy=f.d(1);  // Value of df/dy (index 1 of 2)

    std::cout << "    f(x,y) = " << fval << std::endl;
    std::cout << "df/dx(x,y) = " << dfdx << std::endl;
    std::cout << "df/dy(x,y) = " << dfdy << std::endl;

    return 0;
}

输出为：

    f(x,y) = 2
df/dx(x,y) = 2
df/dy(x,y) = 2

另一个例子，假设我们对 sin() 的一阶导数感兴趣。从解析上来看，它是 cos。这很好，因为我们需要比较给定函数的真实导数和其数值近似值以计算真实误差。

#include <iostream>
#include "fadiff.h"

using namespace fadbad;

F<double> func(const F<double>& x)
{
    return sin(x);
}



int main()
{
    F<double> f,x;
    double dfdx;
    x = 0.0;
    x.diff(0,1);
    f = func(x);
    dfdx=f.d(0);


    for (int i(0); i < 8; ++i ){
        std::cout << "       x: " << x.val()        << "\n"
                  << "    f(x): " << f.x()          << "\n" 
                  << " fadDfdx: " << dfdx           << "\n"
                  << "trueDfdx: " << cos(x.val())   << std::endl;
        std::cout << "=========================="   << std::endl;

        x += 0.1;
        f = func(x);
        dfdx=f.d(0);
    }


    return 0;
}

结果是

       x: 0
    f(x): 0
 fadDfdx: 1
trueDfdx: 1
==========================
       x: 0.1
    f(x): 0.0998334
 fadDfdx: 0.995004
trueDfdx: 0.995004
==========================
       x: 0.2
    f(x): 0.198669
 fadDfdx: 0.980067
trueDfdx: 0.980067
==========================
       x: 0.3
    f(x): 0.29552
 fadDfdx: 0.955336
trueDfdx: 0.955336
==========================
       x: 0.4
    f(x): 0.389418
 fadDfdx: 0.921061
trueDfdx: 0.921061
==========================
       x: 0.5
    f(x): 0.479426
 fadDfdx: 0.877583
trueDfdx: 0.877583
==========================
       x: 0.6
    f(x): 0.564642
 fadDfdx: 0.825336
trueDfdx: 0.825336
==========================
       x: 0.7
    f(x): 0.644218
 fadDfdx: 0.764842
trueDfdx: 0.764842
==========================

- CroCo

Answer 8

这其实是一个和Lisp有关而非C/C++的旁注，但它或许能帮助其他寻找类似任务的人，或者你可以从中得到一些在C/C++中实现类似功能的想法。SICP关于这个主题有一些关于lisp的讲座：

1. 导数规则 3b
2. 代数规则 4a

在Lisp中，实现相对简单（以及其他拥有强大模式匹配和多态类型的函数式语言）。在C中，你可能需要大量使用枚举和结构体以获得同样的功能（更不用说分配/释放内存）。在ocaml中，你可以在不到一个小时的时间内编写出所需的代码——我认为打字速度是限制因素。如果你需要C，你实际上可以从C中调用ocaml（反之亦然）。

- nlucaroni

---

没错。在 C 语言中进行符号数学运算将是一场噩梦。仅仅是正确地使用 malloc 和 free 就需要比使用高级函数式语言更长的时间。 - Jules

我最初甚至没有考虑到这一点（可能是因为我太习惯使用带有GC的函数式语言了）。但你也是对的。那将是一场噩梦。 - nlucaroni

Mathematica的核心是用C语言编写的。因此，这是可能的。 - alfC

Answer 9

看看Theano，它支持符号微分（在神经网络的上下文中）。该项目是开源的，因此您应该能够看到他们如何实现它。

- Andrei Pokrovsky

Answer 10

计算一阶导数非常容易实现。 但是让它快速运行则需要一定技巧。 您需要一个包含以下内容的类：

• 值
• 独立变量的导数值数组

然后，您可以编写加法和减法等操作符以及像sin()这样实现此操作的基本和众所周知规则的函数。

计算更高阶导数应该使用截断的泰勒级数。 您还可以将上述类应用于自身--值和导数值的类型应该是模板参数。 但这意味着需要计算和存储多次导数。

截断的泰勒级数--有两个库可用于此：

http://code.google.com/p/libtaylor/

http://www.foelsche.com/ctaylor

- user678269