Automatic differentiation – Forward Mode

By: Random blog posts about machine learning in Julia

Re-posted from: https://rkube.github.io/julia/autodiff/2021/05/10/autodiff1.html

Automatic differentiation allows one to take derivatives of computer programs which
calculate numerical values. The name refers to a set of techniques that produce a
transformed computer program which in turn calculates various derivatives of these values.

To see why this is extremely useful, consider a traditional data fitting problem.
We have a model \(f_\theta: \mathbb{R}^m \rightarrow \mathbb{R}^n\), where
\(\theta\) are the model parameters and observations \(y_i \in \mathbb{R}^n\)
taken at \(x_i \in \mathbb{R}^{m}\). To fit the model on the data we can
use a loss function, for example the mean-squared error

\[\begin{align}
\mathcal{L} = \sum_{i} \left( f_\theta(x_i) – y_i \right)^2
\end{align}\]

and tune the parameters \(\theta\) in order to minimize \(\mathcal{L}\). This process
is also called learning or solving an inverse problem.

A common method to tune the parameters \(\theta\) is gradient descent. Given a learning
rate \(\eta\) one iterates

\[\theta \leftarrow \theta – \eta \frac{\partial \mathcal{L}}{\partial \theta}\]

until \(\mathcal{L}\) is approximately at a local minimum. Applying the chain rule, we
can expand the derivative as

\[\frac{\partial \mathcal{L}}{\partial \theta} = \sum_i 2 \left( f_\theta(x_i) – y_i \right)
\frac{\partial f_\theta(x_i)}{\partial \theta}.\]

To perform the optization procedure we need to calculate the derivatives described
by the last term in the equation above. This is where automatic differentiation
comes into play.

A Toy Example

To better understand what automatic differentiation does and how it works let us consider
the toy example

\[f(x_1, x_2) = x_1 x_2 \cos(x_2) \left[ \exp(x_1 x_2) – 1 \right] = y.\]

In the following we will investigate how to calculate the derivative of \(f\) with respect
to its inputs \(x_1\) and \(x_2\). In the context of the previous section, we can think of
\(x_1\) as a parameter \(\theta\) that we wish to tone. While the method we present here
is as simple as possible – real-world examples of tunable models often have thousands,
millions, or even billions of parameters – it serves well to introduce the fundamental way
in which automatic differentiation works. For the work below we are following [1]
closely.

This mathematical formula for \(f\) can be decomposed into a succession of elementary functions,
like \(+\), \(\sin\), \(\exp\), \(*\). That is, we can define intermediate variables like
\(v_1(a,b) = a \cdot b\), \(v_2(a) = \cos a\) etc, which are only defined once. Each of these
variables is also initialized by applying a simple expression to previous variables.
Applying this procedure to our toy example gives us the following sequence:

A simple program broken up into intermediate values
\(v_{-1} = x_1\)
\(v_{0} = x_2\)
\(v_{1} = v_{-1} \cdot v_{0}\)
\(v_{2} = \exp v_1\)
\(v_{3} = \cos v_0\)
\(v_{4} = v_1 \cdot v_3\)
\(v_{5} = v_2 – v_4\)
\(v_6 = v_4 v_5\)
\(y = v_6\)

We now define the so-called Tangential Derivative \(\dot{v_i} = \partial v_i / \partial x_1\).
For the tangential derivative, we are fixing the denominator to \(x_1\) so that it describes the sensitivity of the functions
output with respect to the given input. Let us calculate this derivative for the v’s now:

\[\begin{align}
\dot{v}_1 & = \frac{\partial v_1}{\partial x_1} = v_0 \frac{\partial v_{-1}}{\partial v_{-1}} = v_0 \\
\dot{v}_2 & = \frac{\partial v_2}{\partial v_1} \dot{v}_1 = v_0 \exp v_1\\
\dot{v}_3 & = \frac{\partial v_3}{\partial v_0} \dot{v}_0 = 0 \\
\dot{v}_4 & = \frac{\partial v_4}{\partial v_1} \dot{v}_1 + \frac{\partial v_4}{\partial v_3} \dot{v}_3 = v_3 v_0 \\
\dot{v}_5 & = \frac{\partial v_5}{\partial v_2} \dot{v}_2 + \frac{\partial v_5}{\partial v_4} \dot{v}_4 = \dot{v}_2 – \dot{v}_4 \\
\dot{v}_6 & = \frac{\partial v_6}{\partial v_5} \dot{v}_5 + \frac{\partial v_6}{\partial v_4} \dot{v}_4 = \dot{v}_5 v_4 + \dot{v}_4 v_5 \\
\end{align}\]

