# Automatic Differentiation: Theory and Practice

Mario Lezcano-Casado Senior Software Engineer at Quansight
Correspondence<EMAIL_ADDRESS>
###### Abstract

We present the classical coordinate-free formalism for forward and backward
mode ad in the real and complex setting. We show how to formally derive the
forward and backward formulae for a number of matrix functions starting from
basic principles.

Contents

toc

This document is aimed at a reader who has some experience with auto-
differentiation (ad) frameworks (PyTorch, Tensorflow, JAX…) and would want to
get a better understanding of the mathematical ideas behind their forward and
backward ad engines. We expect that the reader will be able to derive from
scratch forward and backward formulae for virtually any function in PyTorch,
Tensorflow, or JAX after reading these notes.

We do not assume any previous mathematical knowledge, but some mathematical
maturity, as we will present the necessary results from different areas of
mathematics in the text. That said, some familiarity with linear algebra and
multivariate calculus would be helpful.

There will be a plethora of examples throughout the text to make the abstract
definitions easier to digest. We will omit some computations in the examples,
but these will always be simple enough so that the reader should be able to
deduce them themselves after some thinking.

Sections 1 and 2 contain a review of a number of definitions and results from
linear algebra and multivariate calculus. We encourage the reader that already
feels comfortable with this material to have a quick read through these, as we
will use them to set the notation.

Sections 3 and 4 show how to use the definitions from the first two sections
to derive the forward and backward (ad) formulae for most functions used in
machine learning or statistics.

Section 5 treats the more advance topic of differentiation of functions of
complex variable. It may be skipped on a first read.

##  1 Coordinate-free Linear Algebra

In this section we answer the question: What is a coordinate-free approach and
why do we care? With this objective in mind, we define abstract vector spaces,
bases, and linear maps. A secondary objective of this section is to show that
thinking of a linear map as something that preserves the operations of a
vector space is easier than thinking about it as a matrix.

###  1.1 Vector spaces

######  Definition 1.1.

A vector space over a field $\mathbb{K}$ (think $\mathbb{K}=\mathbb{R}$ or
$\mathbb{K}=\mathbb{C}$) is a set $V$ with two operations $+\colon V\times
V\to V$ and $\cdot\colon\mathbb{K}\times V\to V$ such that for every $u,v,w\in
V$ and $a,b\in\mathbb{K}$

  * •

Associativity: $(u+v)+w=u+(v+w)$

  * •

Identity: There exists a $0\in V$ such that
Identity: $\forall v\in V$, $v+0=v$

  * •

Inverse: For every $v\in V$ there is a $-v\in V$
Inverse: such that $v+(-v)=0$

  * •

Commutativity: $u+v=v+u$

  * •

Compatibility: $a\cdot(b\cdot v)=(ab)\cdot v$

  * •

Scalar identity: $1\cdot v=v$ where $1\in\mathbb{K}$
Scalar identity: is the multiplicative identity

  * •

Distributivity $1$: $a\cdot(u+v)=a\cdot u+a\cdot v$

  * •

Distributivity $2$: $(a+b)\cdot u=a\cdot u+b\cdot u$

As with any mathematical multiplication, we often omit the symbol and simply
write $av=a\cdot v$.

The elements of $V$ are called vectors and the elements of $\mathbb{K}$ are
called scalars.

Vector spaces over the real (resp. complex) numbers are called real (resp.
complex) vector spaces.

The first time one encounters these axioms, they can be quite overwhelming.
They become a bit easier to digest once one realises that all they are doing
is to model $\mathbb{R}^{n}$ abstractly. Once one sits down and checks that
$\mathbb{R}^{n}$ fulfils all these axioms, we find the first example of a
vector space.

######  Example 1.2.

$\mathbb{R}^{n}$ is a vector space over $\mathbb{R}$.

$\mathbb{R}^{n}$ is the most important example of a vector space. In fact,
virtually all the vector spaces over the real numbers that we will work with
will have an $\mathbb{R}^{n}$ lurking behind one way or another. Whenever we
think about an abstract vector space, it is good to picture in mind something
that looks like $\mathbb{R}^{n}$ to fix the ideas. That being said, there are
other vector spaces.

######  Example 1.3.

$\mathbb{C}^{n}$ is a vector space over $\mathbb{C}$.

$\mathbb{C}^{n}$ is also a real vector space, as addition of vectors is
performed on the real and complex part separately. In fact, $\mathbb{C}^{n}$
is pretty much the same as
$\mathbb{R}^{2n}=\mathbb{R}^{n}\times\mathbb{R}^{n}$ as a real vector space.

We finish with some slightly less conventional examples of vector spaces.

######  Example 1.4.

The set of polynomials in one variable of degree less or equal to $n$ is a
real vector space.

The set of polynomials in one variable of degree exactly equal to $n$ is not a
real vector space.

The set of matrices $\mathbb{R}^{m\times n}$ is a real vector space.

The set of tensors with real entries of a fixed shape is a real vector space.

The set of infinite sequences of real numbers
$(a_{i})_{i=1}^{\infty}=(a_{1},a_{2},\dots)$ is a real vector space.

It is a good exercise to convince oneself that all these examples are what
they claim to be. In general, to check that an object is indeed a vector
space, it tends to be enough to check that multiplying a vector by a scalar
gives a vector in the set and adding two vectors in the set gives a vector in
the set (the second example above is a counterexample of this).

###  1.2 Vector subspaces and bases

######  Definition 1.5.

A subset $U$ of a vector space $V$ over $\mathbb{K}$ is a vector subspace or
linear subspace if it is closed under addition and multiplication by scalars.
In symbols,

$u+v\in U\quad\text{and}\quad av\in U\mathrlap{\qquad\forall u,v\in
U,a\in\mathbb{K}}$

######  Example 1.6.

The vectors of the form $(0,a_{2},\ldots,a_{n})$ with $a_{i}\in\mathbb{R}$
form a linear subspace of $\mathbb{R}^{n}$.

The polynomials in one real variable of degree less than $n$ form a linear
subspace of the vector space of all polynomials in one real variable.

A simple way to define subspaces is to pick a number of vectors and consider
all the possible combinations of additions and multiplication by scalars that
one can form—i.e., consider their closure.

######  Definition 1.7.

The vector space spanned by $v_{1},\dots,v_{n}\in V$ is defined as

$\operatorname{span}(v_{1},\dots,v_{n})=\Big{\\{}\sum_{i=1}^{n}a_{i}v_{i}\>\Big{|}\>a_{i}\in\mathbb{K}\Big{\\}}\subseteq
V.$

The vector space spanned by a number of vectors is always a vector subspace of
$V$ by definition, as it is closed with respect to sums and product by
scalars.

######  Example 1.8.

We compute the span of some sets of vectors in $V=\mathbb{R}^{3}$:

  * •

Let $v_{1}=\lparen 0,1,0\rparen$. $\operatorname{span}\lparen
v_{1}\rparen=\\{\lparen 0,a,0\rparen\>|\>a\in\mathbb{R}\\}$.

  * •

Let $v_{2}=\lparen 0,1,1\rparen$. $\operatorname{span}\lparen
v_{1},v_{2}\rparen=\\{\lparen 0,a+b,b\rparen\>|\>a,b\in\mathbb{R}\\}$. Setting
$a=c-b$ for a new variable $c\in\mathbb{R}$, we have that
$\operatorname{span}\lparen v_{1},v_{2}\rparen=\\{\lparen
0,c,b\rparen\>|\>c,b\in\mathbb{R}\\}$. In plain words, $v_{1}$ and $v_{2}$
generate the subspace of all the vectors in $\mathbb{R}^{3}$ with first
coordinate equal to zero.

  * •

Let $v_{3}=\lparen 0,1,2\rparen$. $\operatorname{span}\lparen
v_{1},v_{2},v_{3}\rparen=\operatorname{span}\lparen
v_{1},v_{2}\rparen=\\{\lparen 0,a,b\rparen\>|\>a,b\in\mathbb{R}\\}$. What
causes this is that $v_{3}=2v_{2}-v_{1}$. As such, any vector that we can form
using $v_{3}$, we could already form as a linear combination of $v_{1}$ and
$v_{2}$. Intuitively, the third vector is “redundant”.

If we have vectors $v_{1},\dots,v_{n}$, and we express another vector as a
linear combination of these, $v=\sum_{i=1}^{n}a_{i}v_{i}$, we would like to be
able to represent $v$ as $(a_{1},\dots,a_{n})$, as we do when we are in
$\mathbb{R}^{n}$. This identification may not be possible because of two
reasons.

  1. 1.

Our set might not span all the vectors in $V$. In the examples in Example 1.8,
none of the sets of vectors that we considered was able to span the vector
$v=(1,0,0)$. Formally, we write $(1,0,0)\not\in\operatorname{span}\lparen
v_{1},v_{2},v_{3}\rparen$. More generally, it might be the case that
$\operatorname{span}\lparen v_{1},\dots,v_{n}\rparen\subsetneq V$.

  2. 2.

The representation $v=\sum_{i=1}^{n}a_{i}v_{i}$ might not be unique. This is
what happened in the last example in Example 1.8. Since $v_{3}=2v_{2}-v_{1}$,
if we choose a vector in $\operatorname{span}(v_{1},v_{2},v_{3})$, for example
$v=(0,2,2)$, we have that we can write

$v=2v_{2}=v_{1}+v_{3}=\tfrac{2}{3}v_{1}+\tfrac{2}{3}v_{2}+\tfrac{2}{3}v_{3}.$

Thus, we cannot identify $v=\sum_{i=1}^{3}a_{i}v_{i}$ with
$(a_{1},a_{2},a_{3})$ because the choice of $a_{i}$ is not unique.

These problems motivate the following definition:

######  Definition 1.9.

A set of vector $B=\\{v_{1},\dots,v_{n}\\}$ is called a basis of a vector
space $V$ if $\operatorname{span}(B)=V$ and any element $v\in V$ may be
represented uniquely in terms of the vectors in $B$. If there is a basis of
$n$ vectors for $V$, we say that $V$ has dimension $n$ or $\dim V=n$.

In plain words, a basis of a vector space is a set of vectors such that

  * •

Spans the whole $V$.

  * •

Has no redundancies: No vector in $B$ can be expressed as a linear combination
of others.

###  1.3 Linear maps

The other object that linear algebra studies are linear maps. These are maps
between vector spaces that respect the vector space structure.

######  Definition 1.10.

A map $T\colon V\to W$ between vector spaces over a field $\mathbb{K}$ is said
to be $\mathbb{K}$-linear—or simply linear—if

$T(u+v)=T(u)+T(v)\qquad T(av)=aT(v)\mathrlap{\qquad\forall u,v\in
V,a\in\mathbb{K}.}$

The following lemma shows that there is a very close connection between linear
maps and bases.

######  Lemma 1.11 (A linear map is defined by its values on a basis).

Let $V$ be a vector space with a basis $B_{V}=\\{v_{1},\dots,v_{n}\\}$. If we
choose $w_{1},\dots,w_{n}\in W$, there exists a unique linear map $T\colon
V\to W$ such that $T(v_{i})=w_{i}$.

###### Proof.

Since $B_{V}$ is a basis, there exists a unique representation of an arbitrary
element $v\in V$ as $v=a_{1}v_{1}+\dots+a_{n}v_{n}$. We can define $T$ on an
arbitrary vector $v$ as

$T(v)=T(a_{1}v_{1}+\dots+a_{n}v_{n})=a_{1}T(v_{1})+\dots+a_{n}T(v_{n}).$

so $T$ exists. To show that it is unique, suppose that there is another linear
map $T^{\prime}$ that takes the same values on $B_{V}$. We then have that they
are equal

$T(v)=a_{1}T(v_{1})+\dots+a_{n}T(v_{n})=a_{1}T^{\prime}(v_{1})+\dots+a_{n}T^{\prime}(v_{n})=T^{\prime}(a_{1}v_{1}+\dots+a_{n}v_{n})=T^{\prime}(v)$

where we have used that $T(v_{i})=T^{\prime}(v_{i})$ and the fact that both
$T$ and $T^{\prime}$ are linear. ∎

The best linear maps are those that map bases to bases.

######  Definition 1.12.

A linear map $T\colon V\to W$ is called a linear isomorphism if for any basis
$B$ of $V$, $T(B)=\\{T(v_{1}),\dots,T(v_{n})\\}$ is a basis of $W$.

This definition packs quite a bit of information. First, it says that $T$ is
surjective. This is because, since $T(B)$ is a basis of $W$,
$\operatorname{span}(T(B))=W$ and any $w\in W$ is in the image of $T$:

$w=a_{1}T(v_{1})+\dots+a_{n}T(v_{n})=T(a_{1}v_{1}+\dots+a_{n}v_{n}).$

It also says that this representation of a vector in terms of $T(B)$ is
unique, as that is the other property that bases have. In summary, we can
represent any vector in $W$ uniquely in terms of $B$ and $T$. Reciprocally, if
we denote basis $T(v_{i})=w_{i}$, by Lemma 1.11, there exists a unique linear
map $T^{\prime}\colon W\to V$ such that $T^{\prime}(w_{i})=v_{i}$. By using
the linearity of $T$ and $T^{\prime}$ It is easy to prove that

$T^{\prime}(T(v))=v\quad\text{and}\quad
T(T^{\prime}(w))=w\mathrlap{\qquad\forall v\in V,w\in W.}$

It should be clear now that a better name for this linear map would be
$T^{-1}$. In fact, we have proved the following:

######  Proposition 1.13.

A linear isomorphism as per Definition 1.12 has a linear inverse.111The
reciprocal also holds making this another possible definition of a linear
isomorphism.

This is another characterisation of linear isomorphisms which is very useful
in practice:

