Incorrect answer for binary division involving two `TaylorScalar`s

[facusapienza21]

Hi!

I encounter this undesirable behaviour as I was doing some TaylorDiff:

derivative(x -> 1 + 1/x, 1.0, Val(1)) 
# returns -1.0, correct answer

derivative(x -> (x+1)/x, 1.0, Val(1))
# returns -2.0, should return -1.0

# If we change the order of the operations, this last example works
derivative(x -> 1/x*(x+1), 1.0, Val(1))
# returns -1.0

It seems that TaylorDiff.jl cannot handle very well the division of TaylorScalars, which is even more evident by looking at the following example:

derivative(x -> x/x, 1.0, Val(1))
# returns -1.0, should return 0.0

This behaviour is the same as when operating directly with TaylorScalars:

x_ty = TaylorScalar((1.0, 1.0))
# TaylorScalar{Float64, 1}(1.0, (1.0,))

x_ty / x_ty
# TaylorScalar{Float64, 1}(1.0, (-1.0,))

I tried to find the origin of this bug but I couldn't find it... I checked on the definition of the ration between TaylorScalars and it seems right, although it is difficult to debug with the @immutable macro in the definition.

Happy to help debugging this! Definitively something that needs to be solved soon!

[tansongchen]

Thanks for mentioning this. Earlier this month, I tried refactoring the code to use a macro to build the instructions to manipulate Taylor polynomials, in order to reduce hard-to-read generated functions everywhere. However, this line https://github.com/JuliaDiff/TaylorDiff.jl/blob/main/src/utils.jl#L94 accidentally introduced the error you mentioned.

I will fix this soon!

[facusapienza21]

Thank you @tansongchen ! Yeah, it would be great to fix it soon, since I don't know how much this is currently affecting the higher-order derivatives I am computing ;)

[tansongchen]

Solved by 29b1312

[facusapienza21]

Thank you @tansongchen ! Can you give me some light on what exactly these part of the code does? Or if you have some external reference that may help me understand this, I will appreciate it!

[tansongchen]

Sure! Here are the breakdowns:

1. Each function marked with @immutable converts a regular function to a @generated function, which makes it possible to utilize the size (P) info at compile time. This is needed because in Julia there is no simple way to generate a fixed-size NTuple that is recursively defined (element $j$ depends on $1,\cdots,j-1$), so I do that by knowing size at compile time and handle it correspondingly;
2. In the @generated function, expressions are transformed with function process. During these transformations, $P + 1$ variables are created, calculated and, at the end, combined into a NTuple and converted to a TaylorScalar finally. Specifically, v[i] are "magic names" that mark variables to be generated. As a mental model, you can imagine that v is a mutable fixed size array (like std::array in C++) initialized to be all zero(T).

The incorrectness happens because I want to replace the first += or -= assignment to these variables to be = (or else they are undefined), however I forgot to handle the case where the first assignment to a variable is =. The diff in the commit above is only 1 line: take all assignments into account and remember the "legal" names to be added into.

These are hacky and might not be very robust, but did make the functions like *, /, exp easy to read and correspond better to mathematical formulas. If you find that it doesn't work in some scenarios, please let me know...

[tansongchen]

In addition, could you let me know what specific problems you are currently solving that require higher-order autodiff?

[facusapienza21]

Thank you for the clarification @tansongchen !

I recently started exploring higher-order AD and nested AD features in Julia. (See for example my recent project SphereUDE.jl, where the goal is to use UDEs to fit directional data with regularization. Sorry the docs are not updated, I will do soon!) I am in a very exploratory phase of applications of higher-order derivatives for other problems, for example to evaluate the smoothness of a given numerical solution.

Hopefully I will keep using the library and we can collaborate!

[tansongchen]

Thank you for introducing your project, it looks really cool! Sure, let me know if there's anything we can collaborate on.