I am trying to calculate higher derivatives with ForwardDiff, but I get DomainError for the third derivative. I definitely need at least order 7, but up to 13 would be great. Does anybody have an idea, if a) this is supposed to work b) if yes what am I doing wrong. Thanks in advance…

The code is

g0(z) = [z[2]^2*z[1]^3]

g1(z) = ForwardDiff.jacobian(g0,z)

g2(z) = ForwardDiff.jacobian(g1,z)

g3(z) = ForwardDiff.jacobian(g2,z)

x=[ 2 2 ]

g1(x)

1×2 Array{Int64,2}:

48 32

g2(x)

2×2 Array{Int64,2}:

48 48

48 16

g3(x)

ERROR: DomainError:

Stacktrace:

[1] ^ at /home/rs1909/.julia/v0.6/ForwardDiff/src/dual.jl:395 [inlined] (repeats 3 times)

etc…

Use `x = [2., 2.]`

instead (giving you a vector of floats instead of integers). The error you get says:

```
Cannot raise an integer x to a negative power -n.
Make x a float by adding a zero decimal (e.g. 2.0^-n instead of 2^-n), or write 1/x^n, float(x)^-n, or (x//1)^-n.
```

Order 13 derivatives will create very complicated nested types, it is possible compilation time / inference time will become an issue but perhaps it will work well.

I suggest you try TaylorSeries.jl if you need high derivatives.

2 Likes

I wondered about that, but it seems to me that I would need to expand about a constant value. I am trying to implement the Faa Di Bruno formula (high order chain rule) and I need to dynamically change where the derivative is calculated.

Thank, this is great help. I should have realised the integer type problem. It works now.

Can you specify a bit more precisely exactly what you need? Here’s a starter example with TaylorSeries.jl (which is a registered package, so just `Pkg.add("TaylorSeries")`

to install):

```
julia> using TaylorSeries
julia> f(x, y) = cos(x*y) + sin(x - y^2);
julia> xx, yy = set_variables("x y", order=15);
julia> f(xx, yy)
1.0 + 1.0 x² - 1.0 y² - 0.5 x² y² - 0.16666666666666666 x⁶ + 0.5 x⁴ y² - 0.5 x² y⁴ + 0.16666666666666666 y⁶ + 0.041666666666666664 x⁴ y⁴ + 0.008333333333333333 x¹⁰ - 0.041666666666666664 x⁸ y² + 0.08333333333333333 x⁶ y⁴ - 0.08333333333333333 x⁴ y⁶ + 0.041666666666666664 x² y⁸ - 0.008333333333333333 y¹⁰ - 0.0013888888888888887 x⁶ y⁶ - 0.00019841269841269839 x¹⁴ + 0.0013888888888888885 x¹² y² - 0.004166666666666666 x¹⁰ y⁴ + 0.006944444444444443 x⁸ y⁶ - 0.006944444444444443 x⁶ y⁸ + 0.004166666666666666 x⁴ y¹⁰ - 0.0013888888888888885 x² y¹² + 0.00019841269841269839 y¹⁴ + 𝒪(‖x‖¹⁶)
```

There are various `derivative`

methods available to extract derivatives, or you can just access the coefficients of this multivariate polynomial.

If you need to expand e.g. around (-3, 4), you can do

```
f(3+xx, -4+yy)
```

Here, `xx`

and `yy`

now represent small perturbations around -3 and 4, respectively.

Thanks for the interest, the full problem is the following. I have a differential equation defined in a rather strange way. Let n > 2, and denote the total differentiation with respect to ‘t’ by capital D. My differential equation is

D^(n+1) u(t) = f(t,u(t)) + D^n { g(t,u(t)) },

where u is a vector and t is a scalar, f and g are also vector valued, the same dimension as u.

My biggest problem is that I need to calculate the right-hand side, which involves calculating the n-th derivative of g(t,u(t)). For this I need to use the generalised chain rule. By rewriting the equation into first order form, the values of D^k u(t), 0 <= k <= n become state variables, so in essence I need to calculate various partial derivatives of g(t,u) and combine them with D^k u(t).

The first issue is that I really need the full tensor form of the derivative, because they are used as multi-linear forms, multiplied by various D^k u(t) values. As far as I can tell Taylor expansion assumes that all arguments of the multi linear forms are the same (hence you get a polynomial)

The second issue is that I could not see an example how to parameterise the function by those shifted values, such that the derivative functions do not need to be recompiled every time a different shift is used (which is every call).

I can now see that ForwardDiff will work (in theory) even if it is very slow to compile. I was also thinking of just using computer algebra to do the derivation, but the function ‘g’ is not a nice algebraic expression.

As you can see this is highly non-standard and at the moment I think Julia is my best bet.

2 Likes

Thanks for the explanation. Is ForwardDiff set up to return the nth derivative as a tensor, though?