######  Proposition 1.14 (Characterisation of Linear Isomorphisms).

Let $T\colon V\to W$ be a linear map. If $T(B)$ is a basis of $W$ for one
basis $B$ of $V$, then $T(B^{\prime})$ is a basis of $W$ for any basis
$B^{\prime}$ of $V$, that is, $T$ is a linear isomorphism.

Linear isomorphisms are the nicest maps in linear algebra, as they allow us to
translate computations on one space to computations on another and back.

######  Example 1.15.

A basis of the set of polynomials $\mathcal{P}_{<n}$ with real coefficients of
degree less than $n$ is given by $B=\\{1,x,x^{2},\dots,x^{n-1}\\}$, as any
polynomial $p\in\mathcal{P}_{<n}$ is represented as a linear combination of
these as $p(x)=\sum a_{i}x^{i}$. Denoting $p_{i}(x)=x^{i}$, we may map this
basis of $\mathcal{P}_{<n}$ into the canonical basis of $\mathbb{R}^{n}$ via

$T(p_{0})=(1,0,\dots,0)\quad T(p_{1})=(0,1,\dots,0)\quad\dots\quad
T(p_{n-1})=(0,0,\dots,1)$

This extends to a linear map $T\colon\mathcal{P}_{<n}\to\mathbb{R}^{n}$ by
Lemma 1.11 and since it maps a basis to a basis, by Proposition 1.14, it is a
linear isomorphism. In other words, $\mathcal{P}_{<n}$ and $\mathbb{R}^{n}$
are _isomorphic_ as real vector spaces.

###  1.4 Differences between a real vector space and $\mathbb{R}^{n}$

For now, we have set up the basic ideas of linear algebra: We have a vector
space, which is an abstract space in which we can add vectors and multiply
vectors by scalars, we can encode most of the information of a vector space
into a basis, and we can map vector spaces to other vector spaces via linear
maps. In particular, if we have a linear isomorphism between vector spaces,
this allows us to translate operations on one to operations on the other and
back, rendering them somewhat equivalent.

On the other hand, while doing all this, we have barely talked about
$\mathbb{R}^{n}$. This is a bit odd given that it is the main space that we
want to study. The only thing that we have mentioned is that “real vector
spaces model $\mathbb{R}^{n}$” and that, as one could expect, $\mathbb{R}^{n}$
is a real vector space. But why would we care about real vector spaces in the
abstract? Why not simply work on $\mathbb{R}^{n}$?

There are two main differences between real vector spaces and
$\mathbb{R}^{n}$. One huge and one subtler.

The huge one is that a general vector space may have infinite dimensions. We
put an example of such a vector space in Example 1.4, when we talked about
infinite sequences of real numbers. Infinite-dimensional vector spaces are
whole different beasts, and we will not talk about them here.

Now, a real vector space $V$ of dimension $n$, it looks much more similar to
$\mathbb{R}^{n}$ as

######  Theorem 1.16.

Let $V$ be a real vector space of dimension $n$, there exists a linear
isomorphism $T\colon V\to\mathbb{R}^{n}$.222The real numbers are not important
here. A vector space of dimension $n$ over $\mathbb{K}$ is isomorphic to
$\mathbb{K}^{n}$ using the same argument.

###### Proof.

Choose a basis $B=\\{v_{1},\dots,v_{n}\\}$ of $V$ and perform the same
construction that we did in Example 1.15 constructing $T$ by sending
$T(v_{i})=e_{i}$. ∎

This isomorphism says that we can identify every element from $V$ with a
vector from $\mathbb{R}^{n}$. As such, if we need to do computations in $V$,
we can map $V$ into $\mathbb{R}^{n}$ by $T$, do the computations on
$\mathbb{R}^{n}$ and map the result back to $V$ via $T^{-1}$, potentially
simplifying the abstract space to working with $\mathbb{R}^{n}$.

In the same way, if we have a linear map between vector spaces, and we have
bases for these vector spaces, we can write this map in coordinates, giving
raise to the following well-known concept.

######  Definition 1.17.

Given a linear map between real vector spaces $T\colon V\to W$ with bases
$B_{V}=\\{v_{1},\dots,v_{n}\\},B_{W}=\\{w_{1},\dots,w_{m}\\}$, we define the
matrix associated to $T$ as the element $A\in\mathbb{R}^{m\times n}$ where the
$i$-th column is given by the coordinates of $T(v_{i})\in W$ in the basis
$B_{W}$. In symbols,

$T(v_{i})=A_{1,i}w_{1}+\dots+A_{m,i}w_{m}\mathrlap{\qquad\text{for
}i=1,\dots,n.}$

Matrices are representations of linear maps, and any $n$-dimensional real
vector space looks exactly the same as $\mathbb{R}^{n}$, so it looks like we
are saying that all these things are the same? Well, not quite, as there is a
subtle yet fundamental difference:

$\mathbb{R}^{n}$ is an abstract vector space $V$ of dimension $n$ together
with a choice of an ordered basis.

The point here is that, given a finite-dimensional abstract vector space,
there is no canonical choice of a basis. It is for this reason that the
approach to linear algebra that simply talks about vector spaces and linear
maps and not about matrices and $\mathbb{R}^{n}$ is often referred to as
coordinate-free. There are times, as it happened in Example 1.15, that there
exists a clear choice, but others there is not distinguished basis.

We will spend the rest of these notes showing how thinking about a map $T$
being linear if

$T(u+v)=T(u)+T(v)\qquad\text{and}\qquad T(au)=aT(u)\mathrlap{\qquad\forall
u,v\in V,a\in\mathbb{R}.}$

makes computations easier than thinking about $T$ as being a matrix. We show
an example as a taster:

######  Example 1.18 (Linear maps are simpler than matrices).

If we think of a linear map as being a matrix, it might not be clear at first
sight that the trace of a matrix $\operatorname{tr}\colon\mathbb{R}^{n\times
n}\to\mathbb{R}$ is a linear map. Now, when we look at it using its abstract
definition (Definition 1.10) this becomes obvious as

$\operatorname{tr}(A+B)=\operatorname{tr}(A)+\operatorname{tr}(B)\qquad\operatorname{tr}(cA)=c\operatorname{tr}(A)\mathrlap{\qquad\forall
A,B\in\mathbb{R}^{n\times n},c\in\mathbb{R}.}$

Transposing a matrix $A\mapsto A^{\intercal}$ or taking the first two columns
of a matrix are other examples of linear maps on matrices. We will see many
more in Section 3.

##  2 Multivariable Calculus

In this section, we go over the standard definition of a differentiable map
from multivariable calculus. We will avoid the language of coordinates—e.g.,
Jacobians and partial derivatives—whenever possible.

######  Definition 2.1.

A map $f\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ is differentiable at a point
$x\in\mathbb{R}^{m}$ if there exists a linear map
$L_{x}\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ such that

$\lim_{h\to 0_{\mathbb{R}^{m}}}\frac{f(x+h)-f(x)-L_{x}(h)}{\lVert
h\rVert_{\mathbb{R}^{m}}}=0_{\mathbb{R}^{n}}.$

In this case, we denote the linear map $\lparen df\rparen_{x}=L_{x}$, and we
call it the differential of $f$ at $x$. We say that $\lparen df\rparen_{x}(v)$
is the directional derivative of $f$ at $x$ in the direction $v$.

######  Remark (The differential is a first order linear approximation).

All this definition is saying is that there exists a linear first order
approximation of $f$ at $x$. In other words, if we subtract the approximation
from $f$, what we have left vanishes at zero slower than linearly—i.e., we
have removed all the “linear terms” of $f$. Another way to look at this is by
thinking that $\lparen\mathrm{d}f\rparen_{x}$ is the first term of the Taylor
expansion of $f$

$f(x+h)=f(x)+\lparen\mathrm{d}f\rparen_{x}(h)+o(h)$ (1)

where $f(x)$ is the $0^{\text{th}}$ order approximation to $f$ at $x$,
$\lparen\mathrm{d}f\rparen_{x}(h)$ is the $1^{\text{st}}$ order approximation
and the $o(h)$ denotes that the difference vanishes at zero slower than
linearly.333Formally, $o(h)$ means that there exists a function $E$ such that
$f(x+h)=f(x)+\lparen\mathrm{d}f\rparen_{x}(h)+E(h)\lVert h\rVert$ with
$\lim_{h\to 0}E(h)=0$. This is clear choosing
$E(h)=\frac{f(x+h)-f(x)-\lparen\mathrm{d}f\rparen_{x}(h)}{\lVert h\rVert}$.

Even though it is possible to generalise some of these concepts to deal with
functions that are not differentiable such as
$\operatorname{relu}(x)=\max(x,0)$,

We will always assume that the maps we consider are differentiable.

######  Remark (Directional derivatives).

We defined the directional derivative of a differentiable function $f$ at $x$
in the direction $v$ as $\lparen\mathrm{d}f\rparen_{x}(v)$. Now, since the
function is differentiable, the limit when $h$ approaches zero always exist,
so we can give a much more reasonable definition of directional derivative. We
can let $h$ tend to zero in the direction of $v\in\mathbb{R}^{m}$ with $\lVert
v\rVert=1$, by letting $h=\varepsilon v$ for $\varepsilon\in\mathbb{R}$,
getting that

$\lparen\mathrm{d}f\rparen_{x}(v)=\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\varepsilon}$}\Big{|}_{\varepsilon=0}f(x+\varepsilon
v)$

where
$\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\varepsilon}$}\Big{|}_{\varepsilon=0}$
differentiates $f(x+\varepsilon v)\in\mathbb{R}^{n}$ coordinate-wise at
$\varepsilon=0$. This shows how
$\lparen\mathrm{d}f\rparen_{x}(v)\in\mathbb{R}^{n}$ represents how $f$ varies
when $x$ is approached in the direction $v$.

######  Remark (The matrix associated to the differential is the Jacobian).

The differential of $f$ is a linear map
$\lparen\mathrm{d}f\rparen_{x}\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ and as
such, it has a matrix representation (Definition 1.17). The columns of this
matrix representation of $\lparen\mathrm{d}f\rparen_{x}$ are given by
evaluating $\lparen\mathrm{d}f\rparen_{x}$ on the vectors $\\{e_{i}\\}$ of the
basis of $\mathbb{R}^{n}$

$\lparen\mathrm{d}f\rparen_{x}(e_{i})=\mbox{{}\small$\displaystyle\frac{\mathrm{d}}{\mathrm{d}\varepsilon}$}\Big{|}_{\varepsilon=0}f(x+\varepsilon
e_{i})=\frac{\partial f}{\partial x_{i}}(x).$

In other words, the matrix associated to the differential is the usual
Jacobian matrix.

It will come to no surprise that we will never talk about partial derivatives
nor use the Jacobian. We will just use the definition of the differential
together with two results. The first one is the chain rule:

######  Theorem 2.2 (Chain rule).

Let $f\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ and
$g\colon\mathbb{R}^{n}\to\mathbb{R}^{p}$ be two differentiable maps. We have
that

$\mathrm{d}\lparen g\circ
f\rparen_{x}=\lparen\mathrm{d}g\rparen_{f(x)}\circ\lparen\mathrm{d}f\rparen_{x}\mathrlap{\qquad\forall
x\in\mathbb{R}^{m}}$

or equivalently

$\mathrm{d}\lparen g\circ f\rparen_{x}\lparen
v\rparen=\lparen\mathrm{d}g\rparen_{f(x)}[\lparen\mathrm{d}f\rparen_{x}\lparen
v\rparen]\mathrlap{\qquad\forall x,v\in\mathbb{R}^{m}.}$

This formula is the crux of all ad engines, as we will see in the sequel. Note
that this formula defines the equality between a linear map and the
composition of two other linear maps. It is also worth noting that the domains
and codomains of these maps are compatible:

Differentiable Maps
$\mathbb{R}^{m}$$\mathbb{R}^{n}$$\mathbb{R}^{p}$$f$$g$$g\circ f$ Linear Maps
$\mathbb{R}^{m}$$\mathbb{R}^{n}$$\mathbb{R}^{p}$$\lparen\mathrm{d}f\rparen_{x}$$\lparen\mathrm{d}g\rparen_{f(x)}$$\mathrm{d}\lparen
g\circ f\rparen_{x}$

The other rule that we will use repeatedly is a more abstract version of the
derivative of the product. This roughly says that “differentiating a function
of two variables accounts for differentiating the first variable fixing the
second one plus differentiating the second variable fixing he first one”.

######  Proposition 2.3 (Leibnitz rule).

Let $f\colon\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{p}$ be a
differentiable function. Define $f_{1,x}\colon\mathbb{R}^{m}\to\mathbb{R}^{p}$
for $x\in\mathbb{R}^{n}$ as $f$ partially evaluated on the first variable on
$x\in\mathbb{R}^{n}$, that is, $f_{1,x}(y)\coloneqq f(x,y)$. Define also
$f_{2,y}(x)=f(x,y)$. We have that

$\lparen\mathrm{d}f\rparen_{(x,y)}\lparen
e_{1},e_{2}\rparen=\lparen\mathrm{d}f_{2,y}\rparen_{x}(e_{1})+\lparen\mathrm{d}f_{1,x}\rparen_{y}(e_{2})$

The formula in this proposition may look quite difficult to parse, but it will
be much easier to understand when we use it in the next section to compute the
differential of some matrix functions.

######  Remark (On gradients).

Note that we have not defined the gradient of a function yet, but just its
differential. The definition of the gradient of a function will have to wait
until Section 4.

##  3 Forward Mode ad

###  3.1 The model and definition

After introducing the necessary concepts from linear algebra and multivariate
calculus, we are ready to put them to good use in the context of ad. Let us
consider a model with two differentiable maps

