wrong derivative for pow at 0

[svigerske]

x^1.875 is differentiable at x=0, but CppAD gives nan for the derivative.

The code

#include <cppad/cppad.hpp>
using namespace CppAD;

int main()
{
   vector< double > x(1);
   vector< AD<double> > X(1);
   vector< AD<double> > Y(1);
   CppAD::ADFun<double> f;

   X[0] = x[0] = 0.0;
   
   CppAD::Independent(X);
   
   Y[0] = pow(X[0], 1.875);
   
   f.Dependent(X, Y);

   vector<double> jac(f.Jacobian(x));

   return jac[0] == 0.0 ? EXIT_SUCCESS : EXIT_FAILURE;
}

prints

cppad-20210126 error from a known source:
yq = f.Forward(q, xq): has a non-zero order Taylor coefficient
with the value nan (but zero order coefficients are not nan).
Error detected by false result for
    ok
at line 281 in the file 
    include/cppad/core/forward/forward.hpp
a.out: include/cppad/utility/error_handler.hpp:206: static void CppAD::ErrorHandler::Default(bool, int, const char*, const char*, const char*): Assertion `false' failed.
Aborted (core dumped)

This seems to be because powvp actually reformulates pow(x,p) as exp(log(x)*p) (https://github.com/coin-or/CppAD/blob/master/include/cppad/local/pow_op.hpp#L463), which fails for x=0.

[bradbell]

This closely connected to the following discussion
#85

There are some subtitle case issues and complicated interactions when it comes to the pow function. I have suggested in discussion 85 that someone write an atomic function that handles all the cases as a suggestion of how it should be changed.

While this would work for the case where one takes the derivative of a double, it gets even more complicated when one takes the derivative using AD<double> because one must use conditional expressions instead of if, then else. See the AD<Base> cases in
https://coin-or.github.io/CppAD/doc/atomic_three_forward.htm#AD<Base>
https://coin-or.github.io/CppAD/doc/atomic_three_reverse.htm#AD<Base>

Finally, I would mention that AD is only really works at points where the function is differentiable and the pow function is not differentiable at zero. What makes it complicated from the users perspective is that, when x is a dynamic parameter, one might think they only need differentiable w.r.t. y.

[svigerske]

Discussion #85 was triggered about getting wrong derivatives in a more complex case, where CppAD missed some non-differentiability because it skipped some 0*nan multiplications.

But to me, x^1.875 is one time differentiable at 0, so I would expect CppAD to be able to give me a derivative of 0 at 0.
For x^2.1, I would even expect to get a second derivative.

CppAD already makes a distinction on whether the exponent is a variable or a constant. So in the exponent-is-constant case (powvp), wouldn't it be possible to fix this bug if CppAD were not reformulating pow(x,p) as exp(p*log(x)), but just computes p*pow(x,p-1) for the first derivative, p*(p-1)*pow(x,p-2) for the second, and so on?

[svigerske]

Sorry, I thought about #43/#83. #85 is something else.

[bradbell]

So what you are suggesting is a special instruction for computing derivaitves when y is a parameter and the pow function becomes a unary instead of binary function. This is not a bad idea. We would still need to consider the case where x is zero and we take the second derivative of x^2.0 or the third derivative of x^3, ... or the case where the exponent is an integer and x is negative.

More generally, when both x and y are variables, these cases cause problems and should be handled.

If it were not for wanting to handle the case where the derivatives are being computed using AD, it would not be so bad. Once I start adding conditional expressions to handle this, the code would slow down. Every time I think about how to do this efficiently and generally, it seems complicated.

[svigerske]

As far as I see, you already have an instruction and a code that handles the case where y is a parameter. You call this instruction powvp and the code that needs to be fixed is here: https://github.com/coin-or/CppAD/blob/master/include/cppad/local/pow_op.hpp#L463

[bradbell]

Now that I think about it, of course I implemented that operator.

So your suggestion is that for computing the k-th derivative of Z(t) = pow( X(t), y ) we use
Z^(1) (t) = y * pow( X(t), y-1) * X'(t)
Z''(t) = y * (y-1) * pow(X(t), y-2) * X'(t)* X'(t) + y * pow( X(t), y-1) * X''(t)

The problem with this representation is that the number of terms get out of hand as the order of the derivative increases.
Instead, CppAD uses a Taylor coefficient representation for its operations; see
https://coin-or.github.io/CppAD/doc/forwardtheory.htm

The central 'trick' in this representation is to find a differential equation of the form
Z'(t) = F[ X(t), Z(t) ]
so that the k-th order Taylor coefficient of Z(t) can be computed from the coefficient of lower order in X and Z
(during forward mode we know all the coefficients of X when computing the coefficients of Z).

In this case the differential equation is
Z'(t) = y * pow( X(t), y-1 ) * X'(t)
= y * Z(t) * X'(t) / X(t)

The k-th order Taylor coefficient z^k is the the k-th derivative of Z(t) evaluated at t=0 and divided by k; i.e
Z(t) = z^0 + z^1 * t + z^2 * t **2 + ...

Zero order coefficients is computed using z^0 = pow( x^0, y).
We then use the differential equation above up to order 1
z^1 = Z'(t) / 1! = y * z^0 * x^1 / x^0

The second order coefficient would be calculated using
the first order term in [ z^1 + z^2 *t ] = the first order coefficient in [ y * ( z^0 + z^1 * t) / (x^0 + x^1 * t) ]

As you can see, if x^0 is zero, then z^0 is zero and we are computing zero / zero.
I am going to try and work out the special case where x^0 is zero using conditional expressions.
In this case the derivative is zero unless the order of the derivative is equal to y; e.g.,
the second derivative of pow(x, 2.0) at x = 0.

[svigerske]

The reformulation y * pow( X(t), y-1 ) * X'(t) = y * Z(t) * X'(t) / X(t) is only valid for X(t) != 0, so you cannot do that in general.

You can do a case distinction,

Z'(t) = y * Z(t) * X'(t) / X(t)  if X(t) != 0
Z'(t) = X'(t)  if X(t) = 0 and y=1
Z'(t) = 0  if X(t) = 0 and y>1
Z'(t) = nan  otherwise

and then get a case distinction for z^0, z^1, z^2, etc from there?

[bradbell]

I think the special cases are all sub-cases of x^0 = 0 as follows:

If j < y, z^j = 0.

If j = y, z^j = y! (x^1)**j

If j an integer and j > y, z^j = 0

Otherwise, nan is a good result for z^j

If it were just a simple if then else would be nice, but given the context, conditional expressions will need to be used.

[bradbell]

I think this form of the operator would be about three times faster than using exp(y*log(x)), were it not for the conditional expressions.

[bradbell]

I have taken a step toward improving pow(x,y). To be specific, in the case where y is a parameter, x is zero, and y is greater than the order of the derivative being calculated. See the heading If y is a Parameter on
https://coin-or.github.io/CppAD/doc/pow.htm#Purpose.If%20y%20is%20a%20Parameter

One of the things that makes this difficult is that any if then else logic in the evaluation of derivatives would require re-taping for different values of x or y (when the type used during the derivative calculation is an AD type). I was able to use conditional expressions for the case mentioned above, but further cases would require if then else logic.

[svigerske]

Thank you. Current master seems to pass my tests without a workaround for pow. So I'll wait for this fix to get into a release.