A Program whose Automatic Derivative is Wrong on All Inputs

Automatic differentiation (AD) is a process by which one program computing on numbers is transformed into a new program computing the derivative (in the sense of Calculus) of the original program. Over the last decade, AD has become central to modern machine learning systems: “backpropagation” through neural networks is simply a form¹ of AD. I have been personally interested in automatic differentiation for years now as a way to simplify the writing of physical simulation and optimization programs.

What is Automatic Differentiation?

As a simple example consider a program for computing the magnitude of a vector \((x,y)\in\mathbb{R}^2\):

def norm(x, y):
  n2 = x*x + y*y
  return sqrt(n2)

The automatic derivative of this program is given by

def dnorm(x, y, dx, dy):
  n2 = x*x + y*y
  dn2 = 2*x*dx + 2*y*dy
  return 0.5 * 1.0/sqrt(n2) * dn2

Why is this the result? If norm represents a function \(f : \mathbb{R}^2 \to \mathbb{R}\), then its automatic derivative represents a function \(Df : \mathbb{R}^2 \times \mathbb{R}^2 \to \mathbb{R}\), known as the “total derivative.” It is the closest linear approximation to \(f\) “at the point” \((x,y)\), and this is why it has two inputs. That is, \(Df(x,y) : \mathbb{R}^2 \to \mathbb{R}\) (partially evaluated at \((x,y)\)) is a linear function approximating \(f\) as \(f(x+dx,y+dy) \approx f(x,y) + Df(x,y,dx,dy)\). Alternatively, we can specify the total derivative explicitly in terms of coordinates:

Secondly, we might wonder “how was the automatic derivative computed?” To make a long story short, it works by transforming the program text. For each assignment \(y = e\), we add a new assignment \(dy = e'\), where \(e'\) is the derivative of the expression \(e\), taken according to the usual rules of differential calculus. However, there is one small caveat²: the derivative of a lone variable \(x\) is recorded as \(dx\), not as \(1\). It turns out that this is just another expression of the chain rule, but we needn’t digress into that for the sake of this blog post.

One nice thing about automatic differentiation is that it still works on programs that have loops, recursion, and branches. For instance, consider the “relu” function from machine learning and its derivative:

def relu(x):
  if x > 0 then return x else return 0

def drelu(x,dx):
  if x > 0 then return dx else return 0

Loops, recursion and branching are all instances of program control flow. The rule for AD is that we differentiate each branch and keep the expressions in conditionals unmodified. This makes some good sense. Above, so long as \(x \neq 0\), the derivative of relu is unambiguously accurate. But what about at \(x=0\)? At that point, the derivative is undefined. So, we’ve done the best³ we can hope for.

The Program

Computers use floating point numbers to approximate real numbers. We don’t need to know much about floating point for our purposes here other than two facts. First, every floating point number \(q\) is a “dyadic” number, which simply means that it can be written as \(q = n\cdot 2^m\), where \(n\) and \(m\) are both integers. Second, because the number of bits used to encode \(n\) and \(m\) are limited, there is a smallest possible “denominator” for any floating point number: approximately, \(-1024 \leq m < 1024\). In order to avoid mistakes, I will loosen this bound to \(-2048 < m < 2048\).

With this in mind, let us now consider the following program:

def f(x):
  if x == 0.0: return 0.0
  else if x < 0.0: return -f(-x)
  else if x >= 2.0: return 2.0 * f(0.5*x)
  else if x < 1.0: return 0.5 * f(2.0*x)
  else if x == 1.0: return 1.0
  else: # 1.0 < x < 2.0
    return 1.0 + f(x - 1.0)

What does this program do? In effect, it deconstructs a floating point number into its binary representation and then reconstructs it. That is, I claim that the output of the program \(f(x)\) is simply \(x\). To see why this is true, we need only perform a case analysis. I have provided the full case analysis below if you are interested, but you can safely skip it. The rest of the blog post does not rely on these details.