$\mathbb{R}^{m}$$\mathbb{R}^{n}$$\mathbb{R}^{p}$$g$$f$ Figure 1: The automatic
differentiation model.

######  Remark (Several inputs and outputs).

The first simplification that we have performed in Figure 1 is that, for a map
$h\colon\mathbb{R}^{n}\times\mathbb{R}^{m}\to\mathbb{R}^{p}$, we can always
see it as a map from $\mathbb{R}^{n+m}$ into $\mathbb{R}^{p}$. As such, it is
enough to talk about maps with just one input.444Sometimes, it may be
beneficial to look at maps of several inputs though, as some properties such
as linearity are not preserved by this transformation. We will see examples of
this later in this section. The same happens with maps with several outputs,
only that in this case we can see them as separate maps of one output, and we
can differentiate them separately.

######  Remark (Reduction from general ad).

When we do automatic differentiation, we first form a directed acyclic graph
of dependencies between maps—the model—ending in a map that outputs a real
number—the loss function. Without entering into the details of how to define
this graph of dependencies, we can think of $f$ as the last step of the graph
(in topological order) that produces the real number—i.e., set $p=1$ above—and
let $g$ be the rest of the graph. Then, if we want to split $g$ further, we
can consider its last step $g_{2}$ and the rest of the steps $g_{1}$, write
$g=g_{1}\circ g_{2}$ and proceed inductively, as we are in the same situation
as above.

######  Remark (Neural networks).

In the case of neural networks, we have a function $F(\theta,x)$ where
$\theta$ are its parameters concatenated and $x$ is an example from our
dataset. We then want to differentiate with respect to the parameters. This is
exactly the same idea as above, where we have the function $\theta\mapsto
F(\theta,x)$ for a fixed $x$.

Of course, this example does not encode _all_ possible neural network
architectures, as we have not mentioned what to do when the functions involved
are not differentiable, or when we use integers or have if-else constructions.
That being said, extending the theory presented here to all these ideas is not
formally challenging; the complexity lies in the implementation of the
resulting algorithm.

Now that we have everything set, we are ready to define forward mode automatic
differentiation.

######  Definition 3.1.

Forward mode ad for the model represented in Figure 1 accounts for computing

$\mathrm{d}\lparen g\circ f\rparen_{x}\mathrlap{\qquad\text{for
}x\in\mathbb{R}^{m}.}$

In plain words, forward mode ad computes the differentials of the model with
respect to the parameters. This may be achieved incrementally via the chain
rule (Theorem 2.2), as it tells us how to put together the differential of two
functions to compute the differential of the composition. These ideas are
often presented in the literature as “forward mode ad accumulates the product
of the Jacobians”, which is the same idea but in coordinates.

######  Remark (Dual numbers).

Forward mode ad is frequently defined in terms of dual numbers. Dual numbers
are defined through an abstract quantity called $\varepsilon$ with the
property that $\varepsilon^{2}=0$. Points are then described as
$v+\dot{v}\varepsilon$ and the following expansion is stated

$f\lparen
v+\dot{v}\varepsilon\rparen\stackrel{{\scriptstyle?}}{{=}}f(v)+f^{\prime}(v)\dot{v}\varepsilon.$
(2)

Now, this second equality is often not justified, and left to the reader to
interpret.555This is not really true. These identities can be formalised via
_perturbation theory_ [1]. It is typically shown by means of the computation
of the derivative of the product of real numbers.

$\lparen v+\dot{v}\varepsilon\rparen\lparen
u+\dot{u}\varepsilon\rparen=uv+(\dot{u}v+u\dot{v})\varepsilon\mathrlap{\qquad
u,v,\dot{u},\dot{v}\in\mathbb{R}.}$

A moment’s reflection shows that all this approach is encoding is the idea
that the differential is a first order approximation to the function! Having
an $\varepsilon$ such that $\varepsilon^{2}=0$ simply says that we just care
about the terms that are linear in epsilon, and we discard any term of order
two or higher. This is exactly what the idea of the differential of a map
formalises. As such, a formal justification of (2) is then given by the order
$1$ Taylor expansion in (1).

The example of the multiplication of real numbers can then be described in the
language of calculus as letting $f\colon\mathbb{R}^{2}\to\mathbb{R}$ defined
by $f(x,y)=xy$ and computing $\lparen\mathrm{d}f\rparen_{\lparen
u,v\rparen}\lparen\dot{u},\dot{v}\rparen=\dot{u}v+u\dot{v}$.666We will prove
this formula and generalise it to vectors and matrices later in this section.

###  3.2 Computing forward mode ad

We will spend the rest of this section showing how the abstract definition of
the differential (Definition 2.1) makes computations surprisingly easy.

######  Proposition 3.2 (Differential of a linear map).

Let $T\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ be a linear map, we have that

$\lparen\mathrm{d}T\rparen_{x}(v)=T(v)\mathrlap{\qquad\forall
x,v\in\mathbb{R}^{m}.}$

###### Proof.

Plugging $\lparen\mathrm{d}T\rparen_{x}(v)=T(v)$ into the definition of the
differential we get

$\lim_{h\to 0_{\mathbb{R}^{m}}}\frac{T(x+h)-T(x)-T(h)}{\lVert
h\rVert_{\mathbb{R}^{m}}}=\lim_{h\to
0_{\mathbb{R}^{m}}}\frac{T(x)+T(h)-T(x)-T(h)}{\lVert
h\rVert_{\mathbb{R}^{m}}}=\lim_{h\to
0_{\mathbb{R}^{m}}}\frac{0_{\mathbb{R}^{n}}}{\lVert
h\rVert_{\mathbb{R}^{m}}}=0_{\mathbb{R}^{n}}.\qed$

We can use this to compute the differential of a number of functions widely
used in machine learning.

######  Example 3.3 (Differential of a linear layer. Trailing batch
dimension).

Fix $X\in\mathbb{R}^{m\times b}$ a batch of $b$ vectors of size $m$ (trailing
batch dimension), and let $A\in\mathbb{R}^{n\times m}$. We can define

$\\!\begin{aligned} f\colon\mathbb{R}^{n\times m}&\to\mathbb{R}^{n\times b}\\\
A&\mapsto AX\end{aligned}\mathrlap{\qquad\qquad X\in\mathbb{R}^{m\times b}.}$

This is just the usual definition of a linear layer depending on the parameter
and with fixed inputs. We do this because we want to differentiate with
respect to the parameters. It is clear that

$f(A+B)=f(A)+f(B)\qquad f(cA)=cf(A)\mathrlap{\qquad\forall
A,B\in\mathbb{R}^{n\times m},c\in\mathbb{R}.}$

so $f$ is linear and by Proposition 3.2

$\lparen\mathrm{d}f\rparen_{A}(E)=EX\mathrlap{\qquad{\text{for
}E\in\mathbb{R}^{n\times m}}.}$

######  Example 3.4 (Differential of a linear layer. Front batch dimension).

The batch dimension in machine learning is often the first dimension of the
tensor due to the layout of matrices in memory. For this reason, it is common
to write a batch $b$ of $m$-dimensional vectors as $X\in\mathbb{R}^{b\times
m}$. In this case, we may write a linear layer as

$\\!\begin{aligned} f\colon\mathbb{R}^{n\times m}&\to\mathbb{R}^{b\times n}\\\
A&\mapsto XA^{\intercal}\end{aligned}\mathrlap{\qquad\qquad
X\in\mathbb{R}^{b\times m}}$

This function has a transpose and a matrix multiplication. In particular, if
we write $g(A)=A^{\intercal}$ and $h(A)=XA$ we have that $f=h\circ g$.
Luckily, both the transpose and the matrix multiplication are linear
functions, so we can use the chain rule (Theorem 2.2) and the formula for the
differential of a linear function twice to compute the differential of $f$:

$\lparen\mathrm{d}f\rparen_{A}(E)=\mathrm{d}\lparen h\circ
g\rparen_{A}(E)=[\lparen\mathrm{d}h\rparen_{g(A)}\circ\lparen\mathrm{d}g\rparen_{A}](E)=\lparen\mathrm{d}h\rparen_{A^{\intercal}}(g(E))=h(E^{\intercal})=XE^{\intercal}.$

A simpler way of performing this computation is by noting that $f$ is linear
itself so

$\lparen\mathrm{d}f\rparen_{A}(E)=f(E)=XE^{\intercal}.$

######  Example 3.5 (More linear maps).

Linear maps come in different shapes and forms

  * •

Inner product of vectors. Let $f_{2,y}\colon\mathbb{R}^{n}\to\mathbb{R}$ with
$f_{2,y}(x)=x^{\intercal}y=\sum_{i=1}^{n}x_{i}y_{i}$ for a fixed
$y\in\mathbb{R}^{n}$, then
$\lparen\mathrm{d}f_{2,y}\rparen_{x}(v)=v^{\intercal}y$. Fixing the first
variable, $f_{1,x}(y)=x^{\intercal}y$,
$\lparen\mathrm{d}f_{1,x}\rparen_{y}(v)=x^{\intercal}v$.

  * •

Trace of a matrix. Let $f\colon\mathbb{R}^{n\times n}\to\mathbb{R}$ with
$f(A)=\operatorname{tr}\lparen A\rparen$, then
$\lparen\mathrm{d}f\rparen_{A}(E)=\operatorname{tr}\lparen E\rparen$.

  * •

Inner product of matrices. Let $f_{1,X}\colon\mathbb{R}^{m\times
n}\to\mathbb{R}$ with $f_{1,X}(A)=\operatorname{tr}\lparen
A^{\intercal}X\rparen=\sum_{i=1}^{m}\sum_{j=1}^{n}A_{ij}X_{ij}$ for a fixed
$X\in\mathbb{R}^{m\times n}$. 777Note that $\operatorname{tr}\lparen
A^{\intercal}X\rparen$ is just a convenient way to represent the inner product
of matrices as seen as vectors of size $\mathbb{R}^{N}$ with $N=mn$. We then
have $\lparen\mathrm{d}f_{1,X}\rparen_{A}(E)=\operatorname{tr}\lparen
E^{\intercal}X\rparen$ and an analogous formula for the second variable.

######  Example 3.6 (Several inputs).

Consider the inner product of vectors as a function of two arguments

$\\!\begin{aligned} f\colon\mathbb{R}^{n}\times\mathbb{R}^{n}&\to\mathbb{R}\\\
x,y&\mapsto x^{\intercal}y;\end{aligned}$

and define the function partially evaluated in its first and second argument
as $f_{i,x}\colon\mathbb{R}^{n}\to\mathbb{R}$ for $i=1,2$ so that
$f_{1,x}(y)=f_{2,y}(x)=f(x,y)$. We have that $f_{i,x}$ are linear, as for
$i=1,2$

$f_{i,u}\lparen v+w\rparen=f_{i,u}\lparen v\rparen+f_{i,u}\lparen
w\rparen\qquad f_{i,u}\lparen cv\rparen=cf_{i,u}\lparen
v\rparen\qquad\qquad\forall u,v,w\in\mathbb{R}^{n},c\in\mathbb{R}.$

For this reason, we can compute its differential using Leibnitz rule
(Proposition 2.3)

$\lparen\mathrm{d}f\rparen_{(x,y)}\lparen
e_{1},e_{2}\rparen=\lparen\mathrm{d}f_{2,y}\rparen_{x}(e_{1})+\lparen\mathrm{d}f_{1,x}\rparen_{y}(e_{2})=f_{2,y}(e_{1})+f_{1,x}(e_{2})=e_{1}^{\intercal}y+x^{\intercal}e_{2}\qquad
x,y,e_{1},e_{2}\in\mathbb{R}^{n}.$

In contrast, note that $f$ itself is not linear as a function from
$\mathbb{R}^{2n}$ to $\mathbb{R}$ as $f(ax,ay)=a^{2}f(x,y)$.

######  Example 3.7 (Powers of a matrix).

Consider the map that multiplies a matrix with itself $k$ times

$\\!\begin{aligned} f\colon\mathbb{R}^{n\times n}&\to\mathbb{R}^{n\times n}\\\
A&\mapsto A^{k}=A\stackrel{{\scriptstyle k)}}{{\cdots}}A\end{aligned}$

this is the same as evaluating the map $g(A_{1},\dots,A_{k})=A_{1}\dots A_{k}$
for $A_{i}\in\mathbb{R}^{n\times n}$ at $(A,A,\dots,A)$. We can compute the
differential of $g$ using Proposition 2.3 since $g$ is linear in every entry:

$\lparen\mathrm{d}g\rparen_{\lparen A_{1},\dots,A_{k}\rparen}\lparen
E_{1},\dots,E_{k}\rparen=E_{1}A_{2}\cdots A_{k}+A_{1}E_{2}\cdots
A_{k}+\dots+A_{1}A_{2}\cdots E_{k}.$

So the differential of $f$ is given by

$\lparen\mathrm{d}f\rparen_{A}(E)=\lparen\mathrm{d}g\rparen_{\lparen
A,\dots,A\rparen}\lparen
E,\dots,E\rparen=EA^{k-1}+AEA^{k-2}+\dots+A^{k-1}E=\sum_{i=0}^{k-1}A^{i}EA^{k-i-1}.$

This example shows that, morally, if we can write a map $f$ as a map $g$ on
more variables such that $g$ is linear in each of its variables, all we need
to do to compute the differential of $f$ is to substitute each appearance of
$A$ by $E$ on $f$ and add them all together. More generally, if the function
is not linear in some of the variables, we substitute every appearance of
$g(A)$ by $\lparen\mathrm{d}g\rparen_{A}(E)$ as described in Proposition 2.3.
We show this idea in the following example.