Note that we have to use \(\dot{v}_{-1} = 1\) and \(\dot{v}_{0} = 0\) in order to calculate \(v\dot{v}_1\). That
is, the seed determines which part of the sum is non-zero. From the expressions above we can alos see that we
will need to evaluate the derivatives starting at \(v_1\). That is, we begin at the start of the function,
calculate simple derivatives and push this calculation through. Therefore this is also called forward mode automatic differentiation.

In practice, forward mode AD is often implemented by operator overloading. An implementation
would need seed an initial derivative \(x_i = 1\) and then calculate the program and the derivatives simultaneously.
This approach works for our toy example as we see from the evaluation trace shown below:

Forward-mode AD evaluation trace  
\(v_{-1} = x_1 = 1.0\) \(\dot{v}_{-1} = 1.0\)
\(v_{0} = x_2 = 1.1\) \(\dot{v}_{0} = 0.0\)
\(v_1 = v_0 v_{-1} = (1.1) (1.0) = 1.1\) \(\dot{v}_{1} = v_{0} = (1.1)\)
\(v_{4} = v_1 v_3 = (1.1) (0.4536) = 0.4990\) \(\dot{v}_4 = v_3 \dot{v}_1 = (0.4536) * (1.1) = 0.499\)
\(v_{5} = v_2 – v_4 = (3.004) – (0.4990) = 2.5052\) \(\dot{v}_5 = \dot{v}_2 – \dot{v}_4 = 3.305 – 0.499 = 2.851\)
\(v_6 = v_4 v_5 = (0.4990) (2.7180) = 1.2500\) \(\dot{v}_6 = v_4 \dot{v}_5 + \dot{v}_4 v_5 = (0.499) (2.851) + (0.499) (2.505) = 2.650\)

It is illustrative to visualize how the values and derivatives are propagated in the computational graph.
Coputational graph for forward-mode autodiff

Here the intermediate values are shown in yellow and the gradient values are shown in turquoise. Note that
both, values and derivatives are propagated together from the left, the input, to the output on the right.
The implication of this is that forward-mode autodiff can calculate the sensitivities of all outputs y
with respect to one input x in a multidimensional setting.

Summary

Automatic differentiation is a set of tools that allows to calculate the derivatives of
computer programs that return numerical values. Here we looked at the so-called forward mode of
automatic differentiation. It works by de-composing a program into simple functions for which
we know the derivative. We can initialize a program and set the input variable for which we want
to know the output’s derivatives for. Then we can calculate the function value and the derivative
of the output with respect to the chosen value in one sweep.

This works well for situations where we have a small number of inputs and a large number of outputs.
But in data-fitting problems we often end up with the reverse situation: We have a large number
of outputs and few, or maybe only a single output. In these situations one may want to use
reverse-mode automatic differentiation.

References

[1]
A. Griewank, A. Walther – Evaluating Derivatives – SIAM(2008)

Why Julia – Meet Dispatch

By: Miguel Raz Guzmán Macedo

Re-posted from: https://miguelraz.github.io/blog/dispatch/index.html

Hello, World! 👋 I'm Dispatch

My name's Dispatch, your friendly neighbourhood walkie-talkie, and I'm here to talk to you all about Julia – a cool programming language and why I think they made the right calls when they built it.

Reading time: 20 minutes

Summary:

This post is about Julia and multiple dispatch as its foundation. Julia is not

  • the first language to implement multiple dispatch (see the lecture link below)

  • the only language to use multiple dispatch (Dylan does use it, also available as a lib in Python and others)

but it definitely feels like the first to stick the landing. My claim is this: Julia is unique in implementing multiple dispatch as a fundamental design feature and building a language around it. Don't take it too seriously; it's more an invitation to learn than anything else. If you like this topic, check out "The Unreasonable effectiveness of Multiple Dispatch" JuliaCon by Stefan Karpinski, and this article on "The Expression problem" which inspired that talk. If I was being uncharitable, this post is the watered-down, kawaii-ified version of that talk with some Lisp memorabilia sprinkled on top. Your mileage may vary.

