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...
Program:
y = g(x) = x^3
z = f(y) = cuberoot(y)
Want: dz/dxReverse-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 = NaNIt 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.
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.
> a complicated representation of) y=x.
The problem started as that. So no surprise.
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.
#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 >>> 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.