######  Example 3.8 (Matrix inverse).

Let $\operatorname{GL}\lparen n\rparen\subseteq\mathbb{R}^{n\times n}$ be the
set of invertible matrices. Define

$\\!\begin{aligned} f\colon\operatorname{GL}\lparen
n\rparen&\to\operatorname{GL}\lparen n\rparen\\\ A&\mapsto
A^{-1}\end{aligned}$

We have that, by definition of the matrix inverse

$Af(A)=\mathrm{I}_{n}\mathrlap{\qquad\forall A\in\operatorname{GL}\lparen
n\rparen.}$

Defining $g(A)=Af(A)$ this identity can be rewritten as $g(A)=\mathrm{I}_{n}$
for $A\in\operatorname{GL}\lparen n\rparen$. This is an equality between
functions—one of them constant—so we may differentiate them.888Formally, we
would first need to define what does it mean to differentiate over
$\operatorname{GL}\lparen n\rparen$. Luckily, $\operatorname{GL}\lparen
n\rparen$ is an open subset of $\mathbb{R}^{n\times n}$, and since the
definition of differential is local, we can always define the differential at
any matrix $A\in\operatorname{GL}\lparen n\rparen$ by restricting the limit in
Definition 2.1 to a neighbourhood of $A$. It is direct to see from Definition
2.1 that the differential of a constant map is the function that maps any $E$
to the zero matrix. On the left-hand side we apply Proposition 2.3 to get

$EA^{-1}+A\lparen\mathrm{d}f\rparen_{A}(E)=0_{n\times n}$

and solving for $\lparen\mathrm{d}f\rparen_{A}(E)$ we get

$\lparen\mathrm{d}f\rparen_{A}(E)=-A^{-1}EA^{-1}.$

Note that this is a far-reaching generalisation of the result
$(1/x)^{\prime}=-1/x^{2}$ for
$x\in\mathbb{R}\backslash\\{0\\}=\operatorname{GL}\lparen 1\rparen$.

Before giving the last result, we show how to extend functions on the real
numbers to matrix functions.

######  Definition 3.9.

Let $f_{\mathbb{R}}(x)=\sum_{k=0}^{\infty}c_{k}x^{k}$ be an analytic
function—i.e., a function equal to its Taylor series. We define its associated
matrix function as999A matrix function is defined at $A\in\mathbb{R}^{n\times
n}$ if and only if all the eigenvalues of $A$ lie in the domain of definition
of $f$ when seen as a function from $\mathbb{C}$ to $\mathbb{C}$.

$\\!\begin{aligned} f\colon\mathbb{R}^{n\times n}&\to\mathbb{R}^{n\times n}\\\
A&\mapsto\sum_{k=0}^{\infty}c_{k}A^{k}.\end{aligned}$

######  Example 3.10.

Any function with a Taylor series can be turned into a matrix function:

  * •

Exponential: $\exp(A)=\sum_{k=0}^{\infty}\frac{1}{k!}A^{k}$

  * •

Logarithm: $\log\lparen
I_{n}+A\rparen=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{k+1}A^{k}$

  * •

Sine: $\sin(A)=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k+1)!}A^{2k+1}$

  * •

Cosine: $\cos(A)=\sum_{k=0}^{\infty}\frac{(-1)^{k}}{(2k)!}A^{2k}$

We present the last and most general result of this section, which can be
roughly summarised as:

If we know how to approximate a matrix function, we know how to approximate
its differential.

######  Theorem 3.11 (Differential a Matrix Function [2]).

Let $f\colon\mathbb{R}^{n\times n}\to\mathbb{R}^{n\times n}$ be a matrix
function (Definition 3.9). Applying $f$ on a matrix of size $2n\times 2n$, we
get the following result by blocks:

$f\begin{pmatrix}A&E\\\
0&A\end{pmatrix}=\begin{pmatrix}f(A)&\lparen\mathrm{d}f\rparen_{A}(E)\\\
0&f(A)\end{pmatrix}\mathrlap{\qquad\forall A,E\in\mathbb{R}^{n\times n}.}$

###### Proof.

Differentiating the series term by term and using Example 3.7 we have that

$\lparen\mathrm{d}f\rparen_{A}(E)=\sum_{k=0}^{\infty}c_{k}\sum_{i=0}^{k-1}A^{i}EA^{k-i-1}.$
(3)

We can also compute the powers of the block matrix

$\begin{pmatrix}A&E\\\
0&A\end{pmatrix}^{k}=\begin{pmatrix}A^{k}&\sum_{i=0}^{k-1}A^{i}EA^{k-i-1}\\\
0&A^{k}\end{pmatrix}$

so

$f\begin{pmatrix}A&E\\\
0&A\end{pmatrix}=\sum_{k=0}^{\infty}c_{k}\begin{pmatrix}A&E\\\
0&A\end{pmatrix}^{k}=\sum_{k=0}^{\infty}\begin{pmatrix}c_{k}A^{k}&c_{k}\sum_{i=0}^{k-1}A^{i}EA^{k-i-1}\\\
0&c_{k}A^{k}\end{pmatrix}=\begin{pmatrix}f(A)&\lparen\mathrm{d}f\rparen_{A}(E)\\\
0&f(A)\end{pmatrix}.\qed$

##  4 Backward Mode ad

In this section, we go over the most popular method of automatic
differentiation: Backward ad. This method has the advantage that, in order to
compute the backward pass for a model, one does not have to deal with
matrices—the Jacobians associated to the differential of the model—but just
with vectors of the same size of the parameters.

###  4.1 Inner products, gradients, and adjoints

Before introducing the definition of backward mode ad, we need to define one
more mathematical concept from linear algebra.

######  Definition 4.1.

Let $V$ be a real vector space, a (real) inner product is a map

$\\!\begin{aligned} \langle-,-\rangle\colon V\times V&\to\mathbb{R}\\\
x,y&\mapsto\langle x,y\rangle\end{aligned}$

such that it is

  * •

Bilinear: It is linear in each variable.

  * •

Symmetric: $\langle x,y\rangle=\langle y,x\rangle$ for every $x,y\in V$.

  * •

Positive definite: $\langle x,x\rangle>0$ for every $x\in V$, $x\neq 0$.

We will write spd as short for symmetric positive definite matrix, as we will
use them in examples.

######  Example 4.2.

The following are examples of inner products

  * •

Canonical inner product on $\mathbb{R}^{n}$: $\langle
x,y\rangle=x^{\intercal}y$ for $x,y\in\mathbb{R}^{n}$.

  * •

Other inner products on $\mathbb{R}^{n}$: $\langle x,y\rangle=x^{\intercal}Hy$
for $x,y\in\mathbb{R}^{n}$ and a fixed $H\in\mathbb{R}^{n\times n}$
spd.101010To prove that this is positive definite, consider the Cholesky
decomposition of $H=U^{\intercal}U$ with $U$ upper-triangular.

  * •

Canonical inner product on $\mathbb{R}^{m\times n}$: $\langle
A,B\rangle=\operatorname{tr}\lparen A^{\intercal}B\rparen$ for
$A,B\in\mathbb{R}^{m\times n}$.

  * •

Other inner products on $\mathbb{R}^{m\times n}$: $\langle
A,B\rangle=\operatorname{tr}\lparen A^{\intercal}HB\rparen$ for
$A,B\in\mathbb{R}^{m\times n}$ and $H\in\mathbb{R}^{m\times m}$ spd.

Inner products allow us to measure norms of vectors $\lVert
x\rVert=\sqrt{\langle x,x\rangle}$, angles between vectors $\angle\lparen
x,y\rparen=\arccos\frac{\langle x,y\rangle}{\lVert x\rVert\lVert y\rVert}$,
distances $d(x,y)=\lVert x-y\rVert$, and many other metric properties. As
such, it will come to no surprise the fact that inner products are very
important in machine learning and optimisation. For one, we need them to talk
about the distance from a point to the optimum and rates of convergence.
Perhaps less known is the fact that we also require them to talk about
gradients.

######  Remark (Motivating the concept of gradient).

For a function $f\colon\mathbb{R}^{n}\to\mathbb{R}$ and an
$x\in\mathbb{R}^{n}$, the map $v\mapsto\lparen\mathrm{d}f\rparen_{x}(v)$ is a
linear function from $\mathbb{R}^{n}$ to $\mathbb{R}$. Now, if we have an
inner product $\langle-,-\rangle$ on $\mathbb{R}^{n}$, for a fixed
$g\in\mathbb{R}^{n}$, the function $v\mapsto\langle g,v\rangle$ is also a
linear function from $\mathbb{R}^{n}$ to $\mathbb{R}$. The question now is,
given an inner product on $\mathbb{R}^{n}$ and a function $f$, can we always
represent the differential of $f$ as a vector $g_{x}\in\mathbb{R}^{n}$ such
that $\lparen\mathrm{d}f\rparen_{x}(v)=\langle g_{x},v\rangle$? This is, in
fact, the case, and it is the definition of a well-known concept.

######  Definition 4.3.

Let $f\colon\mathbb{R}^{n}\to\mathbb{R}$, and let $\langle-,-\rangle$ be an
inner product on $\mathbb{R}^{n}$. We define the gradient of $f$ at
$x\in\mathbb{R}^{n}$ as the vector $\nabla f(x)\in\mathbb{R}^{n}$ such that

$\lparen\mathrm{d}f\rparen_{x}(v)=\langle\nabla
f(x),v\rangle\mathrlap{\qquad\forall v\in\mathbb{R}^{n}.}$

######  Remark.

A number of remarks are in order.

  * •

As $x\mapsto v^{\intercal}x$ for $x,v\in\mathbb{R}^{n}$ is a linear function,
some people like to think informally of vectors as “column vectors” and linear
functions as “row vectors”. This way, the operation of going from a
differential $v^{\intercal}\cdot-$ to a gradient $v$ for the canonical inner
product looks like “transposing” $v^{\intercal}$.

  * •

Important. The gradient and the differential of a function are not the same
thing. The first one is a function into the real numbers, while the latter one
is a vector.

  * •

The gradient of a function depends on the choice of inner product, the
differential on finite-dimensional spaces does not, since all the norms are
equivalent.

######  Example 4.4.

We compute the gradient of some functions building on results from Section 3.

  * •

Consider $\mathbb{R}^{n}$ with the canonical inner product, and let
$f(x)=\langle g,x\rangle=g^{\intercal}x$ for a fixed $g\in\mathbb{R}^{n}$.
Since $f$ is linear $\lparen\mathrm{d}f\rparen_{x}(v)=\langle g,v\rangle$, and
by definition of a gradient, $\nabla f(x)=g$ for every $x\in\mathbb{R}^{n}$.

  * •

Consider $\mathbb{R}^{n}$ with an arbitrary inner product, and let
$f(x)=\langle g,x\rangle$ then $\nabla f(x)=g$ for every $x\in\mathbb{R}^{n}$.

  * •

Consider $\mathbb{R}^{n}$ with the inner product $\langle
u,v\rangle=u^{\intercal}Hv$ for $H$ spd (see Example 4.2) and let
$f(x)=x^{\intercal}y$ for a fixed $y\in\mathbb{R}^{n}$. As in the first
example, $\lparen\mathrm{d}f\rparen_{x}(v)=v^{\intercal}y$ since $f$ is
linear. On the other hand

$\lparen\mathrm{d}f\rparen_{x}(v)=v^{\intercal}y=v^{\intercal}HH^{-1}y=\langle
v,H^{-1}y\rangle=\langle H^{-1}y,v\rangle$

and so $\nabla f(x)=H^{-1}y$ for every $x\in\mathbb{R}^{n}$.

  * •

Consider $\mathbb{R}^{n\times n}$ with the canonical inner product $\langle
A,B\rangle=\operatorname{tr}\lparen A^{\intercal}B\rparen$, and let
$f(A)=\operatorname{tr}(A)$. Since $f$ is linear
$\lparen\mathrm{d}f\rparen_{A}(E)=\operatorname{tr}(E)=\operatorname{tr}\lparen\lparen\mathrm{I}_{n}\rparen^{\intercal}E\rparen=\langle\mathrm{I}_{n},E\rangle$.
Thus, $\nabla f(A)=\mathrm{I}_{n}$ for every $A\in\mathbb{R}^{n\times n}$.

  * •

Let $f\colon\mathbb{R}^{n}\to\mathbb{R}$, and consider the canonical inner
product on $\mathbb{R}^{n}$. The $i$-th coordinate of $\nabla f(x)$ is equal
to $\frac{\partial f}{\partial x_{i}}(x)$.

  * •

Let $f\colon\mathbb{R}^{n}\to\mathbb{R}$, and consider the inner product on
$\mathbb{R}^{n}$ given by $\langle x,y\rangle=x^{\intercal}Hy$ for
$H\in\mathbb{R}^{n\times n}$ spd. Denote by $g_{x}\in\mathbb{R}^{n}$ the
vector with $i$-th coordinate equal to $\frac{\partial f}{\partial
x_{i}}(x)$—i.e., the gradient of $f$ with respect to the canonical inner
product. Then $\nabla f(x)=H^{-1}g_{x}$, while
$\lparen\mathrm{d}f\rparen_{x}(v)=g_{x}^{\intercal}v$ regardless of the inner
product.

######  Remark (Gradient of a composition).

Consider a function $h\colon\mathbb{R}^{m}\to\mathbb{R}$ defined as a
composition $h=f\circ g$ with $g\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ and
$f\colon\mathbb{R}^{n}\to\mathbb{R}$ and fix inner products on
$\mathbb{R}^{m}$ and $\mathbb{R}^{n}$. How do we compute the gradient of $h$
in terms of $g$ and $f$? By the chain rule (Theorem 2.2) and the definition of
the gradient of $f$ we have that for $x\in\mathbb{R}^{m}$, denoting
$y=g(x)\in\mathbb{R}^{n}$