The audience of this post are programmers who like Python/C++/80s MIT Lisp courses who have heard of Julia but don't get what new things it brings to the table. If you want to dive deeper into emergent Julia features, check out Frame's 2 part posts.

Sure, there's a fancy subtyping algorithm that Jeff Bezanson (Julia co-creator) pulled out of Turing-knows-where and that Swift re-adapted with some unholy term rewriting shenanigans, but that's not the "new contribution" I mean – careful design iteration around multiple dispatch with measured trade-offs is.

The rest of this post are code snippets in different programming languages and how they look like next to a Julian re-implementation.

To talk about all these things, we'll be hearing from our newest battery-powered friend, Dispatch.

Dispatch:

Hey peeps, that's me!

… interspersed with a comments of our not-yet-super-proficient-but-keeps-working-at-it-student, miguelito:

miguelito:

¡Hola Dispatch! Nice to meet you!

Off we go! 🚀

What is dispatch?

To begin our journey, we should define some terms. Per Merriam-Webster, to dispatch means:

to send off or away with promptness or speed

to dispatch a letter

dispatch an ambulance to the scene

Dispatch:

And if you have a dispatcher, which is what I do, then you call what you need at the right time. There's no need for extra fussing about! Sometimes you need to call a cab, sometimes its a delivery driver, sometimes its a helicopter. Just dispatch on what you need.

miguelito:

Ok, call what you nee– Aaaaaah, now I get why you're a walkie-talkie. Clever! I'll have to tell my friend Cool Bear about this one…

At any rate – dispatch sounds simple I guess… but all these programmers from the Object Oriented Land keep asking me about where I stand in inheritance vs composition, and my takeaway was that there's a lot of terms I thought knew what they meant, but I have to re-learn for all this programming stuff. I guess I'm just not that into Rectangle vs Square debates…

Dispatch:

Yup, there's a lot of fancy words. It can be hard to describe where you are standing to other people because the Kingdoms of Nouns monopolized the programming language map-making business decades ago.

Now it sucks to get directions anywhere.

But it's cool – just take a deep breath, and, like a well oiled StepMaster, we'll do this one step at a time.

On namespaces as well, those debates are just not going to be as big a problem in Julia Land – we're not too far from them in the theory, but the language doesn't actively prevent you from doing what you already know. In fact, if you know addition, you already have a natural notion of what multiple dispatch is about.

miguelito:

What? No way.

Dispatch:

