Derivatives in C/C++?

Question

I have some expressions such as x^2+y^2 that I'd like to use for some math calculations. One of the things I'd like to do is to take partial derivatives of the expressions.

So if f(x,y) = x^2 + y^2 then the partial of f with respect to x would be 2x, the partial with respect to y would be 2y.

I wrote a dinky function using a finite differences method but I'm running into lots of problems with floating point precision. For example, I end up with 1.99234 instead of 2. Are there any libraries that support symbolic differentiation? Any other suggestions?

Conceit

Rolling your own would be a bad idea. I've taken a cross-listed grad course in Computer Algebra, and things are not even quite feature complete. Simplistic differentiation could be done by applying maty rules and then try to substitute and evaluate... hard to do in C/C++. – Glochidium

See also #627555 and maybe tag your question derivative – Couchman

Answer 1

I've implemented such libraries in several different languages, but unfortunately not C. If you are dealing only with polynomials (sums, products, variables, constants, and powers) it is pretty easy. Trig functions are also not too bad. Anything more complicated and you will probably be better off taking the time to master somebody else's library.

If you decide to roll your own, I have a few suggestions that will simplify your life:

• Use immutable data structures (purely functional data structures) to represent expressions.

• Use Hans Boehm's garbage collector to manage the memory for you.

• To represent a linear sum, use a finite map (e.g., a binary search tree) to map each variable to its coefficient.

If you're willing to embed Lua into your C code and do your computations there, I have put my Lua code at http://www.cs.tufts.edu/~nr/drop/lua. One of the nicer features is that it can take a symbolic expression, differentiate it, and compile the results into Lua. You will of course find no documentation whatever :-(

Miscegenation

Answer 2

If you're doing numerical differentiation ("evaluate the derivative of f(x) at x = x0") and you know you're equations in advance (ie, not user input), then I'd recommend FADBAD++. It's a C++ template library for solving numeric derivatives using Automatic differentiation. It's very quick, and accurate.

Brittney

Thank you for suggesting this library. I've double checked and posted an answer for the OP. It seems very good library. Thanks again. – Meghannmegiddo

Answer 3

Getting numerical differentiation "right" (in the sense of minimizing errors) can be quite tricky. To get started, you may want to take a look at the Numerical Recipies' section on numerical derivatives.

For free symbolic math packages, you should look at GiNaC. You could also check out SymPy, a self-contained, pure-python symbolic math package. You will find that SymPy is much easier to explore, since you can use it interactively from the Python command line.

On the commercial end, both Mathematica and Maple have C APIs. You need an installed/licensed version of the program to use the libraries, but both can easily do the type of symbolic differentiation you're after.

Assignee

Answer 4

You can improve the accuracy of your numerical differentiation in two simple ways

1. Use a smaller delta. You appear to have used a value of around 1e-2. Start with 1e-8, and test if getting any smaller hurts or helps. Obviously you can't get too close to the machine precision - about 1e-16 for double.

2. Use central differences rather than forward (or backwards) differences. i.e. df_dx =(f(x+delta) - f(x-delta)) / (2.0*delta) For reasons to do with cancellation of higher truncation terms, the error in the central differences estimate is of the order delta^2 rather than the delta of forward differences. See http://en.wikipedia.org/wiki/Finite_difference

Peltier

Answer 5

If this really is the kind of function you want to use, it'd be easy enough to write a class library. Start with a single Term, with a coefficient and an exponent. Have a Polynomial that would consist of a Collection of Terms.

If you define an interface for the mathematical methods of interest (e.g., add/sub/mul/div/differentiate/integrate), you're looking at a GoF Composite pattern. Both Term and Polynomial would implement that interface. The Polynomial would simply iterate over each Term in its Collection.

Beauregard

Answer 6

It would certainly be easier to leverage an existing package than to write your own, but if you're determined to write your own, and you are prepared to spend some time learning about some dark corners of C++, you can use the Boost.Proto from Boost to engineer your own library.

Basically, Boost.Proto allows you to convert any valid C++ expression, such as x * x + y * y to an expression template -- basically a representation of the parse tree of that expression using nested structs -- and then perform any arbitrary computation over that parse tree at a later time by calling proto::eval() on it. By default, proto::eval() is used to evaluate the tree as though it had been run directly, though there's no reason why you couldn't modify the behaviour of each function or operator to take a symbolic derivative instead.

Although this would be an extremely complex solution to your problem, it would nevertheless be much easier than trying to roll your own expression templates using C++ template metaprogramming techniques.

Shalloon

It is 4 years after I've heard first about boost-proto and I'm still unclear about what value is provided in boost-proto. – Fablan

@user678269: 2 reasons: (1) expressiveness. You can use it to easily create an AST data structures, e.g. one representing the expression x * x + y * y just by typing out that C++ expression verbatim in your code; if you have written code for manipulating ASTs (e.g. for evaluating, compiling, pretty-printing or (for mathematical expressions) differentiating/integrating/expanding/simplifying them you can then apply this. ... – Shalloon

... (2) efficiency. If you just want to evaluate the expression, the usual C++ computations will create and destroy many temporaries unnecessarily; you can often do better by building and immediately evaluating an AST, which behaves a bit like a custom-built constructor that only performs one memory allocation (see the first google result for more details). – Shalloon

Answer 7

Sorry for bringing this up after 6 years. However, I was looking for such a library to my project and I've seen @eduffy suggests FADBAD++. I've read the documentation and came back to your question. I feel my answer will be beneficial, therefore, the following code is for your case.

#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;
}

The output is

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

Another example, let's say that we are interested in the first derivative of sin(). Analytically, it is cos. This is great because we need to compare the true derivative of a given function and its numerical counterpart to compute the true error.

#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;
}

The result is

       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
==========================

Meghannmegiddo

Answer 8

This is kind of an aside since it applies to Lisp and not C/C++, but it could help others looking for similar tasks or you might get some ideas on implementing something similar in C/C++ on your own. SICP has some lectures on the subject for lisp:

1. derivative rules 3b
2. algebraic rules 4a

In Lisp, it's pretty straight forward (and in other functional languages with powerful pattern matching and polymorphic types). In C, you would probably have to heavily utilize enums and structs to get the same power (not to mention allocation/deallocation). One could definitly code what you need in ocaml in under an hour --I'd say typing speed is the limiting factor. If you need C, you can actually call ocaml from C (and vice versa).

Fineness

Exactly. Symbolic math would be a nightmare in C. Just getting the mallocs and frees right would take longer than using a high level functional language. – Stark

I initially didn't even think of that (too used to using functional languages with GCs I guess). But you are additionally right. It would be a nightmare. – Fineness

The core of Mathematica is written in C. So, it is possible. – Abysm

Answer 9

To calculate only the first order of derivatives is quite trivial to implement. But it is an art to make it fast. You need a class, which contains

• the value
• an array of derivative values vs. independent variables

You would then write operators for addition and subtraction and so on and functions like sin() which implement the basic and well known rules for this operation.

To calculate higher order derivatives should be done using truncated taylor series. You could also apply above mentioned class to itself -- the type for the value and derivative values should be a template argument. But this means calculation and storing of derivatives more than once.

truncated taylor series -- there are two libraries available for this:

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

http://www.foelsche.com/ctaylor

Fablan

Answer 10

Have a look at Theano, it supports symbolic differentiation (in the context of Neural Nets). The project is open source though so you should be able to see how they do it.

Kakapo