$\lparen\mathrm{d}h\rparen_{x}(v)=\lparen\mathrm{d}f\rparen_{y}[\lparen\mathrm{d}g\rparen_{x}(v)]=\langle\nabla
f(y),\lparen\mathrm{d}g\rparen_{x}(v)\rangle.$

To be able to compute the gradient of $h$ at $x$, we would have to solve for
$v$ on the last equality, sending the linear map
$\lparen\mathrm{d}g\rparen_{x}$ to the left-hand side of the inner product.
This is exactly what the adjoint of a linear map achieves.

######  Definition 4.5.

Let $T\colon V\to W$ be a linear map between real finite-dimensional vector
spaces with inner products $\langle-,-\rangle_{V},\langle-,-\rangle_{W}$. We
define its adjoint as the linear map $T^{\ast}\colon W\to V$ such that

$\langle w,T(v)\rangle_{W}=\langle
T^{\ast}(w),v\rangle_{V}\mathrlap{\qquad\forall v\in V,w\in W.}$

Before giving examples of the adjoint of some linear maps, we formalise the
motivation that led to the definition of the adjoint.

######  Proposition 4.6 (Gradient of a composition).

Let $g\colon\mathbb{R}^{m}\to\mathbb{R}^{n}$ and
$f\colon\mathbb{R}^{n}\to\mathbb{R}$ and fix inner products on
$\mathbb{R}^{m}$ and $\mathbb{R}^{n}$. We have that for every
$x\in\mathbb{R}^{m}$, denoting $y=g(x)\in\mathbb{R}^{n}$,

$\nabla\lparen f\circ
g\rparen(x)=\lparen\mathrm{d}g\rparen^{\ast}_{x}\lparen\nabla f(y)\rparen$

###### Proof.

We finish the computation that we started before

$\mathrm{d}\lparen f\circ
g\rparen_{x}(v)=\lparen\mathrm{d}f\rparen_{g(x)}[\lparen\mathrm{d}g\rparen_{x}(v)]=\langle\nabla
f(y),\lparen\mathrm{d}g\rparen_{x}(v)\rangle_{\mathbb{R}^{n}}=\langle\lparen\mathrm{d}g\rparen^{\ast}_{x}\lparen\nabla
f(y)\rparen,v\rangle_{\mathbb{R}^{m}}.\qed$

######  Example 4.7 (Adjoint of the matrix multiplication).

Consider the linear map of multiplying on the right by a matrix
$R_{X}\colon\mathbb{R}^{m\times n}\to\mathbb{R}^{m\times p}$, $R_{X}(A)=AX$
for a fixed matrix $X\in\mathbb{R}^{n\times p}$. For the canonical inner
products on $\mathbb{R}^{m\times n}$ and $\mathbb{R}^{m\times p}$:

$\langle B,R_{X}(A)\rangle_{\mathbb{R}^{m\times p}}=\operatorname{tr}\lparen
B^{\intercal}AX\rparen=\operatorname{tr}\lparen
XB^{\intercal}A\rparen=\operatorname{tr}\lparen\lparen
BX^{\intercal}\rparen^{\intercal}A\rparen=\langle
BX^{\intercal},A\rangle_{\mathbb{R}^{m\times n}}$

In other words $\lparen
R_{X}\rparen^{\ast}(B)=BX^{\intercal}=R_{X^{\intercal}}(B)$, or simply
$\lparen R_{X}\rparen^{\ast}=R_{X^{\intercal}}$.

An analogous computation gives that, for the left multiplication $L_{X}(A)=XA$
with respect to the canonical inner products, $\lparen
L_{X}\rparen^{\ast}=L_{X^{\intercal}}$.

######  Example 4.8 (The adjoint depends on the choice of inner product).

Consider $R_{X}(A)=AX$ for $A,X\in\mathbb{R}^{n\times n}$ as defined in
Example 4.7 and consider the inner product $\langle
A,B\rangle=\operatorname{tr}\lparen A^{\intercal}HB\rparen$ for a fixed
$H\in\mathbb{R}^{n\times n}$ spd (cf., Example 4.2). We have

$\langle B,R_{X}(A)\rangle=\operatorname{tr}\lparen
B^{\intercal}HAX\rparen=\operatorname{tr}\lparen\lparen
BX^{\intercal}\rparen^{\intercal}HA\rparen=\langle
R_{X^{\intercal}}(B),A\rangle$

so $R^{\ast}_{X}=R_{X^{\intercal}}$, as before. On the other hand, for the
left multiplication $L_{X}(A)=XA$,

$\langle B,L_{X}(A)\rangle=\operatorname{tr}\lparen
B^{\intercal}HXA\rparen=\operatorname{tr}\lparen
B^{\intercal}HXH^{-1}HA\rparen=\langle L_{H^{-1}X^{\intercal}H}(B),A\rangle$

so $L_{X}^{\ast}=L_{H^{-1}X^{\intercal}H}$, where we have used that the
inverse of an spd matrix is symmetric.

We finish this section enumerating two properties that will be particularly
useful in Section 4.3.

######  Proposition 4.9.

Let $S\colon U\to V$ and $T,T_{1},T_{2}\colon V\to W$ be linear maps between
(finite-dimensional real) vector spaces with inner products, then

  * •

The adjoint is linear. Defining $(aT)(v)=aT(v)$ and
$(T_{1}+T_{2})(v)=T_{1}(v)+T_{2}(v)$ for $v\in V$, $a\in\mathbb{R}$, then
$(aT)^{\ast}=aT^{\ast}$ and $\lparen
T_{1}+T_{2}\rparen^{\ast}=T_{1}^{\ast}+T_{2}^{\ast}$.

  * •

The adjoint reverses the order of the composition: $(T\circ
S)^{\ast}=S^{\ast}\circ T^{\ast}$.

###  4.2 The model and definition

As in the case of forward ad, we have a model described by a composition of
functions. The difference is that, in this case, the last map will be a
function mapping the result into the real numbers. In machine learning this is
called the _loss function_.

$\mathbb{R}^{k}$$\mathbb{R}^{m}$$\mathbb{R}^{n}$$\mathbb{R}$$h$$g$$f$ Figure
2: The automatic differentiation model.

######  Definition 4.10.

Backward mode ad for the model represented in Figure 2 with respect to the
canonical inner product on $\mathbb{R}^{k}$ accounts for computing the
gradient

$\nabla\lparen f\circ g\circ h\rparen(x)\mathrlap{\qquad\text{for
}x\in\mathbb{R}^{k}.}$

All the theory explained before about linear maps, differentials, inner
products, gradients, and adjoints comes together to give this particularly
simple definition. Even better, we have laid out the theory in such a way that
we have all the tools to compute this quantity. In Section 3, we saw how to
compute the differential of different maps and how to compose them together.
Using Proposition 4.6 we can compute the gradient of the model in terms of the
gradient of $f$ and the adjoints of the differentials of $g$ and $h$—for
example by choosing the canonical inner product on $\mathbb{R}^{n}$ and
$\mathbb{R}^{m}$.111111It is possible to prove that the result does not depend
on the choice of inner product in the intermediate spaces, the only thing that
changes is the matrix representation of the functions. Finally, we use
Proposition 4.9 to compute the adjoint of the composition as the reversed
composition of the adjoints.121212It should be clear why in some fields in
applied mathematics backpropagation is called the adjoint method. All this
together gives, denoting $y=h(x)$ and $z=g(y)=g(h(x))$,

$\nabla\lparen f\circ g\circ
h\rparen(x)=\lparen\mathrm{d}h\rparen^{\ast}_{x}[\lparen\mathrm{d}g\rparen^{\ast}_{y}\lparen\nabla
f(z)\rparen]\mathrlap{\qquad\forall x\in\mathbb{R}^{k}.}$

Note that the last function being applied, $f$, is the first one that we
differentiate. We then pass its gradient “backwards” to $g$ and then $h$,
hence the name of the method.

###  4.3 Computing backward mode ad

We already computed the adjoint of some linear functions in Example 4.7. We
now show how these formulae together with the properties from Proposition 4.9
are enough to compute the adjoint of the differentials of the maps we
considered in Section 3.

######  Example 4.11 (Adjoint of the powers of a matrix).

In Example 3.7 we showed that for $f(A)=A^{k}$

$\lparen\mathrm{d}f\rparen_{A}(E)=\sum_{i=0}^{k-1}A^{i}EA^{k-i-1}=\sum_{i=0}^{k-1}L_{A^{i}}\lparen
R_{A^{k-i-1}}\lparen E\rparen\rparen\mathrlap{\qquad A,E\in\mathbb{R}^{n\times
n}}$

with $L_{X}(A)=XA$ and $R_{X}(A)=AX$ being the left and right multiplication.
For the canonical inner product on $\mathbb{R}^{n\times n}$, using that the
adjoint is linear (Proposition 4.9), and the formulae for the adjoint of $L$
and $R$ (Example 4.7) we get

$\lparen\mathrm{d}f\rparen^{\ast}_{A}=\Big{\lparen}\sum_{i=0}^{k-1}L_{A^{i}}\circ
R_{A^{k-i-1}}\Big{\rparen}^{\ast}=\sum_{i=0}^{k-1}\lparen L_{A^{i}}\circ
R_{A^{k-i-1}}\rparen^{\ast}=\sum_{i=0}^{k-1}R_{\lparen
A^{\intercal}\rparen^{k-i-1}}\circ L_{\lparen A^{\intercal}\rparen^{i}}.$

or more explicitly

$\lparen\mathrm{d}f\rparen^{\ast}_{A}(E)=\sum_{i=0}^{k-1}\lparen
A^{\intercal}\rparen^{i}E\lparen A^{\intercal}\rparen^{k-i-1}.$

######  Example 4.12 (Adjoint of matrix function).

For a matrix function $f\colon\mathbb{R}^{n\times n}\to\mathbb{R}^{n\times
n}$, $f(A)=\sum_{k=0}^{\infty}c_{k}A^{k}$ (cf., Definition 3.9) and the
canonical inner product on $\mathbb{R}^{n\times n}$ we have that

$\lparen\mathrm{d}f\rparen^{\ast}_{A}=\lparen\mathrm{d}f\rparen_{A^{\intercal}}\mathrlap{\qquad\forall
A\in\mathbb{R}^{n\times n}.}$

This follows from the formula for the differential of $f$ computed in Equation
3 and the properties of the adjoint (Proposition 4.9). As a corollary, we get
that the adjoint of the differential of an analytic function on matrices can
be computed by applying $f$ to a larger function, using the formula in Theorem
3.11.

###  4.4 Exercises

We leave here a number of exercises to help the reader wrapping their head
around the material. If you want to try just two, have a look at Exercises
4.13 and 4.16.

######  Exercise 4.13.

Consider a feed-forward network $F_{A,b}(x)=\ell\lparen\sigma\lparen
Ax+b\rparen\rparen$ for a function $\ell\colon\mathbb{R}^{m}\to\mathbb{R}$,
fixed $x\in\mathbb{R}^{n},A\in\mathbb{R}^{m\times n},b\in\mathbb{R}^{m}$ and
an element-wise function $\sigma$. Show that the gradient with respect to $A$
has rank $1$. In other words, show that if
$\tilde{F}(A)=F_{A,b}(x)=\ell\lparen\sigma\lparen Ax+b\rparen\rparen$ for
fixed $x,b$, $\nabla\tilde{F}(A)=uv^{\intercal}$ for two vectors
$u\in\mathbb{R}^{m}$ and $v\in\mathbb{R}^{n}$

Hint. What is the adjoint of the map $A\mapsto Ax+b$ for fixed $x,b$ evaluated
on a vector $g\in\mathbb{R}^{m}$?

The next two exercises look more difficult, but they follow via the same
argument as the one above.

######  Exercise 4.14.

Same as above, but with a feed-forward network of depth $d$.

######  Exercise 4.15.

Same as above, but with a feed-forward network of depth $d$ and in the
stochastic setting, where we define the total loss as
$\widehat{F}(A)=\frac{1}{r}\sum_{i=1}^{r}F_{A,b}(x_{i})$ for input vectors
$\\{x_{i}\\}_{i=1}^{r}$.

######  Exercise 4.16.

Implement the gradient for the two layer feedforward network in LABEL:lst:ffn.

Hint. Name more intermediate variables in FFN.forward to be able to store
them.

Hint. Reverse the order of the arguments in FFN.forward for it to be easier to
debug. The gradient with respect to b1 is not going to be correct if the
gradient with respect to b2 is not correct.

⬇

1import torch

2

3class FFN(torch.autograd.Function):

4 staticmethod def forward(ctx, x, y, A1, b1, A2, b2): x = (A1 x +
b1).sigmoid()

5 x = (A2 x + b2).sigmoid() loss = (x - y).pow(2).sum()
ctx.save_for_backward(...) return loss staticmethod

6 def backward(ctx, g_l):

7 t1, t2, ... = ctx.saved_tensors

8 ...

9 return None, None, g_A1, g_b1, g_A2, g_b2

10

11class Model(torch.nn.Module):

12 def __init__(self, in_features, hidden_features, out_features):

13 super().__init__()

14 def make_param(*size):

15 return torch.nn.Parameter(torch.empty(*size, dtype=torch.double))

16 self.register_parameter("A1", make_param(hidden_features, in_features))

17 self.register_parameter("b1", make_param(hidden_features))

18 self.register_parameter("A2", make_param(out_features, hidden_features))

19 self.register_parameter("b2", make_param(out_features))

20 torch.nn.init.xavier_normal_(self.A1)

21 torch.nn.init.xavier_normal_(self.A2)