• First, either \(x\) is \(0\) or not. If \(x\) is \(0\), then \(f(0) = 0\) as claimed. Otherwise, none of the other cases will make a recursive call with argument \(0\), so we won’t need to reconsider this case.
• Next, if \(x < 0\), we recursively call \(f\) with \(-x\), which must now be positive. If the recursive call computes \(f(-x) = -x\) as claimed, then we just return \(-f(-x) = x\), as expected. Furthermore, so long as \(x\) is positive, all further cases will only pass positive arguments to recursive calls. From here on out, we can simply assume that \(x > 0\).
• Next, if \(x \geq 2\) we will make recursive calls to \(f(x / 2)\), until the argument lies in the range \(0 < x < 2\). For each of these calls we can again confirm that so long as the recursive calls guarantee \(f(x) = x\), then \(2\cdot f(x/2) = x\). After this point, \(x\) must remain in the range \(0 < x < 2\).
• To understand the last three cases, let us first observe that any number \(0 < x < 2\) can either be written in the form \(1.\cdots\) or \(0.\cdots\), both being binary decimal expansions of the number. And so long as \(x\) is a floating point number, it can only have \(1\)s in the first \(k\) digits of this binary expansion. (where \(k is bounded by say\)4096$$)
• \(x=1\) exactly when this expansion is \(1.000\cdots\). At this point we terminate and return \(f(1) = 1\) as expected.
• When \(0 < x < 1\), it has a binary expansion \(0.\cdots\). The recursive call to \(f(2x)\) is equivalent to left-shifting this expansion, which will eventually move a \(1\) digit into this high-order bit position. As with other cases, \(f(2x)/2 = 2x/2 = x\).
• Finally, if \(1 < x < 2\), then it has a binary expansion \(1.\cdots\). Taking \(x-1\) consumes this bit and reduces to the left-shifting case for \(0 < x < 1\). As with other cases, \(1 + f(x-1) = 1 + x - 1 = x\).

To recap, we have analyzed our program to show that (1) it computes the identity function and (2) terminates for arbitrary floating point number inputs⁴.

Now, let us take the automatic derivative of this program. We simply follow the rules we already specified for how to take derivatives of expressions and through branches.

def Df(x,dx):
  if x == 0.0: return 0.0
  else if x < 0.0: return -Df(-x, -dx)
  else if x >= 2.0: return 2.0 * Df(0.5*x, 0.5*dx)
  else if x < 1.0: return 0.5 * Df(2.0*x, 2.0*dx)
  else if x == 1.0: return 0.0
  else: # 1.0 < x < 2.0
    return Df(x - 1, dx)

A quick examination of this function will reveal that every base case of the recursion returns \(0\), and that if the result of a recursive call is \(0\), then \(0\) will again be returned in every other case. Therefore, we can conclude that the above program computes \(Df(x,dx) = 0\) for all floating point numbers \(x\).

However, \(f(x)\) computes \(x\), so its total derivative ought to compute \(dx\). The \(Df(x,dx)\) program given above computes \(0\) for every floating point number. That is to say, \(Df\) is “wrong on all inputs.”

Let us now reflect on our sins.

Programs have branches, and we had no good way to define the derivative at a branching point. Therefore, we waved our hands and declared aloud before all who would hear, “I have done the best I possibly could! The derivative at this point is undefined!” But the derivative at a branch point was not always undefined. In our pathological example, the function happened to be continuously differentiable at every branch point, and every input was a branch point! Did we really do our best?

If only a want of effort was our worst sin, we might be forgiven. But lo and behold how deeply we have profaned against common sense! Suppose we replace the line

  else if x == 1.0: return 1.0

in our original, undifferentiated program f with the following alternate, and semantically equivalent line

  else if x == 1.0: return x

Now the derivative of this line would be

  else if x == 1.0: return dx

which gives the correct result. Then the result of automatic differentiation would be correct everywhere other than \(x=0\) (which could be similarly fixed).

