Jump into Julia

Today we build two new kinds of numbers in a few lines each: integers mod $p$, and dual numbers that compute derivatives. Then we use them with Julia's own code, such as its matrix solver and determinant, without changing that code. All you need is the Julia REPL and one small package.

0. Getting started

Install Julia

Go to julialang.org/downloads for full instructions. In short:

• Windows: in a Command Prompt, run winget install --name Julia --id 9NJNWW8PVKMN -e -s msstore
• macOS or Linux: in a terminal, run curl -fsSL https://install.julialang.org | sh

This installs the latest stable version of Julia. Then type julia in a new terminal window to start it. If you have time before the talk, also run ] add Primes (see section 2), since the first package install can take a few minutes.

REPL tips

• Start Julia by typing julia in a terminal, or open the Julia app.
• Type ? for help mode, ] for package mode, and Backspace to go back.
• Unicode: type \varepsilon then press Tab to get ε.
• methods(f) lists all the methods of a function f. @which f(x) shows which method Julia picks.

1. Multiple dispatch in one minute

A Julia function can have many methods. Julia picks the most specific method based on the types of all the arguments.

describe(x)          = "something"
describe(x::Number)  = "a number"
describe(x::Integer) = "an integer"
describe(x::String)  = "a string"

describe(1), describe(1.5), describe("hi"), describe([1, 2])
# ("an integer", "a number", "a string", "something")

methods(+)           # + itself has around 200 methods

Arithmetic such as + and * is just ordinary functions with many methods. We can add our own.

2. Integers mod $p$

"Mod 7" means we keep only the remainder after dividing by 7, as with days of the week: 100 days after a Tuesday is a Thursday, since dividing 100 by 7 leaves a remainder of 2. The integers mod a prime $p$ form a field: we can add, subtract, multiply, and divide by anything nonzero. We write $a \equiv b \pmod 7$ when $a$ and $b$ leave the same remainder after dividing by 7. Division means undoing multiplication: $3/5$ is the $x$ with $5x \equiv 3$, which is $x = 2$, since $5 \times 2 = 10$, and dividing 10 by 7 leaves a remainder of 3. (Mod 12 this fails. To compute $1/3$ we need an $x$ with $3x \equiv 1$. But multiples of 3 divided by 12 only ever leave remainders 0, 3, 6, or 9, never 1. So $1/3$ does not exist mod 12. When the number we divide by is prime, this never happens.)

We put $p$ in the type, so Mod{7} and Mod{11} are different types.

Installing a package

Julia has thousands of free packages. You can browse them at juliapackages.com. We need one, Primes.jl, for its isprime function. At the julia> prompt, press ] to enter package mode. The prompt changes to (@v1.13) pkg>. Type add Primes, then press Backspace to go back.

] add Primes              # one time only; the first add can take a minute
using Primes              # load it in each new session
isprime(7), isprime(91)   # (true, false)   91 = 7 × 13

The type

struct Mod{p} <: Integer
    k::Int
    function Mod{p}(k) where p
        @assert isprime(p)
        new(mod(k, p))
    end
end

import Base: +, -, *, /, inv, show, promote_rule, oneunit, abs

show(io::IO, a::Mod{p}) where p =
    print(io, get(io, :compact, false) ? a.k : "$(a.k) (mod $p)")

+(a::Mod{p}, b::Mod{p}) where p = Mod{p}(a.k + b.k)
-(a::Mod{p}) where p = Mod{p}(-a.k)
-(a::Mod{p}, b::Mod{p}) where p = a + -b
*(a::Mod{p}, b::Mod{p}) where p = Mod{p}(a.k * b.k)
inv(a::Mod{p}) where p = Mod{p}(invmod(a.k, p))
/(a::Mod{p}, b::Mod{p}) where p = a * inv(b)

promote_rule(::Type{Mod{p}}, ::Type{Int}) where p = Mod{p}   # Mod + Int gives Mod
oneunit(::Type{Mod{p}}) where p = Mod{p}(1)
abs(a::Mod{p}) where p = a.k        # used by \ to choose pivots

Try it