22

23 def forward(self, x, y):

24 return FFN.apply(x, y, self.A1, self.b1, self.A2, self.b2)

25

26x = torch.rand(32, dtype=torch.double) # Batch size 1

27y = torch.rand(8, dtype=torch.double)

28model = Model(32, 16, 8)

29

30args = (x, y, model.A1, model.b1, model.A2, model.b2)

31torch.autograd.gradcheck(FFN.apply, args, atol=0.01)

Listing 1: Modify this PyTorch 1.9 code so that autograd passes.

######  Exercise 4.17.

Generalise your code in FFN.backwards to handle batches of arbitrary size.
Then, use the code you have implemented to fit mnist and feel good about
yourself.

######  Exercise 4.18.

Compute the gradient for a recurrent neural network (rnn) with respect to the
recurrent kernel on PyTorch or just on paper.

######  Exercise 4.19.

Find where the adjoint for matrix_exp is implemented in PyTorch and make sure
you understand its code. Hint. Look for the function matrix_exp_backward.

##  5 Complex Maps

###  5.1 Forward mode ad

When we derived the formulae for forward mode ad, they all followed from the
definition of differential (Definition 2.1) and the formula differential of a
linear map (Proposition 3.2). As such, if we can generalise these two to
complex maps, we should be able to generalise all the forward mode ad to
complex numbers.

In order to do this, we recall the point that we made in Example 1.3,
$\mathbb{C}^{n}$ is a real vector space of dimension $2n$. This means that for
$a\in\mathbb{R}$,

$u+v\in\mathbb{C}^{n}\qquad au\in\mathbb{C}^{n}\mathrlap{\qquad\qquad\forall
u,v\in\mathbb{C}^{n}.}$

Furthermore, it means that these operations—again, with $a\in\mathbb{R}$, not
$a\in\mathbb{C}$—satisfy all the axioms of a real vector space in Definition
1.1.

This real vector space structure treats the $n$ real components and $n$
imaginary components as independent, as if they were two parts of a vector of
size $2n$ in $\mathbb{R}^{2n}=\mathbb{R}^{n}\times\mathbb{R}^{n}$. As such,
the norm of a vector in $\mathbb{C}^{n}$ as a real vector space is given by

$\lVert
v\rVert^{2}_{\mathbb{C}^{n}}=\sum_{k=1}^{n}a_{k}^{2}+b_{k}^{2}\mathrlap{\qquad\text{for
}v=\lparen a_{1}+ib_{1},\dots,a_{n}+ib_{n}\rparen.}$

Using this norm, we can extend the definition of a differential of a real map
to complex maps.

######  Definition 5.1.

A map $f\colon\mathbb{C}^{m}\to\mathbb{C}^{n}$ is real differentiable at a
point $x\in\mathbb{C}^{m}$ if there exists a map
$\lparen\mathrm{d}f\rparen_{x}\colon\mathbb{C}^{m}\to\mathbb{C}^{n}$ which is
linear over the real numbers (cf., Definition 1.10) such that

$\lim_{h\to
0_{\mathbb{C}^{m}}}\frac{f(x+h)-f(x)-\lparen\mathrm{d}f\rparen_{x}(h)}{\lVert
h\rVert_{\mathbb{C}^{m}}}=0_{\mathbb{C}^{n}}.$

######  Remark (Real differentiable vs. complex differentiable).

Here we have defined the real differential as an $\mathbb{R}$-linear map, that
is, a map such that
$\lparen\mathrm{d}f\rparen_{x}(av)=a\lparen\mathrm{d}f\rparen_{x}(v)$ for
$a\in\mathbb{R}$. If we require the differential to be
$\mathbb{C}$-linear—that is,
$\lparen\mathrm{d}f\rparen_{x}(av)=a\lparen\mathrm{d}f\rparen_{x}(v)$ for
$a\in\mathbb{C}$—we get the definition of a complex differentiable map, often
called holomorphic map.

It should be clear that, if a complex map is complex differentiable, it is
also real differentiable, but the opposite is not true. Consider for example
$f(z)=\overline{z}$ for $z\in\mathbb{C}$. We have that
$f(az)=\overline{a}f(z)$ for $a\in\mathbb{C}$, so it is not
$\mathbb{C}$-linear, but it is $\mathbb{R}$-linear as $f(az)=af(z)$ for
$a\in\mathbb{R}$. Luckily, we will not need to use holomorphic maps, as real
differentiable maps will be enough to compute differentials and gradients.

The chain rule (Theorem 2.2) and the Leibnitz rule (Proposition 2.3) also hold
verbatim for real differentiable maps. We also have the following equivalent
to Proposition 3.2.

######  Proposition 5.2 (Differential of a linear map).

Let $T\colon\mathbb{C}^{m}\to\mathbb{C}^{n}$ be an $\mathbb{R}$-linear map, we
have that

$\lparen\mathrm{d}T\rparen_{x}(v)=T(v)\mathrlap{\qquad\forall
x,v\in\mathbb{C}^{m}.}$

###### Proof.

The proof is the same as in the real case. ∎

Having this, we can compute the differential of many maps, as we did in the
real case.

######  Example 5.3.

We compute the differential of some linear maps from $\mathbb{C}^{m}$ to
$\mathbb{C}^{n}$ or to $\mathbb{R}^{n}\subseteq\mathbb{C}^{n}$.

  * •

Let $f\colon\mathbb{C}^{n}\to\mathbb{R}^{n}$, $f(x)=\operatorname{Im}\lparen
x\rparen$ be the imaginary part of a vector. Since

$\operatorname{Im}(x+y)=\operatorname{Im}(x)+\operatorname{Im}(y)\qquad\operatorname{Im}\lparen
ax\rparen=a\operatorname{Im}\lparen x\rparen\mathrlap{\qquad\forall
x,y\in\mathbb{C}^{n},a\in\mathbb{R},}$

$f$ is $\mathbb{R}$-linear and
$\lparen\mathrm{d}f\rparen_{x}(v)=\operatorname{Im}\lparen v\rparen$ for
$v\in\mathbb{C}^{n}$. Note that $f$ is not $\mathbb{C}$-linear.

  * •

Analogously, if $f\colon\mathbb{C}^{n}\to\mathbb{R}^{n}$,
$f(x)=\operatorname{Re}\lparen x\rparen$, $f$ is $\mathbb{R}$-linear and
$\lparen\mathrm{d}f\rparen_{x}(v)=\operatorname{Re}\lparen v\rparen$ for
$v\in\mathbb{C}^{n}$.

  * •

Let $L_{X}(A)=XA$ for $X\in\mathbb{C}^{m\times n}$, $A\in\mathbb{C}^{n\times
k}$. Since $L_{X}$ is $\mathbb{C}$-linear, it is in particular
$\mathbb{R}$-linear, so $\lparen\mathrm{d}L_{X}\rparen_{A}(E)=XE$ for
$E\in\mathbb{C}^{n\times k}$.

  * •

Let $f_{X}(A)=A^{\textup{H}}X$ for $X\in\mathbb{C}^{m\times n}$,
$A\in\mathbb{C}^{k\times n}$ where $A^{\textup{H}}=\overline{A}^{\intercal}$.
Note that this is not a $\mathbb{C}$-linear map as
$f_{X}(cA)=\overline{c}f_{X}(A)$ for $c\in\mathbb{C}$, but it is an
$\mathbb{R}$-linear map, and as such,
$\lparen\mathrm{d}f_{X}\rparen_{A}(E)=E^{\textup{H}}X$.

  * •

The formulae for the differential of the powers of a matrix (Example 3.7),
inverse of a matrix (Example 3.8) and differential of a matrix function
(Theorem 3.11) are also valid for complex matrices.

These examples show that formulae for forward ad for complex maps are the same
as their real counterparts, as the basic formulae (Proposition 5.2, chain
rule, and Leibnitz rule) are the same.

###  5.2 Backward mode ad

For backward mode ad, all we need is a real inner product (cf., Definition
4.1). To do that all we need to do is to consider $\mathbb{C}^{n}$ as a real
vector space, as we did in the previous section.

######  Proposition 5.4.

The canonical real inner product on $\mathbb{C}^{n}$ as a real vector space
for $x,y\in\mathbb{C}^{n}$ can be written as

$\langle
x,y\rangle_{\mathbb{C}^{n}}=\sum_{k=1}^{n}a_{k}c_{k}+b_{k}d_{k}=\operatorname{Re}x^{\textup{H}}y\mathrlap{\qquad\text{for}\quad\begin{matrix}x=\sum
a_{k}+ib_{k}\\\ y=\sum c_{k}+id_{k}\end{matrix}.}$

The canonical real inner product on $\mathbb{C}^{m\times n}$ as a real vector
space for $X,Y\in\mathbb{C}^{m\times n}$ can be written as

$\langle X,Y\rangle_{\mathbb{C}^{m\times
n}}=\sum_{j=1,k=1}^{n}A_{jk}C_{jk}+B_{jk}D_{jk}=\operatorname{Re}\operatorname{tr}\lparen
X^{\textup{H}}Y\rparen\mathrlap{\qquad\text{for}\quad\begin{matrix}X=A+iB\\\
Y=C+iD\end{matrix}.}$

###### Proof.

Note that the first equality in the vector case comes the definition of the
canonical real inner product on $\mathbb{C}^{n}$, which is just the inner
product on $\mathbb{R}^{2n}$ (cf., Example 4.2). Same happens for the
$\mathbb{C}^{m\times n}$ case.

We prove this proposition for the matrix case, as vectors can be seen as the
case $\mathbb{C}^{n\times 1}=\mathbb{C}^{n}$.

We start by rewriting the left-hand side in a coordinate-free way

$\sum_{j,k}A_{jk}C_{jk}+B_{jk}D_{jk}=\langle A,C\rangle_{\mathbb{R}^{m\times
n}}+\langle B,D\rangle_{\mathbb{R}^{m\times n}}.$

Thus, we just need to prove that $\operatorname{Re}\operatorname{tr}\lparen
X^{\textup{H}}Y\rparen=\langle A,C\rangle_{\mathbb{R}^{m\times n}}+\langle
B,D\rangle_{\mathbb{R}^{m\times n}}$, but this is direct as

$\operatorname{tr}\lparen
X^{\textup{H}}Y\rparen=\operatorname{tr}\lparen\lparen
A^{\intercal}-iB^{\intercal}\rparen\lparen
C+iD\rparen\rparen=\underbrace{\operatorname{tr}\lparen
A^{\intercal}C\rparen+\operatorname{tr}\lparen
B^{\intercal}D\rparen}_{\text{real
part}}+i\underbrace{[\operatorname{tr}\lparen
A^{\intercal}D\rparen-\operatorname{tr}\lparen
C^{\intercal}B\rparen]}_{\text{imaginary part}}.$

Note that since $\langle A,C\rangle_{\mathbb{R}^{m\times n}}+\langle
B,D\rangle_{\mathbb{R}^{m\times n}}$ is a (real) inner product,131313Note that
in mathematics we also find complex inner products. These are sesquilinear
maps rather than bilinear. These products are more general than the real inner
products, as their real part is always a real inner product, while their
imaginary part is a non-degenerate symplectic (i.e., skew-symmetric) bilinear
form. Luckily, we do not need these to compute gradients. so is
$\operatorname{Re}\operatorname{tr}\lparen X^{\textup{H}}Y\rparen$. In other
words, it is a symmetric positive definite (real) bilinear map (Definition
4.1). ∎

All this proposition says is that $\operatorname{Re}\operatorname{tr}\lparen
X^{\textup{H}}Y\rparen$ is a convenient way to write the canonical real inner
product on $\mathbb{C}^{m\times n}$. Now, since $\mathbb{C}^{m\times n}$ is a
real vector space, and we have a real inner product on it, all the definitions
and general results in Section 4 translate to this setting. Note that the
gradients are just defined for functions with values in $\mathbb{R}$, not
$\mathbb{C}$, while the adjoints are defined for arbitrary maps.

######  Example 5.5.

Consider the canonical real inner products in each of the spaces

  * •

Let $f\colon\mathbb{C}^{n}\to\mathbb{R}$,
$f(x)=\operatorname{Re}g^{\intercal}x$ for a fixed $g\in\mathbb{C}^{n}$. Since
$f$ is $\mathbb{R}$-linear
$\lparen\mathrm{d}f\rparen_{x}(v)=\operatorname{Re}g^{\intercal}x=\langle\overline{g},v\rangle_{\mathbb{C}^{n}}$,
and by definition of a gradient, $\nabla f(x)=\overline{g}$ for every
$x\in\mathbb{C}^{n}$.

  * •

For $L_{X}(A)=XA$, $X\in\mathbb{C}^{m\times n},A\in\mathbb{C}^{n\times k}$, we
have that $L_{X}^{\ast}\colon\mathbb{C}^{n\times k}\to\mathbb{C}^{m\times n}$
is given by $L_{X}^{\ast}(E)=X^{\textup{H}}E$, so that
$L_{X}^{\ast}=L_{X^{\textup{H}}}$.141414This shows why some areas of
mathematics abuse the notation and write $X^{\textup{H}}$ as $X^{\ast}$.

  * •

Similarly, if $R_{X}(A)=AX$ for complex $A,X$,
$R_{X}^{\ast}=R_{X^{\textup{H}}}$.

  * •

Example 4.12 translates to
$\lparen\mathrm{d}f\rparen^{\ast}_{A}=\lparen\mathrm{d}f\rparen_{A^{\textup{H}}}$
for an analytic function $f\colon\mathbb{C}^{n\times n}\to\mathbb{C}^{n\times
n}$ and $A\in\mathbb{C}^{n\times n}$.