Let \(g\) be this “fixed” program. The two programs \(f\) and \(g\) compute/represent the same function, but their automatic derivatives are different. Sacre bleu! We cannot even say that automatic differentiation preserves the meaning of our programs!

Theory: The Impossibility of a Systematic Fix

Automatic differentiation seemed so beautiful and alluring. What went wrong?

Permit me a brief theoretical digression. Let \(P\) be a program, meaning the literal syntactic text of the program, not the function it computes. Then, we say that \(\mathcal{E}(P)\) is the function that \(P\) denotes. In programming language theory, \(\mathcal{E}\) represents the semantics of our programming language. It tells us what programs mean. This distinction between syntax and semantics is foundational to the study of programming languages.

With these ideas in mind, we can now provide a formal specification for what we want automatic differentiation to do. That is, we seek a program transformation \(AD\) that transforms a program \(P\) (i.e. its syntactic form) into a new program \(AD(P)\) (also literal program text). In particular, we want the resulting program to mean the derivative of the original program, in the sense that

holds for all \(x,dx\) inputs such that the derivative of \(\mathcal{E}(P)\) is well-defined for \(x\). Surely, this is what automatic differentiation should do, but as we just saw it is not what automatic differentiation actually does. In truth, automatic differentiation is a syntactic manipulation with much more limited semantic guarantees. The result depends not only on what result a program computes, but on the internal details of how it chooses to compute that result. AD is fundamentally and inextricably syntactic.

The natural instinct of many readers at this point is to rebel. Perhaps we can create a better version of automatic differentiation. If we could simply analyze the function at each branch point, then we could atone and build the one, true, correct way of differentiating computer programs.

Unfortunately, this ideal is impossible—provably so. In a classic paper⁵, John Myhill gives a construction for a computable function \(f\) for which \(\frac{df}{dx}\) exists but is uncomputable. To be precise, the notion of computable here is that there is a program \(P_f\) which can produce arbitrarily close approximations to \(f\) so long as you give it enough time to run. From this citation we can deduce our claim via contradiction. If there was a way to compute automatic differentiation of programs that satisfied our correctness criterion, then we could apply it to \(P_f\) to get computable approximations of the derivative. Since we provably can’t do this, the program transformation function \(AD\) (whatever it is, and however it works) must at least sometimes fail to “compute the derivative.”

Silver Linings

Despite these shortcomings, automatic differentiation is incredibly useful and important. One might be inclined to say that code riddled with branches doesn’t occur in practice. However, neural networks full of relu activation functions are the backbone of modern machine learning. It might be better to say that many applications are inherently tolerant to the shortcomings of AD, or that they can be made tolerant.

When there exists a correct answer (but AD does not compute it) a clever programmer can rewrite the original program to avoid the issue. Even in cases where there doesn’t exist a derivative to correctly compute, programmers can often patch around the issue with application specific fixes. For instance, the vector-norm function at the start of this post is pervasive in computer graphics. In order to avoid division by zero exceptions, one can replace \(\sqrt{x^2 + y^2} < \epsilon\) in denominators with \(\max(\epsilon,\sqrt{x^2 + y^2})\). As \((x,y)\) tends towards \(0\), this hack will ensure that dnorm also decays continuously to \(0\). This is the wrong answer, but is much more likely to produce unobjectionable rather than catastrophic results in the context of a larger application. For instance, when performing shading calculations during 3d rendering, this hack might produce some abnormally dark pixels for very bad geometry. Such pixels are rare, and their influence is further reduced by anti-aliasing.

For those of us who work on programming languages (especially domain-specific languages) this is also a fantastic opportunity. Even though there cannot be a universally correct automatic differentiation system, there are many ways to improve the practical reliability and robustness of existing AD systems. This is especially true when derivatives are combined with statistical estimation and integration (e.g. via Monte-Carlo).

Regardless, I hope this example program helps make the (necessary) theoretical deficiencies of automatic differentiation more tangible for you.

Footnotes