a = Mod{7}(3)
b = Mod{7}(5)
a + b        # 1 (mod 7)   3 + 5 = 8, and dividing 8 by 7 leaves a remainder of 1
a * b        # 1 (mod 7)   3 × 5 = 15, and dividing 15 by 7 leaves a remainder of 1
inv(a)       # 5 (mod 7)   since 3 × 5 ≡ 1
a / b        # 2 (mod 7)   since 5 × 2 = 10 ≡ 3
a + 10       # 6 (mod 7)   3 + 10 = 13 ≡ 6; Julia turns 10 into a Mod{7} for us
a^100        # 4 (mod 7)   we never wrote ^ ; Julia's generic power code works

Solving linear equations mod 7, for free

A matrix is a compact way to write a system of linear equations:

The answer is $x = 2$, $y = 3$ (check it). If we know the inverse matrix $A^{-1}$, then $\mathbf{x} = A^{-1}\mathbf{b}$, just as $5x = 10$ gives $x = 5^{-1} \cdot 10 = 2$. Julia can solve $A\mathbf{x} = \mathbf{b}$ either way:

A = [2 1; 1 3]
b = [7, 11]
x = A \ b    # [2.0, 3.0]
A * x        # [7.0, 11.0], which is b
inv(A)       # the inverse of A
inv(A) * b   # also [2.0, 3.0]

Both give the same answer. In practice A \ b is the better choice: it is faster and more accurate, because it never forms the inverse.

Now keep the same left-hand sides, but work mod 7:

Julia solves it with the same A \ b:

A = Mod{7}.([2 1; 1 3])
b = Mod{7}.([1, 2])
x = A \ b        # the solution, mod 7
A * x            # [1, 2], which is b

inv(A)           # the inverse matrix, mod 7
inv(A) * b       # the same solution (but in practice, use A \ b)

Check Julia's answer by hand: put its $x$ and $y$ into both equations, then find the remainders after dividing by 7.

Nobody who wrote Julia's \ knew about our Mod type. It works because the solver only needs +, -, *, /, and a way to pick pivots, and we supplied those.

3. Automatic differentiation with dual numbers

You may know $i$, a new number with $i^2 = -1$. A dual number uses a new number $\varepsilon$ that is not zero but has $\varepsilon^2 = 0$. No real number does this, just as no real number has square $-1$. Then

exactly, and $6$ is the derivative of $x^2$ at $x = 3$. In general,

So if we evaluate $f$ at $a + \varepsilon$, the $\varepsilon$ part of the answer is $f'(a)$. The rules below are the sum, product, and quotient rules from calculus.

struct Dual{T<:Number} <: Number
    p::T     # value
    d::T     # derivative
end
Dual(p, d) = Dual(promote(p, d)...)

import Base: +, -, *, /, inv, sin, cos, exp, show, promote_rule, convert

+(x::Dual, y::Dual) = Dual(x.p + y.p, x.d + y.d)              # sum rule
-(y::Dual)          = Dual(-y.p, -y.d)
-(x::Dual, y::Dual) = x + -y
*(x::Dual, y::Dual) = Dual(x.p * y.p, x.d * y.p + x.p * y.d)  # product rule
inv(y::Dual)        = Dual(1 / y.p, -y.d / y.p^2)
/(x::Dual, y::Dual) = x * inv(y)                               # quotient rule

sin(y::Dual) = Dual(sin(y.p),  cos(y.p) * y.d)                 # chain rule
cos(y::Dual) = Dual(cos(y.p), -sin(y.p) * y.d)
exp(y::Dual) = Dual(exp(y.p),  exp(y.p) * y.d)

promote_rule(::Type{Dual{T}}, ::Type{S}) where {T<:Number, S<:Real} =
    Dual{promote_type(T, S)}
convert(::Type{Dual{T}}, x::Real) where T<:Number = Dual{T}(x, zero(T))
convert(::Type{Dual{T}}, x::Dual) where T<:Number = Dual{T}(x.p, x.d)

show(io::IO, x::Dual) = print(io, x.p, " + ", x.d, "*ε")

ε = Dual(0, 1)

Try it

(3 + ε)^2                    # 9 + 6*ε
ε * ε                        # 0 + 0*ε
f(x)  = sin(exp(x + 1))
fp(x) = cos(exp(x + 1)) * exp(x + 1)    # derivative by hand

f(4.0), fp(4.0)              # (-0.6876914116945115, -107.74841574308525)
f(4.0 + ε)                   # -0.6876914116945115 + -107.74841574308525*ε