In this case, the formulae are almost the same, but it tends to happen that
when a matrix or a vector is transposed in the real case, it is transposed and
conjugated in the complex case.

## References

  * [1] Tosio Kato  “Perturbation theory for linear operators” Reprint of the 1980 edition  Springer, 1995
  * [2] Roy Mathias  “A chain rule for matrix functions and applications”  In _SIAM Journal on Matrix Analysis and Applications_ 17.3, 1996, pp. 610–620  DOI: 10.1137/S0895479895283409


# Binarized Weight Error Networks With a Transition Regularization Term

###### Abstract

This paper proposes a novel binarized weight network (BT) for a resource-
efficient neural structure. The proposed model estimates a binary
representation of weights by taking into account the approximation error with
an additional term. This model increases representation capacity and
stability, particularly for shallow networks, while the computation load is
theoretically reduced. In addition, a novel regularization term is introduced
that is suitable for all threshold-based binary precision networks. This term
penalizes the trainable parameters that are far from the thresholds at which
binary transitions occur. This step promotes a swift modification for binary-
precision responses at train time. The experimental results are carried out
for two sets of tasks: visual classification and visual inverse problems.
Benchmarks for Cifar10, SVHN, Fashion, ImageNet2012, Set5, Set14, Urban and
BSD100 datasets show that our method outperforms all counterparts with binary
precision.

Index Terms—  Binarized Neural Networks, Efficient Neural Networks, Binary
Networks, Ternary Networks

##  1 Introduction

Over the past decade, neural networks have shown remarkable performance in
various fields [1, 2]. In particular, this structure transforms high-
dimensional data into low-dimensional latent space by taking advantage of non-
linear transformations and vast parameters. However, in addition to its
popularity, this competence also multiplies the computation load. Increasing
computational and storage demands, therefore, inevitably reduce its
applicability on portable devices.

In recent years, several critical studies aim to reduce the need for these
limitations with a slight decrease in performance than full precision
counterparts. In general, compressing full precision parameters into binarized
ones plays a prominent role because of the simplicity and efficiency [3, 4, 5,
6, 7, 8]. Mathematically, this type of methods minimizes the approximation
error between full and binary weights as follows:

$\displaystyle
J(\alpha,\mathbf{W}^{b})=\underset{\alpha,\mathbf{W}^{b}}{\mathrm{argmin}}{||\mathbf{W}-\alpha\mathbf{W}^{b}||}$
(1)

Here, $\mathbf{W}$ and $\mathbf{W}^{b}$ indicate full and binary precision
trainable parameters, respectively. Moreover, $\alpha$ denotes a scale factor
where its optimal solution is equal to mean of the absolute value of full
precision weights $\mathbb{E}[|\mathbf{W}|]$. In the end, this approximation
leads to $\mathbf{W}\approx\alpha\mathbf{W}^{b}$ with a significant
quantization error. Observe that this quantization error is even worse for
shallow networks since the allowable range for parameters is large. This issue
eventually weakens the stability of solutions. Additionally, error propagation
for binary precision parameters can be severe due to limited operations in
back-propagation steps (i.e., multiplication with binarized weights instead of
double precision). This limitation causes additional problems for converging
optimal solutions.

In this work, the quantization error in Eq. 1 is reduced by an iterative
quantization step. In particular, full weights are iteratively quantized in
several steps, and a binary representation is computed with a negligible
increase in the computational load. As noticed, this schema intuitively shares
an assumption similar to hierarchical clustering, which tries to create a
hierarchy of clusters. Similarly, this step helps improve the stability of
binarized weights because of the better partitions. Indeed, our bitwise
operations induce a fast inference similar to other quantized neural network
by consuming less power and memory (appx. $16\times$) on the devices. However,
performance can slightly reduce due to the bit-length compared to full-
precision counterparts.

In addition, a regularization term is introduced that is feasible for all
threshold-based binarized networks. This term penalizes parameters that are
far from the thresholds at which binary transitions occur. Hence the
parameters are concentrated close to the binary transitions. This term favors
changing weights quickly in binary precision to observe their effects to the
performance and avoid overfitting for parameters.

The rest of the paper is formed as follows. First, the related work for
binarized weight networks is summarized. Later, the concept of binarized
weight networks and our contributions are described. Lastly, we explain the
test results and final remarks.

##  2 Related Work

Although there is a large body of studies related to the compression of
computation and memory requirements for deep neural networks [9], binary-
precision networks play a crucial role due to their efficiency and simplicity
[3, 4].

In baseline methods, all trainable parameters are binarized with a quantizer.
Later, the backpropagation step is also realized on these binary-precision
weights. [6] shows that the performance for binarized-parameters can be
further improved by adding a single scale factor. Moreover, activation and
parameters are also converted to binary-precision to achieve a binarized
inference [6, 7, 8].

Note that these methods limit weights to {-1, 1} (with some scale factor) and
decrease the sparsity property of neural networks. [10, 11] explain that the
use of zero weights {-1,0,1} contributes to performance by promoting sparsity
in their solutions. In this way, robust representations can be computed. [11]
also shows that multiple scale factors obtain better results compared to the
cases where a single scale factor is utilized. Recently, SBNN prunes some of
NN layers in a binarized network and refines the network with a softsign
function [12]. The most similar work [13, 14] to our method is to quantize
continuous parameters to a set of integer numbers with multiple binary
quantizers. However, this can lead to overly quantized representations where
instability occurs. We recommend [15, 16] to read for the survey in binary
neural networks.

Lastly, the applications of this terminology are proliferated on different
tasks [17, 18].

##  3 Binarized Weight Error Networks

Fig. 1: The flow of conventional (first) and our method (second) for the
binarized network. Differently, our method quantizes full-precision parameter
$W$ by counting approximation error $E$.

The motivation of binarized neural networks is to minimize the error between
full-precision $\mathbf{W}$ and binary-precision $\mathbf{W}^{b}$ weights as
formulated in Eq. 1. Generally, two binarized network structures come into
prominence in the literature.

In binary-weight networks, an approximation problem is solved by assigning a
binary weight as:

$\displaystyle\mathbf{W}^{b}=\begin{cases}+1,&\text{if }\mathbf{W}\geq 0\\\
-1,&\text{if }\mathbf{W}<0\end{cases}$ (2)

Implicitly, this overall solution corresponds to the operation
$\mathbf{W}^{b}=\text{sign}(\mathbf{W})$. Note that $\mathbf{W}^{b}$ takes
only $\\{-1,1\\}$, and the responses of weights are always active since there
is no zero operation. This assumption conflicts with the conclusions reached
about filter characteristics, as explained in [19] (i.e., Gabor-type filtering
and sparsity).

On the contrary, ternary-weight networks compute a representation with a
threshold function in Eq. 3, since there is no deterministic solution for
binary-weights as claimed in [20]:

$\displaystyle\mathbf{W}^{b}=\begin{cases}+1,&\text{if }\mathbf{W}>\Delta\\\
0,&\text{if }|\mathbf{W}|<\Delta\\\ -1,&\text{if
}\mathbf{W}<-\Delta\end{cases}$ (3)

Here, this assumption is exploited to approximate the threshold value $\Delta$
as $0.66\cdot\mathbb{E}[|\mathbf{W}|]$. This model introduces more
representation capacity because of the sparsity, while stability is also
increased.

###  3.1 Binarized Weight Error Networks

As excepted, quantization errors for binary weights are moderately high
compared to ternary weights. However, we believe there is room to take
additional steps to reduce this error.

In particular, the idea emphasized in our work is based on the fact that an
iterative approach can be more suitable for binarization process to increase
the capacity of representations than the quantization of full-precision
parameters in a single step. For this purpose, the problem definition in Eq. 1
is rewritten as:

$\displaystyle
J_{i}(\alpha,\mathbf{W}_{i}^{b})=\underset{\alpha,\mathbf{W}_{i}^{b}}{\mathrm{argmin}}{||\mathbf{E}_{i}-\alpha_{i}\mathbf{W}_{i}^{b}||}$
(4)

Here, $\mathbf{E}_{i}$ indicates the quantization error that occurs in the
previous iteration step $i-1$. This is equal to $\mathbf{E}_{0}=\mathbf{W}$ in
the first iteration. When the total iteration is set to 2, the final
representation can be approximated as:

$\displaystyle\mathbf{W}^{b}=\alpha_{opt}\cdot(\mathbf{W}_{0}^{b}+\mathbf{W}_{1}^{b})$
(5)

Here, $\alpha_{opt}$ indicates the optimum scale value for this conversion
that needs to be recalculated. The flows of conventional and our method for
binarized networks are illustrated in Fig. 1.

An observation for Eq. 5 is that the use of two binary-weight models for both
$\mathbf{W}_{0}^{b}$ and $\mathbf{W}_{1}^{b}$ implicitly reproduces a ternary-
weight model. This result is critical because it confirms the mathematical
proof of our method. Note that this step also shows that binary-representation
can be improved further with simple replacements in the formulation.

Second, we compute the optimum scale value $\alpha_{opt}$ as:

$\displaystyle\alpha_{opt}=0.75\cdot\alpha_{1}-0.25\cdot\alpha_{2}$ (6)

where $\alpha_{1}$ and $\alpha_{2}$ are equal to
$\mathbb{E}[|\mathbf{E_{0}}|]$ and $\mathbb{E}[|\mathbf{E_{1}}|]$,
respectively. If binary and ternary models are sequentially applied by using
Eq. 4 (that we call BT), the representation will take one of the values
$\\{-2,-1,0,1,2\\}$. Hence, an ideal scale value must be between the lower
bound (i.e., $(\alpha_{1}-\alpha_{2})$) and upper bound (i.e.,
$0.5\cdot(\alpha_{1}+\alpha_{2})$, $0.5$ is because of the maximum binary
range which is $2$). Intuitively, the mean value of these scales yields the
optimum scale value as given in Eq. 6.

###  3.2 Transition Regularization

As explained in [10], even if full-precision parameters are sampled from
either a normal or uniform distribution, the allowable parameter range can
still be large. A high approximation error is therefore expected. Furthermore,
parameter updates with binary-precision are problematic as accurate error
propagation is nearly impossible and greatly overfits the parameters.

We propose a novel regularization term that forces parameters to concentrate
around binary transitions and penalizes those that are far from those points.
For clarity, transition points correspond to $0$ for binary-weight models,
$\\{-\Delta,\Delta\\}$ for ternary-weight models and $\alpha_{opt}$s for our
models. Hence, the effects of parameter changes for binary-precision can be
learned more quickly with this regularization term at train time.

To favor these changes, a simple and effective trick is used. Mathematically,
full-precision parameters $\mathbf{W}$ are corrupted with a random noise
$\mathbf{\tilde{W}}=\mathbf{W}+0.1\cdot\mathcal{N}(0,\alpha)$. Later, these
noisy weights are also quantized $\mathbf{\tilde{W}}^{b}$ similarly with Eq.
1. Note that at the boundaries of transition points, binarized-representations
for corrupted parameters differ from the original versions:

$\displaystyle\mathcal{L}=\mathcal{L}_{pb}-\lambda||\mathbf{W}^{b}-\mathbf{\tilde{W}}^{b}||_{1}$
(7)

Fig. 2: Visual samples from Fashion (Upper-Left), ImageNet2012 (Upper-Right),
Cifar10 (Lower-Left) and SVHN (Lower-Right) datasets.

Here, $\lambda$ is a coefficient that scales the contribution of
regularization term to the overall loss function with a task-specific loss
function $\mathcal{L}_{pb}$. Note that this coefficient is empirically set to
$0.1$. From this equation, the absolute distance between actual and corrupted
binary-precision weights is expected to be maximized. Hence, full-parameters
are pushed close to these transition points. This notion is similar to Dropout
[21] regularization where parameter space is intentionally corrupted.

In the experiments, we will show that this term is particularly useful for
classification tasks than inverse problems. The main reason is that as
explained in [22], regression tasks can lead to holes due to continuous
solution spaces, so that sparsity cannot be maintained correctly in such
problems.

###  3.3 Implementation Detail

For the implementation detail111https://github.com/savasozkan/BofT, each model
is trained with Adam stochastic optimizer at train time. Mini-batch size is
set 128, and learning rate scaling is adopted. In the classification task,
VGG-6 (K) with batch normalization is used for acceleration and parameter
regularization (i.e., 2×(K-C3) + MP2 +2×(2×K-C3) + MP2 + 2×(4×K-C3) + MP2 +
8×K-FC + Softmax and K is the number of filters). Here, ”K-C3” is convolution
layer with K filters of a $3\times 3$ kernel. Also, ”MP2” and ”K-FC” denote a
max pooling with stride $2$ and a fully-connected layer with K filters
respectively. Similarly, espcn architecture with skip-connection is utilized
for super-resolution and denoising tasks.

##  4 Experiments

Experiments are conducted on two sets of tasks as visual classification and
visual inverse problems to show the superiority of our method. The performance
of our model (BT) is compared with binary-weight (BWN) [6], ternary-weight
(TWN) [10], trained ternary quantization (TTQ) [11] and LQ-Nets (LQ) [14]
models that are selected as baselines. Note that all results are obtained by
ours running on their publicly available codes. For classification tasks, the
mean average precision (MAP)@1 metric is used, while PSNR scores are
calculated for super-resolution and denoising tasks. For clarity, two problem
sets are explained in two individual sections.

###  4.1 Visual Classification

Table 1: Comparisons with baselines on ImageNet2012 (INet), Fashion, Cifar10
and SVHN datasets. MAP@1 scores are reported. Best results for binary
precision are bold.