Of course! What's more natural than addition? You just haven't used the name multiple dispatch for it, which is fine. Bear with me for some simple arithmetic: Can you describe the procedure you know for adding integers? Say you want to add 123 + 890. (This is a good moment to install and open julia from julialang.org if you haven't already done so.)

miguelito:

Sure, you line up the digits, add the columns right to left, carry the ones…

julia> 123 + 890
1013

Dispatch:

Yup. What about trying to add 1//2 + 1//3. Those "fractions" are what we call Rationals in Julia. What procedure do you follow then?

miguelito:

Oh right, the famous regla de tres – find the common denominator, cancel, add up apples to apples…

julia> 1//2 + 1//3
5//6

Dispatch:

Swell. And .25 + 5.2? You're still "adding" things right?

miguelito:

Yup, like adding natural numbers – line 'em up, carry the one.

julia> .25 + 5.5
5.75

Dispatch:

Excelente, ignore the decimals for now. Did you ever face matrices? Can you try to add 2x2 matrix of 1s and a 2x2 matrix of 3s?

miguelito:

Yup – I think we defined it as element to element. So you end up with a 2x2 matrix of 4s.

julia> [1 1; 1 1] + [3 3; 3 3]
2×2 Matrix{Int64}:
 4  4
 4  4

Dispatch:

Alright – thanks for following along. Now here's the tickler question – Who does the + belong to?

miguelito:

…qué.

Dispatch:

Exactly. If you know you're adding 2 Natural numbers, you just use, (or when coding, call) the right procedure. If you + 2 decimal numbers, you call the right thing. It's the addition you've always known! There's this notion of what + should do in different cases, and we wrap them all up under the same + umbrella – even though addition of matrices, decimals and natural numbers mean different procedures. You can check all the ones that come out of the box if you try this:

julia> methods(+)
# 195 methods for generic function "+":
[1] +(x::T, y::T) where T<:Union{Int128, Int16, Int32, Int64, Int8, UInt128, UInt16, UInt32, UInt64, UInt8} in Base at int.jl:87
[2] +(c::Union{UInt16, UInt32, UInt64, UInt8}, x::BigInt) in Base.GMP at gmp.jl:528
[3] +(c::Union{Int16, Int32, Int64, Int8}, x::BigInt) in Base.GMP at gmp.jl:534
...

miguelito:

Oh! But you stated that + "belongs" to someone, so what does + have to do with property?

Dispatch:

Nothing – that's the point! At least not in Julia. It doesn't make sense to say that the + belongs to the 1 or the 2 in the statement 1 + 2, and that's where the headaches come from: when you bind identity to objects, you're gonna have a bad time, as 80s lisp hackers and philosophers alike have struggled with that question for a long, long time – it's just devilishly hard to reason about identity when objects change. Avoid worrying about that if you can. Just worry about calling the right procedure in the right case. In other words, "just dispatch and carry on."

Single dispatch

This post promised comparing different languages. Now that we've discussed what we mean by dispatching, let's see how it's implemented. Astute readers have picked up that dispatching already exists in other languages, albeit in limited form: single dispatch.

I'm going to do some serious hand-waving here, so strap in: in Python and C++, and other object oriented languages, you have the barebones version of dispatch. If a function takes an argument f(a) you can only dispatch on the type of the first argument.

Specifically, we want the function f to behave differently when the types of a are different. Recall our previous examples for addition: if a is an Int, we want it to do one thing – if it's an Float64, another thing.

Whipping up a Python REPL:

>>> class Foo:
...     x = 1
...
>>> class Bar:
...     x = 1
...
>>> def f(a: Foo):
...     print("Just Foo")
...
>>> a = Foo()
>>> b = Bar()
>>> f(a)
Just Foo
>>> def f(b: Bar):
...     print("Just Bar")
...
>>> f(b)
"Just Bar"
>>> f(a)
"Just Bar"

💥WHOOPS!💥 – something went wrong, we should have gotten a "Just Foo" in that last line!

Multiple Dispatch

In the previous section, we wanted f to have two different behaviours depending on the types (and call it polymorphic operator overloading, if we want to bait some Dunning-Krugers on the god-forsaken orange-site). But we won't worry about fancy terms – we just want to our programming language tools to able to behave like the + we know from primary school.

miguelito:

OK – I'll take a crack at this. I've read and reread the Julia Manual page for Methods, and I think I have a better idea of this now. Here's what my Julia code looks like:

abstract type Things end # We'll come back to this line
struct Foo <: Things end
struct Bar <: Things end
f(::Things) = "Just a Thing"
f(x::Foo) = "Just Foo"
f(x::Bar) = "Just Bar"
x = Foo()
y = Bar()
f(x)
f(y)

And it works! I like thinking about this as a table, just like what we talked about with +: I just check what types the arguments I'm applying + to are, and apply to proper procedure. Adding integers means line them up and carry the digits. Fractions, common denominators, etc. For f, if I apply it to an Foo, the procedure is to return the string "Just Foo". If I apply f to an object Bar, it returns the string "Just Bar.". Dispatch and carry on, ach so…

Dispatch:

Great! You're on your way to learn the Zen of Julia!. It usually looks like

  1. Setup an abstract type like Things.

  2. Make some structs that are subtypes of Things with <:

  3. Define functions f that act on specific subtypes of Things – aka dispatching.

  4. Create your structs/objects with constructors like x = Foo() and then call f(x) to handle different behaviors.

That's it!

Small nitpick – You actually can dispatch in Python, but it requires a bit of boilerplate, but that's a secondary concern. Since it wasn't a foundational design of the language, people didn't build a vocabulary or an ecosystem around it. In the words of Turing award winner Barbara Liskov,

The quality of software depends primarily on the programming methodology in use. […] A methodology can be easy or difficult to apply in a given language, depending on how well the language constructs match the structures that the methodology deems desirable. […] – Barbara Liskov, Abstractions Mechanisms in CLU (1977)

(Shoutout to Quinn Wilton and her great Gleam presentation were I took this quote from). Basically, not all tools can fit nicely into a particular niche, and insisting otherwise is a recipe for frustration, but I guess some companies just have billions in cash for R+D to burn. It's not only honest, but necessary, to know the limitations of your own tools.

miguelito:

Huh – sounds like Julia got really… lucky (?) in that it didn't need to be the first to run up against these problems? That knowledge accrues over decades by loads of smart people by trying, failing, and figuring things out.

Dispatch:

🎉 Correct! 🎉 Julia has benefitted immensely from the efforts of others. We gain nothing from being smug about recent successes – there's still lots of problems to solve and it's in our best interests that we nurture a diverse community of people that we can cross-pollinate ideas with. Maybe someone implements multiple dispatch with some different tradeoffs in Python (like the Plum library!), or type class resolution in Lean or whatever they're building with F* that shows us a new way of thinking. We lose nothing by encouraging people to experiment, far and wide.

miguelito:

Hold up, you had mentioned that Julia's not the first to get multiple dispatch. Why didn't it pick up in the other languages?

Dispatch:

Hmmm, hard to say, I think we'd need to reach out to a legit PL historian for that. However, looking at some of the other key components that coalesce together helps suss some of it out:

  1. Common Lisp had an easy opt-in multiple dispatch system, but it was slow. People didn't buy in because it cost performance. There's a social factor to this – if your paradigm takes more effort to use, it's less likely to be grow.

  2. Performance was not an afterthought. Look at the graveyard of attempts to speed up Python, all mutually incompatible. The Julia devs designed the abstractions to match LLVM semantics to optimize for performance. At some point, you have to ask if you're standing in the right place to begin with, like the London cabby:

Excuse me, what's the best way to get to Manchester?

Ah, Manchester? I wouldn't start from here…

I don't think Python for performance has failed, but wouldn't start with a language that takes 28 bytes to store an integer:

Python 3.9.3 (default, Apr  8 2021, 23:35:02)
[GCC 10.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from sys import getsizeof
>>> getsizeof(1)
28

whereas

julia> sizeof(1) # This is in bytes, so an Int on my system is 64 bits
8

  1. Speaking of LLVM – JIT compilers are a recent invention, and if LLVM didn't have an open source implementation, who knows if Julia would have picked up. Caching the computation really helps to overcome the bureaucratic overhead of the dispatching system, and it's not at all trivial to see those two things and put them together. But don't take my word for it, hear them say why it wasn't practical then:

Lisp hackers preaching the Julia gospel

  1. That video tells precisely the pain-point that Stefan talked about for functional languages in his talk: In Functional Land, it's easy to add new procedures to old types, but unwieldy to add new types to old functions.

Concretely, the code that stumps object oriented PLs is

f(x::Foo) = "Just Foo"
f(x::Bar) = "Just Bar"

and the one that stumps functional language oriented PLs without multiple dispatch is

struct Qux <: Things end
t = Qux()
f(t)

miguelito:

Alright so what I'm really getting here is

  1. I should really, really just go watch Stefan's talk.

  2. A big "the Romans had a steam engine but didn't industrialize because they didn't have carbon deposits" vibe from all the convergence of PL ideas and techniques.

There's a time and a place I guess.

Dispatch:

Yeesh.

Remember the starting claim for this discussion? It felt so long ago… but the gist was that for all the bells and whistles that Julia has, they needed time and effort to figure out some hard problems that other people had come up against (and whose expertise they drew from!). Julia is the place to park a decision until it gets done right, with oodles of discussions from experts back and forth. That's not a linear process, but I can't complain, we're still increasing the SmileFactor of all the things that feel like they should work, and do. Like the REPL.

Until next time. Toodles. 👋

If you want to see more posts like this, consider chucking a buck or two on my GitHub sponsors, or, you know, hire me as a grad student.

Note:

I created Dispatch by copy/pasting the icon from flaticon.com under the terms of their Flaticon License. It is free for personal and commercial purpose with attribution. I changed the colors to match the Julia dot logo colors. If you plan to use it for commercial purposes, please donate a non-trivial part of your profits from the Dispatch merch to Doctors without Borders.

Thanks a lot to the Julia community for helping with this post, but especially to Thiebaut Lienart and the Franklin.jl team, Stefan Karpinski for his talk and Frames for her blog posts diving into these similar materials.