C / C ++中的衍生物？

我有一些表达式，如x^2+y^2 ，我想用它来进行一些数学计算。 我想做的一件事就是对表达式进行偏导。

因此，如果f(x,y) = x^2 + y^2那么f相对于x的部分将是2x ，相对于y的部分将是2y 。

我使用有限差分方法编写了一个很小的函数，但是我遇到了很多浮点精度问题。 例如，我最终得到1.99234而不是2 。 有没有支持符号差异化的图书馆？ 还有其他建议吗？

Answer 1

我已经用几种不同的语言实现了这样的库，但不幸的是没有C.如果你只处理多项式（求和，乘积，变量，常数和幂），这很容易。 Trigfunction也不算太差。 任何更复杂的事情，你可能会更好地花时间去掌握别人的图书馆。

如果您决定推出自己的产品，我会提出一些简化您生活的建议：

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

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

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

如果您愿意将Lua嵌入到C代码中并在那里进行计算，我将我的Lua代码放在http://www.cs.tufts.edu/~nr/drop/lua上 。 其中一个更好的function是它可以采用符号表达，区分它，并将结果编译成Lua。 你当然会找不到任何文件:-(

Answer 2

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

Answer 3

将数值区分“正确”（在最小化错误的意义上）可能非常棘手。 首先，您可能需要查看关于数值导数的 “数字关联”部分。

对于免费的符号数学包，你应该看看GiNaC 。 您还可以查看SymPy ，这是一个独立的纯python符号数学包。 您会发现SymPy更容易探索，因为您可以从Python命令行以交互方式使用它。

在商业端，Mathematica和Maple都有C API。 您需要已安装/许可的程序版本才能使用这些库，但两者都可以轻松实现您所追求的符号差异化。

Answer 4

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

1. 使用较小的delta。 您似乎使用了大约1e-2的值。 从1e-8开始，测试是否有任何较小的伤害或帮助。 显然你不能太接近机器精度 – 大约1e-16双倍。

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

Answer 5

如果这确实是您想要使用的function，那么编写类库就很容易了。 从单个Term开始，带有系数和指数。 有一个包含术语集合的多项式。

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

Answer 6

利用现有的软件包肯定比编写自己的软件包更容易，但是如果你决定编写自己的软件包，并且你准备花一些时间来学习C ++的一些黑暗角落，你可以使用Boost.Proto从Boost到工程师自己的图书馆。

基本上，Boost.Proto允许您将任何有效的C ++表达式（例如x * x + y * y转换为表达式模板 – 基本上是使用嵌套struct s表示该表达式的解析树 – 然后执行任意任意稍后通过调用proto::eval()对该解析树进行计算。 默认情况下， proto::eval()用于评估树，就像它已经直接运行一样，尽管没有理由不修改每个函数或运算符的行为来取代符号导数。

虽然这对于您的问题来说是一个非常复杂的解决方案，但它比使用C ++模板元编程技术尝试滚动自己的表达式模板要容易得多。

Answer 7

我用C ++创建了这样一个库，你可以在这里看到样本： https ： //github.com/MartinKosicky/auto_diff/blob/master/sample/main.cpp

用法很简单，它也支持矩阵。 我用它来创建递归神经网络…我需要添加GPU矩阵乘法

Answer 8

这是一种旁观，因为它适用于Lisp而不是C / C ++，但它可以帮助其他人寻找类似的任务，或者你可能会得到一些关于在你自己的C / C ++中实现类似的东西的想法。 SICP就lisp主题进行了一些讲座：

1. 衍生规则3b
2. 代数规则4a

在Lisp中，它非常直接（在其他function语言中具有强大的模式匹配和多态类型）。 在C中，您可能不得不大量使用枚举和结构来获得相同的function（更不用说分配/释放）。 人们可以在一小时内明确地编码你需要的东西 – 我会说打字速度是限制因素。 如果你需要C，你实际上可以从C调用ocaml（反之亦然）。

Answer 9

仅计算衍生物的第一顺序是非常微不足道的。 但要快速实现它是一门艺术。 你需要一个包含的类

• 价值
• 一系列导数值与自变量

然后，您将为加法和减法等编写运算符，并执行sin（）等函数，这些函数实现了此操作的基本和众所周知的规则。

要计算高阶导数，应使用截断的泰勒级数。 您还可以将上面提到的类应用于自身 – 值和派生值的类型应该是模板参数。 但这意味着不止一次地计算和存储衍生物。

截断的泰勒系列 – 有两个可用的库：

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

http://www.foelsche.com/ctaylor

Answer 10

看看Theano，它支持符号区分（在神经网络的背景下）。 该项目是开源的，所以你应该能够看到他们是如何做到的。

Answer 11

很抱歉在6年后提出这个问题。 但是，我正在为我的项目寻找这样一个库，我看到@eduffy建议使用FADBAD ++ 。 我已阅读文档并回到您的问题。 我觉得我的答案是有益的，因此，以下代码适用于您的情况。

 #include  #include "fadiff.h" using namespace fadbad; F func(const F& x, const F& y) { return x*x + y*y; } int main() { F 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=fx(); // Value of function double dfdx=fd(0); // Value of df/dx (index 0 of 2) double dfdy=fd(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  #include "fadiff.h" using namespace fadbad; F func(const F& x) { return sin(x); } int main() { F f,x; double dfdx; x = 0.0; x.diff(0,1); f = func(x); dfdx=fd(0); for (int i(0); i < 8; ++i ){ std::cout << " x: " << x.val() << "\n" << " f(x): " << fx() << "\n" << " fadDfdx: " << dfdx << "\n" << "trueDfdx: " << cos(x.val()) << std::endl; std::cout << "==========================" << std::endl; x += 0.1; f = func(x); dfdx=fd(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 ==========================