lars/FYS-STK4155
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

added articlels

Browse Source
This commit is contained in:
Morten Hjorth-Jensen
2023-12-01 15:21:57 +01:00
parent e2d1a64bc9
commit 48bfb84b91
21 changed files with 12771 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
doc/Articles/Autodiff/17-468.bbl
+1475
View File
File diff suppressed because it is too large Load Diff
doc/Articles/Autodiff/17-468.tex
+694
View File
@@ -0,0 +1,694 @@
\documentclass[twoside,11pt]{article}
% Any additional packages needed should be included after jmlr2e.
% Note that jmlr2e.sty includes epsfig, amssymb, natbib and graphicx,
% and defines many common macros, such as 'proof' and 'example'.
%
% It also sets the bibliographystyle to plainnat; for more information on
% natbib citation styles, see the natbib documentation, a copy of which
% is archived at http://www.jmlr.org/format/natbib.pdf
\PassOptionsToPackage{hyphens}{url}
\usepackage{jmlr2e_mod}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage[inline]{enumitem}
\usepackage{paralist}
\usepackage{adjustbox}
\usepackage{tikz}
\usetikzlibrary{decorations.pathmorphing}
\usetikzlibrary{arrows}
\usetikzlibrary{positioning}
\usepackage{standalone}
\usepackage{rotating}
\usepackage{tabularx}
\usepackage{booktabs}
\usepackage{multirow}
\usepackage{todonotes}
\usepackage{microtype}
\hyphenation{DiffSharp}
\hyphenation{TensorFlow}
\hyphenation{auto-diff}
\hyphenation{ADIFOR}
% Definitions of handy macros can go here
\newcommand{\dataset}{{\cal D}}
\newcommand{\fracpartial}[2]{\frac{\partial #1}{\partial #2}}
\DeclareMathOperator*{\argmin}{arg\,min}
\newcommand{\Del}{\mathrm{\Delta}}
\newcommand{\RedMark}{{\color[rgb]{1,0,0}$X$}}
\newcommand{\GreenMark}{{\color[rgb]{0,0.69,0}$Y$}}
\newcommand{\BlueMark}{{\color[rgb]{0,0,1}$Z$}}
% Heading arguments are {volume}{year}{pages}{date submitted}{date published}{paper id}{author-full-names}
% \jmlrheading{1}{2000}{1-48}{4/00}{10/00}{meila00a}{Marina Meil\u{a} and Michael I. Jordan}
\jmlrheading{18}{2018}{1-43}{8/17}{2/18}{17-468}{Atılım Güneş Baydin, Barak~A.~Pearlmutter, Alexey~Andreyevich~Radul, and Jeffrey~Mark~Siskind}
% Short headings should be running head and authors last names
\ShortHeadings{Automatic Differentiation in Machine Learning: a Survey}{Baydin, Pearlmutter, Radul, and Siskind}
\firstpageno{1}
\begin{document}
\title{Automatic Differentiation\\in Machine Learning: a Survey}
\author{\name Atılım Güneş Baydin \email gunes@robots.ox.ac.uk \\
\addr Department of Engineering Science\\
University of Oxford\\
Oxford OX1 3PJ, United Kingdom
\AND
\name Barak~A.~Pearlmutter \email barak@pearlmutter.net \\
\addr Department of Computer Science\\
National University of Ireland Maynooth\\
Maynooth, Co. Kildare, Ireland
\AND
\name Alexey~Andreyevich~Radul \email axch@mit.edu \\
\addr Department of Brain and Cognitive Sciences\\
Massachusetts Institute of Technology\\
Cambridge, MA 02139, United States
\AND
\name Jeffrey~Mark~Siskind \email qobi@purdue.edu \\
\addr School of Electrical and Computer Engineering\\
Purdue University\\
West Lafayette, IN 47907, United States}
\editor{Léon Bottou}
\maketitle
\begin{abstract}% <- trailing '%' for backward compatibility of .sty file
Derivatives, mostly in the form of gradients and Hessians, are ubiquitous in machine learning. Automatic differentiation (AD), also called algorithmic differentiation or simply ``autodiff'', is a family of techniques similar to but more general than backpropagation for efficiently and accurately evaluating derivatives of numeric functions expressed as computer programs. AD is a small but established field with applications in areas including computational fluid dynamics, atmospheric sciences, and engineering design optimization. Until very recently, the fields of machine learning and AD have largely been unaware of each other and, in some cases, have independently discovered each other's results. Despite its relevance, general-purpose AD has been missing from the machine learning toolbox, a situation slowly changing with its ongoing adoption under the names ``dynamic computational graphs'' and ``differentiable programming''. We survey the intersection of AD and machine learning, cover applications where AD has direct relevance, and address the main implementation techniques. By precisely defining the main differentiation techniques and their interrelationships, we aim to bring clarity to the usage of the terms ``autodiff'', ``automatic differentiation'', and ``symbolic differentiation'' as these are encountered more and more in machine learning settings.
\end{abstract}
\begin{keywords}
Backpropagation, Differentiable Programming
\end{keywords}
\section{Introduction}
Methods for the computation of derivatives in computer programs can be classified into four categories:
\begin{inparaenum}[(1)]
\item manually working out derivatives and coding them;
\item \emph{numerical differentiation} using finite difference approximations;
\item \emph{symbolic differentiation} using expression manipulation in computer algebra systems such as Mathematica, Maxima, and Maple; and
\item \emph{automatic differentiation}, also called \emph{algorithmic differentiation}, which is the subject matter of this paper.
\end{inparaenum}
Conventionally, many methods in machine learning have required the evaluation of derivatives and most of the traditional learning algorithms have relied on the computation of gradients and Hessians of an objective function \citep{Sra2011}. When introducing new models, machine learning researchers have spent considerable effort on the manual derivation of analytical derivatives to subsequently plug these into standard optimization procedures such as L-BFGS \citep{Zhu1997} or stochastic gradient descent \citep{Bottou1998}. Manual differentiation is time consuming and prone to error. Of the other alternatives, numerical differentiation is simple to implement but can be highly inaccurate due to round-off and truncation errors \citep{Jerrell1997}; more importantly, it scales poorly for gradients, rendering it inappropriate for machine learning where gradients with respect to millions of parameters are commonly needed. Symbolic differentiation addresses the weaknesses of both the manual and numerical methods, but often results in complex and cryptic expressions plagued with the problem of ``expression swell'' \citep{Corliss1988}. Furthermore, manual and symbolic methods require models to be defined as closed-form expressions, ruling out or severely limiting algorithmic control flow and expressivity.
We are concerned with the powerful fourth technique, automatic differentiation (AD). AD performs a non-standard interpretation of a given computer program by replacing the domain of the variables to incorporate derivative values and redefining the semantics of the operators to propagate derivatives per the chain rule of differential calculus. Despite its widespread use in other fields, general-purpose AD has been underused by the machine learning community until very recently.\footnote{See, e.g., \url{https://justindomke.wordpress.com/2009/02/17/automatic-differentiation-the-most-criminally-underused-tool-in-the-potential-machine-learning-toolbox/}} Following the emergence of deep learning \citep{lecun2015deep,goodfellow2016deep} as the state-of-the-art in many machine learning tasks and the modern workflow based on rapid prototyping and code reuse in frameworks such as Theano \citep{Bastien2012}, Torch \citep{collobert2011torch7}, and TensorFlow \citep{abadi2016tensorflow}, the situation is slowly changing where projects such as autograd\footnote{\url{https://github.com/HIPS/autograd}} \citep{maclaurin2016modeling}, Chainer\footnote{\url{https://chainer.org/}} \citep{tokui2015chainer}, and PyTorch\footnote{\url{http://pytorch.org/}} \citep{paszke2017automatic} are leading the way in bringing general-purpose AD to the mainstream.
The term ``automatic'' in AD can be a source of confusion, causing machine learning practitioners to put the label ``automatic differentiation'', or just ``autodiff'', on any method or tool that does not involve manual differentiation, without giving due attention to the underlying mechanism. We would like to stress that AD as a technical term refers to a specific family of techniques that compute derivatives through accumulation of values during code execution to generate numerical derivative evaluations rather than derivative expressions. This allows accurate evaluation of derivatives at machine precision with only a small constant factor of overhead and ideal asymptotic efficiency. In contrast with the effort involved in arranging code as closed-form expressions under the syntactic and semantic constraints of symbolic differentiation, AD can be applied to regular code with minimal change, allowing branching, loops, and recursion. Because of this generality, AD has been applied to computer simulations in industry and academia and found applications in fields including engineering design optimization \citep{forth2002aerofoil,casanova2002application}, computational fluid dynamics \citep{Muller2005,thomas2006using,Bischof2006}, physical modeling \citep{Ekstrom2010}, optimal control \citep{Walther2007}, structural mechanics \citep{haase2002optimal}, atmospheric sciences \citep{Carmichael1997,Charpentier2000}, and computational finance \citep{Bischof2002,Capriotti2011}.
In machine learning, a specialized counterpart of AD known as the backpropagation algorithm has been the mainstay for training neural networks, with a colorful history of having been reinvented at various times by independent researchers \citep{Griewank2012,schmidhuber2015deep}. It has been one of the most studied and used training algorithms since the day it became popular mainly through the work of \citet{rumelhart1986learning}. In simplest terms, backpropagation models learning as gradient descent in neural network weight space, looking for the minima of an objective function. The required gradient is obtained by the backward propagation of the sensitivity of the objective value at the output (Figure~\ref{FigureBackpropagation}), utilizing the chain rule to compute partial derivatives of the objective with respect to each weight. The resulting algorithm is essentially equivalent to transforming the network evaluation function composed with the objective function under reverse mode AD, which, as we shall see, actually generalizes the backpropagation idea. Thus, a modest understanding of the mathematics underlying backpropagation provides one with sufficient background for grasping AD techniques.
\begin{figure}
\centering
\trimbox{0cm -0.4cm}{\resizebox{0.75\textwidth}{!}{\includegraphics{figures/backprop/backprop}}}
\caption{Overview of backpropagation. (a) Training inputs $x_i$ are fed forward, generating corresponding activations $y_i$. An error $E$ between the actual output $y_3$ and the target output $t$ is computed. (b) The error adjoint is propagated backward, giving the gradient with respect to the weights $\nabla_{w_i}E = \left(\frac{\partial E}{\partial w_1},\dots,\frac{\partial E}{\partial w_6}\right)$, which is subsequently used in a gradient-descent procedure. The gradient with respect to inputs $\nabla_{x_i}E$ can be also computed in the same backward pass.}
\label{FigureBackpropagation}
\end{figure}
In this paper we review AD from a machine learning perspective, covering its origins, applications in machine learning, and methods of implementation. Along the way, we also aim to dispel some misconceptions that we believe have impeded wider recognition of AD by the machine learning community. In Section~\ref{SectionWhatADIsNot} we start by explicating how AD differs from numerical and symbolic differentiation. Section~\ref{SectionPreliminaries} gives an introduction to the AD technique and its forward and reverse accumulation modes. Section~\ref{SectionDerivativesAndMachineLearning} discusses the role of derivatives in machine learning and examines cases where AD has relevance. Section~\ref{SectionImplementations} covers various implementation approaches and general-purpose AD tools, followed by Section~\ref{SectionConclusions} where we discuss future directions.
\section{What AD Is Not}
\label{SectionWhatADIsNot}
Without proper introduction, one might assume that AD is either a type of numerical or symbolic differentiation. Confusion can arise because AD does in fact provide numerical values of derivatives (as opposed to derivative expressions) and it does so by using symbolic rules of differentiation (but keeping track of derivative values as opposed to the resulting expressions), giving it a two-sided nature that is partly symbolic and partly numerical \citep{Griewank2003}. We start by emphasizing how AD is different from, and in several aspects superior to, these two commonly encountered techniques of computing derivatives.
\begin{figure*}
\centering
\trimbox{0cm -0.4cm}{\resizebox{\textwidth}{!}{\small\input{figures/differentiation/differentiation.tex}}}
\caption{The range of approaches for differentiating mathematical expressions and computer code, looking at the example of a truncated logistic map (upper left). Symbolic differentiation (center right) gives exact results but requires closed-form input and suffers from expression swell; numerical differentiation (lower right) has problems of accuracy due to round-off and truncation errors; automatic differentiation (lower left) is as accurate as symbolic differentiation with only a constant factor of overhead and support for control flow.}
\label{FigureDifferentiation}
\end{figure*}
\subsection{AD Is Not Numerical Differentiation}
Numerical differentiation is the finite difference approximation of derivatives using values of the original function evaluated at some sample points \citep{Burden2001} (Figure~\ref{FigureDifferentiation}, lower right). In its simplest form, it is based on the limit definition of a derivative. For example, for a multivariate function $f:\mathbb{R}^n \to \mathbb{R}$, one can approximate the gradient $\nabla f=\left(\frac{\partial f}{\partial x_1},\dots,\frac{\partial f}{\partial x_n}\right)$ using
\begin{equation}
\label{EquationForwardDifference}
\frac{\partial f(\mathbf{x})}{\partial x_i} \approx \frac{f(\mathbf{x} + h \mathbf{e}_i) - f(\mathbf{x})}{h}\;,
\end{equation}
where $\mathbf{e}_i$ is the $i$-th unit vector and $h > 0$ is a small step size. This has the advantage of being uncomplicated to implement, but the disadvantages of performing $O(n)$ evaluations of $f$ for a gradient in $n$ dimensions and requiring careful consideration in selecting the step size $h$.
Numerical approximations of derivatives are inherently ill-conditioned and unstable,\footnote{Using the limit definition of the derivative for finite difference approximation commits both cardinal sins of numerical analysis: \emph{``thou shalt not add small numbers to big numbers''}, and \emph{``thou shalt not subtract numbers which are approximately equal''}.} with the exception of complex variable methods that are applicable to a limited set of holomorphic functions \citep{Fornberg1981}. This is due to the introduction of truncation\footnote{Truncation error is the error of approximation, or inaccuracy, one gets from $h$ not actually being zero. It is proportional to a power of $h$.} and round-off\footnote{Round-off error is the inaccuracy one gets from valuable low-order bits of the final answer having to compete for machine-word space with high-order bits of $f(\mathbf{x} + h \mathbf{e}_i)$ and $f(\mathbf{x})$ (Eq.~\ref{EquationForwardDifference}), which the computer has to store just until they cancel in the subtraction at the end. Round-off error is inversely proportional to a power of $h$.} errors inflicted by the limited precision of computations and the chosen value of the step size $h$. Truncation error tends to zero as $h \to 0$. However, as $h$ is decreased, round-off error increases and becomes dominant (Figure~\ref{FigureApproximationError}).
\begin{figure*}
\centering
\resizebox{0.82\textwidth}{!}{\small\input{figures/approx-error/approx-error.tex}}
\caption{Error in the forward (Eq.~\ref{EquationForwardDifference}) and center difference (Eq.~\ref{EquationCenterDifference}) approximations as a function of step size $h$, for the derivative of the truncated logistic map ${f(x)=64x(1 - x)(1 - 2 x)^2 (1 - 8 x + 8 x^2)^2}$. Plotted errors are computed using ${E_{\mathrm{forward}}(h,x_0)=\left|\frac{f(x_0+h)-f(x_0)}{h} - \frac{d}{dx}f(x)\big|_{x_0}\right|}$ and ${E_{\mathrm{center}}(h,x_0)=\left|\frac{f(x_0+h)-f(x_0-h)}{2h} - \frac{d}{dx}f(x)\big|_{x_0}\right|}$ at $x_0=0.2$\;.}
\label{FigureApproximationError}
\end{figure*}
Various techniques have been developed to mitigate approximation errors in numerical differentiation, such as using a center difference approximation
\begin{equation}
\label{EquationCenterDifference}
\frac{\partial f(\mathbf{x})}{\partial x_i} = \frac{f(\mathbf{x} + h \mathbf{e}_i) - f(\mathbf{x} - h \mathbf{e}_i)}{2h} + O(h^{2})\;,
\end{equation}
where the first-order errors cancel and one effectively moves the truncation error from first-order to second-order in $h$.\footnote{This does not avoid either of the cardinal sins, and is still highly inaccurate due to truncation.} For the one-dimensional case, it is just as costly to compute the forward difference (Eq.~\ref{EquationForwardDifference}) and the center difference (Eq.~\ref{EquationCenterDifference}), requiring only two evaluations of $f$. However, with increasing dimensionality, a trade-off between accuracy and performance is faced, where computing a Jacobian matrix of a function $f: \mathbb{R}^n \to \mathbb{R}^m$ requires $2mn$ evaluations.
Other techniques for improving numerical differentiation, including higher-order finite differences, Richardson extrapolation to the limit \citep{Brezinski1991}, and differential quadrature methods using weighted sums \citep{Bert1996}, have increased computational complexity, do not completely eliminate approximation errors, and remain highly susceptible to floating point truncation.
The $O(n)$ complexity of numerical differentiation for a gradient in $n$ dimensions is the main obstacle to its usefulness in machine learning, where $n$ can be as large as millions or billions in state-of-the-art deep learning models \citep{shazeer2017outrageously}. In contrast, approximation errors would be tolerated in a deep learning setting thanks to the well-documented error resiliency of neural network architectures \citep{gupta2015deep}.
\subsection{AD Is Not Symbolic Differentiation}
Symbolic differentiation is the automatic manipulation of expressions for obtaining derivative expressions \citep{Grabmeier2003} (Figure~\ref{FigureDifferentiation}, center right), carried out by applying transformations representing rules of differentiation such as
\begin{equation}
\begin{aligned}
\frac{d}{dx} \left(f(x) + g(x)\right) &\leadsto \frac{d}{dx} f(x) + \frac{d}{dx} g(x)\\
\frac{d}{dx} \left(f(x)\,g(x)\right) &\leadsto \left(\frac{d}{dx} f(x)\right) g(x) + f(x) \left(\frac{d}{dx} g(x)\right)\; .
\end{aligned}
\label{EquationMultiplicationRule}
\end{equation}
When formulae are represented as data structures, symbolically differentiating an expression tree is a perfectly mechanistic process, considered subject to mechanical automation even at the very inception of calculus \citep{Leibniz1685}. This is realized in modern computer algebra systems such as Mathematica, Maxima, and Maple and machine learning frameworks such as Theano.
In optimization, symbolic derivatives can give valuable insight into the structure of the problem domain and, in some cases, produce analytical solutions of extrema (e.g., solving for $\frac{d}{dx}f(x)=0$) that can eliminate the need for derivative calculation altogether. On the other hand, symbolic derivatives do not lend themselves to efficient runtime calculation of derivative values, as they can get exponentially larger than the expression whose derivative they represent.
Consider a function $h(x)=f(x)g(x)$ and the multiplication rule in Eq.~\ref{EquationMultiplicationRule}. Since $h$ is a product, $h(x)$ and $\frac{d}{dx}h(x)$ have some common components, namely $f(x)$ and $g(x)$. Note also that on the right hand side, $f(x)$ and $\frac{d}{dx}f(x)$ appear separately. If we just proceeded to symbolically differentiate $f(x)$ and plugged its derivative into the appropriate place, we would have nested duplications of any computation that appears in common between $f(x)$ and $\frac{d}{dx}f(x)$. Hence, careless symbolic differentiation can easily produce exponentially large symbolic expressions which take correspondingly long to evaluate. This problem is known as \emph{expression swell} (Table~\ref{TableExpressionSwell}).
\begin{table}
\centering
\renewcommand{\arraystretch}{1.2}
\caption{Iterations of the logistic map $l_{n+1}=4l_n (1-l_n)$, $l_1=x$ and the corresponding derivatives of $l_n$ with respect to $x$, illustrating expression swell.}
\label{TableExpressionSwell}
{\small
\begin{tabularx}{\columnwidth}{@{}lp{2.8cm}XX@{}}
\toprule
$n$ & $l_n$ & $\frac{d}{dx}l_n$ & $\frac{d}{dx}l_n$ (Simplified form)\\
\addlinespace
\midrule
1 & $x$ & $1$ & $1$\\
\addlinespace
2 & $4x(1 - x)$ & $4(1 - x) -4x$ & $4 - 8x$\\
\addlinespace
3 & $16x(1 - x)(1 - 2 x)^2$ & $16(1 - x)(1 - 2 x)^2 - 16x(1 - 2 x)^2 - 64x(1 - x)(1 - 2 x)$ & $16 (1 - 10 x + 24 x^2 - 16 x^3)$\\
\addlinespace
4 & $64x(1 - x)(1 - 2 x)^2$ $(1 - 8 x + 8 x^2)^2$ & $128x(1 - x)(-8 + 16 x)(1 - 2 x)^2 (1 - 8 x + 8 x^2) + 64 (1 - x)(1 - 2 x)^2 (1 - 8 x + 8 x^2)^2 - 64x(1 - 2 x)^2 (1 - 8 x + 8 x^2)^2 - 256x(1 - x)(1 - 2 x)(1 - 8 x + 8 x^2)^2$ & $64 (1 - 42 x + 504 x^2 - 2640 x^3 + 7040 x^4 - 9984 x^5 + 7168 x^6 - 2048 x^7)$\\
\bottomrule
\end{tabularx}}
\end{table}
When we are concerned with the accurate numerical evaluation of derivatives and not so much with their actual symbolic form, it is in principle possible to significantly simplify computations by storing only the values of intermediate sub-expressions in memory. Moreover, for further efficiency, we can interleave as much as possible the differentiation and simplification steps. This interleaving idea forms the basis of AD and provides an account of its simplest form: \emph{apply symbolic differentiation at the elementary operation level and keep intermediate numerical results, in lockstep with the evaluation of the main function.} This is AD in the forward accumulation mode, which we shall introduce in the following section.
\section{AD and Its Main Modes}
\label{SectionPreliminaries}
AD can be thought of as performing a non-standard interpretation of a computer program where this interpretation involves augmenting the standard computation with the calculation of various derivatives. All numerical computations are ultimately compositions of a finite set of elementary operations for which derivatives are known \citep{Verma2000,Griewank2008}, and combining the derivatives of the constituent operations through the chain rule gives the derivative of the overall composition. Usually these elementary operations include the binary arithmetic operations, the unary sign switch, and transcendental functions such as the exponential, the logarithm, and the trigonometric functions.
On the left hand side of Table~\ref{TableForwardADExample} we see the representation of the computation $y = f(x_1, x_2) = \ln(x_1) + x_1 x_2 - \sin(x_2)$ as an \emph{evaluation trace} of elementary operations---also called a Wengert list \citep{Wengert1964}. We adopt the three-part notation used by \citet{Griewank2008}, where a function $f: \mathbb{R}^n \to \mathbb{R}^m$ is constructed using intermediate variables $v_i$ such that
\begin{compactitem}
\item variables $v_{i-n} = x_i,\;i = 1, \dotsc, n$ are the input variables,
\item variables $v_i\;i = 1, \dotsc, l$ are the working (intermediate) variables, and
\item variables $y_{m-i} = v_{l-i},\;i = m - 1, \dotsc, 0$ are the output variables.
\end{compactitem}
Figure~\ref{FigureComputationalGraph} shows the given trace of elementary operations represented as a computational graph \citep{Bauer1974}, useful in visualizing dependency relations between intermediate variables.
\begin{figure}
\centering
\trimbox{0cm -0.5cm}{\resizebox{0.8\textwidth}{!}{\normalsize\input{figures/comp-graph/comp-graph.tex}}}
\caption{Computational graph of the example $f(x_1, x_2) = \ln(x_1) + x_1 x_2 - \sin(x_2)$. See the primal trace in Tables \ref{TableForwardADExample} or \ref{TableReverseADExample} for the definitions of the intermediate variables $v_{-1} \dots v_5$\;.}
\label{FigureComputationalGraph}
\end{figure}
Evaluation traces form the basis of the AD techniques. An important point to note here is that AD can differentiate not only closed-form expressions in the classical sense, but also algorithms making use of control flow such as branching, loops, recursion, and procedure calls, giving it an important advantage over symbolic differentiation which severely limits such expressivity. This is thanks to the fact that any numeric code will eventually result in a numeric evaluation trace with particular values of the input, intermediate, and output variables, which are the only things one needs to know for computing derivatives using chain rule composition, regardless of the specific control flow path that was taken during execution. Another way of expressing this is that AD is blind with respect to any operation, including control flow statements, which do not directly alter numeric values.
\subsection{Forward Mode}
AD in forward accumulation mode\footnote{Also called \emph{tangent linear} mode.} is the conceptually most simple type. Consider the evaluation trace of the function $f(x_1, x_2) = \ln(x_1) + x_1 x_2 - \sin(x_2)$ given on the left-hand side in Table~\ref{TableForwardADExample} and in graph form in Figure~\ref{FigureComputationalGraph}. For computing the derivative of $f$ with respect to $x_1$, we start by associating with each intermediate variable $v_i$ a derivative
\begin{equation*}
\dot{v}_i = \frac{\partial v_i}{\partial x_1}\; .
\end{equation*}
Applying the chain rule to each elementary operation in the forward primal trace, we generate the corresponding tangent (derivative) trace, given on the right-hand side in Table~\ref{TableForwardADExample}. Evaluating the primals $v_i$ in lockstep with their corresponding tangents $\dot{v}_i$ gives us the required derivative in the final variable $\dot{v}_5=\frac{\partial y}{\partial x_1}$\;.
%\footnote{In an implementation, the order in which $v_i$ and its corresponding $\dot{v}_i$ are computed can make a difference.}
\begin{table}
\centering
\renewcommand{\arraystretch}{1.2}
\caption{Forward mode AD example, with $y = f(x_1, x_2) = \ln(x_1) + x_1 x_2 - \sin(x_2)$ evaluated at $(x_1, x_2) = (2, 5)$ and setting $\dot{x}_1 = 1$ to compute $\frac{\partial y}{\partial x_1}$. The original forward evaluation of the primals on the left is augmented by the tangent operations on the right, where each line complements the original directly to its left.}
\label{TableForwardADExample}
\begin{minipage}[c]{0.47\textwidth}
{\footnotesize
\begin{tabularx}{\textwidth}{p{0.2mm}p{2mm}p{16mm}X}
\toprule
\multicolumn{4}{l}{Forward Primal Trace}\\
\multirow{9}{2mm}{\begin{tikzpicture}\draw[->,>=triangle 60,thick](0,0)--(0,-3.8);\end{tikzpicture}} & $v_{-1}$ & $=x_1$ & $=2$\\
& $v_0$ & $=x_2$ & $=5$\\
\cmidrule{2-4}
& $v_1$ & $=\ln{v_{-1}}$ & $=\ln{2}$\\
& $v_2$ & $=v_{-1} \times v_0$ & $=2 \times 5$\\
& $v_3$ & $=\sin{v_0}$ & $=\sin{5}$\\
& $v_4$ & $=v_1+v_2$ & $=0.693+10$\\
& $v_5$ & $=v_4-v_3$ & $=10.693+0.959$\\
\cmidrule{2-4}
& $y$ & $=v_5$ & $=11.652$\\
\bottomrule
\end{tabularx}}
\end{minipage}
\begin{minipage}[c]{0.52\textwidth}
\setlength{\fboxsep}{0pt}\colorbox{gray!20}
{\footnotesize
\begin{tabularx}{\textwidth}{p{0.2mm}p{2mm}p{28mm}X}
\toprule
\multicolumn{4}{l}{Forward Tangent (Derivative) Trace}\\
\multirow{9}{2mm}{\begin{tikzpicture}\draw[->,>=triangle 60,thick](0,0)--(0,-3.8);\end{tikzpicture}} & $\dot{v}_{-1}$ & $=\dot{x}_1$ & $=1$\\
& $\dot{v}_0$ & $=\dot{x}_2$ & $=0$\\
\cmidrule{2-4}
& $\dot{v}_1$ & $=\dot{v}_{-1}/v_{-1}$ & $=1/2$\\
& $\dot{v}_2$ & $=\dot{v}_{-1} \times v_0 + \dot{v}_0 \times v_{-1}$ & $=1 \times 5 + 0 \times 2$\\
& $\dot{v}_3$ & $=\dot{v}_0 \times \cos{v_0}$ & $=0 \times \cos{5}$\\
& $\dot{v}_4$ & $=\dot{v}_1+\dot{v}_2$ & $=0.5+5$\\
& $\dot{v}_5$ & $=\dot{v}_4-\dot{v}_3$ & $=5.5-0$\\
\cmidrule{2-4}
& \boldmath$\dot{y}$ & \boldmath$=\dot{v}_5$ & \boldmath$=5.5$\\
\bottomrule
\end{tabularx}}
\end{minipage}
\end{table}
This generalizes naturally to computing the Jacobian of a function $f : \mathbb{R}^n \to \mathbb{R}^m$ with $n$ independent (input) variables $x_i$ and $m$ dependent (output) variables $y_j$. In this case, each forward pass of AD is initialized by setting only one of the variables $\dot{x}_i=1$ and setting the rest to zero (in other words, setting $\dot{\mathbf{x}} = \mathbf{e}_i$, where $\mathbf{e}_i$ is the $i$-th unit vector). A run of the code with specific input values $\mathbf{x}=\mathbf{a}$ then computes
\begin{equation*}
\dot{y}_j = \left.\frac{\partial y_j}{\partial x_i}\right|_{\mathbf{x}=\mathbf{a}},\;j = 1, \dotsc, m\;,
\end{equation*}
giving us one column of the Jacobian matrix
\begin{equation*}
\mathbf{J}_f = \left. \begin{bmatrix}
\frac{\partial y_1}{\partial x_1} & \cdots & \frac{\partial y_1}{\partial x_n} \\
\vdots & \ddots & \vdots \\
\frac{\partial y_m}{\partial x_1} & \cdots & \frac{\partial y_m}{\partial x_n}
\end{bmatrix} \right|_{\mathbf{x}\; = \; \mathbf{a}}
\end{equation*}
evaluated at point $\mathbf{a}$. Thus, the full Jacobian can be computed in $n$ evaluations.
Furthermore, forward mode AD provides a very efficient and matrix-free way of computing Jacobian--vector products
\begin{equation}
\mathbf{J}_f\,\mathbf{r} = \begin{bmatrix}
\frac{\partial y_1}{\partial x_1} & \cdots & \frac{\partial y_1}{\partial x_n} \\
\vdots & \ddots & \vdots \\
\frac{\partial y_m}{\partial x_1} & \cdots & \frac{\partial y_m}{\partial x_n}
\end{bmatrix}
\begin{bmatrix}
r_1 \\
\vdots \\
r_n
\end{bmatrix}\; ,
\label{EquationJacobianVectorProduct}
\end{equation}
simply by initializing with $\dot{\mathbf{x}}=\mathbf{r}$. Thus, we can compute the Jacobian--vector product in just one forward pass. As a special case, when $f: \mathbb{R}^n \to \mathbb{R}$, we can obtain the directional derivative along a given vector $\mathbf{r}$ as a linear combination of the partial derivatives
\begin{equation*}
\nabla f \cdot \mathbf{r}
\end{equation*}
by starting the AD computation with the values $\dot{\mathbf{x}}=\mathbf{r}$.
Forward mode AD is efficient and straightforward for functions $f: \mathbb{R} \to \mathbb{R}^m$, as all the derivatives $\frac{d y_i}{d x}$ can be computed with just one forward pass. Conversely, in the other extreme of $f: \mathbb{R}^n \to \mathbb{R}$, forward mode AD requires $n$ evaluations to compute the gradient
\begin{equation*}
\nabla f = \left( \frac{\partial y}{\partial x_1}, \dots, \frac{\partial y}{\partial x_n}\right)\; ,
\end{equation*}
which also corresponds to a $1 \times n$ Jacobian matrix that is built one column at a time with the forward mode in $n$ evaluations.
In general, for cases $f: \mathbb{R}^n \to \mathbb{R}^m$ where $n \gg m$, a different technique is often preferred.
We will describe AD in \emph{reverse accumulation mode} in Section~\ref{sec:reverse-mode}.
\subsubsection{Dual Numbers}
\label{SectionDualNumbers}
Mathematically, forward mode AD (represented by the left- and right-hand sides in Table~\ref{TableForwardADExample}) can be viewed as evaluating a function using dual numbers,\footnote{First introduced by \citet{Clifford1873}, with important uses in linear algebra and physics.} which can be defined as truncated Taylor series of the form
\begin{equation*}
v + \dot{v}\epsilon \;,
\end{equation*}
where $v, \dot{v} \in \mathbb{R}$ and $\epsilon$ is a nilpotent number such
that $\epsilon^2 = 0$ and $\epsilon \neq 0$. Observe, for example, that
\begin{align*}
(v + \dot{v}\epsilon) + (u + \dot{u}\epsilon) &= (v + u) + (\dot{v} + \dot{u})\epsilon\\
(v + \dot{v}\epsilon)(u + \dot{u}\epsilon) &= (vu) + (v\dot{u} + \dot{v}u)\epsilon\;,
\end{align*}
in which the coefficients of $\epsilon$ conveniently mirror symbolic differentiation rules (e.g., Eq.~\ref{EquationMultiplicationRule}). We can utilize this by setting up a regime where
\begin{equation}\label{EquationDualRule}
f(v + \dot{v}\epsilon) = f(v) + f'(v)\dot{v}\epsilon
\end{equation}
and using dual numbers as data structures for carrying the tangent value together with the primal.\footnote{Just as the complex number written $x + y i$ is represented in the computer as a pair in memory $(x, y)$ whose two slots are reals, the dual number written $x + \dot{x}\epsilon$ is represented as the pair $(x, \dot{x})$. Such pairs are sometimes called Argand pairs \citep[][p107 Eqs.~(157) and (158)]{Hamilton1837}.} The chain rule works as expected on this representation: two applications of Eq.~\ref{EquationDualRule} give
\begin{align*}
f(g(v + \dot{v}\epsilon)) &= f(g(v) + g'(v)\dot{v}\epsilon)\\
&= f(g(v)) + f'(g(v))g'(v)\dot{v}\epsilon\;.
\end{align*}
The coefficient of $\epsilon$ on the right-hand side is exactly the derivative of the composition of $f$ and $g$. This means that since we implement elementary operations to respect the invariant Eq.~\ref{EquationDualRule}, all compositions of them will also do so. This, in turn, means that we can extract the derivative of a function by interpreting any non-dual number $v$ as $v + 0 \epsilon$ and evaluating the function in this non-standard way on an initial input with a coefficient $1$ for $\epsilon$:
\begin{align*}
\left.\frac{df(x)}{dx}\right|_{x=v} = \textrm{epsilon-coefficient}(\textrm{dual-version}(f)(v + 1\epsilon))\;.
\end{align*}
This also extends to arbitrary program constructs, since dual numbers, as data types, can be contained in any data structure. As long as a dual number remains in a data structure with no arithmetic operations being performed on it, it will just remain a dual number; and if it is taken out of the data structure and operated on again, then the differentiation will continue.
In practice, a function $f$ coded in a programming language of choice would be fed into an AD tool, which would then augment it with corresponding extra code to handle the dual operations so that the function and its derivative are simultaneously computed. This can be implemented through calls to a specific library, in the form of source code transformation where a given source code will be automatically modified, or through operator overloading, making the process transparent to the user. We discuss these implementation techniques in Section~\ref{SectionImplementations}.
\subsection{Reverse Mode}
\label{sec:reverse-mode}
AD in the reverse accumulation mode\footnote{Also called \emph{adjoint} or \emph{cotangent linear} mode.} corresponds to a generalized backpropagation algorithm, in that it propagates derivatives backward from a given output. This is done by complementing each intermediate variable $v_i$ with an adjoint
\begin{equation*}
\bar{v}_i = \frac{\partial y_j}{\partial v_i}\; ,
\end{equation*}
which represents the sensitivity of a considered output $y_j$ with respect to changes in $v_i$. In the case of backpropagation, $y$ would be a scalar corresponding to the error $E$ (Figure~\ref{FigureBackpropagation}).
In reverse mode AD, derivatives are computed in the second phase of a two-phase process. In the first phase, the original function code is run \emph{forward}, populating intermediate variables $v_i$ and recording the dependencies in the computational graph through a bookkeeping procedure. In the second phase, derivatives are calculated by propagating adjoints $\bar{v}_i$ in \emph{reverse}, from the outputs to the inputs.
Returning to the example $y = f(x_1, x_2) = \ln(x_1) + x_1 x_2 - \sin(x_2)$, in Table~\ref{TableReverseADExample} we see the adjoint statements on the right-hand side, corresponding to each original elementary operation on the left-hand side. In simple terms, we are interested in computing the contribution $\bar{v}_i = \frac{\partial y}{\partial v_i}$ of the change in each variable $v_i$ to the change in the output $y$. Taking the variable $v_0$ as an example, we see in Figure~\ref{FigureComputationalGraph} that the only way it can affect $y$ is through affecting $v_2$ and $v_3$, so its contribution to the change in $y$ is given by
\begin{align*}
\frac{\partial y}{\partial v_0} &= \frac{\partial y}{\partial v_2}\frac{\partial v_2}{\partial v_0} + \frac{\partial y}{\partial v_3}\frac{\partial v_3}{\partial v_0}&
\text{or}&&
\bar{v}_0 &= \bar{v}_2\frac{\partial v_2}{\partial v_0} + \bar{v}_3\frac{\partial v_3}{\partial v_0}\;.
\end{align*}
In Table~\ref{TableReverseADExample}, this contribution is computed in two incremental steps
\begin{align*}
\bar{v}_0 &= \bar{v}_3\frac{\partial v_3}{\partial v_0} &
\text{and} &&
\bar{v}_0 &= \bar{v}_0 + \bar{v}_2\frac{\partial v_2}{\partial v_0}\;,
\end{align*}
lined up with the lines in the forward trace from which these expressions originate.
After the forward pass on the left-hand side, we run the reverse pass of the adjoints on the right-hand side, starting with $\bar{v}_5 = \bar{y} = \frac{\partial y}{\partial y} = 1$. In the end we get the derivatives $\frac{\partial y}{\partial x_1} = \bar{x}_1$ and $\frac{\partial y}{\partial x_2} = \bar{x}_2$ in just one reverse pass.
\begin{table}
\centering
\renewcommand{\arraystretch}{1.2}
\caption{Reverse mode AD example, with $y = f(x_1, x_2) = \ln(x_1) + x_1 x_2 - \sin(x_2)$ evaluated at $(x_1, x_2) = (2, 5)$. After the forward evaluation of the primals on the left, the adjoint operations on the right are evaluated in reverse (cf.\ Figure~\ref{FigureBackpropagation}). Note that both $\frac{\partial y}{\partial x_1}$ and $\frac{\partial y}{\partial x_2}$ are computed in the same reverse pass, starting from the adjoint $\bar{v}_5 = \bar{y} = \frac{\partial y}{\partial y} = 1$.}
\label{TableReverseADExample}
\begin{minipage}[t]{0.41\textwidth}
{\footnotesize
\begin{tabularx}{\textwidth}[t]{p{0.5mm}p{0.8mm}p{18mm}@{}X}
\toprule
\multicolumn{4}{l}{Forward Primal Trace}\\
\multirow{9}{1mm}{\begin{tikzpicture}\draw[->,>=triangle 60,thick](0,0)--(0,-5.4);\end{tikzpicture}} & $v_{-1}$ & $=x_1$ & $=2$\\
& $v_0$ & $=x_2$ & $=5$\\
\cmidrule{2-4}
& $v_1$ & $=\ln{v_{-1}}$ & $=\ln{2}$\vspace{0.25mm}\\
& $v_2$ & $=v_{-1} \times v_0$ & $=2 \times 5$\vspace{0.25mm}\\
&\vspace{0.25mm}\\
& $v_3$ & $=\sin{v_0}$ & $=\sin{5}$\vspace{0.25mm}\\
& $v_4$ & $=v_1+v_2$ & $=0.693+10$\vspace{0.25mm}\\
&\vspace{0.25mm}\\
& $v_5$ & $=v_4-v_3$ & $=10.693+0.959$\vspace{0.25mm}\\
&\vspace{0.25mm}\\
\cmidrule{2-4}
& $y$ & $=v_5$ & $=11.652$\\
\bottomrule
\end{tabularx}}\vspace{1mm}
\end{minipage}
\begin{minipage}[t]{0.58\textwidth}
\setlength{\fboxsep}{0pt}\colorbox{gray!20}
{\footnotesize
\begin{tabularx}{\textwidth}[t]{p{0.5mm}p{1mm}p{23mm}@{\hspace{1mm}}p{24mm}@{}X}
\toprule
\multicolumn{5}{l}{Reverse Adjoint (Derivative) Trace}\\
\multirow{9}{1mm}{\begin{tikzpicture}\draw[<-,>=triangle 60,thick](0,0)--(0,-5.4);\end{tikzpicture}} & \boldmath$\bar{x}_1$ & \boldmath$=\bar{v}_{-1}$ & & \boldmath$=5.5$\\
& \boldmath$\bar{x}_2$ & \boldmath$=\bar{v}_0$ & & \boldmath$=1.716$\\
\cmidrule{2-5}
& $\bar{v}_{-1}$ & $=\bar{v}_{-1} + \bar{v}_1 \frac{\partial v_1}{\partial v_{-1}}$ & $=\bar{v}_{-1} + \bar{v}_1 / v_{-1}$ & $=5.5$\\
& $\bar{v}_0$ & $=\bar{v}_0 + \bar{v}_2 \frac{\partial v_2}{\partial v_0}$ & $=\bar{v}_0 + \bar{v}_2 \times v_{-1}$ & $=1.716$\\
& $\bar{v}_{-1}$ & $=\bar{v}_2 \frac{\partial v_2}{\partial v_{-1}}$ & $=\bar{v}_2 \times v_0$ & $=5$\\
& $\bar{v}_0$ & $=\bar{v}_3 \frac{\partial v_3}{\partial v_0}$ & $=\bar{v}_3 \times \cos{v_0}$ & $=-0.284$\\
& $\bar{v}_2$ & $=\bar{v}_4 \frac{\partial v_4}{\partial v_2}$ & $=\bar{v}_4 \times 1$ & $=1$\\
& $\bar{v}_1$ & $=\bar{v}_4 \frac{\partial v_4}{\partial v_1}$ & $=\bar{v}_4 \times 1$ & $=1$\\
& $\bar{v}_3$ & $=\bar{v}_5 \frac{\partial v_5}{\partial v_3}$ & $=\bar{v}_5 \times (-1)$ & $=-1$\\
& $\bar{v}_4$ & $=\bar{v}_5 \frac{\partial v_5}{\partial v_4}$ & $=\bar{v}_5 \times 1$ & $=1$\\
\cmidrule{2-5}
& $\bar{v}_5$ & $=\bar{y}$ & $=1$\\
\bottomrule
\end{tabularx}}
\end{minipage}
\end{table}
Compared with the straightforwardness of forward accumulation mode, reverse mode AD can, at first, appear somewhat ``mysterious'' \citep{Dennis1996}. \citet{Griewank2008} argue that this is in part because of the common acquaintance with the chain rule as a mechanistic procedure propagating derivatives forward.
An important advantage of the reverse mode is that it is significantly less costly to evaluate (in terms of operation count) than the forward mode for functions with a large number of inputs. In the extreme case of $f: \mathbb{R}^n \to \mathbb{R}$, only one application of the reverse mode is sufficient to compute the full gradient $\nabla f = \left(\frac{\partial y}{\partial x_1},\dots,\frac{\partial y}{\partial x_n}\right)$, compared with the $n$ passes of the forward mode needed for populating the same. Because machine learning practice principally involves the gradient of a scalar-valued objective with respect to a large number of parameters, this establishes the reverse mode, as opposed to the forward mode, as the mainstay technique in the form of the backpropagation algorithm.
In general, for a function $f: \mathbb{R}^n \to \mathbb{R}^m$, if we denote the operation count to evaluate the original function by $\textrm{ops}(f)$, the time it takes to calculate the $m \times n$ Jacobian by the forward mode is $n\;c\;\textrm{ops}(f)$, whereas the same computation can be done via reverse mode in $m\;c\;\textrm{ops}(f)$, where $c$ is a constant guaranteed to be $c<6$ and typically $c \sim [2,3]$ \citep{Griewank2008}. That is to say, reverse mode AD performs better when $m \ll n$.
Similar to the matrix-free computation of Jacobian--vector products with forward mode (Eq.~\ref{EquationJacobianVectorProduct}), reverse mode can be used for computing the transposed Jacobian--vector product
\begin{equation*}
\mathbf{J}^{\intercal}_f\,\mathbf{r} = \begin{bmatrix}
\frac{\partial y_1}{\partial x_1} & \cdots & \frac{\partial y_m}{\partial x_1} \\
\vdots & \ddots & \vdots \\
\frac{\partial y_1}{\partial x_n} & \cdots & \frac{\partial y_m}{\partial x_n}
\end{bmatrix}
\begin{bmatrix}
r_1 \\
\vdots \\
r_m
\end{bmatrix}\;,
\end{equation*}
by initializing the reverse phase with $\bar{\mathbf{y}}=\mathbf{r}$.
The advantages of reverse mode AD, however, come with the cost of increased storage requirements growing (in the worst case) in proportion to the number of operations in the evaluated function. It is an active area of research to improve storage requirements in implementations by using advanced methods such as checkpointing strategies and data-flow analysis \citep{Dauvergne2006,siskind2017divide}.
\subsection{Origins of AD and Backpropagation}
Ideas underlying AD date back to the 1950s \citep{Nolan1953,Beda1959}. Forward mode AD as a general method for evaluating partial derivatives was essentially discovered by \citet{Wengert1964}. It was followed by a period of relatively low activity, until interest in the field was revived in the 1980s mostly through the work of \citet{Griewank1989}, also supported by improvements in modern programming languages and the feasibility of an efficient reverse mode AD.
Reverse mode AD and backpropagation have an intertwined history. The essence of the reverse mode, cast in a continuous-time formalism, is the Pontryagin maximum principle \citep{Rozonoer-Pontryagin-1959a, Boltyanskii-Gamkrelidze-Pontryagin-1960a}. This method was understood in the control theory community \citep{Bryson-1962a, Bryson-Ho-1969a} and cast in more formal terms with discrete-time variables topologically sorted in terms of dependency by \citet{Werbos-1974a}. Prior to Werbos, the work by \citet{linnainmaa1970representation,linnainmaa1976taylor} is often cited as the first published description of the reverse mode. \citet{Speelpenning80} subsequently introduced reverse mode AD as we know it, in the sense that he gave the first implementation that was actually automatic, accepting a specification of a computational process written in a general-purpose programming language and automatically performing the reverse mode transformation.
Incidentally, \citet{Hecht1989} cites the work of \citet{Bryson-Ho-1969a} and \citet{Werbos-1974a} as the two earliest known instances of backpropagation. Within the machine learning community, the method has been reinvented several times, such as by \citet{Parker1985}, until it was eventually brought to fame by \citet{rumelhart1986learning} and the Parallel Distributed Processing (PDP) group. The PDP group became aware of Parker's work only after their own discovery; similarly, Werbos' work was not appreciated until it was found by Parker \citep{Hecht1989}. This tells us an interesting story of two highly interconnected research communities that have somehow also managed to stay detached during this foundational period.
For a thorough review of the development of AD, we advise readers to refer to \citet{Rall2006}. Interested readers are highly recommended to read \citet{Griewank2012} for an investigation of the origins of the reverse mode and \citet{schmidhuber2015deep} for the same for backpropagation.
\section{AD and Machine Learning}
\label{SectionDerivativesAndMachineLearning}
In the following, we examine the main uses of derivatives in machine learning and report on a selection of works where general-purpose AD, as opposed to just backpropagation, has been successfully applied in a machine learning context. Areas where AD has seen use include optimization, neural networks, computer vision, natural language processing, and probabilistic inference.
\subsection{Gradient-Based Optimization}
Gradient-based optimization is one of the pillars of machine learning \citep{bottou2016optimization}. Given an objective function $f: \mathbb{R}^n \to \mathbb{R}$, classical gradient descent has the goal of finding (local) minima $\mathbf{w}^* = \argmin_{\mathbf{w}} f(\mathbf{w})$ via updates of the form $\Del \mathbf{w} = -\eta \nabla f$, where $\eta>0$ is a step size. Gradient-based methods make use of the fact that $f$ decreases steepest if one goes in the direction of the negative gradient. The convergence rate of gradient-based methods is usually improved by adaptive step-size techniques that adjust the step size $\eta$ on every iteration \citep{duchi2011adaptive,schaul2013no,kingma2015adam}.
As we have seen, for large $n$, reverse mode AD provides a highly efficient method for computing gradients.\footnote{See \url{http://DiffSharp.github.io/DiffSharp/examples-gradientdescent.html} for an example of a general-purpose AD-based gradient descent routine using DiffSharp.} Figure~\ref{FigureHelmholtz} and Table~\ref{TableHelmholtz} demonstrate how gradient computation scales differently for forward and reverse mode AD and numerical differentiation, looking at the Helmholtz free energy function that has been used in AD literature for benchmarking gradient calculations \citep{Griewank1989,Griewank2008,griewank2012numerical}.
\begin{figure}
\centering
\trimbox{0cm -1cm}{\resizebox{0.65\textwidth}{!}{\normalsize\input{figures/helmholtz/helmholtz.tex}}}
\caption{Evaluation time of the Helmholtz free energy function of a mixed fluid, based on the Peng-Robinson equation of state \citep{Peng1976}, ${f(\mathbf{x}) = R \, T \sum_{i = 0}^{n} \log \frac{x_i}{1 - \mathbf{b^T} \mathbf{x}} - \frac{\mathbf{x^T} \mathbf{A} \mathbf{x}}{\sqrt{8} \mathbf{b^T} \mathbf{x}} \log \frac{1 + (1 + \sqrt{2}) \mathbf{b^T} \mathbf{x}}{1 + (1 - \sqrt{2}) \mathbf{b^T} \mathbf{x}}}$, where $R$ is the universal gas constant, $T$ is the absolute temperature, $\mathbf{b} \in \mathbb{R}^n$ is a vector of constants, $\mathbf{A} \in \mathbb{R}^{n \times n}$ is a symmetric matrix of constants, and $\mathbf{x} \in \mathbb{R}^n$ is the vector of independent variables describing the system. The plots show the evaluation time of $f$ and the gradient $\nabla f$ with numerical differentiation (central difference), forward mode AD, and reverse mode AD, as a function of the number of variables $n$. Reported times are relative to the evaluation time of $f$ with $n=1$. The lower plot uses logarithmic scale for illustrating the behavior for small $n$. Numerical results are given in Table~\ref{TableHelmholtz}. (Code: {\small\url{http://DiffSharp.github.io/DiffSharp/misc/Benchmarks-h-grad-v0.5.7.fsx}})}
\label{FigureHelmholtz}
\end{figure}
\begin{table}
\centering
\renewcommand{\arraystretch}{1.2}
\setlength{\tabcolsep}{1.85mm}
\caption{Evaluation times of the Helmholtz free energy function and its gradient (Figure~\ref{FigureHelmholtz}). Times are given relative to that of the original function with both (1) $n=1$ and (2) $n$ corresponding to each column. (For instance, reverse mode AD with $n=43$ takes approximately twice the time to evaluate relative to the original function with $n=43$.) Times are measured by averaging a thousand runs on a machine with Intel Core i7-4785T 2.20 GHz CPU and 16 GB RAM, using DiffSharp 0.5.7. The evaluation time for the original function with $n=1$ is 0.0023 ms.}
\label{TableHelmholtz}
{\small
\begin{tabularx}{\columnwidth}{@{}p{42.8mm}rrrrrrrr@{}}
\toprule
& \multicolumn{4}{l}{$n$, number of variables}\\
\cmidrule(l){2-9}
& 1 & 8 & 15 & 22 & 29 & 36 & 43 & 50 \\
\midrule
$f$, original\\
\hspace{2mm} Relative $n=1$ & 1 & 5.12 & 14.51 & 29.11 & 52.58 & 84.00 & 127.33 & 174.44 \\
$\nabla f$, numerical diff.\\
\hspace{2mm} Relative $n=1$ & 1.08 & 35.55 & 176.79 & 499.43 & 1045.29 & 1986.70 & 3269.36 & 4995.96 \\
\hspace{2mm} Relative $n$ in column & 1.08 & 6.93 & 12.17 & 17.15 & 19.87 & 23.64 & 25.67 & 28.63 \\
$\nabla f$, forward AD\\
\hspace{2mm} Relative $n=1$ & 1.34 & 13.69 & 51.54 & 132.33 & 251.32 & 469.84 & 815.55 & 1342.07\\
\hspace{2mm} Relative $n$ in column & 1.34 & 2.66 & 3.55 & 4.54 & 4.77 & 5.59 & 6.40 & 7.69 \\
$\nabla f$, reverse AD\\
\hspace{2mm} Relative $n=1$ & 1.52 & 11.12 & 31.37 & 67.27 & 113.99 & 174.62 & 254.15 & 342.33 \\
\hspace{2mm} Relative $n$ in column & 1.52 & 2.16 & 2.16 & 2.31 & 2.16 & 2.07 & 1.99 & 1.96 \\
\bottomrule
\end{tabularx}}
\end{table}
%Second order descent methods, Newton, Gauss-Newton, Levenberg-Marquardt
Second-order methods based on Newton's method make use of both the gradient $\nabla f$ and the Hessian $\mathbf{H}_f$, working via updates of the form $\Del \mathbf{w} = -\eta\,\mathbf{H}^{-1}_f \nabla f$ and providing significantly faster convergence \citep{Press2007}. AD provides a way of automatically computing the exact Hessian, enabling succinct and convenient general-purpose implementations.\footnote{See \url{http://DiffSharp.github.io/DiffSharp/examples-newtonsmethod.html} for an implementation of Newton's method with the full Hessian.} Newton's method converges in fewer iterations, but this comes at the cost of having to compute $\mathbf{H}_f$ in each iteration. In large-scale problems, the Hessian is usually replaced by a numerical approximation using first-order updates from gradient evaluations, giving rise to quasi-Newton methods. A highly popular such method is the BFGS\footnote{After Broyden–Fletcher–Goldfarb–Shanno, who independently discovered the method in the 1970s.} algorithm, together with its limited-memory variant L-BFGS \citep{Dennis1996}. On the other hand, Hessians arising in large-scale applications are typically sparse. This sparsity along with symmetry can be readily exploited by AD techniques such as computational graph elimination \citep{Dixon1991}, partial separability \citep{Gay1996}, and matrix coloring and compression \citep{Gebremedhin2009}.
In many cases one does not need the full Hessian but only a Hessian--vector product $\mathbf{H} \mathbf{v}$, which can be computed efficiently using a reverse-on-forward configuration of AD by applying the reverse mode to take the gradient of code produced by the forward mode.\footnote{\citet{Christianson2012ALN} demonstrates that the second derivative can be computed with the same arithmetic operation sequence using forward-on-reverse, reverse-on-forward, and reverse-on-reverse. The taping overheads of these methods may differ in implementation-dependent ways.} Given the function $f:\mathbb{R}^n\to\mathbb{R}$, the evaluation point $\mathbf{x}$, and the vector $\mathbf{v}$, one can accomplish this by first computing the directional derivative $\nabla f \cdot \mathbf{v}$ through the forward mode via setting $\dot{\mathbf{x}}=\mathbf{v}$ and then applying the reverse mode on this result to get $\nabla^2 f \cdot \mathbf{v}=\mathbf{H}_f \mathbf{v}$ \citep{Pearlmutter1994}. This computes $\mathbf{H} \mathbf{v}$ with $O(n)$ complexity, even though $\mathbf{H}$ is a $n \times n$ matrix. Availability of robust AD tools may make more sophisticated optimization methods applicable to large-scale machine-learning problems. For instance, when fast stochastic Hessian--vector products are available, these can be used as the basis of stochastic Newton's methods \citep{agarwal2016second}, which have the potential to endow stochastic optimization with quadratic convergence.
Another approach for improving the rate of convergence of gradient-based methods is to use gain adaptation methods such as stochastic meta-descent (SMD) \citep{Schraudolph1999}, where stochastic sampling is introduced to avoid local minima and reduce the computational expense. An example using SMD with AD Hessian--vector products is given by \citet{Vishwanathan2006} on conditional random fields (CRF). Similarly, \citet{Schraudolph2003} use Hessian--vector products in their model combining conjugate gradient techniques with stochastic gradient descent.
\subsection{Neural Networks, Deep Learning, Differentiable Programming}
Training of a neural network is an optimization problem with respect to its set of weights, which can in principle be addressed by using any method ranging from evolutionary algorithms \citep{such2017deep} to gradient-based methods such as BFGS \citep{Apostolopoulou2009} or the mainstay stochastic gradient descent \citep{bottou2010large} and its many variants \citep{kingma2015adam,tieleman2012lecture,duchi2011adaptive}. As we have seen, the backpropagation algorithm is only a special case of AD: by applying reverse mode AD to an objective function evaluating a network's error as a function of its weights, we can readily compute the partial derivatives needed for performing weight updates.\footnote{See \url{http://DiffSharp.github.io/DiffSharp/examples-neuralnetworks.html} for an implementation of backpropagation with reverse mode AD.}
The LUSH system \citep{LUSH2002}, and its predecessor SN \citep{bottou-lecun-88}, were the first production systems that targeted efficient neural network simulation while incorporating both a general-purpose programming language and AD. Modern deep learning frameworks provide differentiation capability in one way or another, but the underlying mechanism is not always made clear and confusion abounds regarding the use of the terms ``autodiff'', ``automatic differentiation'', and ``symbolic differentiation'', which are sometimes even used interchangeably. In mainstream frameworks including Theano\footnote{Theano is a computational graph optimizer and compiler with GPU support and it currently handles derivatives in a highly optimized form of symbolic differentiation. The result can be interpreted as a hybrid of symbolic differentiation and reverse mode AD, but Theano does not use the general-purpose reverse accumulation as we describe in this paper. (Personal communication with the authors.)} \citep{Bastien2012}, TensorFlow \citep{abadi2016tensorflow}, Caffe \citep{jia2014caffe}, and CNTK \citep{seide2016cntk} the user first constructs a model as a computational graph using a domain-specific mini language, which then gets interpreted by the framework during execution. This approach has the advantage of enabling optimizations of the computational graph structure (e.g., as in Theano), but the disadvantages of having limited and unintuitive control flow and being difficult to debug. In contrast, the lineage of recent frameworks led by autograd \citep{maclaurin2016modeling}, Chainer \citep{tokui2015chainer}, and PyTorch \citep{paszke2017automatic} provide truly general-purpose reverse mode AD of the type we outline in Section~\ref{SectionPreliminaries}, where the user directly uses the host programming language to define the model as a regular program of the forward computation. This eliminates the need for an interpreter, allows arbitrary control flow statements, and makes debugging simple and intuitive.
Simultaneously with the ongoing adoption of general-purpose AD in machine learning, we are witnessing a modeling-centric terminology emerge within the deep learning community. The terms \emph{define-and-run} and \emph{static computational graph} refer to Theano-like systems where a model is constructed, before execution, as a computational graph structure, which later gets executed with different inputs while remaining fixed. In contrast, the terms \emph{define-by-run} and \emph{dynamic computational graph} refer to the general-purpose AD capability available in newer PyTorch-like systems where a model is a regular program in the host programming language, whose execution dynamically constructs a computational graph on-the-fly that can freely change in each iteration.\footnote{Note that the terms ``static'' and ``dynamic'' here are used in the sense of having a fixed versus non-fixed computational graph topology and not in the sense of data flow architectures.}
\emph{Differentiable programming}\footnote{A term advocated by Christopher Olah (\url{http://colah.github.io/posts/2015-09-NN-Types-FP/}), David Dalrymple (\url{https://www.edge.org/response-detail/26794}), and Yann LeCun (\url{https://www.facebook.com/yann.lecun/posts/10155003011462143}) from a deep learning point of view. Note the difference from \emph{differential} dynamic programming \citep{jacobson1970differential} in optimal control.} is another emerging term referring to the realization that deep learning practice essentially amounts to writing program templates of potential solutions to a problem, which are constructed as differentiable directed graphs assembled from functional blocks whose parameters are learned from examples using gradient-based optimization. Expressed in this paradigm, neural networks are just a class of parameterized differentiable programs composed of building blocks such as feed-forward, convolutional, and recurrent elements. We are increasingly seeing these traditional building blocks freely composed in arbitrary algorithmic structures using control flow, as well as the introduction of novel differentiable architectures such as the neural Turing machine \citep{graves2014neural}, a range of controller--interface abstractions \citep{graves2016hybrid,zaremba2016learning,joulin2015inferring,sukhbaatar2015end}, and differentiable versions of data structures such as stacks, queues, deques \citep{grefenstette2015learning}. Availability of general-purpose AD greatly simplifies the implementation of such architectures by enabling their expression as regular programs that rely on the differentiation infrastructure. Although the differentiable programming perspective on deep learning is new, we note that programming with differentiable functions and having differentiation as a language infrastructure has been the main research subject of the AD community for many decades and realized in a wide range of systems and languages as we shall see in Section~\ref{SectionImplementations}.
There are instances in neural network literature---albeit few---where explicit reference has been made to AD for computing error gradients, such as \citet{Eriksson1998} using AD for large-scale feed-forward networks, and the work by \citet{Yang2008}, where the authors use AD to train a neural-network-based proportional-integral-derivative (PID) controller. Similarly, \citet{Rollins2009} uses reverse mode AD in conjunction with neural networks for the problem of optimal feedback control. Another example is given for continuous time recurrent neural networks (CTRNN) by \citet{AlSeyab2008}, where the authors apply AD for the training of CTRNNs predicting dynamic behavior of nonlinear processes in real time and report significantly reduced training time compared with other methods.
\subsection{Computer Vision}
Since the influential work by \citet{krizhevsky2012imagenet}, computer vision has been dominated by deep learning, specifically, variations of convolutional neural networks \citep{lecun1998gradient}. These models are trained end-to-end, meaning that a mapping from raw input data to corresponding outputs is learned, automatically discovering the representations needed for feature detection in a process called representation learning \citep{bengio2013representation}.
Besides deep learning, an interesting area where AD can be applied to computer vision problems is inverse graphics \citep{horn1977understanding,hinton1997generative}---or analysis-by-synthesis \citep{yildirim2015efficient}---where vision is seen as the inference of parameters for a generative model of a scene. Using gradient-based optimization in inverse graphics requires propagating derivatives through whole image synthesis pipelines including the renderer. \citet{eslami2016attend} use numerical differentiation for this purpose. \citet{loper2014opendr} implement the Open Differentiable Renderer (OpenDR), which is a scene renderer that also supplies derivatives of the image pixels with respect to scene parameters, and demonstrate it in the task of fitting an articulated and deformable 3D model of the human body to image and range data from a Kinect device. Similarly, \citet{kulkarni2015picture} implement a differentiable approximate renderer for the task of inference in probabilistic programs describing scenes.
\citet{srajer2016benchmark} investigate the use of AD for three tasks in computer vision and machine learning, namely bundle adjustment \citep{triggs1999bundle}, Gaussian mixture model fitting, and hand tracking \citep{taylor2014user}, and provide a comprehensive benchmark of various AD tools for the computation of derivatives in these tasks.
\citet{Pock2007} make use of AD in addressing the problems of denoising, segmentation, and recovery of information from stereoscopic image pairs, and note the usefulness of AD in identifying sparsity patterns in large Jacobian and Hessian matrices. In another study, \citet{Grabner2008} use reverse mode AD for GPU-accelerated medical 2D/3D registration, a task involving the alignment of data from different sources such as X-ray images or computed tomography. The authors report a six-fold increase in speed compared with numerical differentiation using center difference (cf.\ our benchmark with the Helmholtz function, Figure~\ref{FigureHelmholtz} and Table~\ref{TableHelmholtz}).
\citet{Barrett2013} present a use of general-purpose AD for the task of video event detection using hidden Markov models (HMMs) and \citet{Dalal2005} object detectors, performing training on a corpus of pre-tracked video using an adaptive step size gradient descent with reverse mode AD. Initially implemented with the R6RS-AD package\footnote{\url{https://github.com/qobi/R6RS-AD}} which provides forward and reverse mode AD in Scheme, the resulting gradient code was later ported to C and highly optimized.\footnote{Personal communication.}
\subsection{Natural Language Processing}
Natural language processing (NLP) constitutes one of the areas where rapid progress is being made by applying deep learning techniques \citep{goldberg2016primer}, with applications in tasks including machine translation \citep{bahdanau2014neural}, language modeling \citep{mikolov2010recurrent}, dependency parsing \citep{chen2014fast}, and question answering \citep{pmlr-v48-kumar16}. Besides deep learning approaches, statistical models in NLP are commonly trained using general purpose or specialized gradient-based methods and mostly remain expensive to train. Improvements in training time can be realized by using online or distributed training algorithms \citep{Gimpel2010}. An example using stochastic gradient descent for NLP is given by \citet{Finkel2008} optimizing conditional random field parsers through an objective function. Related with the work on video event detection in the previous section, \citet{Yu2013} report their work on sentence tracking, representing an instance of grounded language learning paired with computer vision, where the system learns word meanings from short video clips paired with descriptive sentences. The method uses HMMs to represent changes in video frames and meanings of different parts of speech. This work is implemented in C and computes the required gradients using AD through the ADOL-C tool.\footnote{An implementation of the sentence tracker applied to video search using sentence-based queries can be accessed online: \url{http://upplysingaoflun.ecn.purdue.edu/~qobi/cccp/sentence-tracker-video-retrieval.html}}
\subsection{Probabilistic Modeling and Inference}
Inference in probabilistic models can be static, such as compiling a given model to Bayesian networks and using algorithms such as belief propagation for inference; or they can be dynamic, executing a model forward many times and computing statistics on observed values to infer posterior distributions. Markov chain Monte Carlo (MCMC) \citep{Neal1993} methods are often used for dynamic inference, such as the Metropolis--Hastings algorithm based on random sampling \citep{Chib1995}. \citet{Meyer2003} give an example of how AD can be used to speed up Bayesian posterior inference in MCMC, with an application in stochastic volatility. Amortized inference \citep{gershman2014amortized,stuhlmuller2013learning} techniques based on deep learning \citep{le2016inference,ritchie2016deep} work by training neural networks for performing approximate inference in generative models defined as probabilistic programs \citep{gordon2014probabilistic}.
When model parameters are continuous, the Hamiltonian---or, hybrid---Monte Carlo (HMC) algorithm provides improved convergence characteristics avoiding the slow exploration of random sampling, by simulating Hamiltonian dynamics through auxiliary ``momentum variables'' \citep{Duane1987}. The advantages of HMC come at the cost of requiring gradient evaluations of complicated probability models. AD is highly suitable here for complementing probabilistic modeling, because it relieves the user from the manual derivation of gradients for each model.\footnote{See \url{http://diffsharp.github.io/DiffSharp/examples-hamiltonianmontecarlo.html} for an implementation of HMC with reverse mode AD.} For instance, the probabilistic programming language Stan \citep{carpenter2016stan} implements automatic Bayesian inference based on HMC and the No-U-Turn sampler (NUTS) \citep{Hoffman2014} and uses reverse mode AD for the calculation of gradients for both HMC and NUTS \citep{carpenter2015stan}. Similarly, \citet{Wingate2011} demonstrate the use of AD as a non-standard interpretation of probabilistic programs enabling efficient inference algorithms. \citet{kucukelbir2017automatic} present an AD-based method for deriving variational inference (VI) algorithms.
PyMC3 \citep{salvatier2016probabilistic} allows fitting of Bayesian models using MCMC and VI, for which it uses gradients supplied by Theano. Edward \citep{tran2016edward} is a library for deep probabilistic modeling, inference, and criticism \citep{tran2017deep} that supports VI using TensorFlow. Availability of general-purpose AD in this area has enabled new libraries such as Pyro\footnote{\url{http://pyro.ai/}} and ProbTorch \citep{siddharth2017learning} for deep \emph{universal} probabilistic programming with support for recursion and control flow, relying, in both instances, on VI using gradients supplied by PyTorch's reverse mode AD infrastructure.
When working with probabilistic models, one often needs to backpropagate derivatives through sampling operations of random variables in order to achieve stochastic optimization of model parameters. The score-function estimator, or REINFORCE \citep{williams1992simple}, method provides a generally applicable unbiased gradient estimate, albeit with high variance. When working with continuous random variables, one can substitute a random variable by a deterministic and differentiable transformation of a simpler random variable, a method known as the ``reparameterization trick'' \citep{williams1992simple,kingma2014auto,rezende2014stochastic}. For discrete variables, the REBAR \citep{tucker2017rebar} method provides a lower-variance unbiased gradient estimator by using continuous relaxation. A generalization of REBAR called RELAX \citep{grathwohl2017backpropagation} works by learning a free-form control variate parameterized by a neural network and is applicable in both discrete and continuous settings.
\section{Implementations}
\label{SectionImplementations}
It is useful to have an understanding of the different ways in which AD can be implemented. Here we cover major implementation strategies and provide a survey of existing tools.
A principal consideration in any AD implementation is the performance overhead introduced by the AD arithmetic and bookkeeping. In terms of computational complexity, AD guarantees that the amount of arithmetic goes up by no more than a small constant factor \citep{Griewank2008}. On the other hand, managing this arithmetic can introduce a significant overhead if done carelessly. For instance, naïvely allocating data structures for holding dual numbers will involve memory access and allocation for every arithmetic operation, which are usually more expensive than arithmetic operations on modern computers.\footnote{The implementation of forward mode in Julia \citep{revels2016forward} attempts to avoid this, and some current compilers can avoid this expense by unboxing dual numbers \citep{leroy1997effectiveness, jones1993glasgow, jones1991unboxed, siskind2016efficient}. This method is also used to reduce the memory-access overhead in the implementations of forward mode in Stalingrad and the Haskell \emph{ad} library.} Likewise, using operator overloading may introduce method dispatches with attendant costs, which, compared to raw numerical computation of the original function, can easily amount to a slowdown of an order of magnitude.\footnote{Flow analysis \citep{shivers1991control} and/or partial evaluation \citep{jones1993partial}, together with tag stripping \citep{appel1989runtime, peterson1989untagged}, can remove this method dispatch. These, together with unboxing, can often make it possible to completely eliminate the memory access, memory allocation, memory reclamation, and method dispatch overhead of dual numbers \citep{siskind2016efficient}.}
Another major issue is the risk of hitting a class of bugs called ``perturbation confusion'' \citep{SiskindPearlmutter2005a,manzyuk2012confusion}. This essentially means that if two ongoing differentiations affect the same piece of code, the two formal epsilons they introduce (Section~\ref{SectionDualNumbers}) need to be kept distinct. It is very easy to have bugs---particularly in performance-oriented AD implementations---that confuse these in various ways. Such situations can also arise when AD is nested, that is, derivatives are computed for functions that internally compute derivatives.
Translation of mathematics into computer code often requires attention to numeric issues. For instance, the mathematical expressions $\log(1 + x)$ or $\sqrt{x^2+y^2+z^2}$ or $\tan^{-1}(y/x)$ should not be naïvely translated, but rather expressed as \texttt{log1p(x)}, \texttt{hypot(x,hypot(y,z))}, and \texttt{atan2(y,x)}. In machine learning, the most prominent example of this is probably the so-called log-sum-exp trick to improve the numerics of calculations of the form $\log\sum_i \exp x_i$. AD is not immune to such numeric considerations. For example, code calculating $E=\sum_i E_i$, processed by AD, will calculate $\nabla_w E=\sum_i\nabla_w E_i$. If the system is seeking a local minimum of $E$ then $\nabla_w E = \sum_i\nabla_w E_i \rightarrow_t 0$, and naïvely adding a set of large numbers whose sum is near zero is numerically fraught. This is to say that AD is not immune to the perils of floating point arithmetic, and can sometimes introduce numeric issues which were not present in the primal calculation. Issues of numeric analysis are outside our present scope, but there is a robust literature on the numerics of AD (e.g., \citet{griewank2012numerical}) involving using subgradients to allow optimization to proceed despite non-differentiability of the objective, appropriate subgradients and approximations for functions like $\lvert \cdot \rvert$ and $\lVert\cdot\rVert_2$ and $\sqrt{\cdot}$ near zero, and a spate of related issues.
One should also be cautious about approximated functions and AD \citep{sirkes-tziperman-1997a}. In this case, if one has a procedure \emph{approximating} an ideal function, AD always gives the derivative of the procedure that was actually programmed, which may not be a good approximation of the derivative of the ideal function that the procedure was approximating. For instance, consider $e^x$ computed by a piecewise-rational approximation routine. Using AD on this routine would produce an approximated derivative in which each piece of the piecewise formula will get differentiated. Even if this would remain an approximation of the derivative of $e^x$, we know that $\frac{de^x}{dx} = e^x$ and the original approximation itself was already a better approximation for the derivative of $e^x$.\footnote{In modern systems this is not an issue, because $e^x$ is a primitive implemented in hardware.} Users of AD implementations must be therefore cautious to \emph{approximate the derivative, not differentiate the approximation}. This would require explicitly approximating a known derivative, in cases where a mathematical function can only be computed approximately but has a well-defined mathematical derivative.
We note that there are similarities as well as differences between machine learning workloads and those studied in the traditional AD literature \citep{baydin2016tricks}. Deep learning systems are generally compute-bound and spend a considerable amount of computation time in highly-optimized numerical kernels for matrix operations \citep{hadjis2015caffe,chetlur2014cudnn}. This is a situation which is arguably amenable to operator-overloading-based AD implementations on high-level operations, as is commonly found in current machine learning frameworks. In contrast, numerical simulation workloads in traditional AD applications can be bandwidth-bound, making source code transformation and compiler optimization approaches more relevant. Another difference worth noting is that whereas high numerical precision is desirable in traditional application domains of AD such as computational fluid dynamics \citep{cohen2009fast}, in deep learning lower-precision is sufficient and even desirable in improving computational efficiency, thanks to the error resiliency of neural networks \citep{gupta2015deep,courbariaux2015binaryconnect}.
There are instances in recent literature where implementation-related experience from the AD field has been put to use in machine learning settings. One particular area of recent interest is implicit and iterative AD techniques \citep{Griewank2008}, which has found use in work incorporating constrained optimization within deep learning \citep{amos2017optnet} and probabilistic graphical models and neural networks \citep{johnson2016composing}. Another example is checkpointing strategies \citep{Dauvergne2006,siskind2017divide}, which allow balancing of application-specific trade-offs between time and space complexities of reverse mode AD by not storing the full tape of intermediate variables in memory and reconstructing these as needed by re-running parts of the forward computation from intermediate checkpoints. This is highly relevant in deep learning workloads running on GPUs with limited memory budgets. A recent example in this area is the work by \citet{gruslys2016memory}, where the authors construct a checkpointing variety of the backpropagation through time (BPTT) algorithm for recurrent neural networks and demonstrate it saving up to 95\% memory usage at the cost of a 33\% increase in computation time in one instance.
In Table~\ref{TableADImplementations} we present a review of notable general-purpose AD implementations.\footnote{Also see the website \url{http://www.autodiff.org/} for a list of tools maintained by the AD community.} A thorough taxonomy of implementation techniques was introduced by \citet{Juedes1991}, which was later revisited by \citet{Bischof2008} and simplified into \emph{elemental}, \emph{operator overloading}, \emph{compiler-based}, and \emph{hybrid} methods. We adopt a similar classification for the following part of this section.
% Remarks concerning further implementation methods in \citet{Gay2006}.
\addtolength{\tabcolsep}{-3pt}
\begin{sidewaystable}
\centering
\renewcommand{\arraystretch}{1.5}
\caption{Survey of AD implementations. Tools developed primarily for machine learning are highlighted in bold.}
\label{TableADImplementations}
{\tiny
\begin{tabularx}{\textwidth}{@{}p{12mm}p{20mm}p{5mm}p{6mm}p{62mm}p{38mm}p{65mm}@{}}
\toprule
Language & Tool & Type & Mode & Institution / Project & Reference & URL\\
\midrule
AMPL & AMPL & INT & F, R & Bell Laboratories & \citet{Fourer2002} & \tiny\url{http://www.ampl.com/}\\
C, C++ & ADIC & ST & F, R & Argonne National Laboratory & \citet{Bischof1997} & \tiny\url{http://www.mcs.anl.gov/research/projects/adic/}\\
& ADOL-C & OO & F, R & Computational Infrastructure for Operations Research & \citet{Walther2012} & \tiny\url{https://projects.coin-or.org/ADOL-C}\\
C++ & Ceres Solver & LIB & F & Google & & \tiny\url{http://ceres-solver.org/}\\
& CppAD & OO & F, R & Computational Infrastructure for Operations Research & \citet{Bell2008} & \tiny\url{http://www.coin-or.org/CppAD/}\\
& FADBAD++ & OO & F, R & Technical University of Denmark & \citet{Bendtsen1996} & \tiny\url{http://www.fadbad.com/fadbad.html}\\
& Mxyzptlk & OO & F & Fermi National Accelerator Laboratory & \citet{Ostiguy2007} & \\
C\# & AutoDiff & LIB & R & George Mason Univ., Dept. of Computer Science & \citet{Shtof2013} & \tiny\url{http://autodiff.codeplex.com/}\\
F\#, C\# & \textbf{DiffSharp} & OO & F, R & Maynooth University, Microsoft Research Cambridge & \citet{baydin2016diffsharp} & \tiny\url{http://diffsharp.github.io}\\
Fortran & ADIFOR & ST & F, R & Argonne National Laboratory & \citet{Bischof1996} & \tiny\url{http://www.mcs.anl.gov/research/projects/adifor/}\\
& NAGWare & COM & F, R & Numerical Algorithms Group & \citet{Naumann2005} & \tiny\url{http://www.nag.co.uk/nagware/Research/ad_overview.asp}\\
& TAMC & ST & R & Max Planck Institute for Meteorology & \citet{Giering1998} & \tiny\url{http://autodiff.com/tamc/}\\
Fortran, C & COSY & INT & F & Michigan State Univ., Biomedical and Physical Sci. & \citet{Berz1996} & \tiny\url{http://www.bt.pa.msu.edu/index_cosy.htm}\\
& Tapenade & ST & F, R & INRIA Sophia-Antipolis & \citet{Hascoet2013} & \tiny\url{http://www-sop.inria.fr/tropics/tapenade.html}\\
Haskell & ad & OO & F, R & Haskell package & & \tiny\url{http://hackage.haskell.org/package/ad}\\
Java & ADiJaC & ST & F, R & University Politehnica of Bucharest & \citet{slusanschi2016adijac} & \tiny\url{http://adijac.cs.pub.ro}\\
& Deriva & LIB & R & Java \& Clojure library & & \tiny\url{https://github.com/lambder/Deriva}\\
Julia & JuliaDiff & OO & F, R & Julia packages & \citet{RevelsLubinPapamarkou2016} & \tiny\url{http://www.juliadiff.org/}\\
Lua & \textbf{torch-autograd} & OO & R & Twitter Cortex & & \tiny\url{https://github.com/twitter/torch-autograd}\\
MATLAB & ADiMat & ST & F, R & Technical University of Darmstadt, Scientific Comp. & \citet{Willkomm2013} & \tiny\url{http://adimat.sc.informatik.tu-darmstadt.de/}\\
& INTLab & OO & F & Hamburg Univ. of Technology, Inst. for Reliable Comp. & \citet{Rump1999} & \tiny\url{http://www.ti3.tu-harburg.de/rump/intlab/}\\
& TOMLAB/MAD & OO & F & Cranfield University \& Tomlab Optimization Inc. & \citet{Forth2006} & \tiny\url{http://tomlab.biz/products/mad}\\
Python & ad & OO & R & Python package & & \tiny\url{https://pypi.python.org/pypi/ad}\\
& \textbf{autograd} & OO & F, R & Harvard Intelligent Probabilistic Systems Group & \citet{maclaurin2016modeling} & \tiny\url{https://github.com/HIPS/autograd}\\
& \textbf{Chainer} & OO & R & Preferred Networks & \citet{tokui2015chainer} & \tiny\url{https://chainer.org/}\\
& \textbf{PyTorch} & OO & R & PyTorch core team & \citet{paszke2017automatic} & \tiny\url{http://pytorch.org/}\\
& \textbf{Tangent} & ST & F, R & Google Brain & \citet{van2017tangent} & \tiny\url{https://github.com/google/tangent}\\
Scheme & R6RS-AD & OO & F, R & Purdue Univ., School of Electrical and Computer Eng. & & \tiny\url{https://github.com/qobi/R6RS-AD}\\
& Scmutils & OO & F & MIT Computer Science and Artificial Intelligence Lab. & \citet{Sussman2001} & \tiny\url{http://groups.csail.mit.edu/mac/users/gjs/6946/refman.txt}\\
& Stalingrad & COM & F, R & Purdue Univ., School of Electrical and Computer Eng. & \citet{pearlmutter2008reverse} & \tiny\url{http://www.bcl.hamilton.ie/~qobi/stalingrad/}\\
% Odyssee reverse mode source transformation (predecessor of Tapenade)
% TAMC reverse mode source transformation
\bottomrule
\addlinespace
\multicolumn{7}{l}{F: Forward, R: Reverse; COM: Compiler, INT: Interpreter, LIB: Library, OO: Operator overloading, ST: Source transformation}
\end{tabularx}
}
\end{sidewaystable}
\addtolength{\tabcolsep}{3pt}
\subsection{Elemental Libraries}
These implementations form the most basic category and work by replacing mathematical operations with calls to an AD-enabled library. Methods exposed by the library are then used in function definitions, meaning that the decomposition of any function into elementary operations is done manually when writing the code.
The approach has been utilized since the early days of AD, with prototypical examples being the WCOMP and UCOMP packages of \citet{Lawson1971}, the APL package of \citet{Neidinger1989}, and the work by \citet{Hinkins1994}. Likewise, \citet{Hill1992} formulate their implementation of AD in MATLAB using elemental methods.
Elemental libraries still constitute the simplest strategy to implement AD for languages without operator overloading.
\subsection{Compilers and Source Code Transformation}
These implementations provide extensions to programming languages that automate the decomposition of algorithms into AD-enabled elementary operations. They are typically executed as preprocessors\footnote{Preprocessors transform program source code before it is given as an input to a compiler.} to transform the input in the extended language into the original language.
Classical instances of source code transformation include the Fortran preprocessors GRESS \citep{Horwedel1988} and PADRE2 \citep{Kubo1990}, which transform AD-enabled variants of Fortran into standard Fortran 77 before compiling. Similarly, the ADIFOR tool \citep{Bischof1996}, given a Fortran source code, generates an augmented code in which all specified partial derivatives are computed in addition to the original result. For procedures coded in ANSI C, the ADIC tool \citep{Bischof1997} implements AD as a source code transformation after the specification of dependent and independent variables. A recent and popular tool also utilizing this approach is Tapenade \citep{Pascual2008,Hascoet2013}, implementing forward and reverse mode AD for Fortran and C programs. Tapenade itself is implemented in Java and can be run locally or as an online service.\footnote{\url{http://www-tapenade.inria.fr:8080/tapenade/index.jsp}}
In addition to language extensions through source code transformation, there are implementations introducing new languages with tightly integrated AD capabilities through special-purpose compilers or interpreters. Some of the earliest AD tools such as SLANG \citep{Adamson1969} and PROSE \citep{Pfeiffer1987} belong to this category. The NAGWare Fortran 95 compiler \citep{Naumann2005} is a more recent example, where the use of AD-related extensions triggers automatic generation of derivative code at compile time.
As an example of interpreter-based implementation, the algebraic modeling language AMPL \citep{Fourer2002} enables objectives and constraints to be expressed in mathematical notation, from which the system deduces active variables and arranges the necessary AD computations. Other examples in this category include the FM/FAD package \citep{Mazourik1991}, based on the Algol-like DIFALG language, and the object-oriented COSY language \citep{Berz1996} similar to Pascal.
The Stalingrad compiler \citep{pearlmutter2008reverse,Siskind2008}, working on the Scheme-based AD-aware VLAD language, also falls under this category. The newer DVL compiler\footnote{\url{https://github.com/axch/dysvunctional-language}} is based on Stalingrad and uses a reimplementation of portions of the VLAD language.
Motivated by machine learning applications, the Tangent library \citep{van2017tangent} implements AD using source code transformation, and accepts numeric functions written in a syntactic subset of Python and Numpy.
\subsection{Operator Overloading}
In modern programming languages with polymorphic features, operator overloading provides the most straightforward way of implementing AD, exploiting the capability of redefining elementary operation semantics.
A popular tool implemented with operator overloading in C++ is ADOL-C \citep{Walther2012}. ADOL-C requires the use of AD-enabled types for variables, and records arithmetic operations on variables in tape data structures, which can subsequently be ``played back'' during reverse mode AD computations. The Mxyzptlk package \citep{Michelotti1990} is another example for C++ capable of computing arbitrary-order partial derivatives via forward propagation. The FADBAD++ library \citep{Bendtsen1996} implements AD for C++ using templates and operator overloading. For Python, the \emph{ad} package\footnote{\url{http://pythonhosted.org/ad/}} uses operator overloading to compute first- and second-order derivatives, while the newer autograd package\footnote{\url{https://github.com/HIPS/autograd}} provides forward and reverse mode AD with support for higher-order derivatives.
For functional languages, examples include R6RS-AD\footnote{\url{https://github.com/NUIM-BCL/R6RS-AD}} and the AD routines within the Scmutils library\footnote{\url{http://groups.csail.mit.edu/mac/users/gjs/6946/refman.txt}} for Scheme, the \emph{ad} library\footnote{\url{http://hackage.haskell.org/package/ad}} for Haskell, and DiffSharp\footnote{\url{http://diffsharp.github.io}} for F\# and C\#.
%\subsection{Parallelization}
% Discuss advantages of AD for parallelization.
%Parallel implementation of Hessian calculations, automated generation of parallel code \citep{Bucker2008}
%\citet{Bischof2008} give an example of parallel reverse mode AD for a plasma simulation code
\section{Conclusions}
\label{SectionConclusions}
Backpropagation and gradient-based optimization are behind virtually all recent successes in machine learning, yielding state-of-the-art results in computer vision, speech recognition and synthesis, and machine translation. We expect these techniques to remain at the core of machine learning for the foreseeable future. Research in the field involves a rapid prototyping and development cycle for testing new models and ideas, using a collection of increasingly higher-quality machine learning frameworks. These frameworks are in the process of transition from coarse-grained (module level) backpropagation towards fine-grained, general-purpose AD, allowing models to be implemented as regular programs in general-purpose programming languages with differentiation as an integral part of the infrastructure. We strongly believe that general-purpose AD is the future of gradient-based machine learning and we expect it to become an indispensable tool in the machine learning toolbox.
It is an exciting time for working at the intersection of AD and machine learning, and there are many opportunities for bringing advanced techniques and expertise from AD literature to bear on machine learning problems. Techniques that have been developed by the AD community such as tape reduction and elimination \citep{naumann2004optimal}, fixed-point iterations \citep{christianson1994reverse}, utilizing sparsity by matrix coloring \citep{Gebremedhin2009,gebremedhin2013colpack}, and reverse AD checkpointing \citep{Dauvergne2006} are just a few examples that can find potential use in machine learning for increasing performance, improving convergence of optimization, using hardware more efficiently, and even enabling new types of machine learning models to be implemented. Similarly, exciting new AD modes like direct propagation of the inverse Jacobian \citep{srinivasan-todorov-2015a} have emerged from the machine learning community, but have yet to be examined and formalized by the AD community.
An important direction for future work is to make use of nested AD techniques in machine learning, allowing differentiation to be nested arbitrarily deep with referential transparency \citep{Siskind2008b,pearlmutter2008reverse}. Nested AD is highly relevant in hyperparameter optimization as it can effortlessly provide exact hypergradients, that is, derivatives of a training objective with respect to the hyperparameters of an optimization routine \citep{Maclaurin2015,baydin2017online}. Potential applications include Bayesian model selection \citep{Rasmussen2006} and gradient-based tuning of Hamiltonian Monte Carlo step sizes and mass matrices \citep{Salimans2014}. Besides hyperparameters, models internally using higher-order derivatives constitute a straightforward usage case for nested AD. The Riemannian manifold Langevin and Hamiltonian Monte Carlo methods \citep{Girolami2011} use higher-order derivative information to more closely track the information geometry of the sampled distribution for faster convergence and exploration. In neural networks, it is very natural to use nested derivatives in defining objective functions that take input transformations into account, such as the Tangent Prop method \citep{Simard1998} for imposing invariance under a set of chosen transformations.
% Acknowledgements should go at the end, before appendices and references
\acks{We thank the anonymous reviewers whose comments helped improve this manuscript. This work was supported, in part, by Science Foundation Ireland grant 09/IN.1/I2637, by the Army Research Laboratory, accomplished under Cooperative Agreement Number W911NF-10-2-0060, by the National Science Foundation under Grants 1522954-IIS and 1734938-IIS, and by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/Interior Business Center (DOI/IBC) contract number D17PC00341. Any opinions, findings, views, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views, official policies, or endorsements, either expressed or implied, of the sponsors. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes, notwithstanding any copyright notation herein.
}
% Manual newpage inserted to improve layout of sample file - not
% needed in general before appendices/bibliography.
%
% \newpage
%
% \appendix
% \section*{Appendix A.}
% \label{app:theorem}
% Note: in this sample, the section number is hard-coded in. Following
% proper LaTeX conventions, it should properly be coded as a reference:
%In this appendix we prove the following theorem from
%Section~\ref{sec:textree-generalization}:
\vskip 0.2in
\bibliography{17-468}
\end{document}
doc/Articles/AutodiffArticle → doc/Articles/Autodiff/adiff.tar
View File
doc/Articles/Autodiff/figures/approx-error/approx-error.tex
+2189
View File
File diff suppressed because it is too large Load Diff
doc/Articles/Autodiff/figures/backprop/backprop-eps-converted-to.pdf
BIN
View File
Binary file not shown.
doc/Articles/Autodiff/figures/backprop/backprop.eps
+650
View File
@@ -0,0 +1,650 @@
%!PS-Adobe-3.0 EPSF-3.0
%%Creator: cairo 1.8.8 (http://cairographics.org)
%%CreationDate: Tue Jan 23 23:31:51 2018
%%Pages: 1
%%BoundingBox: 1 1 377 199
%%DocumentData: Clean7Bit
%%LanguageLevel: 2
%%EndComments
%%BeginProlog
/cairo_eps_state save def
/dict_count countdictstack def
/op_count count 1 sub def
userdict begin
/q { gsave } bind def
/Q { grestore } bind def
/cm { 6 array astore concat } bind def
/w { setlinewidth } bind def
/J { setlinecap } bind def
/j { setlinejoin } bind def
/M { setmiterlimit } bind def
/d { setdash } bind def
/m { moveto } bind def
/l { lineto } bind def
/c { curveto } bind def
/h { closepath } bind def
/re { exch dup neg 3 1 roll 5 3 roll moveto 0 rlineto
0 exch rlineto 0 rlineto closepath } bind def
/S { stroke } bind def
/f { fill } bind def
/f* { eofill } bind def
/B { fill stroke } bind def
/B* { eofill stroke } bind def
/n { newpath } bind def
/W { clip } bind def
/W* { eoclip } bind def
/BT { } bind def
/ET { } bind def
/pdfmark where { pop globaldict /?pdfmark /exec load put }
{ globaldict begin /?pdfmark /pop load def /pdfmark
/cleartomark load def end } ifelse
/BDC { mark 3 1 roll /BDC pdfmark } bind def
/EMC { mark /EMC pdfmark } bind def
/cairo_store_point { /cairo_point_y exch def /cairo_point_x exch def } def
/Tj { show currentpoint cairo_store_point } bind def
/TJ {
{
dup
type /stringtype eq
{ show } { -0.001 mul 0 cairo_font_matrix dtransform rmoveto } ifelse
} forall
currentpoint cairo_store_point
} bind def
/cairo_selectfont { cairo_font_matrix aload pop pop pop 0 0 6 array astore
cairo_font exch selectfont cairo_point_x cairo_point_y moveto } bind def
/Tf { pop /cairo_font exch def /cairo_font_matrix where
{ pop cairo_selectfont } if } bind def
/Td { matrix translate cairo_font_matrix matrix concatmatrix dup
/cairo_font_matrix exch def dup 4 get exch 5 get cairo_store_point
/cairo_font where { pop cairo_selectfont } if } bind def
/Tm { 2 copy 8 2 roll 6 array astore /cairo_font_matrix exch def
cairo_store_point /cairo_font where { pop cairo_selectfont } if } bind def
/g { setgray } bind def
/rg { setrgbcolor } bind def
/d1 { setcachedevice } bind def
%%EndProlog
%!FontType1-1.1 f-0-0 1.0
11 dict begin
/FontName /f-0-0 def
/PaintType 0 def
/FontType 1 def
/FontMatrix [0.001 0 0 0.001 0 0] readonly def
/FontBBox {0 -250 703 750 } readonly def
/Encoding 256 array
0 1 255 {1 index exch /.notdef put} for
dup 1 /g1 put
dup 2 /uni0061 put
dup 3 /uni0029 put
dup 4 /uni0046 put
dup 5 /uni006F put
dup 6 /uni0072 put
dup 7 /uni0077 put
dup 8 /uni0064 put
dup 9 /uni0070 put
dup 10 /uni0073 put
dup 11 /uni0062 put
dup 12 /uni0042 put
dup 13 /uni0063 put
dup 14 /uni006B put
readonly def
currentdict end
currentfile eexec
f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467
b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de
800a9c2cbaae3e98d910ec768e3ec0a68d6662ed0a4786decd0bc494681b7a8f496910cd7b4630
6dc8ae9298e870deef7a0ea25578226e409aa0b4147f7b82b7598ad503870b1045c168eb0c100e
42b0f810d306f8f7f38011c973c0fb43b3d4511ebd07847dd8587bb33de670b6b0b94ed7772fe2
92f5b1bfe83ad2c9d4f17ed4f34a3ab5b1c268ae43a5947d16525f9ffef6d057d70f07acc67574
27d7d3f6aa307d2bc8f59c82b8579918945de0cbfe985b93304f9c14bf82fa81a93da770200793
2dd1508907c081ece0f1ee6c8b723038fa488f2c81ddbb1ed051535205378c6d03a283ffea585c
6512c1ac0452fb99c4190b42973a1583935f2b9ecfda4f5c6879c310a39088a1c9f6c7d48e44ff
eff56f83cb68d530cc40ea36342fdb74d52256cb883f6759820be5b0f9a34c5da927f964272e52
dc34b64ab876eaebf2e744f2f0d8f5845003e27a3a27e098f035c208639194ae202e355a00fea8
136991d4991f2536e7f6bde1bc7143538d29199f4f20230fd5ecb614ab715150d159a0918fcfb6
358f638698f9442655d673f28397c64f0cad2df281ad1308353e3583992169fe2170384c6d33ce
abb707113d40c13481d0b05be5614b687e83d58bfdd52339f971d1edcb77cf95f6520452ec32bf
d5261a0c0e07d7992c376a677fcf40e4fd044ba412e003d0332d9b1f1b0e4456b089ce9837d0f1
15bd55065e60319b0a843800ca3b9d6fdc794179b340d1e77ce705a439ca07daccab58f9186e94
cdaa29beb312cd97993f8a353c705a405f445106127a7f902199a9dd5dcbb2329be7f2f40125e0
3d22c73b271e28479d0194cd949e8848c243ba17c2d37fbb86326e12840dfcaaaf5bdcda5bb181
1800aa2938e98ee8b67d291bd5e5b95d8e7dd36b4a71da9486eacc2e89923792520818cebdb349
ca34479be318df799d2f42fd59836e3422caa3c25a647622e0da55d98ba7591408713aba7b40c5
b8e3cbb63a5552ca72bf360dcd2ae223ef68c37ce2f8bc494e05805e5bacc706190fcf693c80eb
63f77aa65289faa69895f49e491410c6d5f1cc2300a9012dbfd915937f2bc8bdc4ad7a5bdc0a9a
a6dcec317f94a4d22bc1aa1464fa30132bb897136957f69f80e4b9fe051b9552da17dd4d776328
bbe9db2c435d733d0023160559e903dd283d7d834c3a5117b67919e9888d80b4596c1a95ff11b5
8fd6f817680add6d340040e497d0b5a326d4eb4221b8359a72da0bd0849d4aea1c5d66d9ef5c80
a068411472672882862c9cc6e67ec4f3d80d0b4e10645da2edfa586183e1f5884b3413339f7d36
4dfe9771fa7641efe09542369e3d2a67631c16c9db4ad96f58a937e84d390fc3929edb40453fa4
9ce08f6f4f9cc828b510547000d0b1950a9df6893395c4224eb4bb09c9cc952f1951cc7697df04
5f33f0445361bda924fab0c96ef285a87766237f3153020ab09bcfe9d6ab2b17ce3eb94a28598d
c905592fb8ff800170a36cb22e54dba7bf25168084e9914ed68f9d8564876c5eb8c4d1cf696141
720948c7fff51dd431baa98edc00aada78d4a94cef0a01434f3f57724b5cfd426c79feb930a95d
7c886a6e427726d0cfedd5060d6b432c92c1bd752f6ac6f2df5bd07b53cbed80f1413e41e06284
3288de3c08a038d0741b410f60aa4071767129604ff540d876e258afec2d64f00f179e147b0fb9
43101e2391840e8b507a8c4bacd82f2ac158a72eb06aeea7523eaf286baaf0e9558f13886c69b0
de61bb64126986e1a3c3b20337a806144051f9c0aa1a42a05b2c599ca827d687f4133e1e7d7776
8d1014b8c07ab69b56bfb9e7edb5d7ab66fd739077f910fec38e75af74302053bbf174be79f719
9bf0185ce294137a2804878c20c647246c37b367a0a96171e0a417866cab70bf94da9e0eac092e
169d44bd37d22b47a69de9ed3c0948051a8aa4de3ebf65b2cfa78361502261e98be16281c2541e
58d2a2624603b25bb01a237207092c0b2fe1bcd98f921b4ab15e72b2c08385d9b7375de6c82f7c
0ce46d643c85f229100162735f51ecf749bb51039fadf1b8d166a0cff86aae824ba1fec3845eac
1b6f50c3ec9bd0719a18651372a05e618e5a7d11bd8f359b0705ecf396d8856b3cced22f428ba6
78b3dea13299710ff7028785510f5519acf2cfa2b9eacd17963afd391c011d3710848500b284b3
ca2f0dcc9cb7e54434d088900be49c396f0b1db2a625fcd16fe4f9c28147ededf18508ded01053
e8830f6a6c9102bb5aa91ec9a47cf59d00e4abc987b5560f43f5ebc098757b3eb010d0775dd761
802bbabe16934ced432a9113dfbc89cde064557eecaaad8b65fd7390e455699cf0a66211d5b232
45f13e968775359e777981ccf4c08611f266e9b45e2e897d229e82d61b75481d626c63696108c7
fef9d5a03c9bcbf4f19e68f8185ab67090ea782303f03f5cdc3db65636d347c4ebe13a5e553e54
b0d11e53f9765a0caa21417471d0bb861beb2e57fb62c051d2e241ff524b6de48cffa72dca1868
1173d8b1626e68e9438be808ea73060264a0588320380c53fd4e22aced9b58b73a02c7c7844e39
69bb647a58a484aba8f05658aba1f1ac4f6cd1defd2efb330982db6d2ad76baa3cc500589c2346
6a7e5930e2cf7a312348797ad03fb006d5a956d8aea7648b791c0d1b6e7378822a7d31d9c8985c
72b79d0f4e8da8c4b1085785d1d74962f96406156ce2f8b2ce1e1e9eb0f2842d69a8e9ff0427c3
502adc1c5884f85dbfa44d2f9c0d4e5dd6f3d9997983c8cee03d8c02ce9f0aad219b159113b6ab
01368ac92822dc00af917a4a0514e70228ecc2293f57ba9a3daca3db6997f4d55d7e0dcd1cf52a
c9323b2529b10eaf558bd1271dcf42692778d144bc2df03a52e7bc25d4717dfa93f75bb4a04c3d
fd8a4c35253b6acb0217f7a88036fbb1f2e74b4a6063aee99612adf3066d2ba977eda6c39ff806
a07d639f985d82bd5850ef3b12d8224b95b440179ddab8844290677b425a17e071e0d0cba5de07
1bb6cc188343615cb79370604e57fb04eae8e8d5f325410f809c4bd63df311771fbad98f606e89
a2028e3ab1be4cc0f97c8817cd5b9a7de239b3d9603d34ff03a4a3c86302ab43ad7f583e83e51a
45e6ed849060cd6a442496141d247dde24be1c8c8440c15ebe12de69ef8eb626181537b1286274
b856549bb065d6f092a9823da544cf378726703359e0feda8063cddf141d68bb486c0b5085527c
2deb51c15f2da1d760182ecd326461ea7f206c39debdad7d0cebc97f8f40da2236776b1087258d
d07e3176b0bea9b6b316c510bae30fbedd118b36cbd5d813c5cdb6b67fe6ef28b5488f8dc51309
91216ff807199e36c98b1d1183a36a645b69a89450bf33792cf825b41e8f63bd4a978445b90ba5
5c71b5c631904d2f10ec6090e33fb3fdfee2d940c5940e81f50c9e4923a3af72cdebeab8511f98
7e96279e964df28f4ac9069a423f7f44dbc3b0ac3232298429622bffce0bfccc81a9c1b6934ba5
d094f10255651cfac77028a67e1a6df71c74066468f2b6e318a761ad0ec8e727daf45361c0a22e
1bf1cde946209cc355535b7ad556edb98998800db13ccc7058846cab6062a9ddb5f74724c729cb
5e4ac8f6134bccbfe336125b8e4a78e5bfb7766e0bc13a4bc604688c82323342b118c13bdf08e0
d8ec40f1a4e902157233e95823872f08a114d0bc634e7b9653d2d918c912c549ae0b57f7b40f7c
1fc41b4496151a257b761323ee9b206709adbc58911e5f53dc14eac7d0c4423d10f8b9c1bb275d
239885a0ca26cbee44a461749fc2e92fa953e9c6a4bce9a5768fd636598ecd5985d2adc4df4ba6
52d393a7e007bb1a64aa8ca95ca355bc170000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
cleartomark
%!FontType1-1.1 f-1-0 1.0
11 dict begin
/FontName /f-1-0 def
/PaintType 0 def
/FontType 1 def
/FontMatrix [0.001 0 0 0.001 0 0] readonly def
/FontBBox {0 -250 766 750 } readonly def
/Encoding 256 array
0 1 255 {1 index exch /.notdef put} for
dup 1 /uni0078 put
dup 2 /uni0045 put
dup 3 /uni0079 put
dup 4 /g4 put
dup 5 /uni0074 put
dup 6 /uni0040 put
dup 7 /uni003D put
dup 8 /uni0077 put
readonly def
currentdict end
currentfile eexec
f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467
b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de
800a9c2cbaae3e98d910ec768e3ec0a68d6662ed0a4786decd0bc494681b7a8f496910cd7b4630
6dc8ae9298e870deef7a0ea25578226e409aa0b4147f7b82b7598ad503870b1045c168eb0c100e
42b0f810d306f8f7f38011c973c0fb43b3d4511ebd07847dd8587bb33de670b6b0b94ed7772fe2
92f5b1bfe83ad2c9d4f17ed4f34a3ab5b1c268ae43a5947d16525f9ffef6d057d70f07acc67574
27d7d3f6aa307d2bc8f59c82b8579918945de0cbfe985b93304f9c14b9ebce3f978066a15e8e1c
b8872669361a9912c11ae5a12685e20e3a574fb8f5c503c55bdc4f2f97b63599b38a4b82078158
6824607f86ed6a4e2550c5a80c5c7238ab889d73953a4fbeae2b569260d980e20d6ac48f9813ba
d76f7624f1a9776f46fa475931f4acdb576a47b2db103722545b54f191651f855c67233c270449
688ae651b88d78c8972bf0ff8b6496cad1e36d3e70284b45de7161ae3788d098898aa40327d971
35861e5217018a802384bc6f712d6325f3a1e8c6c1ab7b79abce11cf82d88cd1b5eda2d146f6a5
cc6c2a29dd858b906d9ee39d6a18cf70b4b2b793b9371d594dc81a27ab5023bbf158e07b1ea085
ce28f7b81902277c241ecb16a8930f5d06cedfd905334cde1524e15f619ce117bd3df9f9766404
2370ffc21dd784b576afcb1e362564481c30c72deb8fcda4f0323901eaf210470211fb872e1115
5b233a2a330e13a961528d201baccde13309d62ea565cc874fbd88934c9cbec1499f8b19d65d0e
843ec372fb5369864328f57be3e321d286d1264f0343d1916ba2384d6ee269ae7b6eae2b372103
ca9742c0091d9f0ad3332a19839d2be908816cea939182b79e8268b18b6798af21976d47ded678
7bf13e3b539bc51b7d92792299719de3b8fe3c0816df7db52c04c022561c1d97e83048f450e9ab
db2835f9906edea4f4c47edc57a4fb2f257c17a07754b37a02fb73ba91dad8d4ba8f48c0bb1324
ca13b442cd9614b922415e37267f83466f4f832b56cea8c0b144881b08aca238f23e6d2f892885
3d1786304224452000b8c5bd8b5e9d49d98411a0ca6e1bdb6d362eeba41cb0fa6aca537c2d244d
a048a719446d9562677f47386a45e5ebc9421ad85de3e0473e7a29ea1a423dd4ab10fabb2f95ea
6e371f0e7ed0a66335284f4a25a6842f835662d928c37a5913607ca98ae1c48be483cc31969d36
34a2545101541ab389d52ef8d3324016ffb72f99f4705022eeed8c9560086937d84d57af894f03
7dff92cf9f83c3b67374a05171a141423e8bf5d552a808d0c671fe16bdb64f17d26dce9d9ddb91
cb04863eaaf1390aa5602e1a91532caff41b3d45639964bb765f3d6ac92bb8b9105bb18cb73301
65b4ffa03bcc1217576991d3c7bebfbef10627df4e685309444cd3fa92ca9fe80754863b06cf2a
36df89af19fc3f581f03bbb8c535f978416fd2fc56fcdba67791c5343420d9e85530127f0a5c8b
05c89ed37d97dd69adbf30262ac8b9db752181d9c51a4dfc04b63bbf98c7c4b413fdfbc101550f
af889cc30c445ebbe86737e7257bf4657cbcd0d6828536fc51836ff3eab24c7878867bfe9ad577
c28d4d064d5e1a51934e592a21739a53c137af4c79afd1940a94d004b8f9af2f56921500c0c5e4
37ada19550294f554c520e7d232311854c579c3db90ddd720cb7bfda6495d80040be25ce92d292
b0b72f8267da6e7861b152454c46b9bdf334a5d20d15e270bed5f5920d6085aa85c10b3704255e
fb9d4637c68fe595b0db73e03cd8dff167d9205a6e26275da84b7c5089fe3984be304ee8645df9
519be21afe772252524b32d2cf4585258ffc91f7eb6579478789219263f40c72d798fc1e1a9b8f
c8cf3f05fabaec84a2ae2c8cb596344996abc2707457618cc6d9b5b98090f2d52aa72f5c7493fe
5d0bc62d690cf0cabd64d2c4bb223c2575c93288adf6106f6ca36f574d037a7fa6f0d2dd784191
302f2b16b138f7d67048b3a4eceb007af35c009873a9574e88256e33b85cef4b2cc4d89e0fd74a
6f7a407267bf8a04b3e6867b38d0db2c96ff9ffe6262a2ada612e101a6c314d0955d5f13ba8712
48226fb301483ac93da9cc80e4da8d7d6fe69a41eeb19e91915599799027b7fff39121e591b700
db385f4afdf7e1e6325529552decc00ce96b0c8a7fbe97c55f3c3050b289a76e6233ff4f04f7f3
3894220591150f0bcefc76c0cc46043fc8ed2ae52e1b9ab09f928246f11d63645b01e59705d97c
aa32bf9d5c2979779aeea89b2d45411ceb451a244ce83154a0903729972246982e8830f5f870f6
648eca98f462471f4e964a2cf3bb286bb21dad73d1cd73ecd5435b21da80f112be0c3e70835357
16193f81afa2a6c3f76b5e4932da2dbc851ffcbd44617b7e1d1b38e6989c3dab8f8d98f3b038bc
69a660d2b30d37e636c4b9d5a9e5488fa3c2d52c20d6518e0a78683e0d1001da23e0341001a1e6
a6a704f7f561009a55fda5191b240819503312a100e2597c6d90a0c08ae4168031e4b9e705fd50
0477a439b58c3e814f4faf005612c7ca9b7f9ea629367177da069eab27a75d9390b8296ce3f1d5
6103daf4207001f915c70b25c404df6770f49cd1528fbe4eace30f55e4d9a8a8310dc0afc942a0
78d128d117adeb69d8530c9d8b6109a69469617c2ba0c3289e63f300a41d0c74ba5991491e6286
268e49e5815a7199cd7762303bbd2366c1821e64ce1a1ac6655f50d61cdcc87240ef15d9de6bd7
c8530000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
cleartomark
%!FontType1-1.1 f-2-0 1.0
11 dict begin
/FontName /f-2-0 def
/PaintType 0 def
/FontType 1 def
/FontMatrix [0.001 0 0 0.001 0 0] readonly def
/FontBBox {0 -20 529 674 } readonly def
/Encoding 256 array
0 1 255 {1 index exch /.notdef put} for
dup 1 /g1 put
dup 2 /uni0032 put
dup 3 /uni0033 put
dup 4 /uni0034 put
dup 5 /uni0035 put
dup 6 /uni0036 put
readonly def
currentdict end
currentfile eexec
f983ef0097ece636fb4a96c74d26ab84185f6dfa4a16a7a1c27bbe3f1156aea698df336d20b467
b10e7f33846656653c5ac6962759d3056cbdb3190bac614b984bf5a132dc418192443014ba63de
800a9c2cbaae3e98d910ec768e3ec0a68d6662ed0a4786decd0bc494681b7a8f496910cd7b4630
6dc8ae9298e870deef7a0ea25578226e409aa0b4147f7b82b7598ad503870b1045c168eb0c100e
42b0f810d306f8f7f38011c973c0fb43b3d4511ebd07847dd8587bb33de670b6b0b94ed7772fe2
92f5b1bfe83ad2c9d4f17ed4f34a3ab5b1c268ae43a5947d16525f9ffef6d057d70f07acc67574
27d7d3f6aa307d2bc8f59c82b8579918945de0cbfe985b93304f9c1debee09fd9cb05390ad5555
1edc479557e28c77bf40ae0ca1e76c9e9d5e9fd505cb2d360d436ba51b1eb6a57ed14c172717fa
9bb2c7a10c31769dc26929a055f956f7cf149ab4cba49646d92263d9bfee7cea7d985c42241e6e
c88d5f4d218a47241552eb4d05d25b1cb159184c752d34efb284de3305f87a7be8099e4f7c56f7
ae7e6b54f225ee0622bf95bb82bcc6d918117844dcf54aded8483f8e00fd0468defec73a6e95fc
3db701aafe06cbc8ea4d7bd4bfb424d91dabdca38f7ad5a466cec9679a751e60f09882d0e8bbcb
b2d556f429c1120168384e48feea2883b6c79384234f3048bac469d1f48d65101201e176a818a1
c1730e3843bae40fe84eb6c321482354833dcb3fe914ee9284b94942ee8bce81b5ca69224011d1
a4afd56e3806faa7fda0078a5c112add219abba1a16346ce5cc904dfaca53701e132fef0eb59a8
de7335876f57d22ae2b3ea046513df283556be1a4bc8d08ff28a80d58158e93125060fe79f22fe
dd8d269a5720c0ec14b139b0ff895b79bed2d383e5b7c5ec615b697020e344c60310c80a035cfe
2ffaf97d111e289f4bbd47f2f483dd1fec2a50ba7d809c810e91c3428d0e8194e7f531704eb45a
8fb198386426fa42a9baf9220254609045953f2fdc846b521d8e4c4c66f9dadd67267499b2833c
41d5037f6455b6d99c5be4bfaa07f03be8b476a0ad2bdf1b5a7224f452398a61aa7faaae6e9b71
12257c51bf7274c7a3859f6db6228c0105aee63991e71a7573b8bfe4db74ce1f993fb138481852
176e65a7b900fdd664707d931ad887be0401b89a34f3f049dd1c01877bc862a30b9b5fc78f8514
88c9eadf68ff7195898a9fea872ef6da1c6869db919d3319c875eab952f04dd48115643d9ee173
f76d6571ec9a297ce679dd1465de3d796dfe35a0ef8edb13fae907d4db3af55bd738245db5047f
6cbfdab3ed718cc8af3c907396c0c6b7b3396e70363660322617460417b9ee73ec79b1fe34ab13
48e13b679d43f9ba5e299d877983c8c75e4445628d9af01adbf0cd6b434c83657b97681dd6971a
aaab2d91744913c484a43b4b01aea8b901892752d5df833c367a9e7a475cd2bbd9902d82be506b
2adc3df42ad4538e986630a1964d1aac90cb949f0b50bbf516a8a4500eb0b4b516b826190e2fca
a96713f071d113b4769a25100e9ccbe77fb69d700e957ef5bf8037ce5a5bb6fdf6c5e699921ab3
cf705b02f235944d634af833f1a7f7c2b02657fa5ed68ced32b34bd8d1a65ac2ae71aef98489ee
79a0f205d2f8a683f1bae9cc4d6fe79ac36062dff0e27bb082914e1a762bc6377b5797e7ea2bf2
8a0262632e9a49238a8ce58c222bee6836729041561ff67c310631634eba45e2c5208bef811ba1
cd32e0725a672e9ddc3f9330ce19fef56b7ff88b2262e646c760a017aa87039c5bd16b00c36081
626c5adb38b08d83f6496d968384eaea9e153ce52120b3b2cb7a3d37b6003807e71638be92c3c8
8809a0fd54f53b0202d17232493d783b62b8824ab9df7a4aa536d37669c969a09ede073c66fcd7
498a3ec623ed603c0094862f5a589e7890a875a273394decdd423f0214a8ec6500e04b6a89b8dd
e284f787adec90657bcfedfbd551a348bc200af24ecb2a7f66cb0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
0000000000000000000000000000000000000000000000000000000000000000
cleartomark
%%Page: 1 1
%%BeginPageSetup
%%PageBoundingBox: 1 1 377 199
%%EndPageSetup
q
0 g
0.4 w
0 J
1 j
[] 0.0 d
10 M 112.508 136.895 m 115.473 142.82 113.07 150.031 107.141 152.992 c
101.215 155.957 94.008 153.555 91.043 147.629 c 88.078 141.699 90.48
134.492 96.41 131.527 c 102.336 128.562 109.543 130.965 112.508 136.895
c h
112.508 136.895 m S
BT
9.9626 0 0 9.9626 1 190.528 Tm
/f-0-0 1 Tf
[<010203>-333<04>83<050607>28<020608>-333<09020a0a>]TJ
ET
17 142.527 m 89 142.527 l S
89 142.527 m 84 144.191 l 84 140.863 l h
89 142.527 m f*
89 142.527 m 84 144.191 l 84 140.863 l h
89 142.527 m S
BT
9.9626 0 0 9.9626 1 142.531 Tm
/f-1-0 1 Tf
<01>Tj
6.9738 0 0 6.9738 6.694 141.037 Tm
/f-2-0 1 Tf
<01>Tj
9.9626 0 0 9.9626 1 62.531 Tm
/f-1-0 1 Tf
<01>Tj
6.9738 0 0 6.9738 6.694 61.037 Tm
/f-2-0 1 Tf
<02>Tj
9.9626 0 0 9.9626 289 102.528 Tm
/f-1-0 1 Tf
<02>Tj
/f-0-0 1 Tf
[<>-58<01>]TJ
/f-1-0 1 Tf
<03>Tj
6.9738 0 0 6.9738 305.687 101.033 Tm
/f-2-0 1 Tf
<03>Tj
9.9626 0 0 9.9626 310.157 102.528 Tm
/f-1-0 1 Tf
[<04>-167<05>]TJ
/f-0-0 1 Tf
<03>Tj
ET
17 62.527 m 89 62.527 l S
89 62.527 m 84 64.191 l 84 60.863 l h
89 62.527 m f*
89 62.527 m 84 64.191 l 84 60.863 l h
89 62.527 m S
209 102.527 m 281 102.527 l S
281 102.527 m 276 104.191 l 276 100.863 l h
281 102.527 m f*
281 102.527 m 276 104.191 l 276 100.863 l h
281 102.527 m S
112.121 57.195 m 115.086 63.125 112.684 70.332 106.754 73.297 c 100.828
76.262 93.617 73.855 90.656 67.93 c 87.691 62 90.094 54.793 96.02
51.828 c 101.949 48.867 109.156 51.27 112.121 57.195 c h
112.121 57.195 m S
207.516 96.852 m 210.48 102.777 208.078 109.988 202.148 112.949 c
196.223 115.914 189.016 113.512 186.051 107.586 c 183.086 101.656
185.488 94.449 191.418 91.484 c 197.344 88.52 204.551 90.922 207.516
96.852 c h
207.516 96.852 m S
BT
9.9626 0 0 9.9626 97 62.534 Tm
/f-1-0 1 Tf
<03>Tj
6.9738 0 0 6.9738 101.884 61.04 Tm
/f-2-0 1 Tf
<02>Tj
ET
0.827 g
101.812 34.32 39.926 15.941 re f
0 g
BT
9.9626 0 0 9.9626 104.800078 39.798592 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<03>]TJ
6.9738 0 0 6.9738 134.280078 38.304592 Tm
/f-2-0 1 Tf
<02>Tj
ET
0.827 g
289 1.047 87.062 15.941 re f
0 g
BT
9.9626 0 0 9.9626 291.989 6.526 Tm
/f-0-0 1 Tf
[<010b03>-333<0c020d>27<0e07>28<020608>-333<09020a0a>]TJ
/f-1-0 1 Tf
-25.193122 4.015518 Td
<08>Tj
6.9738 0 0 6.9738 48.132 45.037 Tm
/f-2-0 1 Tf
<04>Tj
ET
0.827 g
40.613 26.668 42.176 15.941 re f
0 g
BT
9.9626 0 0 9.9626 43.601 32.14872 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<08>]TJ
6.9738 0 0 6.9738 75.328 30.65472 Tm
/f-2-0 1 Tf
<04>Tj
9.9626 0 0 9.9626 41 170.531 Tm
/f-1-0 1 Tf
<08>Tj
6.9738 0 0 6.9738 48.132 169.037 Tm
/f-2-0 1 Tf
<01>Tj
ET
0.827 g
40.613 150.668 42.176 15.941 re f
0 g
BT
9.9626 0 0 9.9626 43.601 156.14872 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<08>]TJ
6.9738 0 0 6.9738 75.328 154.65472 Tm
/f-2-0 1 Tf
<01>Tj
9.9626 0 0 9.9626 25 130.531 Tm
/f-1-0 1 Tf
<08>Tj
6.9738 0 0 6.9738 32.132 129.037 Tm
/f-2-0 1 Tf
<02>Tj
ET
114.086 62.16 m 186.211 95.898 l S
186.211 95.898 m 180.977 95.289 l 182.387 92.27 l h
186.211 95.898 m f*
186.211 95.898 m 180.977 95.289 l 182.387 92.27 l h
186.211 95.898 m S
BT
9.9626 0 0 9.9626 25 90.531 Tm
/f-1-0 1 Tf
<08>Tj
6.9738 0 0 6.9738 32.132 89.037 Tm
/f-2-0 1 Tf
<03>Tj
9.9626 0 0 9.9626 97 142.534 Tm
/f-1-0 1 Tf
<03>Tj
6.9738 0 0 6.9738 101.884 141.04 Tm
/f-2-0 1 Tf
<01>Tj
9.9626 0 0 9.9626 193 102.534 Tm
/f-1-0 1 Tf
<03>Tj
6.9738 0 0 6.9738 197.884 101.04 Tm
/f-2-0 1 Tf
<03>Tj
ET
0.827 g
197.406 73.918 39.926 15.941 re f
0 g
BT
9.9626 0 0 9.9626 200.396511 79.395525 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<03>]TJ
6.9738 0 0 6.9738 229.876511 77.901525 Tm
/f-2-0 1 Tf
<03>Tj
9.9626 0 0 9.9626 153 142.531 Tm
/f-1-0 1 Tf
<08>Tj
6.9738 0 0 6.9738 160.132 141.037 Tm
/f-2-0 1 Tf
<05>Tj
9.9626 0 0 9.9626 153 70.531 Tm
/f-1-0 1 Tf
<08>Tj
6.9738 0 0 6.9738 160.132 69.037 Tm
/f-2-0 1 Tf
<06>Tj
ET
0.827 g
152.613 50.668 42.176 15.941 re f
0 g
BT
9.9626 0 0 9.9626 155.6006 56.14872 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<08>]TJ
6.9738 0 0 6.9738 187.3276 54.65472 Tm
/f-2-0 1 Tf
<06>Tj
ET
97 190.527 m 281 190.527 l S
281 190.527 m 274 192.859 l 274 188.195 l h
281 190.527 m f*
281 190.527 m 274 192.859 l 274 188.195 l h
281 190.527 m S
[ 4 4] 0 d
281 6.527 m 97 6.527 l S
97 6.527 m 104 4.195 l 104 8.859 l h
97 6.527 m f*
[] 0.0 d
97 6.527 m 104 4.195 l 104 8.859 l h
97 6.527 m S
0.827 g
288.812 82.32 38.5 15.941 re f
0 g
BT
9.9626 0 0 9.9626 291.799578 87.798592 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<02>]TJ
ET
[ 4 4] 0 d
209.438 98.582 m 281.438 98.582 l S
209.438 98.582 m 214.438 96.918 l 214.438 100.246 l h
209.438 98.582 m f*
[] 0.0 d
209.438 98.582 m 214.438 96.918 l 214.438 100.246 l h
209.438 98.582 m S
[ 4 4] 0 d
114.086 58.16 m 186.211 91.898 l S
114.086 58.16 m 119.32 58.77 l 117.906 61.785 l h
114.086 58.16 m f*
[] 0.0 d
114.086 58.16 m 119.32 58.77 l 117.906 61.785 l h
114.086 58.16 m S
[ 4 4] 0 d
17 58.527 m 89 58.527 l S
17 58.527 m 22 56.863 l 22 60.191 l h
17 58.527 m f*
[] 0.0 d
17 58.527 m 22 56.863 l 22 60.191 l h
17 58.527 m S
[ 4 4] 0 d
113.695 138.047 m 185.438 104.309 l S
113.695 138.047 m 117.512 134.41 l 118.93 137.426 l h
113.695 138.047 m f*
[] 0.0 d
113.695 138.047 m 117.512 134.41 l 118.93 137.426 l h
113.695 138.047 m S
113.695 142.047 m 185.438 108.309 l S
185.438 108.309 m 181.621 111.941 l 180.203 108.93 l h
185.438 108.309 m f*
185.438 108.309 m 181.621 111.941 l 180.203 108.93 l h
185.438 108.309 m S
15.117 137.664 m 91.98 70.691 l S
91.98 70.691 m 89.305 75.23 l 87.117 72.723 l h
91.98 70.691 m f*
91.98 70.691 m 89.305 75.23 l 87.117 72.723 l h
91.98 70.691 m S
15.375 67.957 m 92.754 133.902 l S
92.754 133.902 m 87.871 131.926 l 90.027 129.391 l h
92.754 133.902 m f*
92.754 133.902 m 87.871 131.926 l 90.027 129.391 l h
92.754 133.902 m S
[ 4 4] 0 d
17 138.527 m 89 138.527 l S
17 138.527 m 22 136.863 l 22 140.191 l h
17 138.527 m f*
[] 0.0 d
17 138.527 m 22 136.863 l 22 140.191 l h
17 138.527 m S
[ 4 4] 0 d
15.117 133.664 m 91.98 66.691 l S
15.117 133.664 m 17.793 129.121 l 19.98 131.633 l h
15.117 133.664 m f*
[] 0.0 d
15.117 133.664 m 17.793 129.121 l 19.98 131.633 l h
15.117 133.664 m S
[ 4 4] 0 d
15.375 63.957 m 92.754 129.902 l S
15.375 63.957 m 20.262 65.934 l 18.102 68.469 l h
15.375 63.957 m f*
[] 0.0 d
15.375 63.957 m 20.262 65.934 l 18.102 68.469 l h
15.375 63.957 m S
0.827 g
24.613 70.668 42.176 15.941 re f
0 g
BT
9.9626 0 0 9.9626 27.600643 76.148759 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<08>]TJ
6.9738 0 0 6.9738 59.327643 74.654759 Tm
/f-2-0 1 Tf
<03>Tj
ET
0.827 g
101.461 113.621 39.926 15.941 re f
0 g
BT
9.9626 0 0 9.9626 104.449471 119.098494 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<03>]TJ
6.9738 0 0 6.9738 133.929471 117.604494 Tm
/f-2-0 1 Tf
<01>Tj
ET
0.827 g
152.609 122.672 42.176 15.941 re f
0 g
BT
9.9626 0 0 9.9626 155.600152 128.149072 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<08>]TJ
6.9738 0 0 6.9738 187.327152 126.655072 Tm
/f-2-0 1 Tf
<05>Tj
ET
0.827 g
24.613 110.672 42.176 15.941 re f
0 g
BT
9.9626 0 0 9.9626 27.600803 116.14912 Tm
/f-1-0 1 Tf
[<06>-56<02>-57<0706>-56<08>]TJ
6.9738 0 0 6.9738 59.327803 114.65512 Tm
/f-2-0 1 Tf
<02>Tj
ET
Q
showpage
%%Trailer
count op_count sub {pop} repeat
countdictstack dict_count sub {end} repeat
cairo_eps_state restore
%%EOF
doc/Articles/Autodiff/figures/comp-graph/comp-graph.tex
+43
View File
@@ -0,0 +1,43 @@
\documentclass[tikz]{standalone}
\usepackage{framed}
\usepackage{amsmath}
\usepackage{amsfonts}
\DeclareMathOperator{\f}{f}
\usepackage{xcolor}
\usetikzlibrary{arrows}
\usetikzlibrary{positioning}
\begin{document}
\begin{tikzpicture}[]
\tikzstyle{vnode} = [circle,draw,thick,fill=white,minimum size=9mm]
\tikzstyle{vedge} = [->,>=latex,thick]
\node[vnode] (v-1) at (-8.5,0.5) {$v_{-1}$};
\node[vnode] (v0) at (-8.5,-2.5) {$v_0$};
\node[vnode] (v1) at (-6,0.5) {$v_1$};
\node[vnode] (v2) at (-6,-1) {$v_2$};
\node[vnode] (v3) at (-3.5,-2.5) {$v_3$};
\node[vnode] (v4) at (-3.5,0.5) {$v_4$};
\node[vnode] (v5) at (-1,-1) {$v_5$};
\node[] (x1) at (-11,0.5) {$x_1$};
\node[] (x2) at (-11,-2.5) {$x_2$};
\node[] (f) at (1.5,-1) {$f(x_1,x_2)$};
\draw (v-1) edge [vedge] (v1);
\draw (v-1) edge [vedge] (v2);
\draw (v0) edge [vedge] (v2);
\draw (v0) edge [vedge] (v3);
\draw (v1) edge [vedge] (v4);
\draw (v2) edge [vedge] (v4);
\draw (v3) edge [vedge] (v5);
\draw (v4) edge [vedge] (v5);
\draw (x1) edge [vedge] (v-1);
\draw (x2) edge [vedge] (v0);
\draw (v5) edge [vedge] (f);
\end{tikzpicture}
\end{document}
doc/Articles/Autodiff/figures/differentiation/differentiation.tex
+79
View File
@@ -0,0 +1,79 @@
\documentclass[tikz]{standalone}
\usepackage{framed}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{xcolor}
\usetikzlibrary{decorations.pathmorphing}
\usetikzlibrary{arrows}
\begin{document}
\begin{tikzpicture}[
pencildraw/.style={
decorate,
decoration={random steps,segment length=2pt,amplitude=1pt}
}]
%\draw[help lines] (0,0) grid (20,10);
\tikzstyle{paperbox} = [pencildraw,draw,thick,fill=white,text width=7cm,inner sep=5mm]
\tikzstyle{codebox} = [draw,thick,text width=7cm,inner sep=5mm]
\tikzstyle{edge} = [->,>=triangle 60,thick]
\node[paperbox] (a) at (-5.5,-12.5) {
$l_1=x$\\
$l_{n+1}=4l_n(1-l_n)$\\
\vspace{4mm}
$f(x)=l_4=64x(1 - x)(1 - 2 x)^2 (1 - 8 x + 8 x^2)^2$
};
\node[paperbox,fill=gray!20] (b) at (5.5,-12.5) {$f'(x)=128x(1 - x)(-8 + 16 x)(1 - 2 x)^2(1 - 8 x + 8 x^2) + 64 (1 - x)(1 - 2 x)^2 (1 - 8 x + 8 x^2)^2 - 64x(1 - 2 x)^2 (1 - 8 x + 8 x^2)^2 - 256x(1 - x)(1 - 2 x)(1 - 8 x + 8 x^2)^2$};
\node[codebox,fill=white] (c) at (-5.5,-18.5) {\texttt{\parbox{7cm}{
f(x):\\
\hphantom{tt} v = x\\
\hphantom{tt} for i = 1 to 3\\
\hphantom{tttt} v = 4*v*(1 - v)\\
\hphantom{tt} return v\\
\hphantom{t}\\
\textrm{or, in closed-form,}\\
\hphantom{t}\\
f(x):\\
\hphantom{tt}\parbox{6cm}{\hangindent=0.4cm \hangafter=1 return 64*x*(1-x)*((1-2*x)\^{}2)\\ *(1-8*x+8*x*x)\^{}2}\\
}}};
\node[codebox,black,fill=gray!20] (d) at (5.5,-18.5) {\color{black}\texttt{\parbox{7cm}{\textbf{
f'(x):\\
\hphantom{tt}\parbox{6cm}{\hangindent=0.4cm \hangafter=1 return 128*x*(1 - x)*(-8 + 16*x)\\ *((1 - 2*x)\^{}2)*(1 - 8*x + 8*x*x)\\+ 64*(1 - x)*((1 - 2*x)\^{}2)*((1 - 8*x + 8*x*x)\^{}2) - (64*x*(1 - 2*x)\^{}2)*(1 - 8*x + 8*x*x)\^{}2 - 256*x*(1 - x)*(1 - 2*x)*(1 - 8*x + 8*x*x)\^{}2}\\}
\flushright \color{gray} f'($\mathtt{x_0}$) $= f'(x_0)$\\\textrm{Exact}
}}};
\node[codebox,black,fill=gray!20] (e) at (-5.5,-25.5) {\color{black}\texttt{\parbox{7cm}{\textbf{
f'(x):\\
\hphantom{tt} (v,dv) = (x,1)\\
\hphantom{tt} for i = 1 to 3\\
\hphantom{tttt} (v,dv) = (4*v*(1-v), 4*dv-8*v*dv)\\
\hphantom{tt} return (v,dv)\\}
\flushright \color{gray} f'($\mathtt{x_0}$) $= f'(x_0)$\\\textrm{Exact}
}}};
\node[codebox,black,fill=gray!20] (f) at (5.5,-25.5) {\color{black}\texttt{\parbox{7cm}{\textbf{
f'(x):\\
\hphantom{tt} h = 0.000001\\
\hphantom{tt} return (f(x + h) - f(x)) / h\\}
\flushright \color{gray} f'($\mathtt{x_0}$) $\approx f'(x_0)$\\\textrm{Approximate}
}}};
\draw (a) edge [edge] (b);
\node[align=center,below] at (0,-12.5) {Manual\\Differentiation};
\draw (c) edge [edge] (d);
\node[align=center,below] at (0,-18.5) {Symbolic\\Differentiation\\of the Closed-form};
\draw (a) edge [edge] (c);
\node[align=left,right] at (-5.5,-15) {Coding};
\draw (b) edge [edge] (d);
\node[align=left,right] at (5.5,-15) {Coding};
\draw (c) edge [edge] (f);
\node[align=left,left] at (0.55,-22.25) {Numerical\\Differentiation};
\draw (c) edge [edge] (e);
\node[align=left,right] at (-5.5,-22.25) {Automatic\\Differentiation};
\end{tikzpicture}
\end{document}
doc/Articles/Autodiff/figures/helmholtz/helmholtz.tex
+598
View File
@@ -0,0 +1,598 @@
% Created by tikzDevice version 0.8.1 on 2015-02-19 05:31:36
% !TEX encoding = UTF-8 Unicode
\begin{tikzpicture}[x=1pt,y=1pt]
\definecolor{fillColor}{RGB}{255,255,255}
\path[use as bounding box,fill=fillColor,fill opacity=0.00] (0,0) rectangle (325.21,433.62);
\begin{scope}
\path[clip] ( 54.00,198.30) rectangle (313.21,421.62);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60,206.61) --
( 97.89,206.78) --
(132.18,207.17) --
(166.46,207.78) --
(200.75,208.75) --
(235.04,210.05) --
(269.33,211.84) --
(303.61,213.79);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 61.61,204.62) rectangle ( 65.59,208.61);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 95.89,204.79) rectangle ( 99.88,208.78);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (130.18,205.18) rectangle (134.17,209.17);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (164.47,205.78) rectangle (168.46,209.77);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (198.76,206.75) rectangle (202.75,210.74);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (233.05,208.05) rectangle (237.03,212.04);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (267.33,209.84) rectangle (271.32,213.83);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (301.62,211.79) rectangle (305.61,215.78);
\end{scope}
\begin{scope}
\path[clip] ( 0.00, 0.00) rectangle (325.21,433.62);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,206.57) -- ( 54.00,413.35);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,206.57) -- ( 48.00,206.57);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,247.93) -- ( 48.00,247.93);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,289.28) -- ( 48.00,289.28);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,330.64) -- ( 48.00,330.64);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,371.99) -- ( 48.00,371.99);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,413.35) -- ( 48.00,413.35);
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,206.57) {0};
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,247.93) {1000};
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,289.28) {2000};
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,330.64) {3000};
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,371.99) {4000};
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,413.35) {5000};
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,198.30) --
(313.21,198.30) --
(313.21,421.62) --
( 54.00,421.62) --
( 54.00,198.30);
\end{scope}
\begin{scope}
\path[clip] ( 54.00,198.30) rectangle (313.21,421.62);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 4pt off 4pt ,line join=round,line cap=round] ( 63.60,206.62) --
( 97.89,208.04) --
(132.18,213.88) --
(166.46,227.23) --
(200.75,249.80) --
(235.04,288.73) --
(269.33,341.78) --
(303.61,413.18);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60,203.80) --
( 66.42,206.62) --
( 63.60,209.44) --
( 60.78,206.62) --
( 63.60,203.80);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 97.89,205.22) --
(100.71,208.04) --
( 97.89,210.86) --
( 95.07,208.04) --
( 97.89,205.22);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (132.18,211.06) --
(135.00,213.88) --
(132.18,216.70) --
(129.36,213.88) --
(132.18,211.06);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (166.46,224.41) --
(169.28,227.23) --
(166.46,230.05) --
(163.64,227.23) --
(166.46,224.41);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (200.75,246.98) --
(203.57,249.80) --
(200.75,252.62) --
(197.93,249.80) --
(200.75,246.98);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (235.04,285.91) --
(237.86,288.73) --
(235.04,291.55) --
(232.22,288.73) --
(235.04,285.91);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (269.33,338.96) --
(272.15,341.78) --
(269.33,344.60) --
(266.51,341.78) --
(269.33,338.96);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61,410.36) --
(306.43,413.18) --
(303.61,416.00) --
(300.79,413.18) --
(303.61,410.36);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 1pt off 3pt ,line join=round,line cap=round] ( 63.60,206.63) --
( 97.89,207.14) --
(132.18,208.70) --
(166.46,212.04) --
(200.75,216.96) --
(235.04,226.00) --
(269.33,240.30) --
(303.61,262.07);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60,210.13) --
( 66.63,204.88) --
( 60.57,204.88) --
( 63.60,210.13);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 97.89,210.64) --
(100.92,205.39) --
( 94.86,205.39) --
( 97.89,210.64);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (132.18,212.20) --
(135.21,206.95) --
(129.15,206.95) --
(132.18,212.20);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (166.46,215.54) --
(169.49,210.29) --
(163.43,210.29) --
(166.46,215.54);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (200.75,220.46) --
(203.78,215.22) --
(197.72,215.22) --
(200.75,220.46);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (235.04,229.50) --
(238.07,224.25) --
(232.01,224.25) --
(235.04,229.50);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (269.33,243.80) --
(272.36,238.55) --
(266.30,238.55) --
(269.33,243.80);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61,265.57) --
(306.64,260.32) --
(300.58,260.32) --
(303.61,265.57);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 1pt off 3pt on 4pt off 3pt ,line join=round,line cap=round] ( 63.60,206.63) --
( 97.89,207.03) --
(132.18,207.87) --
(166.46,209.35) --
(200.75,211.29) --
(235.04,213.79) --
(269.33,217.08) --
(303.61,220.73);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60,203.14) --
( 66.63,208.38) --
( 60.57,208.38) --
( 63.60,203.14);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 97.89,203.53) --
(100.92,208.78) --
( 94.86,208.78) --
( 97.89,203.53);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (132.18,204.37) --
(135.21,209.62) --
(129.15,209.62) --
(132.18,204.37);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (166.46,205.85) --
(169.49,211.10) --
(163.43,211.10) --
(166.46,205.85);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (200.75,207.79) --
(203.78,213.03) --
(197.72,213.03) --
(200.75,207.79);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (235.04,210.29) --
(238.07,215.54) --
(232.01,215.54) --
(235.04,210.29);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (269.33,213.58) --
(272.36,218.83) --
(266.30,218.83) --
(269.33,213.58);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61,217.23) --
(306.64,222.48) --
(300.58,222.48) --
(303.61,217.23);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60,411.28) rectangle (207.00,351.28);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 66.30,399.28) -- ( 84.30,399.28);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 4pt off 4pt ,line join=round,line cap=round] ( 66.30,387.28) -- ( 84.30,387.28);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 1pt off 3pt ,line join=round,line cap=round] ( 66.30,375.28) -- ( 84.30,375.28);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 1pt off 3pt on 4pt off 3pt ,line join=round,line cap=round] ( 66.30,363.28) -- ( 84.30,363.28);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 73.31,397.29) rectangle ( 77.29,401.28);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 75.30,384.46) --
( 78.12,387.28) --
( 75.30,390.10) --
( 72.48,387.28) --
( 75.30,384.46);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 75.30,378.78) --
( 78.33,373.53) --
( 72.27,373.53) --
( 75.30,378.78);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 75.30,359.78) --
( 78.33,365.03) --
( 72.27,365.03) --
( 75.30,359.78);
\node[text=drawColor,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 93.30,395.84) {$f$, original function};
\node[text=drawColor,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 93.30,383.84) {$\nabla f$, numerical diff.};
\node[text=drawColor,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 93.30,371.84) {$\nabla f$, forward AD};
\node[text=drawColor,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 93.30,359.84) {$\nabla f$, reverse AD};
\end{scope}
\begin{scope}
\path[clip] ( 54.00, 54.00) rectangle (313.21,185.10);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60, 58.86) --
( 97.89, 82.16) --
(132.18, 96.98) --
(166.46,106.90) --
(200.75,115.33) --
(235.04,122.01) --
(269.33,127.93) --
(303.61,132.42);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 61.61, 56.86) rectangle ( 65.59, 60.85);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 95.89, 80.16) rectangle ( 99.88, 84.15);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (130.18, 94.99) rectangle (134.17, 98.98);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (164.47,104.91) rectangle (168.46,108.90);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (198.76,113.34) rectangle (202.75,117.32);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (233.05,120.01) rectangle (237.03,124.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (267.33,125.94) rectangle (271.32,129.93);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (301.62,130.43) rectangle (305.61,134.41);
\end{scope}
\begin{scope}
\path[clip] ( 0.00, 0.00) rectangle (325.21,433.62);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 58.70, 54.00) -- (303.61, 54.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 58.70, 54.00) -- ( 58.70, 48.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (107.68, 54.00) -- (107.68, 48.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (156.67, 54.00) -- (156.67, 48.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (205.65, 54.00) -- (205.65, 48.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (254.63, 54.00) -- (254.63, 48.00);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61, 54.00) -- (303.61, 48.00);
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 58.70, 32.40) {0};
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at (107.68, 32.40) {10};
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at (156.67, 32.40) {20};
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at (205.65, 32.40) {30};
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at (254.63, 32.40) {40};
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at (303.61, 32.40) {50};
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 54.00) --
(313.21, 54.00) --
(313.21,185.10) --
( 54.00,185.10) --
( 54.00, 54.00);
\end{scope}
\begin{scope}
\path[clip] ( 0.00, 0.00) rectangle (325.21,197.10);
\definecolor{drawColor}{RGB}{0,0,0}
\node[text=drawColor,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at (183.61, 8.40) {$n$};
\end{scope}
\begin{scope}
\path[clip] ( 0.00, 0.00) rectangle (325.21,433.62);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 58.86) -- ( 54.00,157.31);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 58.86) -- ( 48.00, 58.86);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 91.67) -- ( 48.00, 91.67);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,124.49) -- ( 48.00,124.49);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,157.31) -- ( 48.00,157.31);
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60, 52.11) {10};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 0.70] at ( 35.51, 62.10) {0};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60, 84.92) {10};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 0.70] at ( 35.51, 94.92) {1};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,117.74) {10};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 0.70] at ( 35.51,127.74) {2};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 39.60,150.56) {10};
\node[text=drawColor,rotate= 90.00,anchor=base west,inner sep=0pt, outer sep=0pt, scale= 0.70] at ( 35.51,160.56) {3};
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 55.68) -- ( 54.00,185.04);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 55.68) -- ( 50.40, 55.68);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 57.35) -- ( 50.40, 57.35);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 68.73) -- ( 50.40, 68.73);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 74.51) -- ( 50.40, 74.51);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 78.61) -- ( 50.40, 78.61);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 81.79) -- ( 50.40, 81.79);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 84.39) -- ( 50.40, 84.39);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 86.59) -- ( 50.40, 86.59);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 88.49) -- ( 50.40, 88.49);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00, 90.17) -- ( 50.40, 90.17);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,101.55) -- ( 50.40,101.55);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,107.33) -- ( 50.40,107.33);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,111.43) -- ( 50.40,111.43);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,114.61) -- ( 50.40,114.61);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,117.21) -- ( 50.40,117.21);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,119.41) -- ( 50.40,119.41);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,121.31) -- ( 50.40,121.31);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,122.99) -- ( 50.40,122.99);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,134.37) -- ( 50.40,134.37);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,140.15) -- ( 50.40,140.15);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,144.25) -- ( 50.40,144.25);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,147.43) -- ( 50.40,147.43);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,150.03) -- ( 50.40,150.03);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,152.22) -- ( 50.40,152.22);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,154.13) -- ( 50.40,154.13);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,155.80) -- ( 50.40,155.80);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,167.19) -- ( 50.40,167.19);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,172.96) -- ( 50.40,172.96);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,177.06) -- ( 50.40,177.06);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,180.24) -- ( 50.40,180.24);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,182.84) -- ( 50.40,182.84);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 54.00,185.04) -- ( 50.40,185.04);
\end{scope}
\begin{scope}
\path[clip] ( 54.00, 54.00) rectangle (313.21,185.10);
\definecolor{drawColor}{RGB}{0,0,0}
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 4pt off 4pt ,line join=round,line cap=round] ( 63.60, 60.04) --
( 97.89,109.75) --
(132.18,132.61) --
(166.46,147.41) --
(200.75,157.94) --
(235.04,167.09) --
(269.33,174.19) --
(303.61,180.23);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60, 57.22) --
( 66.42, 60.04) --
( 63.60, 62.86) --
( 60.78, 60.04) --
( 63.60, 57.22);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 97.89,106.93) --
(100.71,109.75) --
( 97.89,112.57) --
( 95.07,109.75) --
( 97.89,106.93);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (132.18,129.79) --
(135.00,132.61) --
(132.18,135.43) --
(129.36,132.61) --
(132.18,129.79);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (166.46,144.59) --
(169.28,147.41) --
(166.46,150.23) --
(163.64,147.41) --
(166.46,144.59);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (200.75,155.12) --
(203.57,157.94) --
(200.75,160.76) --
(197.93,157.94) --
(200.75,155.12);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (235.04,164.27) --
(237.86,167.09) --
(235.04,169.91) --
(232.22,167.09) --
(235.04,164.27);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (269.33,171.37) --
(272.15,174.19) --
(269.33,177.01) --
(266.51,174.19) --
(269.33,171.37);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61,177.41) --
(306.43,180.23) --
(303.61,183.05) --
(300.79,180.23) --
(303.61,177.41);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 1pt off 3pt ,line join=round,line cap=round] ( 63.60, 63.11) --
( 97.89, 96.15) --
(132.18,115.04) --
(166.46,128.48) --
(200.75,137.62) --
(235.04,146.54) --
(269.33,154.40) --
(303.61,161.50);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60, 66.61) --
( 66.63, 61.36) --
( 60.57, 61.36) --
( 63.60, 66.61);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 97.89, 99.65) --
(100.92, 94.40) --
( 94.86, 94.40) --
( 97.89, 99.65);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (132.18,118.54) --
(135.21,113.29) --
(129.15,113.29) --
(132.18,118.54);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (166.46,131.98) --
(169.49,126.73) --
(163.43,126.73) --
(166.46,131.98);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (200.75,141.12) --
(203.78,135.87) --
(197.72,135.87) --
(200.75,141.12);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (235.04,150.04) --
(238.07,144.79) --
(232.01,144.79) --
(235.04,150.04);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (269.33,157.90) --
(272.36,152.65) --
(266.30,152.65) --
(269.33,157.90);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61,165.00) --
(306.64,159.75) --
(300.58,159.75) --
(303.61,165.00);
\path[draw=drawColor,line width= 0.4pt,dash pattern=on 1pt off 3pt on 4pt off 3pt ,line join=round,line cap=round] ( 63.60, 64.84) --
( 97.89, 93.19) --
(132.18,107.97) --
(166.46,118.84) --
(200.75,126.36) --
(235.04,132.43) --
(269.33,137.78) --
(303.61,142.03);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 63.60, 61.34) --
( 66.63, 66.59) --
( 60.57, 66.59) --
( 63.60, 61.34);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] ( 97.89, 89.69) --
(100.92, 94.94) --
( 94.86, 94.94) --
( 97.89, 89.69);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (132.18,104.47) --
(135.21,109.72) --
(129.15,109.72) --
(132.18,104.47);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (166.46,115.34) --
(169.49,120.59) --
(163.43,120.59) --
(166.46,115.34);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (200.75,122.86) --
(203.78,128.11) --
(197.72,128.11) --
(200.75,122.86);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (235.04,128.94) --
(238.07,134.18) --
(232.01,134.18) --
(235.04,128.94);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (269.33,134.28) --
(272.36,139.53) --
(266.30,139.53) --
(269.33,134.28);
\path[draw=drawColor,line width= 0.4pt,line join=round,line cap=round] (303.61,138.53) --
(306.64,143.78) --
(300.58,143.78) --
(303.61,138.53);
\end{scope}
\begin{scope}
\path[clip] ( 0.00, 0.00) rectangle (325.21,433.62);
\definecolor{drawColor}{RGB}{0,0,0}
\node[text=drawColor,rotate= 90.00,anchor=base,inner sep=0pt, outer sep=0pt, scale= 1.00] at ( 15.60,219.76) {Time};
\end{scope}
\end{tikzpicture}
doc/Articles/Autodiff/jmlr2e_mod.sty
+353
View File
@@ -0,0 +1,353 @@
%
% File: Macros for Journal of Machine Learning Research
% Very minor modification of macros for Journal of Artificial
% Intelligence Research (jair.sty)
%
% Suggestions: Submit an issue or pull request to
% https://github.com/JournalMLR/jmlr-style-file
%
% Last edited October 9, 2000 by Leslie Pack Kaelbling
% Last edited January 23, 2001 by Alex J. Smola (we should set up RCS or CVS)
% Last edited March 29, 2004 Erik G. Learned-Miller
% Last edited January 17, 2016 Charles Sutton
% Last edited January 9, 2017 Charles Sutton
% (We have now set up GIT, good thing that we waited for it to
% be invented.)
%
% The name of this file should follow the article document
% type, e.g. \documentstyle[jmlr]{article}
% Copied and edited from similar file for Machine Learning Journal.
% Original Author: Jeff Schlimmer
% Edited by: Kevin Thompson, Martha Del Alto, Helen Stewart, Steve Minton \& Pandu Nayak.
% Last edited: Mon May 3 20:40:00 1993 by kthompso (Kevin Thompson) on muir
\typeout{Document Style `jmlr' -- January 2016.}
\newif\if@abbrvbib\@abbrvbibfalse
\DeclareOption{abbrvbib}{\@abbrvbibtrue}
\newif\if@usehyper\@usehypertrue
\DeclareOption{nohyperref}{\@usehyperfalse}
\DeclareOption{hyperref}{\@usehypertrue}
\DeclareOption*{\PackageWarning{jmlr}{Unknown ‘\CurrentOption’}}
\ProcessOptions\relax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% REQUIRED PACKAGES
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\RequirePackage{epsfig}
\RequirePackage{amssymb}
\RequirePackage{natbib}
\RequirePackage{graphicx}
\if@usehyper
\RequirePackage[colorlinks=false,allbordercolors={1 1 1}]{hyperref}
\fi
\if@abbrvbib
\bibliographystyle{abbrvnat}
\else
\bibliographystyle{plainnat}
\fi
\bibpunct{(}{)}{;}{a}{,}{,}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% P A G E S I Z E
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change the overall width of the page. If these parameters are
% changed, they will require corresponding changes in the
% maketitle section.
%
\renewcommand{\topfraction}{0.95} % let figure take up nearly whole page
\renewcommand{\textfraction}{0.05} % let figure take up nearly whole page
% Specify the dimensions of each page
\oddsidemargin .25in % Note \oddsidemargin = \evensidemargin
\evensidemargin .25in
\marginparwidth 0.07 true in
%\marginparwidth 0.75 true in
%\topmargin 0 true pt % Nominal distance from top of page to top of
%\topmargin 0.125in
\topmargin -0.5in
\addtolength{\headsep}{0.25in}
\textheight 8.5 true in % Height of text (including footnotes & figures)
\textwidth 6.0 true in % Width of text line.
\widowpenalty=10000
\clubpenalty=10000
\@twosidetrue \@mparswitchtrue \def\ds@draft{\overfullrule 5pt}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% S E C T I O N S
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Definitions for nicer (?) sections, etc., ideas from Pat Langley.
% Numbering for sections, etc. is taken care of automatically.
\def\@startsiction#1#2#3#4#5#6{\if@noskipsec \leavevmode \fi
\par \@tempskipa #4\relax
\@afterindenttrue
\ifdim \@tempskipa <\z@ \@tempskipa -\@tempskipa \@afterindentfalse\fi
\if@nobreak \everypar{}\else
\addpenalty{\@secpenalty}\addvspace{\@tempskipa}\fi \@ifstar
{\@ssect{#3}{#4}{#5}{#6}}{\@dblarg{\@sict{#1}{#2}{#3}{#4}{#5}{#6}}}}
\def\@sict#1#2#3#4#5#6[#7]#8{\ifnum #2>\c@secnumdepth
\def\@svsec{}\else
\refstepcounter{#1}\edef\@svsec{\csname the#1\endcsname}\fi
\@tempskipa #5\relax
\ifdim \@tempskipa>\z@
\begingroup #6\relax
\@hangfrom{\hskip #3\relax\@svsec.\hskip 0.1em}
{\interlinepenalty \@M #8\par}
\endgroup
\csname #1mark\endcsname{#7}\addcontentsline
{toc}{#1}{\ifnum #2>\c@secnumdepth \else
\protect\numberline{\csname the#1\endcsname}\fi
#7}\else
\def\@svsechd{#6\hskip #3\@svsec #8\csname #1mark\endcsname
{#7}\addcontentsline
{toc}{#1}{\ifnum #2>\c@secnumdepth \else
\protect\numberline{\csname the#1\endcsname}\fi
#7}}\fi
\@xsect{#5}}
\def\@sect#1#2#3#4#5#6[#7]#8{\ifnum #2>\c@secnumdepth
\def\@svsec{}\else
\refstepcounter{#1}\edef\@svsec{\csname the#1\endcsname\hskip 0.5em }\fi
\@tempskipa #5\relax
\ifdim \@tempskipa>\z@
\begingroup #6\relax
\@hangfrom{\hskip #3\relax\@svsec}{\interlinepenalty \@M #8\par}
\endgroup
\csname #1mark\endcsname{#7}\addcontentsline
{toc}{#1}{\ifnum #2>\c@secnumdepth \else
\protect\numberline{\csname the#1\endcsname}\fi
#7}\else
\def\@svsechd{#6\hskip #3\@svsec #8\csname #1mark\endcsname
{#7}\addcontentsline
{toc}{#1}{\ifnum #2>\c@secnumdepth \else
\protect\numberline{\csname the#1\endcsname}\fi
#7}}\fi
\@xsect{#5}}
\def\thesection {\arabic{section}}
\def\thesubsection {\thesection.\arabic{subsection}}
\def\section{\@startsiction{section}{1}{\z@}{-0.24in}{0.10in}
{\large\bf\raggedright}}
\def\subsection{\@startsection{subsection}{2}{\z@}{-0.20in}{0.08in}
{\normalsize\bf\raggedright}}
\def\subsubsection{\@startsection{subsubsection}{3}{\z@}{-0.18in}{0.08in}
{\normalsize\sc\raggedright}}
\def\paragraph{\@startsiction{paragraph}{4}{\z@}{1.5ex plus
0.5ex minus .2ex}{-1em}{\normalsize\bf}}
\def\subparagraph{\@startsiction{subparagraph}{5}{\z@}{1.5ex plus
0.5ex minus .2ex}{-1em}{\normalsize\bf}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% F O O T N O T E S
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change the size of the footnote rule
%
% \renewcommand{\footnoterule}{\vspace{10pt}\hrule width 0mm}
\long\def\@makefntext#1{\@setpar{\@@par\@tempdima \hsize
\advance\@tempdima-15pt\parshape \@ne 15pt \@tempdima}\par
\parindent 2em\noindent \hbox to \z@{\hss{\@thefnmark}. \hfil}#1}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% A B S T R A C T
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use \begin{abstract} .. \end{abstract} for abstracts.
\renewenvironment{abstract}
{\centerline{\large\bf Abstract}\vspace{0.7ex}%
\bgroup\leftskip 20pt\rightskip 20pt\small\noindent\ignorespaces}%
{\par\egroup\vskip 0.25ex}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% KEYWORDS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% use \begin{keywords} .. \end{keywords} for keywordss.
\newenvironment{keywords}
{\bgroup\leftskip 20pt\rightskip 20pt \small\noindent{\bf Keywords:} }%
{\par\egroup\vskip 0.25ex}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIRST PAGE, TITLE, AUTHOR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author information can be set in various styles:
% For several authors from the same institution:
% \author{Author 1 \and ... \and Author n \\
% \addr{Address line} \\ ... \\ \addr{Address line}}
% if the names do not fit well on one line use
% Author 1 \\ {\bf Author 2} \\ ... \\ {\bf Author n} \\
% To start a seperate ``row'' of authors use \AND, as in
% \author{Author 1 \\ \addr{Address line} \\ ... \\ \addr{Address line}
% \AND
% Author 2 \\ \addr{Address line} \\ ... \\ \addr{Address line} \And
% Author 3 \\ \addr{Address line} \\ ... \\ \addr{Address line}}
% Title stuff, borrowed in part from aaai92.sty
\newlength\aftertitskip \newlength\beforetitskip
\newlength\interauthorskip \newlength\aftermaketitskip
%% Changeable parameters.
\setlength\aftertitskip{0.1in plus 0.2in minus 0.2in}
\setlength\beforetitskip{0.05in plus 0.08in minus 0.08in}
\setlength\interauthorskip{0.08in plus 0.1in minus 0.1in}
\setlength\aftermaketitskip{0.3in plus 0.1in minus 0.1in}
%% overall definition of maketitle, @maketitle does the real work
\def\maketitle{\par
\begingroup
\def\thefootnote{\fnsymbol{footnote}}
\def\@makefnmark{\hbox to 0pt{$^{\@thefnmark}$\hss}}
\@maketitle \@thanks
\endgroup
\setcounter{footnote}{0}
\let\maketitle\relax \let\@maketitle\relax
\gdef\@thanks{}\gdef\@author{}\gdef\@title{}\let\thanks\relax}
\def\@startauthor{\noindent \normalsize\bf}
\def\@endauthor{}
\def\@starteditor{\noindent \small {\bf ~}}
\def\@endeditor{\normalsize}
\def\@maketitle{\vbox{\hsize\textwidth
\linewidth\hsize \vskip \beforetitskip
{\begin{center} \Large\bf \@title \par \end{center}} \vskip \aftertitskip
{\def\and{\unskip\enspace{\rm and}\enspace}%
\def\addr{\small\it}%
\def\email{\hfill\small\sc}%
\def\name{\normalsize\bf}%
\def\AND{\@endauthor\rm\hss \vskip \interauthorskip \@startauthor}
\@startauthor \@author \@endauthor}
\vskip \aftermaketitskip
\noindent \@starteditor \@editor \@endeditor
\vskip \aftermaketitskip
}}
\def\editor#1{\gdef\@editor{}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%
%%% Pagestyle
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Defines the pagestyle for the title page.
%% Usage: \jmlrheading{1}{1993}{1-15}{8/93}{9/93}{14-115}{Jane Q. Public and A. U. Thor}
%% \jmlrheading{vol}{year}{pages}{Submitted date}{published date}{paper id}{authors}
%%
%% If your paper required revisions that were reviewed by the action editor, then indicate
%% this by, e.g.
%% \jmlrheading{1}{1993}{1-15}{8/93; Revised 10/93}{12/93}{14-115}{Jane Q. Public and A. U. Thor}
\def\firstpageno#1{\setcounter{page}{#1}}
\def\jmlrheading#1#2#3#4#5#6#7{\def\ps@jmlrtps{\let\@mkboth\@gobbletwo%
\def\@oddhead{\scriptsize }%
\def\@oddfoot{\parbox[t]{\textwidth}{\raggedright \scriptsize \hfill}}%
\def\@evenhead{}\def\@evenfoot{}}%
\thispagestyle{jmlrtps}}
%% Defines the pagestyle for the rest of the pages
%% Usage: \ShortHeadings{Minimizing Conflicts}{Minton et al}
%% \ShortHeadings{short title}{short authors}
\def\ShortHeadings#1#2{\def\ps@jmlrps{\let\@mkboth\@gobbletwo%
\def\@oddhead{\hfill {\small\sc #1} \hfill}%
\def\@oddfoot{\hfill \small\rm \thepage \hfill}%
\def\@evenhead{\hfill {\small\sc #2} \hfill}%
\def\@evenfoot{\hfill \small\rm \thepage \hfill}}%
\pagestyle{jmlrps}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MISCELLANY
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define macros for figure captions and table titles
% Figurecaption prints the caption title flush left.
% \def\figurecaption#1#2{\noindent\hangindent 42pt
% \hbox to 36pt {\sl #1 \hfil}
% \ignorespaces #2}
% \def\figurecaption#1#2{\noindent\hangindent 46pt
% \hbox to 41pt {\small\sl #1 \hfil}
% \ignorespaces {\small #2}}
\def\figurecaption#1#2{\noindent\hangindent 40pt
\hbox to 36pt {\small\sl #1 \hfil}
\ignorespaces {\small #2}}
% Figurecenter prints the caption title centered.
\def\figurecenter#1#2{\centerline{{\sl #1} #2}}
\def\figurecenter#1#2{\centerline{{\small\sl #1} {\small #2}}}
%
% Allow ``hanging indents'' in long captions
%
\long\def\@makecaption#1#2{
\vskip 10pt
\setbox\@tempboxa\hbox{#1: #2}
\ifdim \wd\@tempboxa >\hsize % IF longer than one line:
\begin{list}{#1:}{
\settowidth{\labelwidth}{#1:}
\setlength{\leftmargin}{\labelwidth}
\addtolength{\leftmargin}{\labelsep}
}\item #2 \end{list}\par % Output in quote mode
\else % ELSE center.
\hbox to\hsize{\hfil\box\@tempboxa\hfil}
\fi}
% Define strut macros for skipping spaces above and below text in a
% tabular environment.
\def\abovestrut#1{\rule[0in]{0in}{#1}\ignorespaces}
\def\belowstrut#1{\rule[-#1]{0in}{#1}\ignorespaces}
% Acknowledgments
\long\def\acks#1{\vskip 0.3in\noindent{\large\bf Acknowledgments}\vskip 0.2in
\noindent #1}
% Research Note
\long\def\researchnote#1{\noindent {\LARGE\it Research Note} #1}
\renewcommand{\appendix}{\par
\setcounter{section}{0}
\setcounter{subsection}{0}
\def\thesection{\Alph{section}}
\def\section{\@ifnextchar*{\@startsiction{section}{1}{\z@}{-0.24in}{0.10in}%
{\large\bf\raggedright}}%
{\@startsiction{section}{1}{\z@}{-0.24in}{0.10in}
{\large\bf\raggedright Appendix\ }}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PROOF, THEOREM, and FRIENDS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newcommand{\BlackBox}{\rule{1.5ex}{1.5ex}} % end of proof
\newenvironment{proof}{\par\noindent{\bf Proof\ }}{\hfill\BlackBox\\[2mm]}
\newtheorem{example}{Example}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{conjecture}[theorem]{Conjecture}
\newtheorem{axiom}[theorem]{Axiom}
doc/Articles/CNN.pdf
BIN
View File
Binary file not shown.
doc/Articles/Fig1.pdf
BIN
View File
Binary file not shown.
doc/Articles/bibliography.bib
+3268
View File
File diff suppressed because it is too large Load Diff
doc/Articles/fig_qst_beryllium.pdf
BIN
View File
Binary file not shown.
doc/Articles/fig_qst_rydberg.pdf
BIN
View File
Binary file not shown.
doc/Articles/fig_supervised.pdf
BIN
View File
Binary file not shown.
doc/Articles/fig_vmc.pdf
BIN
View File
Binary file not shown.
doc/Articles/main.bbl
+2688
View File
File diff suppressed because it is too large Load Diff
doc/Articles/main.tex
+734
View File
@@ -0,0 +1,734 @@
\documentclass[twocolumn,english,reprint,superscriptaddress,longbibliography,pra]{revtex4-1}
\renewcommand{\familydefault}{\rmdefault}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\setcounter{secnumdepth}{3}
\usepackage{color}
\usepackage{bm}
\usepackage{amstext}
\usepackage{amssymb}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{tabularx}
\usepackage{amsmath,mathrsfs}
\usepackage{enumitem}
\makeatletter
\usepackage{multibib}
\usepackage{listings}
\usepackage{mathtools}
\usepackage{braket}
\DeclareMathOperator*{\argmin}{argmin}
\usepackage{babel}
\usepackage{bbold}
\usepackage{mathrsfs}
\usepackage{hyperref}
\hypersetup{
linktocpage = true,
colorlinks,
citecolor=blue,
filecolor=black,
linkcolor=blue,
urlcolor=blue
}
\newcommand{\BeH}{BeH$_2$}
\newcommand{\tr}{\textrm{Tr}}
\newcommand{\FIG}{\textcolor{blue}{FIGURE\:}}
\newcommand{\REF}{[\textcolor{red}{ref}]}
\usepackage{xcolor}
\definecolor{codegreen}{rgb}{0,0.6,0}
\definecolor{codegray}{rgb}{0.5,0.5,0.5}
\definecolor{codepurple}{rgb}{0.58,0,0.82}
\definecolor{backcolour}{rgb}{0.95,0.95,0.92}
\lstdefinestyle{mystyle}{
backgroundcolor=\color{backcolour},
commentstyle=\color{codegreen},
keywordstyle=\color{magenta},
numberstyle=\tiny\color{codegray},
stringstyle=\color{codepurple},
basicstyle=\ttfamily\footnotesize,
breakatwhitespace=false,
breaklines=true,
captionpos=b,
keepspaces=true,
numbers=left,
numbersep=5pt,
showspaces=false,
showstringspaces=false,
showtabs=false,
tabsize=2
}
\lstset{style=mystyle}
\usepackage{xpatch}
\makeatletter
\patchcmd{\@ssect@ltx}
{\addcontentsline{toc}{#1}{\protect\numberline{}#8}}
{}
{}
{}
\makeatother
\begin{document}
\title{Neural networks in quantum many-body physics: a hands-on tutorial}
\author{Juan Carrasquilla}
\affiliation{Vector Institute, MaRS Centre, Toronto, Ontario, M5G 1M1, Canada}
\author{Giacomo Torlai}
\email{gttorlai@amazon.com}
\thanks{\\This work has been done before Giacomo Torlai joined Amazon.}
\affiliation{AWS Center for Quantum Computing, Pasadena, CA 91125, USA}
\affiliation{Center for Computational Quantum Physics, Flatiron Institute, New York, NY 10010, USA}
\begin{abstract}
Over the past years, machine learning has emerged as a powerful computational tool to tackle complex problems over a broad range of scientific disciplines. In particular, artificial neural networks have been successfully deployed to mitigate the exponential complexity often encountered in quantum many-body physics, the study of properties of quantum systems built out of a large number of interacting particles. In this Article, we overview some applications of machine learning in condensed matter physics and quantum information, with particular emphasis on hands-on tutorials serving as a quick-start for a newcomer to the field. We present supervised machine learning with convolutional neural networks to learn a phase transition, unsupervised learning with restricted Boltzmann machines to perform quantum tomography, and variational Monte Carlo with recurrent neural-networks for approximating the ground state of a many-body Hamiltonian. We briefly review the key ingredients of each algorithm and their corresponding neural-network implementation, and show numerical experiments for a system of interacting Rydberg atoms in two dimensions.
\end{abstract}
\maketitle
\section{Introduction}
Quantum many-body physics refers to the mathematical framework to study the collective behavior of large numbers of interacting particles. The emerging cooperative phenomena that result from seemingly simple interactions can produce an astounding variety of phases of matter such as conventional metals and magnetically ordered states, as well as unanticipated states including high-temperature superconductivity, strange metals, and spin liquids~\cite{Xiao:803748}. In addition to naturally occurring quantum systems, many-body physics studies synthetic quantum matter, e.g., ultracold atoms, superconducting qubits, and trapped ions, which simultaneously reveal new phenomena in highly controlled laboratory settings and advances the development of quantum computers and other quantum information processing devices.
In spite of the simplicity of the physical laws that govern such multi-particle quantum objects, the theoretical and experimental analysis of these systems confront us with complexities which are ultimately rooted in the "curse of dimensionality" associated with the exponential explosion of the size of the space where quantum many-body states live in. Traditionally, the study of many-body systems is performed with the help of tools designed to circumvent this dimensionality explosion and produce a succinct, low-dimensional description that captures the essential aspects of a quantum system. Such descriptions arise from the analysis of data generated in a wide range of theoretical, computational, and experimental devices. These include numerical simulations of model Hamiltonians based on quantum Monte Carlo or variational algorithms, but also experimental arrays of complex electronic-structure images obtained from spectroscopic imaging scanning tunnelling microscopy, or measurements of quantum states prepared on a physical quantum computing platform.
Machine learning, already explored as a tool in several research areas in physics~\cite{RevModPhys.91.045002}, offers a set of alternative approaches to the study of quantum many-body systems in experiments and numerical simulations~\cite{doi:10.1080/23746149.2020.1797528,annurev-conmatphys-031119-050651}. The resurgence of activity at the intersection between physics and machine learning is in part due to a series of scientific breakthroughs in computer vision and natural language processing. Such progress has led to a burst of research where neural networks have been repurposed to tackle fundamental questions in condensed matter physics, quantum computing, statistical physics, and atomic, molecular and optical physics. Machine learning, and in particular deep neural networks, have been used to identify phases of matter in numerical simulations and experiments~\cite{carrasquilla2017nature, evert2017nature, torlai_learning_2016, leiwang2016, chng2017, broecker2017, eun-ah2017, dassarma2017, neupeurt2017, yi-ting2018, huembeli2018,PhysRevB.99.060404,PhysRevB.99.104410,Zhang_MLcuprates,Bohrdt2018,Rem2018,PhysRevLett.122.210503},
to increase the performance of Monte Carlo simulations~\cite{huang2017,junwei2017, xiao_yan2017, inack2018, parolini2019, pilati2019, mcnaughton2020, albergo2019}, to accurately describe the state of classical~\cite{Wu_2019} and quantum systems~\cite{androsiuk1993,LAGARIS19971, Carleo_2017, zi2018, Di_Luo, pfau2019abinitio, hermann2019deep, PhysRevLett.122.250502, PhysRevLett.122.250501,PhysRevLett.122.250503,PhysRevB.99.214306, choo_fermionicnqs2020, PhysRevLett.124.020503, RNNWF_2020, roth2020iterative}, to develop novel quantum control strategies~\cite{PhysRevX.8.031086,PhysRevX.8.031084,PhysRevLett.122.020601,niu_universal_2019,2020arXiv201003655Y,coopmans2020}, to perform quantum tomography~\cite{torlai_Tomo,rocchetto,Torlai_latent,2018arXiv181206693Q,carrasquilla_povm,biamonte_qst,torlai_rydberg19,xin_local-measurement-based_2019,Sehayek2019,torlai_chemistry,PhysRevA.102.022412,NoriGAN,Tiunov:20,Cha2020,PhysRevA.102.042604,DeVlugt2020,2020arXiv200907601S,torlai_QPT,morawetz2020,Nori2020}, to accelerate density functional theory calculations~\cite{PhysRevLett.108.253002,doi:10.1063/1.4834075,PhysRevB.94.245129,brockherde_bypassing_2017,PhysRevA.100.022512,PhysRevLett.125.076402,PhysRevResearch.2.033388}, to develop and elucidate renormalization group analyses~\cite{2014arXiv1410.3831M,koch-janusz_mutual_2018,PhysRevE.97.053304,PhysRevLett.121.260601,PhysRevResearch.2.023369,2020arXiv201005703C}, to devise quantum error correction protocols~\cite{torlai_neural_2016,krastanov_deep_2017,Varsamopoulos_2017,Baireuther2018machinelearning,Chamberland_2018,Breuckmann2018scalableneural,Nautrup2019optimizingquantum,PhysRevLett.122.200501,PhysRevA.99.052351,Andreasson2019quantumerror,Evert2020QEC,Ni2020neuralnetwork,PhysRevResearch.2.033399,PhysRevResearch.2.023230}, among many other examples~\cite{PhysRevB.97.045153,Seif_2018,Melnikov1221,Dunjko_2018,PhysRevE.99.062106, 2019arXiv191211052C,PhysRevB.99.075113,PhysRevLett.124.010508, PhysRevX.10.011006,2020arXiv200600712H,2020arXiv200905580L,2020arXiv201014510L}.
Such an explosion of activity indicates that machine learning techniques may soon become commonplace in quantum many-body physics research, both in experiments and numerical simulation. These clear trends call for the development of resources to stimulate researchers to familiarize with the wealth of concepts, intuition, algorithms, hardware, software, and research culture entailed by the adoption of machine learning and neural networks in physics research. Here, we take a step forward in this direction and develop a set of hands-on tutorials focused on a set of recent prototypical examples of applications of neural network technology to problems in statistical physics, condensed matter and quantum computing.
\subsection*{Outline}
The Article is organized as follows. Starting with a preliminary discussion, we introduce in Sec~\ref{preliminariesA} some fundamental concepts in machine learning and neural networks. In Sec~\ref{preliminariesB} we present a concise description of the physical system studied in our numerical experiments, a two-dimensional array of interacting Rydberg atoms. In Sec~\ref{supervised} we discuss our first application, the classification of phases of matter with supervised machine learning of projective measurement data using a convolutional neural network, and demonstrate it on the quantum phase transition in the Rydberg atoms. In Sec~\ref{qst} we introduce quantum state tomography, and show how this problem can be phrased as an unsupervised machine learning task. Using the restricted Boltzmann machine, we show quantum tomography of the Rydberg ground states, as well as of the ground state of a small molecule from qubit measurement data. In Sec~\ref{vmc} we present the simulation of the ground state of a many-body Hamiltonian using variational Monte Carlo with a recurrent neural network wavefunction. For each of these applications, we also show the key components of the underlying software, with full code tutorials available in an external repository~\cite{coderepo}.
\section{Preliminaries}
\subsection{Machine learning with neural networks}
\label{preliminariesA}
Artificial intelligence, the scientific discipline that deals with the theory and development of computer programs with the ability to perform complex tasks, saw early success solving problems which are relatively straightforward to formalize in an abstract way. The solutions to this breed of problems are typically described by a list of very precise formal rules that computers can process efficiently. As remarkable example, computers have been beating humans at playing chess since 1997, due in part to the fact that chess involves a large set of formal rules.
Modern machine learning, instead, deals with the challenge of automatizing the solution of real world tasks that may be easy for humans to process but that are hard to formally describe by simple rules. These techniques have spurred a recent revolution where algorithms trained using data have started to match humans' ability to recognize objects in an image, decipher speech or translate text to multiple languages, which are tasks that are difficult to formalize and articulate through simple rules.
A key element behind these recent developments can be largely traced back to a series of breakthroughs in the development of powerful neural network models, where data is processed through the sequential combination of multiple nonlinear layers~\cite{Goodfellow-et-al-2016}. Such models solve a fundamental problem in learning real world tasks, namely the problem of automatically extracting knowledge from raw noisy data, rather than relying on hard-coded knowledge directly inscribed in the algorithms by a human. Neural networks automatize the construction of sets of increasingly complex representations of the data, which can be understood as the computational disentangling of complex concepts (e.g. an object in a cluttered image) out of simpler concepts (e.g. pixel values and basic shapes like edges). These representations, in turn, lead to solutions to learning tasks with unprecedented success.
% ML tasks
\begin{figure*}[t]
\noindent \centering{}\includegraphics[width=2.05\columnwidth]{Fig1}
\caption{Rydberg atoms in a two-dimensional square array. ({\bf a}) Schematic representation of the phase diagram at a fixed value of the interaction $V=3$ MHz and $\Omega = 1$ MHz. At large and negative detuning, the system is in a disordered (paramagnetic) phase with all atoms in the ground state. At large and positive detuning, the atoms are found in a checkerboard pattern with N\'eel order. ({\bf b}) Snake-like geometry of the MPS path along the square lattice, used for the DMRG simulations. Ground state energy ({\bf c}) and the staggered magnetization (N\'eel order) ({\bf d}) as a function of the detuning $\delta$, for a $8\times8$ array ($V=3$ MHz). ({\bf e}) Absolute value of the average occupation number in momentum space $|n(\bm{k})|$ deep into the $Z_2$ ordered phase ($\delta=4$ MHz), showing a peak at $\bm{k}=(\pi,\pi)$, a signature of anti-ferromagnetic order. ({\bf f}) Energy gap $\Delta$ between the ground state and first excited state, detecting a quantum phase transition at detuning $\delta\approx1.3$.}
\label{Fig::1}
\end{figure*}
For practical purposes, machine learning algorithms can be divided into the categories of supervised, unsupervised, and reinforcement learning, all of which have found applications to quantum many-body systems~\cite{doi:10.1080/23746149.2020.1797528}. While there is no formal difference between some of the algorithms in these categories when expressed in the language of probability~\cite{Goodfellow-et-al-2016,10.5555/1162264}, such a division is often used as a way to specify the details of the algorithms, the training setup, and the structure of the data sets involved.
Supervised learning tasks aim at predicting a target output vector $\bm y$ associated with input vector $\bm x$, both of which can be discrete or continuous. The training data is thus a list of pairs of input/output tuples $\{ \bm x_i,\bm y_i\}_{i=1}^{M}$, where target output conveys that such a vector corresponds to the ideal output given the input vector~\cite{10.5555/1162264}. Starting with a training data set with $M$ entries, the learning algorithm outputs a function $\hat{\bm y} = f(\bm x)$ which estimates the output values for unseen input vectors $\bm x$. Examples of supervised learning include classification, where the objective is to assign each input vector to one of a set of discrete categories, and the task of
regression, where the output is a vector with continuous entries. Examples for classification and regression are respectively the problem of recognizing images of handwritten digits and the problem of determining the orbits of bodies around the sun from astronomical data.
Unsupervised learning deals with the learning tasks where the training data is composed of a set of input vectors without a corresponding target output~\cite{10.5555/1162264}. These algorithms are typically used to discover hidden structure in the data sets. Examples of tasks in unsupervised learning problems include clustering, where the objective is to discover of groups of similar examples within the data, density estimation, where the objective is to estimate the underlying probability distribution associated with the data, as well as low-dimensional visualization of high-dimensional data algorithms, which depict complex data in two or three dimensions while trying to retain key spatial characteristics in the original data.
Finally, reinforcement learning, although not discussed in this Article, develops algorithms dealing with the problem of discovering actions that maximize a numerical reward signal~\cite{10.5555/551283}. The learning algorithms are not necessarily directly exposed to examples of optimal actions. Instead, it must discover them by a process similar to a guided trial and error. Reinforcement learning augmented by deep neural networks has successfully learned policies from high-dimensional sensory input for game playing achieving human-level performance in several challenging games including Atari 2600~\cite{mnih_human-level_2015} as well as the board game Go~\cite{silver2016}. Likewise, reinforcement learning has been applied to the control of quantum systems~\cite{bukov2018,niu_universal_2019} as well as to the optimization of quantum error correction codes~\cite{Nautrup2019optimizingquantum,Andreasson2019quantumerror,Evert2020QEC}, one key ingredient in the development of fault-tolerant quantum computers.
\subsection{Rydberg atoms in two dimensions}
\label{preliminariesB}
We demonstrate the machine learning algorithms discussed in this Article for a many-body system composed by interacting Rydberg atoms. Engineered arrays of cold Rydberg atoms are increasingly used for highly-controlled quantum simulations of strongly-interacting matter~\cite{Schauss1455,Endres2016,Labuhn,Bernien2017,keesling_quantum_2019,2020arXiv201212281E,2020arXiv201212268S}, as well as for quantum information processing~\cite{Levine2019}. We specifically consider a square array with linear dimension $L$ containing $N = L^2$ atoms. Each atom is described by a local Hilbert space spanned by the states $\{|g\rangle,|e\rangle\}$, referring respectively to the atomic ground and the highly-excited Rydberg states. The atoms are subject to a uniform laser drive with Rabi frequency $\Omega$ and detuning $\delta$, and they interact with one another via the Van der Waals potential $V(x)\approx r^{-6}$ at short distances. The resulting many-body Hamiltonian is
\begin{equation}
\hat{H} = -\Omega \sum_{\bm{r}}\hat{S}^x(\bm{r}) -\delta\sum_{\bm{r}}^N\hat{\Pi}(\bm{r}) +\frac{1}{2}\sum_{\bm{r},\bm{r^\prime}}V(\bm{r}-\bm{r^\prime})\hat{\Pi}(\bm{r}) \hat{\Pi}(\bm{r^\prime})
\label{Eq::RydbergHamiltonian}
\end{equation}
where $\hat{\Pi}(\bm{r}) =|e\rangle\!\langle e|_{\bm{r}} $ is the projector onto the Rydberg state at position $\bm{r}$, $\hat{S}^x(\bm{r})=\frac{1}{2}\hat\sigma^x(\bm{r})$ are spin-$\frac{1}{2}$ operators, and $V(\bm{r}-\bm{r^\prime})=V_0/\|\bm{r}-\bm{r^\prime}\|^6$ is the Van der Waals potential between atoms at position $\bm{r}$ and $\bm{r^{\prime}}$. In the following, we assume $\Omega=1$ MHz.
The phase diagram for the ground state of the Rydberg Hamiltonian is dictated by the mechanism of Rydberg blockade, a constraint that prevents two atoms at sufficiently small distances to be simultaneously excited to the Rydberg states. We can characterize the phase diagram in terms of the detuning $\delta$ and the interaction strength $V_0$. On the square lattice, several different orders have been detected by numerical simulations~\cite{Samajdar2020}. Here, we specifically focus on the $Z_2$ transition between a disordered phase at large and negative detuning, where all atoms are found in the ground state, and an ordered phase at large and positive detuning, where the system is found in one of the two symmetry-broken N\'eel states characterized by a checkerboard pattern in the atomic occupation number (Fig.~\ref{Fig::1}(a)).
We perform numerical simulations of the ground state of Hamiltonian~(\ref{Eq::RydbergHamiltonian}) using the density matrix renormalization group (DMRG)~\cite{PhysRevLett.69.2863,PhysRevB.48.10345,SCHOLLWOCK201196} implemented using the software package ITensor~\cite{itensor}. We adopt a matrix product state (MPS) variational wavefunction $|\Psi\rangle$ with a snake-like geometry shown in Fig.~\ref{Fig::1}(b). We fix the interaction strength to $V_0=3$ MHz, and retain up to the third-nearest-neighbor interactions. For a several values of the detuning $\delta\in\{-5,5\}$ MHz, we run DMRG to find an approximation of the ground state, using a singular value decomposition cutoff of $10^{-10}$ and a target energy accuracy of $10^{-5}$. To certify convergence to the ground state, each run is repeated for different initialization of the starting MPS.
We show the results of the simulations for a $8\times8$ array with open boundary conditions in Fig.~\ref{Fig::1}(c-f). We plot, as a function of the detuning, the ground state energy per site $E_0/N=\langle\Psi_0|\hat{H}|\Psi_0\rangle/N$, and the staggered magnetization $\langle\mathcal{N}\rangle=N^{-1}\sum_{\bm{r}}(-1)^{x+y}\langle\hat{S}^z(\bm{r})\rangle$, which can be used to detect N\'eel order. Whenever all atoms are in the ground state, $\langle\mathcal{N}\rangle\approx0$, while for an ordered state with a checkerboard pattern one has $\langle\mathcal{N}\rangle\approx0.5$. We also show the average occupation number in momentum space,
\begin{equation}
n(\bm{k})=\frac{1}{\sqrt{N}}\sum_{\bm{r}}e^{i\bm{k}\cdot\bm{r}}\langle\hat{n}(\bm{r})\rangle
\end{equation}
where $\hat{n}(\bm{r}) = \frac{1}{2}(1-2\hat{S}^z(\bm{r}))$. We observe a peak at $\linebreak\bm{k}=(\pi,\pi)$ for large detuning $\delta=4$ MHz (Fig.~\ref{Fig::1}(e)), and a featureless state at negative detuning (not shown).
The two phases of the Rydberg atoms are separated by a second-order quantum phase transition at a critical point $\delta_c$. We can extract an approximation of $\delta_c$ by measuring the energy gap $\Delta = |E_0-E_1|$ between the ground state and the first excited state $|\Psi_1\rangle$. We compute $E_1$ by running DMRG on the Hamiltonian $\hat{H}^\prime=\hat{H}+\omega|\Psi\rangle\langle\Psi|$ where $\omega$ is an energy penalty. From the energy gap curve, we estimate the detuning where $\Delta\approx0$ to be $\delta_c\approx 1.3$ MHz. This approximate value will be sufficient for the purpose of this Article, though a more systematic scaling study with appropriate boundary conditions (to minimize finite-size effects) should be performed to accurately determine the critical point and critical exponents of the transition.
Once we have solved for the ground states of the Rydberg Hamiltonian, the corresponding MPSs can be used to generate data to train the neural networks for the different applications. In this case, the data consists of projective measurements in the atomic occupation number basis $|\bm{\sigma}\rangle=|\sigma_1,\dots,\sigma_N\rangle$, where $\sigma_j=0$ and $\sigma_j=1$ refers respectively to the $j$-th atom being in the ground and Rydberg state. Given a wavefunction $|\Psi\rangle$, the probability to observe an atomic pattern $\bm{\sigma}$ following a measurement is simply given by the Born rule $P(\bm{\sigma})=|\langle\bm{\sigma}|\Psi\rangle|^2$. Because of the intrinsic one-dimensional geometry of an MPS, it is possible to efficiently sample the probability distribution $P(\bm{\sigma})$ by exploiting the chain rule of probabilities. Moreover, the sampling is exact in the sense that each sample is completely independent of one another~\cite{PhysRevB.85.165146}. %The sampling procedure consists of iteratively building one-site reduced density matrices conditional on the previous measurement outcomes~\cite{PhysRevB.85.165146}.
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
% SUPERVISED LEARNING
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
\section{Learning a quantum phase transition}
\label{supervised}
An important task in condensed matter and statistical physics is to characterize different phases of matter and the associated phase transitions between them.
Typically, phases of matter are described in terms of simple real-space patterns and their associated order parameters, which are theoretically understood using Landau symmetry-breaking paradigm~\cite{Xiao:803748}. While a wide array of theoretical and experimental tools to study interacting quantum systems have been constructed in relation to these patterns, there is an increasing set of states of matter whose theoretical and experimental understanding eludes the Landau symmetry-breaking paradigm. The characterization of these phases may rely on, e.g., out-of-equilibrium properties of the system as in the many-body localized phase~\cite{basko2006,RevModPhys.91.021001}, or on topological invariants in topological phases and spin liquids~\cite{Xiao:803748,savaryQuantumSpinLiquids2016, doi:10.1146/annurev-conmatphys-031218-013401}.
Machine learning provides an alternative route to the characterization of phases of matter and their associated phase transitions in a semi-automated fashion without a direct use of manually specified real-space patterns and/or other signatures, provided that a sufficiently large training set is available. In its simplest form~\cite{carrasquilla2017nature}, given the existence of a classical or quantum phase transition between two phases in a physical system, one can use supervised learning to attempt to classify experimental or numerical snapshots of the phases of matter separated by the transition. This task can be achieved using most classification algorithms, e.g., those based on a neural network or a support vector machine~\cite{10.5555/1162264}, trained on snapshots of two phases of matter labelled according to the corresponding phase out of which the snapshot originated. Although here we only explore this simple strategy, we stress that machine learning approaches to studying phases and phase transitions have been significantly expanded and they no longer require the precise knowledge of the location of the critical point~\cite{evert2017nature,broecker2017b}, can be fully automatized, and can discover ordered phases~\cite{carrasquilla2017nature,leiwang2016,PhysRevE.96.022140}, topological phases~\cite{evert2017nature,rodriguez-nieva2019,PhysRevLett.125.127401}, and phases such as the many-body localized phase which is characterized by its dynamical properties~\cite{yi-ting2018,rao2018,PhysRevLett.120.257204}.
The nature of the snapshots used to train the learning algorithms is vastly flexible, hence these strategies are of wide applicability, and can include numerically generated configurations visited during a classical or quantum Monte Carlo simulation of the physical system~\cite{carrasquilla2017nature,leiwang2016,evert2017nature,broecker2017,chng2017,PhysRevB.99.060404,PhysRevB.99.121104}, entanglement spectra~\cite{evert2017nature,yi-ting2018}, correlation matrices~\cite{PhysRevB.102.054512,PhysRevLett.125.170603}, tensors in an MPS~\cite{PhysRevLett.125.170603}, numerically generated projective measurements~\cite{berezutskii2020}, high-resolution real-space snapshots of complex many-body systems obtained with quantum gas microscopes for ultracold atoms~\cite{Bohrdt2018,PhysRevA.102.033326}, single-shot experimental momentum-space density images of ultracold quantum gases~\cite{Rem2018}, spectroscopic imaging scanning tunnelling microscopy data~\cite{Zhang_MLcuprates}, among many others.
Below we explore learning a quantum phase in an array of interacting Rydberg atoms using projective measurements. As a classification algorithm, we make use a of a convolutional neural network~\cite{Goodfellow-et-al-2016} that readily takes advantage of the two-dimensional (2D) spatial arrangement of the the Rydberg atoms and their locality, as well as the approximate translation invariance of the system.
\begin{figure}[t]
\noindent \centering{}\includegraphics[width=\columnwidth]{CNN.pdf}
\caption{A schematic representation of a convolutional neural network. The elements of the input $\mathsf{h}^{(0)}_{l,j,k}$ corresponds to the outcome of a projective measurement on the Rydberg system. The first operation is a convolutional layer with $ M_y \times M_x = 3 \times 3$ kernels with $I_{\text{out}} = 32$ output channels and $L_{\text{input}}=1$. This kernel is convolved with an input configuration with $N = 8 \times 8$ Rydberg atoms. Likewise, the second operation corresponds to a convolutional layer with $ M_y \times M_x = 3 \times 3$ kernels with $I_{\text{out}} = 32$ output channels and $L_{\text{input}}=32$. The output of the second convolutional layer is flattened and fed to an FC layer with a ReLU activation, followed by another FC layer with a softmax activation which produces the prediction outcome. }
\label{Fig::CNN}
\end{figure}
\subsection{Convolutional neural networks and their training}
Convolutional neural networks (CNN) employ a mathematical operation called convolution to process information for data that has a natural grid-like topology~\cite{Goodfellow-et-al-2016}. A 2D convolutional layer implements the operation
\begin{align*}
\mathsf{h}^{(q)}_{i,j,k} & = F\left(\sum_{l,m_y,m_x} \mathsf{h}^{(q-1)}_{l,j+m_y,k+m_x} \mathsf{K}^{(q)}_{i,l,m_y,m_x} \right) \\
&\coloneqq F\left( \mathsf{K}^{(q)} *\mathsf{h}^{(q-1)} \right) \nonumber %_{i,j,k} \nonumber
\end{align*}
where the trainable kernel $\mathsf{K}^{(q)}_{i,l,m_y,m_x}$ at layer $q$ specifies the connection strength between a unit in channel $i$ of the output and a unit in channel $l$ of the input, with a spatial offsets of $m_y$ rows (labeled y direction) and $m_x$ columns (labeled x direction) between the output and the input variables. The dimensions of the array $\mathsf{K}^{(q)}_{i,l,m,n}$ are $I_{\text{out}}$, $L_{\text{input}}$, $M_{y}$, $M_{x}$, which corresponds to the number of output channels, input channels, dimension of the filter along the vertical and horizontal directions, respectively. The activation at layer $q$ consists of elements $\mathsf{h}^{(q)}_{l,j,k}$, where $j$ and $k$ label vertical and horizontal directions, respectively, and $l$ specifies the channel. The activation units are labelled by $q$ where $q=0$ corresponds to the raw projective measurement data. Finally, the non-linear function $F(x)$, which in our examples is typically a rectified linear unit (ReLU) $F(x)=\text{max}(0,x)$, is applied element-wise to each of the components of its input. A convolutional neural network equipped with two convolutional layers is schematically shown in Fig.~\ref{Fig::CNN}.
Followed by the convolutional layers, a CNN typically processes information using sets of fully connected (FC) layers which implement a matrix-vector operation followed by a non-linearity $F$ as
\begin{equation}
\mathsf{h}^{(q)}_{i} = F\left(\sum_{l} \mathsf{h}^{(q-1)}_{l} \mathsf{K}^{(q)}_{i,l} + \mathsf{b}^{(q)}_{i}\right),
\end{equation}
where the trainable parameters of the FC layer are the kernel $\mathsf{K}^{(q)}_{i,l}$ and the bias vector $ \mathsf{b}^{(q)}_{i}$.
To feed the output of a convolutional layer $\mathsf{h}^{(q-1)}_{l,j,k}$ to an FC layer, the array is reshaped or ``flattened'' to $\mathsf{h^{\prime}}^{(q-1)}_{l}$ so that all the original components packed into a one-dimensional array with dimension $L_{\text{FC}}$. The last two layers of the CNN in Fig.~\ref{Fig::CNN} correspond to two fully connected layers with a ReLU and a softmax non-linearities, respectively. The softmax function $S$ is given by
\begin{equation}
\text{S}(\bm{v}) = \frac{\exp(\bm{v})}{\sum_i \exp(v_i)}.
\end{equation}
where $v_i$ are the components of a vector $\bm{v}$ and the $\text{exp}$ function acts element-wise on the components of the vector. We note that the input to the CNN and its trainable parameters are real, so that the outcome of the softmax layer can be interpreted as a probability distribution since $0 \le \text{S}(v_i)\le 1$ and $\sum_{i}\text{S}(v_i)=1$.
Finally, we mention that we interpret our CNN as a model for the conditional probability of assigning a phase of matter $y=0,1$ to a projective measurement outcome $\bm{\sigma}=\mathsf{h}^{(0)}$, i.e. $P_{\bm{\theta}}(y|\bm{\sigma})$, where $\bm{\theta}$ encompasses all the trainable parameters of the CNN. The conditional is given by
\begin{widetext}
\begin{equation}
P_{\bm{\theta}}(y|\bm{\sigma}) = S\left( b^{(4)}_{y} + \sum_{m}\mathsf{K}^{(4)}_{y,m} F\left(b^{(3)}_{m} + \sum_{l} \mathsf{K}^{(3)}_{m,l}\,\text{Flatten}\left(F\left((\mathsf{K}^{(2)}*F\left(\mathsf{K}^{(1)}*\bm{\sigma} \right)\right)\right)_{l}\right) \right),
\label{Eq::CNNdist}
\end{equation}
\end{widetext}
where the function $\text{Flatten}()_{l}$ is the $l$-th component of a vector that arises from reshaping the incoming argument of the function to a one-dimensional array.
%\subsection{Phase transition in the Rydberg array}
To estimate the parameters of the CNN we use the maximum likelihood principle, where the parameters of a statistical model are selected by assigning high probability to the observed data. For a dataset with observations $\{\bm{\sigma}_n, y_n \}_{n=1}^{M}$, where $y_n=0,1$ label the phase of matter out of which a projective measurement $\bm{\sigma}_n$ was taken from, the likelihood assigned by the model to the dataset can be written as
\begin{equation}
p(\bm{y}|\bm{\theta}) = \prod_{n=1}^{M}
P_{\bm{\theta}}(y_n|\bm{\sigma}_n)^{y_n}(1-P_{\bm{\theta}}(y_n|\bm{\sigma}_n))^{1-y_n}
\end{equation}
where $\bm{y}=(y_1,...,y_{M})$. Instead of attempting to maximize the likelihood, it is convenient to define a loss function by taking the negative logarithm of the likelihood, which gives the cross-entropy
\begin{widetext}
\begin{equation}
E(\bm{\theta}) = -\ln( p(\bm{y}|\bm{\theta})) - = \sum_{n=1}^{M} \{ y_n \ln \left(P_{\bm{\theta}}(y_n|\bm{\sigma}_n))\right) + (1-y_n)\ln(1 - P_{\bm{\theta}}(y_n|\bm{\sigma}_n)) \}.
\label{Eq::NLL}
\end{equation}
\end{widetext}
To train the model, we minimize $E(\bm{\theta})$ using gradient descent techniques~\cite{10.5555/1162264}. While it is possible to evaluate the gradients of $E(\bm{\theta})$ with respect to the parameters $\bm{\theta}$ in the CNN analytically using the chain rule, a more convenient and less error-prone approach is to use automatic differentiation (AD), which is a set of techniques to numerically evaluate the derivative of a function specified by a computer program. A complete survery detailing AD can be found in Ref.~\cite{AD_2017}.
In addition, instead of using the entire dataset in the calculation of $E(\bm{\theta})$ and its gradients, we use smaller batches of data of size $M_{\text{batch}}<M$, which means that the gradients used during optimization become stochastic since they fluctuate from batch to batch. In the examples below we use $N_{\text{batch}}=32$. The gradient update rule used in our examples is a modified version of the usual gradient descent called Adam~\cite{2014arXiv1412.6980K}.
\subsubsection*{Code walkthrough}
We demonstrate supervised learning with a CNN to learn the quantum phase transition in the Rydberg array, using the machine learning software library TensorFlow~\cite{tensorflow}. We first generate training data by sampling the MPS wavefunctions obtained from DMRG at different detunings $\delta$. This data is then divided into a training set (used to update the neural network parameters) and a test set (used to validate the performance of the model). Each data set consists of a list of atomic occupation patterns $\bm\sigma$ and their ``phase label'' $y$.
We begin by importing the required functionalities and loading the data. Since we are using a CNN with a two-dimensional geometry, the atomic configurations in the training and test data sets need to be appropriately reshaped from the one-dimensional MPS structure.
\begin{lstlisting}[language=Python,numbers=none]
import numpy as np
import tensorflow as tf
from tensorflow.keras import layers, models
# linear dimensions of the system
Lx = 8
Ly = 8
# training set
train_config = np.loadtxt("xtrain.txt")
train_label = np.loadtxt("ytrain.txt",dtype=np.uint8)
# test set
test_config = np.loadtxt("xtest.txt")
test_label = np.loadtxt("ytest.txt",dtype=np.uint8)
# reshaping the training/test configurations
train_config = np.reshape(train_config,(train_config.shape[0],Lx,Ly,1))
test_config = np.reshape(test_config,(test_config.shape[0],Lx,Ly,1))
# reshaping the training/test labels
test_label = np.reshape(test_label,(test_label.shape[0],1))
train_label = np.reshape(train_label,(train_label.shape[0],1))
# Names for the phases
class_names = ["disordered", "ordered"]
\end{lstlisting}
Next, we define the neural network architecture by combining layers pre-defined in TensorFlow. After initialization, we proceed to implement the model of Eq.~\ref{Eq::CNNdist}. First, the raw input data is processed by two stacked convolutional layers, each one with $3\times3$ filters and 32 channels using rectified linear units. The output of the second CNN layer is fed to two stacked fully-connected layers with $64$ hidden units, and two output units. These last units correspond to model output for the disordered and ordered phase.
\begin{lstlisting}[language=Python,numbers=none]
# initialize the model
model = models.Sequential()
# first convolutional layer
model.add(layers.Conv2D(32, (3, 3),
activation = "relu",
input_shape = (Lx, Ly, 1)))
# second convolutional layer
model.add(layers.Conv2D(32, (3, 3),
activation="relu"))
# flatten the output
model.add(layers.Flatten())
# dense layer with 64 output units
model.add(layers.Dense(64, activation='relu'))
# dense layer with two output units
model.add(layers.Dense(2))
\end{lstlisting}
Once the architecture is defined, the model can be compiled by adding the cost function (i.e. the cross-entropy in Eq.~\ref{Eq::NLL} ), and the optimizer to update the model parameters, the Adam optimizer~\cite{2014arXiv1412.6980K}. Internally, TensorFlow builds a computational graph containing each operation being executed from the input state to the final output of the architecture. The gradients of the cost function with respect to each network parameter are then evaluated using AD. At compilation time, we can also add a metric to be monitored during training, in this case being the classification accuracy, i.e. the fraction of the test set samples that are being classified correctly.
\begin{lstlisting}[language=Python,numbers=none]
# Compiling the model
model.compile(optimizer='adam',
loss=tf.keras.losses.SparseCategoricalCrossentropy(from_logits=True),
metrics=['accuracy'])
\end{lstlisting}
The model is now ready to be trained using the Rydberg data for a set number of epochs, which is the number of passes of the entire training data set the machine learning algorithm has completed. After training is complete, we can evaluate the model on the held-out test data to quantify its accuracy.
\begin{lstlisting}[language=Python,numbers=none]
# training the model
history = model.fit(
train_config,
train_label,
epochs=5,
validation_data = (test_config, test_label))
# evaluate the overall model
test_loss,test_acc = model.evaluate(test_config, test_label, verbose=2)
\end{lstlisting}
We show the results of the training in Fig.~\ref{Fig::supervised}, evaluated on various test sets at different detuning values $\delta$. We plot the average output signal for the two output neurons, i.e. an estimate of $f(y,\delta) = \sum_{\bm \sigma} |\Psi(\bm{\sigma})|^2 P_{\bm{\theta}}(y|\bm{\sigma})$, on the top-most dense layer. When the detuning is large and negative, the disordered neuron saturates to one, and the ordered neuron is nearly zero, while the signals reverse at large and positive detuning. We can use the crossing point between the two curves to detect the critical point. We also show the accuracy in the test set, which shows a dip near the critical point. This is when the neural network is most uncertain about assigning a phase label to any given atomic configuration.
\begin{figure}[t]
\noindent \centering{}\includegraphics[width=\columnwidth]{fig_supervised}
\caption{Learning the quantum phase transition in the Rydberg-atom array. ({\bf a}) Output signal for the two units in the top-most dense layer of the neural-network architecture as a function of the detuning. ({\bf b}) Accuracy on the test data set, showing a dip near the critical point. This signal can be used to detect the quantum phase transition.}
\label{Fig::supervised}
\end{figure}
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
% QUANTUM STATE TOMOGRAPHY
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
\section{Quantum state tomography}
\label{qst}
Quantum characterization, verification and validation is a framework for algorithms and routines used to assess the quality and the performance of experimental quantum hardware and characterize its components~\cite{eisert_quantum_2020}. The workflow underlying these algorithms is inherently data-driven: appropriate measurement data is first collected from the quantum device under examination, and then processed by an algorithm running on a classical computer. Depending on the degree of complexity of the algorithm, different amounts of information can be gained. This could be a single figure of merit, such as the average error rate for a set of quantum gates~\cite{PhysRevA.77.012307,PhysRevA.85.042311,PhysRevA.99.052350}, or the fidelity (or a proxy thereof) between a quantum state prepared by the hardware and a desired reference state~\cite{PhysRevLett.106.230501,aolita_reliable_2015,PhysRevLett.120.190501}. One may be also interested in retrieving the full quantum state generated by a device. This procedure -- the reconstruction of an unknown quantum state from measurement data -- is called {\it quantum state tomography} (QST)~\cite{vogel89,PhysRevA.55.R1561,PhysRevA.63.040303,PhysRevA.64.052312,Jezek2003,Blume_Kohout_2010,PhysRevLett.108.070502,Granade_2017}.
%, and it is a fundamental ingredient of any infra-structure for quantum benchmarking and hardware calibration.
There are two assumptions in QST: the ability to prepare many identical copies of the quantum state $\bm{\varrho}$ of interest, and to repeatedly perform measurements on it. The set of measurements is in general described by positive-operator value measures (POVM) $\linebreak\bm{\mathcal{M}}=\{\Pi_k\}$~\cite{nielsen_chuang_2010}. The set $\bm{\mathcal{M}}$ is said to be {\it informationally complete} (IC) if it spans the full Hilbert space. An example of informationally (over)-complete measurements are Von Neumann measurements in the Pauli bases, where for a single qubit $\bm{\mathcal{M}}$ contains the six rank-1 projectors into the eigenstates of the Pauli matrices. For an IC set, any quantum state can be uniquely identified by the probabilities of the measurements in $\bm{\mathcal{M}}$, as specified by the Born rule $p(k) = \tr[\bm{\varrho}\:\Pi_k]$.
The simplest method to perform QST is {\it linear inversion}, which reconstructs the quantum state simply by inverting the Born rule using an empirical approximation of the measurement probabilities. One issue of linear inversion is that the reconstructed density operator $\bm{\rho}$ is not necessarily positive, although negative eigenvalues can be appropriately removed to produce a positive state that is closest to the output of linear inversion~\cite{PhysRevLett.108.070502}. A more powerful approach, but also more computationally intensive, is {\it maximum likelihood estimation}~\cite{PhysRevA.55.R1561,PhysRevA.63.040303,PhysRevA.64.052312}, where the state $\bm{\rho}$ is found by minimizing the likelihood function for the observed data under the constraint $\bm{\rho}\ge0$.
Traditional QST algorithms based on linear inversion or maximum likelihood suffer a complexity that scales exponentially with the number of qubits or particles involved. This exponential scaling stems from two reasons. First, the representation of the quantum state $\bm{\rho}$, which is inevitably exponential in the system size. Second, the {\it sample complexity}, that is the number of measurements that needs to be collected in an experiment to achieve a faithful reconstruction of the quantum state. Typically, statistics from a IC set is required in order to get a good fit, and the size thereof scales exponentially with the number of qubits. For these reasons, traditional QST has remained limited to quantum systems containing only a small number of particles~\cite{haffner_scalable_2005}.
Several algorithms to overcome this severe complexity have been proposed over the last decade. Notable examples are compressed sensing tomography~\cite{PhysRevLett.105.150401,Flammia_2012,PhysRevLett.106.100401,riofrio_experimental_2017}, permutationally invariant tomography~\cite{PhysRevLett.105.250403,Moroder_2012} and tensor-network tomography~\cite{cramer2009efficient,MPOtomo,Lanyon2017,LeiWang2020}, which rely respectively on the sparsity, translational invariance and low-entanglement of the target quantum state. More recently, a new framework built on neural networks and unsupervised learning has been put forward~\cite{torlai_Tomo}, based on the assumption that most physical states of interest typically contains some degree of {\it structure} (i.e. correlations, symmetries, etc), in the sense that they can be described using a reduced number of parameters (much smaller than the dimension of the Hilbert space). The general idea is to leverage the capability of unsupervised machine learning to autonomously identify such structure in raw data, and compress it using a neural-network representation of the quantum state.
In what follows, we will focus on pure quantum states, and discuss the extension to mixed states at the end of this Section. We consider a system of $N$ qubits (or any other two-level system) described by a wavefunction $|\Phi\rangle$ with amplitudes $\Phi(\bm{\sigma}) = \langle\bm{\sigma}|\Phi\rangle$ in an appropriate reference basis $|\bm{\sigma}\rangle=|\sigma_1,\dots,\sigma_N\rangle$ ($\sigma_j\in\{0,1\}$). In order to circumvent the scalability issue of standard QST, we adopt a compact representation of a wavefunction expressed in terms of a neural network~\cite{Carleo_2017}. The resulting {\it neural-network wavefunction} is simply a highly non-linear parametric function of the basis states $\psi_{\bm{\theta}}(\bm{\sigma})$, where $\bm\theta$ is a set of parameters (e.g. weights and biases). Several types of neural networks have been successfully implemented to perform QST, including feed-forward neural networks~\cite{biamonte_qst,xin_local-measurement-based_2019}, variational auto-encoders~\cite{rocchetto}, generative adversarial networks~\cite{NoriGAN}, recurrent neural networks~\cite{carrasquilla_povm,morawetz2020} and transformers~\cite{Cha2020}. Here, we will examine the restricted Boltzmann machine.
\subsection{The restricted Boltzmann machine}
The restricted Boltzmann machine is an energy-based model introduced in the early 1980s for generative modeling~\cite{Ackley85,Smolensky1986}, built on a connection between cognitive science and statistical mechanics~\cite{Little74,Little78,Hopfield82}. The RBM features two layers of stochastic binary units: a visible layer $\bm{\sigma}=(\sigma_1,\sigma_2,\dots)$ and a hidden layer $\linebreak\bm{h}=(h_1,h_2,\dots)$, containing respectively $N$ and $n_h$ neurons (or units). The two layers in the RBM are fully connected by a symmetric weight matrix $\bm{W}$, with no intra-layer connections (hence its {\it restricted} nature). The visible units are used to represent the variables relevant to the specific problem at hand, such as the pixel values in an image (or the computational basis states for qubits). The size of the hidden layer is a natural control parameter for the representational power of the model. Note that, since RBMs are universal function approximators~\cite{LeRoux2008}, they can capture any discrete distribution provided the number of hidden units is sufficiently large (possibly exponential in the number of visible units).
The RBM associates to each configuration of the visible and hidden layer $(\bm \sigma,\bm h)$ the energy
\begin{equation}
E_{\bm{\theta}}(\bm{\sigma},\bm{h})=-\sum_{j}\sum_{i}W_{ij}h_i\sigma_j-\sum_{j}b_j\sigma_j-\sum_{i}c_ih_i\:,
\end{equation}
where $\bm\theta=(\bm{W},\bm{b},\bm{c})$ is the set of parameters, and we also introduced biases $\bm b$ and $\bm c$ for the visible and hidden units respectively. Given this energy functional, the (stochastic) RBM units are distributed according to the Boltzmann distribution at temperature $\beta=1$
\begin{equation}
p_{\bm{\theta}}(\bm{\sigma},\bm{h})=Z_{\bm{\theta}}^{-1}e^{-E_{\bm{\theta}}(\bm{\sigma},\bm{h})}\:,
\end{equation}
where the partition function is
\begin{equation}
Z_{\bm{\theta}}=\sum_{\bm{\sigma},\bm{h}}e^{-E_{\bm{\theta}}(\bm{\sigma}_,\bm{h})}\:.
\end{equation}
Importantly, because the network architecture is restricted, we can trace out the latent space explicitly, obtaining the marginal probability distribution over the visible space
\begin{equation}
p_{\bm{\theta}}(\bm{\sigma})=\sum_{\bm{h}}p_{\bm{\theta}}(\bm{\sigma},\bm{h}) = Z^{-1}_{\bm{\lambda}}e^{\mathcal{E}_{\bm\theta}(\bm\sigma)}\:,
\label{Eq::RBM}
\end{equation}
where we defined an ``effective energy''
\begin{equation}
\mathcal{E}_{\bm\theta}(\bm\sigma)=\sum_jb_j\sigma_j+\sum_{i}\left(1+e^{\:\sum_{j}W_{ij}\sigma_j+c_i}\right)\:.
\end{equation}
The main purpose of the RBM is generative modeling, which is the task of learning a representation of an unknown probability distribution from data, allowing the neural network to produce new data points. In other words, the RBM training attempts to discover low-dimensional features in the data to allow generalization beyond the finite-size data set.
Let us consider a data set $\mathcal{D}=\{\bm{\sigma}_k\}$ with underlying (unknown) probability distribution $q(\bm \sigma)$. The RBM can be trained using unsupervised learning to minimize the distance between the two distributions. Such distance measure is typically expressed in terms of the Kullback-Leibler divergence~\cite{Kullback:1951aa}
\begin{equation}
\text{KL}(q\:|\:p_{\bm{\theta}}) = \sum_{\bm{\sigma}}q(\bm{\sigma})\log\frac{q(\bm{\sigma})}{p_{\bm{\theta}}(\bm{\sigma})}\:,
\end{equation}
with $\text{KL}(q\:|\:p_{\bm{\theta}})>0\:\:\forall q,p_{\bm{\theta}}$ and $\text{KL}(q\:|\:p_{\bm{\theta}})=0$ iff $p_{\bm{\theta}}=q$.
The exponentially large sum over the full configuration space is approximated using the available data, leading to the cost function
\begin{equation}
\mathcal{C}(\bm{\theta})=
-\frac{1}{|\mathcal{D}|}\sum_{\bm{\sigma}\in\mathcal{D}}\log p_{\bm{\theta}}(\bm{\sigma})-H_{\mathcal{D}}\:,
\label{Eq::KL}
\end{equation}
where $|\mathcal{D}|$ is the size of the data set. Note that, up to a constant data set entropy term $H_{\mathcal{D}}$, the KL divergence simply reduces to the negative logarithm of the likelihood function $\mathcal{L}(\mathcal{D}\:|\:p_{\bm\theta})$.
The RBM can be trained using one of the many flavors of gradient descent. It is straightforward to show that the gradients of the cost function are given by\begin{equation}
\nabla_{\bm\theta}\mathcal{C}({\bm\theta}) = \langle\nabla_{\bm\theta} \mathcal{E}_{\bm\theta}(\bm\sigma)\rangle_{p_{\bm\theta}} - \langle\nabla_{\bm\theta} \mathcal{E}_{\bm\theta}(\bm\sigma)\rangle_{\mathcal{D}}
\end{equation}
where the gradients $\nabla_{\bm\theta}\mathcal{E}_{\bm\theta}(\bm\sigma)$ can be computed exactly for any sample $\bm{\sigma}$. We see that the gradients $\nabla_{\bm\theta}\mathcal{C}({\bm\theta})$ have two components. First, there is the average over the data points $\langle\nabla_{\bm\theta} \mathcal{E}_{\bm\theta}(\bm\sigma)\rangle_{\mathcal{D}}$, which is fast to compute. In contrast, the average over the model distribution
\begin{equation}
\Big\langle\frac{\partial \mathcal{E}_{\bm\theta}(\bm\sigma)}{\partial \bm\theta}\Big\rangle_{p_{\bm\theta}}=
\frac{1}{Z_{\bm\theta}}\sum_{\bm\sigma}e^{\mathcal{E}_{\bm\theta}(\bm\sigma)}\nabla_{\bm\theta}\mathcal{E}_{\bm\theta}(\bm\sigma)
\end{equation}
requires the partition function, whose calculation is in general intractable. However, this expectation value can be approximated using Monte Carlo, by drawing $N_S$ samples from the model distribution $\{\bm\sigma_i\}\sim p_{\bm\theta}(\bm\sigma)$:
\begin{equation}
\Big\langle\frac{\partial \mathcal{E}_{\bm\theta}(\bm\sigma)}{\partial \bm\theta}\Big\rangle_{p_{\bm\theta}}\approx
\frac{1}{N_S}\sum_{i=1}^{N_S}\nabla_{\bm\theta}\mathcal{E}_{\bm\theta}(\bm\sigma_i)\:.
\end{equation}
This is the most computationally intensive step of the training, and depending on the specific distribution to be learned, advanced Monte Carlo algorithms may be required to collect sufficiently uncorrelated samples.
\subsection{Reconstruction of Rydberg atoms}
Now that we have introduced the main features of the RBM and its training, we are ready to explore its use for QST. As a first application, we examine the reconstruction of Rydberg-atom wavefunctions. An important property of the ground state wavefunction $|\Phi\rangle$ of the Rydberg Hamiltonian (\ref{Eq::RydbergHamiltonian}) is that it is positive in the occupation number basis $\Phi(\bm\sigma)\ge0$ (where $\sigma_j=0$ and $\sigma_j=1$ refers to the ground and excited states). This property follows directly from the representation of the Hamiltonian in this basis, in which all of its off-diagonal elements can be gauged to be negative (i.e. the Hamiltonian is {\it stoquastic} in this basis~\cite{Bravyi2006}).
The positivity of the target state implies that we may parametrize the neural-network wavefunction simply as $\psi_{\bm\theta}(\bm\sigma) = \sqrt{p_{\bm\theta}(\bm\sigma)}$, for any normalized probability distribution $p_{\bm\theta}(\bm\sigma)$. Here we choose the RBM probability distribution (Eq.~\ref{Eq::RBM}), where the visible units correspond to the atomic occupations. Moreover, because the wavefunction is positive, measurement data from a single basis is sufficient to characterize the state. This means that the QST problem, under these assumptions, is equivalent to unsupervised learning of projective measurement data in the atomic occupation number basis. The tomographic reconstruction of the quantum state is carried out by iteratively changing the RBM parameters to minimize the cost function
\begin{equation}
\mathcal{C}(\bm{\theta})=
%-\frac{1}{|\mathcal{D}|}\sum_{\bm{\sigma}\in\mathcal{D}}\log |\langle \bm \sigma|\psi_{\bm\theta}\rangle|^2=
-\frac{1}{|\mathcal{D}|}\sum_{\bm{\sigma}\in\mathcal{D}}\log p_{\bm\theta}(\bm\sigma)\:.
\label{Eq::KL}
\end{equation}
Upon reaching convergence in the training, the learned RBM can be used to estimate various properties of interest. For a generic observable $\hat{\mathcal{O}}$, the expectation value on the RBM wavefunction reduces to an average over the RBM distribution
\begin{equation}
\langle\hat{\mathcal{O}}\rangle = \frac{\langle\psi_{\bm\theta}|\hat{\mathcal{O}}|\psi_{\bm\theta}\rangle}{\langle\psi_{\bm\theta}|\psi_{\bm\theta}\rangle} =
\frac{1}{Z_{\bm\theta}}\sum_{\bm{\sigma}}p_{\bm\theta}(\bm\sigma)\mathcal{O}_{loc}(\bm\sigma)\:,
\label{Eq::avg_obs}
\end{equation}
where we introduced the so-called ``local observable''
\begin{equation}
\mathcal{O}_{loc}(\bm\sigma) = \frac{\langle\bm\sigma|\hat{\mathcal{O}}|\psi_{\bm\theta}\rangle}{\langle\bm\sigma|\psi_{\bm\theta}\rangle}\:.
\end{equation}
The expectation value in Eq.~\ref{Eq::avg_obs} can then be approximated with a Markov chain using Monte Carlo sampling, similarly to the evaluation of the gradients. We point out that the evaluation of the local observable $\mathcal{O}_{loc}(\bm\sigma)$ remains efficient as long as the matrix representation of $\mathcal{O}$ in the reference basis is sufficiently sparse. This measurement procedure is also useful for monitoring different observables during training, such as average densities and correlation functions. These metrics can be used to assess convergence, since in general the calculation of the KL divergence is intractable as it also requires estimating the partition function $Z_{\bm\theta}$.
\begin{figure}[t]
\noindent \centering{}\includegraphics[width=\columnwidth]{fig_qst_rydberg}
\caption{Quantum state tomography of a $8\times8$ array of Rydberg atoms with unsupervised learning of single-shot atomic occupation data. We compare various observables measured using the MPSs obtained from DMRG (blue line) and the neural-network wavefunctions at different values of the detuning $\delta$. We plot the average energy per spin ({\bf a}), and the average magnetization along the $z$ ({\bf b}) and $x$ ({\bf c}) axes.}
\label{Fig::QST_rydberg}
\end{figure}
\subsubsection*{Code walk-through}
We perform QST on the Rydberg atoms data using the Python package {\it NetKet}~\cite{netket}. First, we import the library and define the relevant parameters for the numerical experiments. For instance, we consider the reconstruction of a square array with linear size $L=8$ (containing $N=64$ spins), with Hamiltonian parameters $V=3.0$, $\Omega=1.0$ and $\delta=2.0$.
\begin{lstlisting}[language=Python,numbers=none]
from mpi4py import MPI
import netket as nk
# Rydberg Hamiltonian parameters
L = 8 # linear size
V = 3.0 # Van der Waals interaction
Omega = 1.0 # Rabi frequency
delta = 2.0 # detuning
\end{lstlisting}
We then define the lattice structure and the Hilbert space for the neural-network wavefunction, and load the training data from file:
\begin{lstlisting}[language=Python,numbers=none]
# define the lattice structure
square_lattice = nk.graph.Hypercube(
length=L,
n_dim=2,
pbc=False)
# build the Hilbert space
hilbert = nk.hilbert.Qubit(graph=square_lattice)
# load the training data
data_path = ".../path_to_data"
rotations, samples, bases = LoadData(hilbert, data_path)
\end{lstlisting}
Note that, since we are training in a single measurement basis, the arrays \texttt{bases} and \texttt{rotations} are ``trivial''. Otherwise, these variables would contain respectively an integer encoding of each distinct basis, and its associated (local) unitary rotations.
Next, we define the main components of the QST algorithm: the neural-network wavefunction, the sampler used to approximate the gradients, and the optimizer for the parameter updates.
\begin{lstlisting}[language=Python,numbers=none]
# Neural-network wavefunction
rbm = nk.machine.RbmSpinReal(hilbert=hilbert, alpha=1)
# Monte Carlo sampler
sa = nk.sampler.MetropolisLocal(machine=rbm)
# Optimizer
op = nk.optimizer.AdaDelta(rho=0.95, epscut = 1.0e-7)
# Initialize tomography object
qst = nk.Qsr(
sampler = sa,
optimizer = op,
n_samples_data = 1000,
n_samples = 2000,
rotations = rotations,
samples = samples,
bases = bases,
sr = None)
\end{lstlisting}
In the definition of the RBM (with real-valued network parameters), the parameter $\alpha = n_h/N$ represents the density of hidden units. We use the AdaDelta optimizer~\cite{2012arXiv1212.5701Z} and a Metropolis sampler using simple single-spin flips. The tomography parameters \texttt{n\_samples\_data} and \texttt{n\_samples} refer respectively to the number of training samples and the number of samples drawn from the model distribution to compute the gradients for a single parameter update (i.e. for batch gradient descent). We refer the reader to the NetKet documentation for additional details. Finally, we generate the Rydberg Hamiltonian to be measured during the learning, and run the QST for a fixed number of training iterations (\texttt{epochs}).
\begin{lstlisting}[language=Python,numbers=none]
# define observable for measurements
H = generatehamiltonian(hilbert, L, L, V, Omega, delta)
qst.add_observable(H, "H")
# run quantum state tomography
for ep in qst.iter(epochs):
obs = qst.get_observable_stats()
\end{lstlisting}
We perform QST on datasets of projective measurements generated using the MPS obtained from DMRG at different detunings. Each data set contains $10^5$ measurements. We train each RBM separately using the hyper-parameters reported above, and measure at each training iteration various observables of interest. We show the results in Fig.~\ref{Fig::QST_rydberg}, where we plot the average energy per spin $\langle\hat{H}\rangle/N$ and the average magnetizations $\langle\hat{S}^{z/x}\rangle = \sum_j\langle\hat{S}_j^{z/x}\rangle/N$ after the training has converged. Each data point is obtained by averaging the expectation values of the observables over the last 100 iterations of the training. The reconstruction show an overall good agreement with the exact values computed using the MPS wavefunctions.
\begin{figure}[t]
\noindent \centering{}\includegraphics[width=\columnwidth]{fig_qst_beryllium}
\caption{Quantum state tomography of the Beryllium Hydride molecule. We show the fidelity $\mathcal{F} = |\langle\Phi|\psi_{\bm\theta}\rangle|^2$ between the neural-network wavefunction $\psi_{\bm\theta}$ and the exact molecular wavefunction $\Phi$ at each iteration during training.}
\label{Fig::QST_beryllium}
\end{figure}
\subsection{Reconstruction of a molecular wavefunction}
We now move to the more general case of a wavefunction that is non-positive or complex-valued. In order to accommodate this different setup, we first need to modify the variational ansatz to allow for complex-valued amplitudes $\psi_{\bm\theta}(\bm\sigma)$. One way to achieve this is to use the RBM to parametrize the probability distribution in the reference basis (as before), and couple it with an additional RBM that parametrizes the phases, $\psi_{\bm\theta\bm\mu}(\bm\sigma)=\sqrt{p_{\bm\theta}(\bm\sigma)}e^{i\log p_{\bm\mu}(\bm\sigma)}$~\cite{torlai_Tomo}. A different strategy, which we adopt here, is to promote the weights and biases to have complex values~\cite{Carleo_2017}. Note that, for the latter representation, we cannot interpret the RBM as a probabilistic graphical model anymore.
The presence of phases in the target wavefunction leads to an additional overhead in the measurement requirements. In fact, projective measurements in the computational basis do not carry enough information to uniquely identify the state. In order to obtain information about the phases, measurement in additional bases are required. To account for this, we can write the training data set as $\mathcal{D}=\{\bm{x}_k\}$, where each single-shot measurement is given by $\bm{x}=(\bm{\tau},\bm{\sigma})$, with $\bm\tau$ and $\bm\sigma$ referring respectively to the measurement basis and the binary measurement outcome. For example, for the case of Pauli measurements, the data point $\bm{x} = (\sigma_1^x=0,\sigma^z_2=1)$ refers to a measurement of qubit 1 (2) in the eigenbasis of the Pauli-X (Z) operator.
The learning algorithm proceeds in a similar fashion to the case of a positive wavefunction, reducing to the optimization of the negative log-likelihood cost function
\begin{equation}
\mathcal{C}(\bm{\theta})=-\frac{1}{|\mathcal{D}|}\sum_{k=1}^{|\mathcal{D}|}\log p_{\bm{\theta}}(\bm{x}_k)\:.
\end{equation}
The important difference is that, in order to evaluate measurement probabilities in bases other than the reference one, we need to appropriately rotate the neural-network wavefunction. We can write this measurement probability as
\begin{equation}
p_{\bm{\theta}}(\bm{x}_k) = \frac{|\langle\bm{\sigma}_k|\bm{U}(\bm{\tau}_k)|\psi_{\bm\theta}\rangle|^2}{\langle\psi_{\bm\theta}|\psi_{\bm\theta}\rangle}\:,
\end{equation}
where $\bm{U}(\bm{\tau})$ refers to the unitary transformation that rotates the reference basis into the measurement basis $\bm{\tau}$. We also assume local measurements, leading to the factorization $\bm{U}(\bm{\tau})=\bigotimes_{j=1}^NU(\tau_j)$ into single-qubit unitary rotations $U(\tau_j)$. For the example described above, with $\linebreak\bm{x} = (\sigma_1^x=0,\sigma^z_2=1)$, the neural-network probability is (up to a normalization factor) $p_{\bm\theta}(\bm{x})\propto|\langle01|H_1\otimes\mathbb{1}_2|\psi_{\bm\theta}\rangle|^2$, where $H$ is the Hadamard gate (i.e. the rotation into the $\sigma^x$ basis).
We show neural-network QST of a non-positive wavefunction for the case of the electronic ground state of a small molecule. The Hamiltonian in second quantization is given by
\begin{equation}
\hat{H} = \sum_{\alpha,\beta}t_{\alpha\beta}\hat{c}^\dagger_\alpha \hat{c}_\beta + \frac{1}{2}\sum_{\alpha,\beta,\gamma,\delta}u_{\alpha\beta\gamma\delta} \hat{c}^\dagger_\alpha \hat{c}^\dagger_\beta \hat{c}_\gamma \hat{c}_\delta
\end{equation}
where $\hat{c}^\dagger$ and $\hat{c}$ are fermionic creation and annihilation operators, and $t_{\alpha\beta}$ and $u_{\alpha\beta\gamma\delta}$ are electronic integrals. This Hamiltonian can be mapped into a qubit Hamiltonian using one of several mappings (i.e. Jordan-Wigner, Bravyi-Kitaev etc)
\begin{equation}
\hat{H} = \sum_k c_k \hat{P}_k\:.
\end{equation}
where $c_k$ are interaction coefficients and $\hat{P}_k$ are operators that belongs to the $N$-qubit Pauli group.
We specifically look at the Beryllium Hydride molecule (\BeH) in the STO-3G basis, mapped to $N=6$ qubits~\cite{kandala_hardware-efficient_2017}. In order to generate the training data, we first obtain the full ground state wavefunction $|\Phi\rangle$ with exact diagonalization of the qubit Hamiltonian. Given the ground state, we can generate measurement data by sampling the full probability distribution obtained from the Born rule. A single measurement is obtained by first selecting a measurement basis $\bm\tau$, then rotating the wavefunction accordingly $|\Phi_{\bm\tau}\rangle=\bm{U}(\bm\tau)|\Phi\rangle$, and finally sampling the measurement outcome $\bm\sigma\sim P_{\bm\tau}(\bm{\sigma})=|\langle\bm\sigma|\Psi_{\bm\tau}\rangle|^2$. We choose the measurement bases according to the Pauli operators $\hat{P}_k$ appearing the Hamiltonian, each one being selected randomly among this set. We show the results of the QST experiment in Fig.~\ref{Fig::QST_beryllium}, where we plot the fidelity between the neural-network wavefunction and the target wavefunction during the training.
\subsection{Discussion}
We have shown that for wavefunctions whose amplitudes can be gauged to be real and positive, QST is equivalent to unsupervised learning of projective measurement data in a single basis. For more general wavefunctions with a sign structure or complex-valued amplitudes, measurements in multiple bases are required to reconstruct the phases. These are processed by the neural network by appropriately rotating the parametrized wavefunction with a unitary $\bm{U}(\bm{\tau})=\bigotimes_{j=1}^NU(\tau_j)$ composed by single-qubit basis rotations $U(\tau_j)$. In practice, this operation entails an exponential cost in the number of non-trivial rotations $U(\tau_j)\ne\mathbb{1}_j$, which means that only a small fraction of any IC set of bases can be use to train the neural network. Nevertheless, this reduced amount of information may be enough to reconstruct sufficiently structured quantum states. For example, a ground state of a local gapped Hamiltonian can be identified by the statistics of projective measurements in a set of bases corresponding to the decomposition of the Hamiltonian in the Pauli group.
An important assumption that was made is the purity of the quantum state under reconstruction, which is violated in any practical setting. There are cases where an approximate pure state reconstruction of an experimental quantum state may be justified, and could still provide valuable insights~\cite{torlai_rydberg19}. However, for benchmarking and noise characterization tasks, one needs to reconstruct the full density matrix. This can be achieved by introducing a neural-network parametrization of a density operator $\rho_{\bm\theta}(\bm\sigma,\bm\sigma^\prime)$, which is trained in an analogous way. The positivity of $\rho_{\bm\theta}$ can be enforced using a purification scheme, where $\rho_{\bm\theta}$ is purified by additional latent units in the neural network~\cite{Torlai_latent}. An iterative procedure to discover the dominant eigenstates of a density operators using RBMs has also been put forward~\cite{PhysRevA.102.022412}.
A different approach for reconstructing generic mixed quantum states consists of directly parametrizing the probability distribution $p(\bm\alpha) = \tr[\bm{\varrho}\:\Pi_{\bm\alpha}]$ of an IC-POVM set $\{\Pi_{\bm\alpha}\}$~\cite{carrasquilla_povm}. Contrary to the purification scheme, this approach does not require unitary rotations, and thus lifts the exponential complexity associated with processing data from arbitrary local bases. This however comes at the cost of possibly violating the positivity of the learned density matrix, since this cannot be enforced at the level of the POVM distribution parametrization.
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
% VARIATIONAL MONTE CARLO
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
%----------------------------------------------------------------------------------------
\section{Variational ground state optimization}
\label{vmc}
The variational principle in quantum mechanics states that the expectation value of the Hamiltonian of a physical system of interest over any valid wavefunction is always greater than or equal to the ground state energy of the system. This principle indicates that a strategy for finding an approximation to the ground state energy of a system is to start from a parameterized wavefunction and vary its parameters until it yields the minimum possible energy.
Historically, the choice of wavefunction, which is critical to the success of the algorithm, has traditionally been made motivated in close connection to a physical understanding of the problem, e.g., from approximate mean-field solutions supplemented with some form of additional correlation such as a Jastrow factor~\cite{becca_sorella_2017}. More recently, motivated by their representation power, efficiency, and generality, neural networks have been explored as trial wavefunctions~\cite{Carleo_2017}. In particular, and as we explore below, recurrent neural networks are naturally well-suited to the study of systems exhibiting strong correlations such as those arising in the study of classical and quantum systems, which are prevalent in condensed matter and statistical physics~\cite{Wu_2019,carrasquilla_povm,RNNWF_2020, roth2020iterative}.
\subsection{Recurrent neural networks}
A recurrent neural network (RNN) models a probability distribution $p(\bm{\sigma})=p(\sigma_1,\dots,\sigma_N)$ using a sequential structure according to the chain rule of probabilities
\begin{equation}
p(\bm{\sigma})=p(\sigma_1)\:p(\sigma_2\,|\,\sigma_1)\:\prod_{j=3}^Np(\sigma_j\,|\bm{\sigma}_{-j}),
\end{equation}
where $\bm{\sigma}_{-j}=(\sigma_1,,\dots,\sigma_{j-2},\sigma_{j-1})$. The recurrent unit in the neural network parametrizes the conditional probability distribution $p(\sigma_j\,|\,\bm{\sigma}_{-j})$ at any {\it time-step} $j$, and it processes the data according to the input $\sigma_{j-1}$ and a recurrent latent vector $\bm{h}_j$ (Fig.~\ref{Fig::rnn}a).
We specifically consider Gated Recurrent Units (GRU)~\cite{cho-etal-2014-learning}, introduced to solve the vanishing gradient issue of vanilla RNNs. The schematic of a GRU units is shown in (Fig.~\ref{Fig::rnn}b). Given a {\it visible} input state $\sigma_{j-1}$ and a latent state $\bm{h}_{j-1}$, the GRU unit at time-step $j$ processes them according to the sequence of operations and outputs an updated latent vector $\bm{h}_j$. The first two operations are the so-called {\it update gate} and {\it reset gate}:
\begin{align}
\bm{z}_j = \text{sig}\Big(\bm{W}_z\:[\bm{h}_{j-1}; \sigma_{j-1}]+\bm{b}_z\Big)\\
\bm{r}_j = \text{sig}\Big(\bm{W}_r\:[\bm{h}_{j-1}; \sigma_{j-1}]+\bm{b}_r\Big)
\end{align}
where $[\bm{h}_{j-1}; \sigma_{j-1}]$ is the concatenation operation, $\text{sig}(x)=(1+e^{-x})^{-1}$ is the sigmoid function, and $\bm{W}_z$, $\bm{W}_r$ and $\bm{b}_z$, $\bm{b}_r$ are respectively weights and biases variational parameters. These gates are used to control how much information about previous time-steps is kept encoded into the latent vector.
Next, given the input visible state $x_j$ at the current time-step, an internal latent state is created according to
\begin{equation}
\bm{\tilde{h}}_j = \tanh\Big(\bm{\tilde{W}}\:[\bm{r}_j\odot\bm{h}_{j-1},x_j]+\bm{\tilde{b}}\Big)
\end{equation}
where $\bm{a}\odot\bm{b}$ is the element-wise multiplication, and $\bm{\tilde{W
}}$, $\bm{\tilde{b}}$ are a new set of parameters. The new latent vector (output of the GRU), is generated as
\begin{equation}
\bm{h}_j=(\bm{1}-\bm{z}_j)\odot\bm{h}_{j-1}+\bm{z}_j\odot\bm{\tilde{h}}_j\:,
\end{equation}
and sent to the GRU unit at time-step $(j+1)$. Finally, the conditionals are computed using a softmax layer
\begin{equation}
p(\sigma_j \,|\,\bm{\sigma}_{-j}) = S\left( U \bm{h}_j + \bm{c}\right)
\end{equation}
where $U$ and $\bm{c}$ are the parameters of the softmax layer.
Given the above parametrization of a probability distributions $p(\bm{\sigma})$, we can now promote RNNs to quantum mechanical wavefunctions $\psi(\bm{\sigma})$. As noted in the QST examples, we stress that {\it stoquastic} many-body
Hamiltonians have ground states $\ket{\Psi}$ with strictly real and positive amplitudes in the standard computational basis~\cite{Bravyi2006}.
This class states can be represented in terms of probability distributions,
\begin{align}
\ket{\Psi} = \sum_{\bm{\sigma}} \psi(\bm{\sigma})\ket{\bm{\sigma}} = \sum_{\bm{\sigma}}\sqrt{p(\bm{\sigma})}\ket{\bm{\sigma}}.
\end{align}
% We note that this property has been exploited in multiple wavefunction representations using generative models such as restricted Boltzmann machines~\cite{Carleo_2017}.
For such family of quantum states, which includes the Rydberg system considered in this work, it is natural to try to approximate $p(\bm{\sigma})$ with an RNN.
\begin{figure}[t]
\noindent \centering{}\includegraphics[width=\columnwidth]{rnn}
\caption{Recurrent neural network (RNN). {\bf a}). A generic RNN cell (left), and its unrolling in time to process sequenced data (right). {\bf b}) The gated recurrent unit (GRU), and the set of operations used to process the input latent state $\bm{h}_{j-1}$ and visible state $x_j$ to generate a new latent state $\bm{h}_j$. Lines joining together means concatenation, while the circles are element-wise operations.}
\label{Fig::rnn}
\end{figure}
\subsection{Variational Monte Carlo simulation of Rydberg atoms}
% \textcolor{red}{JC: Let's write the VMC part here, and just show the results at the end.}
The goal of variational Monte Carlo (VMC) is to iteratively optimize an ansatz wavefunction to approximate ground states of local Hamiltonians~\cite{becca_sorella_2017}. The VMC method make use of trial wavefunction $\ket{\Psi_{\bm{\theta}}}$ endowed with parameters $\bm{\theta}$. Here we consider an GRU-RNN wavefunction. Crucially, we exploit the fact that the RNN wavefunction allows for efficient sampling from the square of the amplitudes of $\ket{\Psi_{\bm{\theta}}}$.
The VMC method iteratively optimizes the expectation value of the
energy
\begin{align}
E \equiv \frac{\braket{\Psi_{\bm{\theta}}|\hat{H}|\Psi_{\bm{\theta}}}}{\braket{\Psi_{\bm{\theta}}|\Psi_{\bm{\theta}}}}.\label{eq:energy}
\end{align}
The minimization is carried out using the gradient descent method or any other variant of it. Since the RNN wavefunction is normalized such that $\braket{\Psi_{\bm{\theta}}|\Psi_{\bm{\theta}}}=1$, the expectation value in Eq.~\ref{eq:energy} can be written as
\begin{align}
E = \braket{\Psi_{\bm{\theta}} |\hat{H}|\Psi_{\bm{\theta}}}
&=\sum_{\bm{\sigma}} |\psi_{\bm{\theta}}(\bm{\sigma})|^2\sum_{\bm{\sigma'}} H_{\bm{\sigma\sigma'}}\frac{\psi_{\bm{\theta}}(\bm{\sigma'})}{\psi_{\bm{\theta}}(\bm{\sigma})} \nonumber \\
&\equiv \sum_{\bm{\sigma}} |\psi_{\bm{\theta}}(\bm{\sigma})|^2E_{loc}(\bm{\sigma}) \nonumber \\
&\approx \frac{1}{N_S}\sum_{\bm{\sigma} \sim |\psi_{\bm{\theta}}(\bm{\sigma})|^2} E_{loc}(\bm{\sigma}),
\label{eq:expectation_value}
\end{align}
which represents a sample average of the local energy $E_{loc}(\bm{\sigma})$. The gradients $\partial_{{\bm{\theta}}}E$ can be similarly written as
\begin{align}
\partial_{{\bm{\theta}}} E =
\sum_{\bm{\sigma}} |\psi_{\bm{\theta}}(\bm{\sigma})|^2\frac{\partial_{{\bm{\theta}}}\psi^{*}_{\bm{\theta}}(\bm{\sigma})}{\psi^{*}_{\bm{\theta}}(\bm{\sigma})} E_{loc}(\bm{\sigma}) + \text{c.c}.
\label{eq:gradient}
\end{align}
An optimization step involves drawing $N_S$ samples $\{ \bm{\sigma}^{(1)}, \bm{\sigma}^{(2)}, \ldots, \bm{\sigma}^{(N_S)}\}$ from $|\psi_{\bm{\theta}}(\bm{\sigma})|^2$, followed by an estimation of the energy gradients
\begin{align}
\partial_{{\bm{\theta}}} E \approx
\frac{2}{N_S} \mathfrak{Re} \left ( \sum_{i=1}^{N_S} \frac{\partial_{{\bm{\theta}}}\psi^{*}_{\bm{\theta}}(\bm{\sigma^{(i)}})}{\psi^{*}_{\bm{\theta}}(\bm{\sigma^{(i)}})} E_{loc}(\bm{\sigma^{(i)}}) \right),
\label{eq:stochastic_gradient1}
\end{align}
using automatic differentiation~\cite{AD_2017} and updating the parameters according to
\begin{align}
{\bm{\theta}} \leftarrow {\bm{\theta}} - \alpha \partial_{{\bm{\theta}}}E
\label{eq:gradient_descent}
\end{align}
with a small learning rate $\alpha$. Instead of this simple gradient descent rule, we can also use the Adam optimizer \cite{2014arXiv1412.6980K} to implement the parameter updates.
We note that the stochastic evaluation of the gradients in Eq.~\eqref{eq:stochastic_gradient1} implies that these may exhibit high variance, which can potentially slow down the convergence of the algorithm. This problem can be alleviated through the introduction of a term in Eq.~\eqref{eq:stochastic_gradient1} that helps reduce the variance of the gradients~\cite{RNNWF_2020}
\begin{align}
\partial_{\bm{\theta}} E &\approx
\frac{2}{N_S} \mathfrak{Re} \left ( \sum_{i=1}^{N_S} \frac{\partial_{\bm{\theta}}\psi^{*}_{\bm{\theta}}(\bm{\sigma^{(i)}})}{\psi^{*}_{\bm{\theta}}(\bm{\sigma^{(i)}})} \left ( E_{loc}(\bm{\sigma^{(i)}}) - E \right ) \right) \nonumber\\
&=
\frac{2}{N_S} \mathfrak{Re} \left ( \sum_{i=1}^{N_S} \partial_{{\bm{\theta}}} \log \psi^{*}_{\bm{\theta}}(\bm{\sigma^{(i)}}) \left ( E_{loc}(\bm{\sigma^{(i)}}) - E \right ) \right).
\label{eq:stochastic_gradient2}
\end{align}
This estimator has improved variance, stabilizes the convergence of the algorithm, and is unbiased~\cite{RNNWF_2020}. We note that in the limit where $E_{loc}(\bm{\sigma^{(i)}}) \approx E$ near convergence, the variance of the gradients $\partial_{{\bm{\theta}}} E$ goes to zero as opposed to the nonzero variance of the gradients in Eq.~\eqref{eq:stochastic_gradient1}.
\subsubsection*{Code walk-through}
We provide an implementation of VMC using a RNN wavefunction based on the TensorFlow library. We first define all the relevant parameters of the Rydberg Hamiltonian and the RNN training.
\begin{lstlisting}[language=Python,numbers=none]
# Hamiltonian parameters
Lx = 4 # linear size in x direction
Ly = 4 # linear size in y direction
N = Lx*Ly # total number of atoms
V = 7.0 # Van der Waals interaction
Omega = 1.0 # Rabi frequency
delta = 1.0 # detuning
# RNN-VMC parameters
lr = 0.001 # learning rate
nh = 32 # number of hidden units
ns = 500 # number of samples
epochs = 1000 # training iterations
seed = 1234 # seed of RNG
# initialize the RNN wavefunction
vmc = VariationalMonteCarlo(
Lx, Ly, V, Omega, delta,
nh, lr, epochs, seed)
\end{lstlisting}
The RNN parameters are the learning rate for the optimizer, the number of hidden units in the GRU cell, the number of samples used to approximate the energy (and its gradients) at each training iteration, and the total number of epochs. The VMC module is then initialized accordingly.
We can now perform the VMC simulation by training the RNN parameters. For a given number of epochs, we first sample the RNN distribution to generate \texttt{ns} samples. These can be used to evaluate the cost function of the optimization and its gradients, computed here using the AD functionalities from TensorFlow. The parameters are then updated according to the gradients.
\begin{lstlisting}[language=Python,numbers=none]
# training loop
for n in range(epochs):
# sample the RNN wavefunction
samples, _ = vmc.sample(ns)
# evaluate the loss function in AD mode
with tf.GradientTape() as tape:
logpsi = vmc.logpsi(samples)
eloc = vmc.localenergy(samples, logpsi)
Eo = tf.stop_gradient(tf.reduce_mean(eloc))
loss = tf.reduce_mean(2.0*tf.multiply(logpsi, tf.stop_gradient(eloc)) - 2.0*Eo*logpsi)
# compute the gradients
gradients = tape.gradient(loss, vmc.trainable_variables)
# update the parameters
vmc.optimizer.apply_gradients(zip(gradients, vmc.trainable_variables))
# get average energy
avg_E = np.mean(eloc.numpy())
\end{lstlisting}
We show in Fig.~\ref{Fig::vmc} the results of VMC simulations for a $8\times8$ array of Rydberg atoms. In Fig.~\ref{Fig::vmc}(a), we plot the average total energy at each training iteration for a RNN wavefunction with a GRU containing $n_h=25$ and $n_h=100$ hidden units. As expected, increasing the number of hidden units in the RNN leads to a improved accuracy in the training. In Fig.~\ref{Fig::vmc}(b) we plot the average energy after convergence as a function of the detuning using $n_h=100$, and compare with the values obtained from DMRG.
\begin{figure}[t]
\noindent \centering{}\includegraphics[width=\columnwidth]{fig_vmc}
\caption{Variational Monte Carlo simulation of the Rydberg Hamiltonian. ({\bf a}) Average energy during as a function of the training step near the phase transition at $\delta=1.3$. We show RNN wavefunctions with $n_h=25$ and $n_h=100$ hidden units, compared to the result obtained from DMRG. ({\bf b}) Average energy over the full phase diagram for a RNN wavefunction with $n_h=100$, in comparison with DMRG.}
\label{Fig::vmc}
\end{figure}
\section{Conclusions}
In this Article, we presented applications of machine learning algorithms based on neural networks to quantum many-body physics. We focused on numerical demonstrations and provided hands-on code tutorials based on open-source software~\cite{tensorflow,netket,itensor}, with the goal of facilitating learning for researcher new to the field, and accelerating the adoption of machine learning in quantum physics.
We showcased distinct machine learning paradigms, implemented with different neural-network architectures. As a testbed for the numerical experiments, we have chosen a system of interacting Rydberg atoms arranged in a two-dimensional square array. First, we demonstrated supervised learning of atomic occupation data, which was generated from ground states of the Rydberg Hamiltonian obtained using the density-matrix renormalization group. Using a convolutional neural network trained on labelled data, we showed how to learn the quantum phase transition between a disordered phase and the anti-ferromagnetically ordered phase.
We presented unsupervised learning of unlabelled data using the restricted Boltzmann machine. We showed that for pure quantum states with real and positive amplitude (e.g. the Rydberg ground states), this procedure is equivalent to quantum state tomography based on a neural-network representation of a quantum state. Using atomic occupation data, we trained Boltzmann machines to learn the ground states of the Rydberg Hamiltonian. We also described a simple extension of this approach to learn quantum states with a sign structure, and showed a demonstration in the context of simulation of chemistry with quantum computers, where we learned the ground state of a molecule using qubit measurements.
The final application we explored is the Monte Carlo optimization of a variational wavefunction to estimate the ground state of a many-body Hamiltonian. We parametrized a wavefunction using a recurrent neural-network and trained its parameters to lower both the expectation value and the variance of the Rydberg Hamiltonian.
It is becoming increasingly evident that machine learning is full of thrilling opportunities and conceptual advances with great potential to energize computational and experimental physics. As these notions continue to spread through the research landscape of strongly-correlated quantum matter and quantum information science, we hope this tutorial will provide a useful first step into the expanding domain of artificial intelligence for the study of quantum many-body systems.
\section*{Acknowledgements}
The numerical simulation were performed on the Simons Foundation Super-Computing Center, using the following software libraries: ITensor~\cite{itensor} (for the simulation of the Rydberg atoms), TensorFlow~\cite{tensorflow} (for supervised learning and variational Monte Carlo), and NetKet~\cite{netket} (for quantum tomography). The Flatiron Institute is supported by the Simons Foundation. JC acknowledges support from Natural Sciences and Engineering Research Council of Canada (NSERC), the Shared Hierarchical Academic Research Computing Network (SHARCNET), Compute Canada, Google Quantum Research Award, and the Canadian Institute for Advanced Research (CIFAR) AI chair program.
\bibliography{bibliography}
\end{document}
doc/Articles/mlquantum.tar
BIN
View File
Binary file not shown.
doc/Articles/rnn.pdf
BIN
View File
Binary file not shown.