The main function : `rdiff()`¶

The differentiation function is called rdiff() and is called with the following parameters:

rdiff( ex::Expr; outsym::Symbol; order::Int, init... )

Arguments¶

ex:	is a Julia Expression containing the code to derive
outsym:	(default = nothing) is the symbol of the variable within `ex` containing the expression output (the result whose derivatives are needed). This variable must evaluate to a `Real`. If not specified, `outsym` defaults to `nothing` which signals to `rdiff` that the last statement is the result of interest for derivation.
order:	(default = 1) is an integer indicating the derivation order (1 for 1st order, etc.). Order 0 is allowed and will produce an expression that is a processed version of `ex` with some variables names rewritten and possibly some optimizations.
init:	(multiple keyword arguments) is one or several symbol / DataType pairs used to indicate for which variable a derivative is needed and how they should be interpreted. By default the generated expression will yield the derivative for each variable given unless the variable is listed in the `ignore` argument.
evalmod:	(default=Main) module where the expression is meant to be evaluated. External variables and functions should be evaluable in this module.
debug:	(default=false) indicates if `rdiff` should dump the graph of the generating expression, instead of returning the expression itself.
allorders:	(default=true) indicates whether to generate the code for all orders up to `order` (true) or only the last order.
ignore:	(default=[]) do not differentiate against the listed variables, useful if you are not interested in having the derivative of one of several variables in `init`.

Output¶

An expression which, when evaluated, will return a tuple containing the expression value and the derivative at first, second , etc.. order.

Usage¶

rdiff takes an expression consisting of a subset of Julia statements ( assignments, getindex, setindex!, for loops, function calls ) and transforms it into a new expression whose evaluation will provide the derivatives at all orders between 0 and the order specified (unless allorders is false).

The generated expression will attempt to remove all unneeded calculations (e.g. x + 0) and factorize repeated function calls as much as possible.

All the variables appearing in the init argument are considered as the expression’s arguments and a derivative is calculated for it (and cross derivatives if order is >= 2), unless they are listed in the ``ignore`` argument. The other variables, if not defined by the expression, are expected to be top level variables in evalmod. If they are not defined there an error will be thrown.

For orders >= 2 only a single variable, of type Real or Vector, is allowed. For orders 0 and 1 variables can be of type Real, Vector or Matrix and can be in an unlimited number:

julia> rdiff( :(x^3) , x=Float64)  # first order
:(begin
    (x^3,3 * x^2.0)
    end)

julia> rdiff( :(x^3) , order=3, x=Float64)  # orders up to 3
:(begin
        (x^3,3 * x^2.0,2.0 * (x * 3),6.0)
    end)

rdiff runs several simplification heuristics on the generated code to remove neutral statements and factorize repeated calculations. For instance calculating the derivatives of sin(x) for large orders will reduce to the calculations of sin(x) and cos(x):

julia> rdiff( :(sin(x)) , order=10, x=Float64)  # derivatives up to order 10
:(begin
        _tmp1 = sin(x)
        _tmp2 = cos(x)
        _tmp3 = -_tmp1
        _tmp4 = -_tmp2
        _tmp5 = -_tmp3
        (_tmp1,_tmp2,_tmp3,_tmp4,_tmp5,_tmp2,_tmp3,_tmp4,_tmp5,_tmp2,_tmp3)
    end)

The expression produced can easily be turned into a function with the @eval macro:

julia> res = rdiff( :(sin(x)) , order=10, x=Float64)
julia> @eval foo(x) = $res
julia> foo(2.)
(0.9092974268256817,-0.4161468365471424,-0.9092974268256817,0.4161468365471424,0.9092974268256817,-0.4161468365471424,-0.9092974268256817,0.4161468365471424,0.9092974268256817,-0.4161468365471424,-0.9092974268256817)

When a second derivative expression is needed, only a single derivation variable is allowed. If you are dealing with a function of several (scalar) variables you will have you aggregate them into a vector:

julia> ex = :( (1 - x[1])^2 + 100(x[2] - x[1]^2)^2 )  # the rosenbrock function
julia> res = rdiff(ex, x=Vector{Float64}, order=2)
:(begin
    _tmp1 = 1
    _tmp2 = 2
    _tmp3 = 100.0
    _tmp4 = _tmp1 - x[_tmp1]
    _tmp5 = length(x)
    _tmp6 = zeros(size(x))
    _tmp7 = x[_tmp2] - x[_tmp1] ^ _tmp2
    _tmp8 = zeros((_tmp5,_tmp5))
    _tmp9 = _tmp2 * (_tmp7 * _tmp3)
    _tmp10 = -_tmp9
    _tmp6[_tmp1] = _tmp6[_tmp1] + (_tmp2 * (x[_tmp1] * _tmp10) + -(_tmp2 * _tmp4))
    _tmp6[_tmp2] = _tmp6[_tmp2] + _tmp9
    for _idx1 = _tmp1:_tmp5
        _tmp11 = zeros(size(_tmp6))
        _tmp12 = zeros(size(x))
        _tmp11[_idx1] = _tmp11[_idx1] + 1.0
        _tmp13 = _tmp11[_tmp2]
        _tmp11[_tmp2] = 0.0
        _tmp11[_tmp2] = _tmp11[_tmp2] + _tmp13
        _tmp14 = _tmp2 * _tmp11[_tmp1]
        _tmp15 = _tmp3 * (_tmp2 * (_tmp13 + -(x[_tmp1] * _tmp14)))
        _tmp12[_tmp1] = _tmp12[_tmp1] + ((_tmp10 * _tmp14 + _tmp2 * (x[_tmp1] * -_tmp15)) + -(_tmp2 * -(_tmp11[_tmp1])))
        _tmp12[_tmp2] = _tmp12[_tmp2] + _tmp15
        _tmp8[(_idx1 - 1) * _tmp5 + 1:_idx1 * _tmp5] = _tmp12
    end
    (_tmp4 ^ _tmp2 + 100 * _tmp7 ^ _tmp2,_tmp6,_tmp8)
    end)
julia> @eval foo(x) = $res
julia> foo([0.5, 2.])
    (306.5,[-351.0,350.0],
    2x2 Array{Float64,2}:
     -498.0  -200.0
     -200.0   200.0)

foo(x) returns a tuple containing respectively the value of the expression at x, the gradient (a 2-vector) and the hessian (a 2x2 matrix)

Limitations¶

The canonical implementation of for loops derivation in reverse accumulation requires the caching of the complete state of each iteration which makes the generated code complex and memory intensive. The current algorithm uses a simpler approach that limits the kind of loops that can be correctly derived : in short, loops should not have any kind of recursivity in them (the calculations of each iteration should not depend on the calculations of previous iterations):
```
# will work
for i in 1:n
    a = f(x[i])
    b = a + g(y[i])
    c[i] = b
end

# will (probably) not work
for i in 1:n
    c[i] = f( c[i-1] )
end
```

However simple accumulations are an instance of recursive calculations that will work:

# will work
for i in 1:n
    a += b[i]    # new a value depends on previous a
end

for loops are limited to a single index. If you have a for i,j in 1:10, 1:10 in your expression you will have to translate it to nested loops as a workaround
All variables should be type-stable (not change from a scalar to a vector for example).
Only a limited set of Julia semantics are supported at this stage. Some frequently used statements such as comprehensions, if else, while loops cannot be used in the expression.
Mutating functions cannot be used (with the exception of setindex! and setfield!).

The main function : rdiff()¶

Arguments¶

Output¶

Usage¶

Limitations¶

The main function : `rdiff()`¶