| Arch. | Full | BWN | TWN | TTQ | LQ | BT(ours)
---|---|---|---|---|---|---|---
INet | K=64 | 42.81 | 35.30 | 38.25 | 37.76 | 39.01 | 40.21
K=128 | 44.13 | 40.08 | 42.59 | 41.18 | 42.65 | 42.80
Fashion | K=16 | 93.48 | 92.51 | 93.20 | 93.32 | 93.34 | 93.41
K=32 | 94.16 | 93.28 | 93.97 | 93.91 | 94.04 | 94.03
K=64 | 94.56 | 94.16 | 94.30 | 94.18 | 94.29 | 94.55
K=128 | 94.58 | 94.39 | 94.32 | 94.34 | 94.45 | 94.51
Cifar10 | K=16 | 87.62 | 78.70 | 82.94 | 83.83 | 84.67 | 84.61
K=32 | 90.64 | 86.45 | 88.37 | 88.53 | 88.93 | 89.49
K=64 | 92.89 | 90.35 | 91.67 | 91.02 | 91.90 | 92.01
K=128 | 93.58 | 92.06 | 92.39 | 92.23 | 92.85 | 92.92
SVHN | K=16 | 94.14 | 93.15 | 93.44 | 93.32 | 94.04 | 94.09
K=32 | 94.93 | 94.32 | 94.84 | 94.69 | 94.83 | 94.88
K=64 | 95.43 | 95.10 | 95.22 | 95.29 | 95.31 | 95.37
K=128 | 95.76 | 95.17 | 95.28 | 95.65 | 95.84 | 95.74

Experimental results are reported on Fashion, ImageNet2012 (64×64), Cifar10
and SVHN datasets. Fashion, Cifar10 and SVHN datasets have 10 object/character
classes, while ImageNet2012 has 1000 object classes. Visual samples for these
datasets are illustrated in Fig. 2. Random crop, padding and cutout policies
are exploited as data augmentation. Moreover, softmax cross-entropy loss is
used as the task-specific loss function $\mathcal{L}_{pb}$. Note that all
trainable parameters, excluding fully-linear layers, are binarized in our
paper.

First, we compare our method with baseline models, and results are illustrated
in Table 1. As observed, our method outperforms all its binary-precision
counterparts. This result validates that it improves the quantization
performance by considering the approximation error. Especially for
ImageNet2012 (INet), the increase is significant, and the proposed model
yields close results with the full-precision architecture. Moreover, the
performance saturates when the number of parameters increases for all-
binarized models as expected.

Table 2: Impact of transition regularization term on ImageNet2012 (INet),
Fashion, Cifar10 and SVHN datasets. MAP@1 scores are reported.

| Arch. | Full | BWN | TWN | BT(ours)
---|---|---|---|---|---
INet | K=64 | 42.81 | 40.29 | 38.45 | 40.97
K=128 | 44.13 | 35.60 | 42.82 | 43.22
Fashion | K=16 | 93.48 | 93.11 | 93.36 | 93.45
K=32 | 94.16 | 93.76 | 94.08 | 94.12
K=16 | 94.56 | 94.34 | 94.42 | 94.55
K=32 | 94.58 | 94.45 | 94.50 | 94.55
Cifar10 | K=16 | 87.62 | 79.42 | 83.42 | 85.03
K=32 | 90.64 | 86.78 | 88.52 | 90.12
K=16 | 92.89 | 90.57 | 91.92 | 92.12
K=32 | 93.58 | 92.26 | 92.57 | 92.99
SVHN | K=16 | 94.14 | 93.67 | 93.82 | 94.12
K=32 | 94.93 | 94.52 | 94.90 | 94.88
K=16 | 95.43 | 95.28 | 95.34 | 95.40
K=32 | 95.76 | 95.48 | 95.54 | 95.72

Finally, the impact of transition regularization term is reported in Table 2.
Results confirm that this term improves all threshold-based binarized weight
models. In particular, these improvements reach up to $1\%$ for the shallow
networks (i.e., K is lower) compared to results in Table 1. The reason is that
when the number of parameters increases, parameters are automatically
concentrated in a narrow range, which has a positive impact on the
quantization.

###  4.2 Visual Inverse Problems

Experiments are carried out on two well-known inverse problem tasks: visual
super-resolution and visual denoise. At train time, images from the COCO
val2017 dataset are used for both tasks. For high resolution, these images are
rescaled by 2, 3, and 4 to simulate changes in resolution. At test time, the
data sets Set5, Set14, Urban, and BSD100 are used. Additionally, different
sigma values (i.e., with a mean of zero) are used to attach images for the
denoise task. Tests are also performed on Set5 and Set14 records.

Table 3: Super-resolution results on different datasets with various scale
factors. PSNR scores are reported.

| Dataset | Full | Bicubic | BWN | TWN | BT(ours)
---|---|---|---|---|---|---
Scale 2 | Set5 | 36.61 | 33.68 | 36.00 | 36.03 | 36.63
Set14 | 32.58 | 30.24 | 32.14 | 32.13 | 32.42
Urban | 29.51 | 26.61 | 28.91 | 28.83 | 29.31
BSD | 31.48 | 29.48 | 31.08 | 31.13 | 31.28
Scale 3 | Set5 | 33.11 | 30.43 | 32.61 | 32.67 | 32.95
Set14 | 29.46 | 27.54 | 29.13 | 29.12 | 29.32
BSD | 28.38 | 27.14 | 28.21 | 28.25 | 28.34
Scale 4 | Set5 | 30.82 | 28.45 | 30.35 | 30.41 | 30.72
Set14 | 27.70 | 26.01 | 27.38 | 27.43 | 27.66
Urban | 24.81 | 23.12 | 24.38 | 24.49 | 26.69
BSD | 27.01 | 25.91 | 26.73 | 26.81 | 26.97

First, an input image is converted to YUV color space for both tasks, and Y
channel is used for both train and test phases. Also, our NN architecture is
fully convolutional. It takes an input image and applies consecutive
convolution layers with residual-based learning, as explained in [23]. At the
output, it is upsampled with a subpixel layer [24] for super-resolution (i.e.,
(64-C3) + 2×(64-C3) + 2×(64-C3) + SubPixel + (1-C3)). Note that the first and
last layers are full-precision networks, while the rest is binarized.

Experimental results for super-resolution are reported in Table 3. Our method
outperforms all other baselines. In particular, for small scale factors,
negligible performance drops can be observed compared to full-precision
networks. Similarly, the reason is that the allowable parameter range is
expected to be small for small scale factors, and the sparsity is decreased.

Table 4: Denoising results on different datasets with various sigma values.
PSNR scores are reported.

| Dataset | Noisy | Full | BWN | TWN | BT(ours)
---|---|---|---|---|---|---
0.05 | Set5 | 26.20 | 33.58 | 33.06 | 33.03 | 33.59
Set14 | 26.14 | 32.41 | 31.62 | 31.64 | 32.37
0.1 | Set5 | 20.26 | 29.73 | 29.63 | 29.42 | 29.71
Set14 | 20.17 | 28.89 | 28.01 | 28.07 | 28.91
0.2 | Set5 | 14.65 | 26.41 | 26.05 | 25.53 | 26.55
Set14 | 14.37 | 25.19 | 24.60 | 24.29 | 25.04

Denoising results are also shown in Table 4. From the results, the proposed
method achieves slightly better performance even for full-precision networks.
Residual-based learning eventually promotes robustness with binarized networks
since parameter range (i.e., additive noise varies $[0.05,0.2]$) is heavily
reduced. Hence, binarized weights can converge better solutions compared to
other tasks.

As indicated, the proposed regularization term yields similar or slightly
worse performance for inverse problems. The main reason is that the learning
process may lead to holes in continuous feature space partitioning.

##  5 Conclusion

This paper proposes a novel model for binarized weight networks that achieves
state-of-the-art performance compared to its counterparts. This model
estimates a binary representation by taking into account the approximation
error with an additional term. In addition, a novel regularization term is
presented. This term is useful for threshold-based binarized models in
particular. This term favors parameters that focus on binary transitions so
that overfitting the parameters for binarized weights is avoided. In the
experiments, we prove that this term is significant for classification-related
tasks and enhances the performance of binarized methods. Experiments with
different benchmarks for different tasks are conducted to confirm the
superiority of our contributions.

## References

  * [1] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun,  “Deep residual learning for image recognition,”  Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 770–778, 2016.
  * [2] Karen Simonyan and Andrew Zisserman,  “Very deep convolutional networks for large-scale image recognition,”  arXiv preprint arXiv:1409.1556, 2014.
  * [3] Daniel Soudry, Itay Hubara, and Ron Meir,  “Expectation backpropagation: Parameter-free training of multilayer neural networks with continuous or discrete weights,”  Advances in Neural Information Processing Systems, pp. 963–971, 2014.
  * [4] Matthieu Courbariaux, Yoshua Bengio, and Jean-Pierre David,  “Binaryconnect: Training deep neural networks with binary weights during propagations,”  Advances in neural information processing systems, pp. 3123–3131, 2015.
  * [5] Matthieu Courbariaux, Itay Hubara, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio,  “Binarized neural networks: Training deep neural networks with weights and activations constrained to+ 1 or-1,”  arXiv preprint arXiv:1602.02830, 2016.
  * [6] Mohammad Rastegari, Vicente Ordonez, Joseph Redmon, and Ali Farhadi,  “Xnor-net: Imagenet classification using binary convolutional neural networks,”  European conference on computer vision, pp. 525–542, 2016.
  * [7] Itay Hubara, Matthieu Courbariaux, Daniel Soudry, Ran El-Yaniv, and Yoshua Bengio,  “Quantized neural networks: Training neural networks with low precision weights and activations,”  The Journal of Machine Learning Research, vol. 18, no. 1, pp. 6869–6898, 2017.
  * [8] Shuchang Zhou, Yuxin Wu, Zekun Ni, Xinyu Zhou, He Wen, and Yuheng Zou,  “Dorefa-net: Training low bitwidth convolutional neural networks with low bitwidth gradients,”  arXiv preprint arXiv:1606.06160, 2016.
  * [9] Forrest N Iandola, Song Han, Matthew W Moskewicz, Khalid Ashraf, William J Dally, and Kurt Keutzer,  “Squeezenet alexnet-level accuracy with 50x fewer parameters and¡ 0.5 mb model size,”  arXiv preprint arXiv:1602.07360, 2016.
  * [10] Fengfu Li, Bo Zhang, and Bin Liu,  “Ternary weight networks,”  arXiv preprint arXiv:1605.04711, 2016.
  * [11] Chenzhuo Zhu, Song Han, Huizi Mao, and William J Dally,  “Trained ternary quantization,”  arXiv preprint arXiv:1612.01064, 2016.
  * [12] Qing Wu, Xiaojin Lu, Shan Xue, Chao Wang, Xundong Wu, and Jin Fan,  “Sbnn: Slimming binarized neural network,”  Neurocomputing, 2020.
  * [13] Jiwei Yang, Xu Shen, Jun Xing, Xinmei Tian, Houqiang Li, Bing Deng, Jianqiang Huang, and Xian-sheng Hua,  “Quantization networks,”  2019, pp. 7308–7316.
  * [14] Dongqing Zhang, Jiaolong Yang, Dongqiangzi Ye, and Gang Hua,  “Lq-nets: Learned quantization for highly accurate and compact deep neural networks,”  2018, pp. 365–382.
  * [15] Haotong Qin, Ruihao Gong, Xianglong Liu, Xiao Bai, Jingkuan Song, and Nicu Sebe,  “Binary neural networks: A survey,”  Pattern Recognition, p. 107281, 2020.
  * [16] Wolfgang Roth, Günther Schindler, Matthias Zöhrer, Lukas Pfeifenberger, Robert Peharz, Sebastian Tschiatschek, Holger Fröning, Franz Pernkopf, and Zoubin Ghahramani,  “Resource-efficient neural networks for embedded systems,”  arXiv preprint arXiv:2001.03048, 2020.
  * [17] Yinglan Ma, Hongyu Xiong, Zhe Hu, and Lizhuang Ma,  “Efficient super resolution using binarized neural network,”  2019, pp. 0–0.
  * [18] Cheng-Ming Chiang, Yu Tseng, Yu-Syuan Xu, Hsien-Kai Kuo, Yi-Min Tsai, Guan-Yu Chen, Koan-Sin Tan, Wei-Ting Wang, Yu-Chieh Lin, Shou-Yao Roy Tseng, et al.,  “Deploying image deblurring across mobile devices: A perspective of quality and latency,”  arXiv preprint arXiv:2004.12599, 2020.
  * [19] Alex Krizhevsky, Ilya Sutskever, and Geoffrey E Hinton,  “Imagenet classification with deep convolutional neural networks,”  Advances in neural information processing systems, pp. 1097–1105, 2012.
  * [20] Kyuyeon Hwang and Wonyong Sung,  “Fixed-point feedforward deep neural network design using weights+ 1, 0, and- 1,”  2014 IEEE Workshop on Signal Processing Systems (SiPS), pp. 1–6, 2014.
  * [21] Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov,  “Dropout: a simple way to prevent neural networks from overfitting,”  The journal of machine learning research, vol. 15, no. 1, pp. 1929–1958, 2014.
  * [22] Venice Erin Liong, Jiwen Lu, Gang Wang, Pierre Moulin, and Jie Zhou,  “Deep hashing for compact binary codes learning,”  Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 2475–2483, 2015.
  * [23] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang,  “Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising,”  IEEE Transactions on Image Processing, pp. 3142–3155, 2017.
  * [24] Wenzhe Shi, Jose Caballero, Ferenc Huszár, Johannes Totz, Andrew P Aitken, Rob Bishop, Daniel Rueckert, and Zehan Wang,  “Real-time single image and video super-resolution using an efficient sub-pixel convolutional neural network,”  Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 1874–1883, 2016.