Second derivatives

A function can return a new function. Here D(g) returns the derivative of g. Apply D twice to get the second derivative. The inner D puts a dual number inside a dual number, so the two $\varepsilon$s stay separate.

D(g) = x -> g(Dual(x, one(x))).d

D(f)(4.0)                    # -107.74841574308525, the same as fp(4.0)
D(D(f))(4.0)                 # 15039.662941328417, which is f''(4.0)

It even works through loops

The Babylonian method computes $\sqrt{x}$ by repeating $t \leftarrow (t + x/t)/2$. We never tell it about derivatives.

function babylonian(x; N = 10)
    t = (1 + x) / 2
    for i in 2:N
        t = (t + x / t) / 2
    end
    return t
end

babylonian(2.0 + ε)          # 1.41421356237... + 0.35355339059...*ε
sqrt(2), 1 / (2sqrt(2))      # check: √2 and d/dx √x at x = 2

Bonus: the derivative of $\log\det X$

Here is a fact most people have not seen. For an invertible matrix $X$,

Why: expanding the determinant along row $i$ shows that $\partial \det X / \partial x_{ij} = C_{ij}$, the $(i,j)$ cofactor. Also $X^{-1} = C^{T} / \det X$. So $\partial \log\det X / \partial x_{ij} = C_{ij} / \det X = (X^{-1})_{ji}$.

Let us check it with dual numbers. We add $\varepsilon$ to one entry of $X$ and call Julia's own det. We need two small additions: a rule for log, and a method that lets Julia's LU code pivot on dual numbers.

using LinearAlgebra

Base.log(x::Dual) = Dual(log(x.p), x.d / x.p)
Base.abs(x::Dual) = abs(x.p)    # LU pivots on the value part

X = [4.0 1 2; 0 3 1; 1 0 5]
Y = Dual.(X, 0.0)        # X with zero ε parts
Y[1, 2] += ε             # add ε to the (1,2) entry
log(det(Y))              # 4.00733... + 0.0181818...*ε
inv(X)[2, 1]             # 0.0181818...  (it is 1/55)

Note the swapped indices: nudging entry $(1,2)$ gives entry $(2,1)$ of the inverse. As a sanity check, add a small real number $t$ instead:

t = 1e-6
E = zeros(3, 3); E[1, 2] = 1
(log(det(X + t*E)) - log(det(X))) / t    # about 0.0181818, close but not exact

Each dual-number pass gives one partial derivative, because $\varepsilon$ sits in one entry. This is forward-mode AD. Getting the whole gradient in one pass takes reverse-mode AD. In machine learning it is called backpropagation, and it is how neural networks are trained.

4. If time allows: why not just use $(f(x+h) - f(x))/h$?

The finite difference $\dfrac{f(x+h) - f(x)}{h}$ gets closer to $f'(x)$ as $h$ shrinks, but only up to a point. When $h$ is tiny, $f(x+h)$ and $f(x)$ agree in almost every digit, and subtracting them throws those digits away. This is called catastrophic cancellation.

fd(f, x, h) = (f(x + h) - f(x)) / h   # f is sin(exp(x + 1)) from section 3

f(4.0 + ε)               # dual numbers: -107.74841574308525 in the ε part
fp(4.0)                  # hand-written derivative: the same number

for k in 1:16
    h = 10.0^(-k)
    err = abs(fd(f, 4.0, h) - fp(4.0))
    println("h = 1e-", k, "    error = ", err)
end
# the error is smallest near h = 1e-9 (about 4e-6), then grows again;
# at h = 1e-16, 4.0 + h == 4.0 in floating point, so the error is the whole
# derivative, about 108

Dual numbers never subtract nearly equal numbers to get the derivative, so there is no $h$ to choose and no cancellation. The answer is as accurate as $f$ itself.

The takeaway

We wrote two small types. We did not write a linear solver, a matrix inverse, a power function, or a determinant rule for derivatives. Julia's generic code and our types work together because methods are chosen by the types of all the arguments. That is multiple dispatch.

Where to go next: juliapackages.com. For real work, ForwardDiff.jl is a full version of our Dual type, and Zygote.jl and Enzyme.jl do reverse-mode AD.