Thank you, this is interesting. I just tried. You need L'Hospital's Rule.

I'd realized this was an issue with forward-mode autodiff, but hadn't thought about reverse-mode to realize it comes up there too.

In general, NaN is a placeholder for "all the derivatives you didn't know you needed in advance". I think the Dual Numbers shed light on this (see "Division": [https://en.wikipedia.org/wiki/Dual_number#Division]).

One could imagine lazily computing higher derivatives as needed for L'Hospital's Rule in either Forward or Reverse mode, however. For Forward Mode, I guess you'd replace your eager "Dual Number" pairs with lazy lists?

For reverse mode it seems you could do extra "passes" on the tape as needed... People who do this stuff for a living must know(?), but a quick Google doesn't turn up an answer for me...

Actually I'm wrong about L'Hospital's Rule in this problem; it solves nothing here. I would insert an "[EDIT ...]" into my original post to avoid misleading people, but it's too late to edit.

It doesn't matter how many times you differentiate x^(1/3); you can not evaluate the result at zero. Which is funny, because you'd think you could differentiate x^3 and use the inverse function theorem. But there's that caveat about it not working when the derivative is zero.

I am now revisiting high school calculus and doubting my own sanity. Will post updates if/when I get my head right.

Here's the (pedantically-worked) example with reverse-mode:

Program:

  y = g(x) = x^3
  z = f(y) = cuberoot(y)

Want: dz/dx

Reverse-mode autodiff:

  // Forward pass
  y = g(0) = 0^3         = 0      g'(0) = 3 x^2                           = 0
  z = f(0) = cuberoot(0) = 0      f'(0) = (1/3) / cuberoot(y^2) = (1/3)/0 = Inf

  // Reverse pass
  dz/dy = f'(0)                   = Inf
  dz/dx = dz/dy * g'(0) = Inf * 0 = NaN

Sorry, I was out.

Ok,

f(x) = cuberoot(x^3)

f(0 + eps) = cuberoot((0 + eps)^3) = cuberoot(eps^3) = eps

e.g. value is 0 and the slope is 1.

    >>> import autograd.numpy as np
    >>> from autograd import grad
    >>> def fn(x):
    ...   return np.power(np.power(x, 3), 1/3)
    ...
    >>> gradfn = grad(fn)
    >>> gradfn(0.0)
    /usr/local/lib/python3.6/dist-packages/autograd/numpy/numpy_vjps.py:59:RuntimeWarning: divide by zero encountered in double_scalars
    lambda ans, x, y : unbroadcast_f(x, lambda g: g * y * x ** anp.where(y, y - 1, 1.)),
    /usr/local/lib/python3.6/dist-packages/autograd/numpy/numpy_vjps.py:59: RuntimeWarning: invalid value encountered in double_scalars
     lambda ans, x, y : unbroadcast_f(x, lambda g: g * y * x ** anp.where(y, y - 1, 1.)),
    nan

… Damn.

    >>> import torch
    >>> x = torch.tensor(0.0, requires_grad=True)
    >>> y = ((x**3) ** (1/3))
    >>> y.backward()
    >>> x.grad
    tensor(nan)

… Damn.

eps^3 is undefined

What is eps^3? ;) You have to turn that into some representation on the computer before you can even pass it to cuberoot.

Theres is no representation of that value.In this case power is defined only in the range 1-2. You need to simplify the roots first.

except by the definition of dual numbers eps^2 is 0, as is eps^3. So the slope you computed is 0, which is rather wrong for (a complicated representation of) y=x.

Eps^3 is undefined. I don't see a problem with simplifying the exponents first.

> a complicated representation of) y=x.

The problem started as that. So no surprise.

> I don't see a problem with simplifying the exponents first.

Sure, we call that computer algebra. As soon as you start doing that, you aren't doing automatic differentiation, you are doing (at least in part) symbolic differentiation.

You could similarly say that a computer can't calculate x^1000/x^998 because you can't fit it in 32/64 bits.

And computers have a hard time with that. Take the following example code:

  #include <math.h>
  #include <stdio.h>
  #include <stdlib.h>
  double f(const double x) {
      return pow(x, 1000) / pow(x, 998);
  }
  int main(int argc, char* argv[]) {
      if(argc != 2) {
          fprintf(stderr, "Usage: %s x\n", argv[0]);
          exit(1);
      }
      const double x = atof(argv[1]); // demo code without error checking
      printf("x = %g, f(x) = %g\n", x, f(x));
      exit(0);
  }

compile with

  gcc -Wall -g -o test test.c -lm

and run it

  petschge@localhost:~$ ./test 0
    x = 0, f(x) = -nan
  petschge@localhost:~$ ./test 1
    x = 1, f(x) = 1
  petschge@localhost:~$ ./test 2
    x = 2, f(x) = 4

The fun thing is, when you compile with

  gcc -Wall -O3 -ffast-math -g -o test2 test.c -lm

you actually get

  petschge@localhost:~$ ./test2 0
    x = 0, f(x) = 0

That's exactly what i'm saying.