arXiv:1708.06799v2 [] 29 Mar 2018

Divide-and-Conquer Checkpointing for Arbitrary Programs
with No User Annotation
Jeffrey Mark Siskinda∗ and Barak A. Pearlmutterb
a

School of Electrical and Computer Engineering, Purdue University, 465 Northwestern Avenue,
West Lafayette, IN, USA; b Dept of Computer Science, Maynooth University, Ireland

Classical reverse-mode automatic differentiation (AD) imposes only a small constant-factor
overhead in operation count over the original computation, but has storage requirements that
grow, in the worst case, in proportion to the time consumed by the original computation.
This storage blowup can be ameliorated by checkpointing, a process that reorders application
of classical reverse-mode AD over an execution interval to tradeoff space vs. time. Application of checkpointing in a divide-and-conquer fashion to strategically chosen nested execution
intervals can break classical reverse-mode AD into stages which can reduce the worst-case
growth in storage from linear to sublinear. Doing this has been fully automated only for computations of particularly simple form, with checkpoints spanning execution intervals resulting
from a limited set of program constructs. Here we show how the technique can be automated
for arbitrary computations. The essential innovation is to apply the technique at the level of
the language implementation itself, thus allowing checkpoints to span any execution interval.
Keywords: Reverse-Mode Automatic Differentiation; Binomial Checkpointing; Treeverse;
Programming Language Theory; Compiler Theory; Lambda Calculus
AMS Subject Classification: 68N20, 68N18, 65F50, 65D25, 46G05, 58C20

1.

Introduction

Reverse-mode automatic differentiation (AD) traverses the run-time dataflow graph of a
calculation in reverse order, in a so-called reverse sweep, so as to calculate a Jacobiantranspose-vector product of the Jacobian of the given original (or primal) calculation
[22]. Although the number of arithmetic operations involved in this process is only a
constant factor greater than that of the primal calculation, some values involved in the
primal dataflow graph must be saved for use in the reverse sweep, thus imposing considerable storage overhead. This is accomplished by replacing the primal computation
with a forward sweep that performs the primal computation while saving the requisite
values on a data structure known as the tape. A technique called checkpointing [25]
reorders portions of the forward and reverse sweeps to reduce the maximal length of the
requisite tape. Doing so, however, requires (re)computation of portions of the primal and
saving the requisite program state to support such as snapshots. Overall space savings
result when the space saved by reducing the maximal length of the requisite tape exceeds the space cost of storing the snapshots. Such space saving incurs a time cost in
(re)computation of portions of the primal. Different checkpointing strategies lead to a
∗

Corresponding author. Email: qobi@purdue.edu

space-time tradeoff.
We introduce some terminology that
(a)
(b)
will be useful in describing checkpointing.
An execution point is a point in time during the execution of a program. A program point is a location in the program
code. Since program fragments might be
invoked zero or more times during the execution of a program, each execution point
(c)
(d)
corresponds to exactly one program point
Figure 1. Checkpointing in reverse-mode AD. See text
but each program point may correspond for description.
to zero or more execution points. An execution interval is a time interval spanning two
execution points. A program interval is a fragment of code spanning two program points.
Program intervals are usually constrained so that they nest, i.e., they do not cross one
boundary of a syntactic program construct without crossing the other. Each program
interval may correspond to zero or more execution intervals, those execution intervals
whose endpoints result from the same invocation of the program interval. Each execution
interval corresponds to at most one program interval. An execution interval might not
correspond to a program interval because the endpoints might not result from the same
invocation of any program interval.
Figs. 1 and 2 illustrate the process of
performing reverse-mode AD with and
(a)
without checkpointing. Control flows from
top to bottom, and along the direction
of the arrow within each row. The symbols u, v, and p0 , . . . , p6 denote execution
(b)
(c)
points in the primal, u being the start of
the computation whose derivative is desired, v being the end of that computation, and each pi being an intermediate
execution point in that computation. Reverse mode involves various sweeps, whose
execution intervals are represented as hor(d)
izontal green, red, and blue lines. Green
lines denote (re)computation of the primal without taping. Red lines denote computation of the primal with taping, i.e.,
the forward sweep of reverse mode. Blue
lines denote computation of the Jacobiantranspose-vector product, i.e., the reverse
(e)
(f)
sweep of reverse mode. The vertical black
2. Divide-and-conquer checkpointing in reverselines denote collections of execution points Figure
mode AD. See text for description.
across the various sweeps that correspond to execution points in the primal, each particular execution point being the intersection of a horizontal line and a vertical line. In
portions of Figs. 1 and 2 other than Fig. 1(a) we refer to execution points for other sweeps
besides the primal in a given collection with the symbols u, v, and p0 , . . . , p6 when the
intent is clear. The vertical violet, gold, pink, and brown lines denote execution intervals
for the lifetimes of various saved values. Violet lines denote the lifetime of a value saved
on the tape during the forward sweep and used during the reverse sweep. The value is
saved at the execution point at the top of the violet line and used once at the execution
u

p0

p1

p2

p3

v

u

p0

p1

p2

p3

v

u

p0

p1

p2

p3

v

u

p0

p1

p2

p3

v

u

p0

v

u

p

v

u

p0

p1

p2

v

u

p0

p1

p2

v

2

u

u

p0

p0

p1

p1

p2

p3

p2

p4

p5

v

p6

v

point at the bottom of that line. Gold and pink lines denote the lifetime of a snapshot.1
The snapshot is saved at the execution point at the top of each gold or pink line and
used at various other execution points during its lifetime. Green lines emanating from a
gold or pink line indicate restarting a portion of the primal computation from a saved
snapshot.
Fig. 1(a) depicts the primal computation, y = f (x), which takes t time steps, with x
being a portion of the program state at execution point u and y being a portion of the
program state at execution point v computed from x. This is performed without taping
(green). Fig. 1(b) depicts classical reverse mode without checkpointing. An uninterrupted
forward sweep (red) is performed for the entire length of the primal, then an uninterrupted reverse sweep (blue) is performed for the entire length. Since the tape values
are consumed in reverse order from which they are saved, the requisite tape length is
O(t). Fig. 1(c) depicts a checkpoint introduced for the execution interval [p0 , p3 ). This
interrupts the forward sweep and delays a portion of that sweep until the reverse sweep.
Execution proceeds by a forward sweep (red) that tapes during the execution interval
[u, p0 ), a primal sweep (green) without taping during the execution interval [p0 , p3 ), a
taping forward sweep (red) during the execution interval [p3 , v ), a reverse sweep (blue)
during the execution interval [v, p3 ), a taping forward sweep (red) during the execution
interval [p0 , p3 ), a reverse sweep (blue) during the execution interval [p3 , p0 ), and then a
reverse sweep (blue) during the execution interval [p0 , u). The forward sweep for the execution interval [p0 , p3 ) is delayed until after the reverse sweep for the execution interval
[v, p3 ). As a result of this reordering, the tapes required for those sweeps are not simultaneously live. Thus the requisite tape length is the maximum of the two tape lengths, not
their sum. This savings comes at a cost. To allow such out-of-order execution, a snapshot
(gold) must be saved at p0 and the portion of the primal during the execution interval
[p0 , p3 ) must be computed twice, first without taping (green) then with (red).
A checkpoint can be introduced into a portion of the forward sweep that has been
delayed, as shown in Fig. 1(d). An additional checkpoint can be introduced for the execution interval [p1 , p2 ). This will delay a portion of the already delayed forward sweep
even further. As a result, the portions of the tape needed for the three execution intervals [p1 , p2 ), [p2 , p3 ), and [p3 , v ) are not simultaneously live, thus further reducing the
requisite tape length, but requiring more (re)computation of the primal (green). The
execution intervals for multiple checkpoints must either be disjoint or must nest; the
execution interval of one checkpoint cannot cross one endpoint of the execution interval
of another checkpoint without crossing the other endpoint.
Execution intervals for checkpoints can be specified in a variety of ways.
program interval
Execution intervals of specified program intervals constitute checkpoints.
subroutine call site
Execution intervals of specified subroutine call sites constitute checkpoints.
subroutine body
Execution intervals of specified subroutine bodies constitute checkpoints [25].
Nominally, these have the same power; with any one, one could achieve the effect of the
other two. Specifying a subroutine body could be accomplished by specifying all call sites
to that subroutine. Specifying some call sites but not others could be accomplished by
having two variants of the subroutine, one whose body is specified and one whose is not,
and calling the appropriate one at each call site. Specifying a program interval could be
1

The distinction between gold and pink lines, the meaning of brown lines, and the meaning of the black tick marks
on the left of the gold and pink lines will be explained in Section 3.10.

3

accomplished by extracting that interval as a subroutine.
Tapenade [12] allows the user to specify program intervals for checkpoints with
the c$ad checkpoint-start and c$ad checkpoint-end pragmas. Tapenade, by default,
checkpoints all subroutine calls [8]. This default can be overridden for named subroutines with the -nocheckpoint command-line option and for both named subroutines and
[u, v )
specific call sites with the c$ad nocheckpoint pragma.
Recursive application of checkpointing in a divide-and-conquer fashion, [u, p) [p, v)
i.e., “treeverse,” can divide the forward and reverse sweeps into stages run
Figure 3.
sequentially [9]. The key idea is that only one stage is live at a time, thus requiring
a shorter tape. However, the state of the primal computation at various intermediate
execution points needs to be saved as snapshots, in order to (re)run the requisite portion
of the primal to allow the forward and reverse sweeps for each stage to run in turn.
This process is illustrated in Fig. 2. Consider a root execution interval [u, v ) of the
derivative calculation. Without checkpointing, the forward and reverse sweeps span the
entire root execution interval, as shown in Fig. 2(a). One can divide the root execution
interval [u, v ) into two subintervals [u, p) and [p, v ) at the split point p and checkpoint
the first subinterval [u, v ). This divides the forward (red) and reverse (blue) sweeps into
two stages. These two stages are not simultaneously live. If the two subintervals are the
same length, this halves the storage needed for the tape at the expense of running the
primal computation for [u, p) twice, first without taping (green), then with taping (red).
This requires a single snapshot (gold) at u. This process can be viewed as constructing
a binary checkpoint tree (Fig. 3) whose nodes are labeled with execution intervals, the
intervals of the children of a node are adjacent, the interval of a node is the disjoint union
of the intervals of its children, and left children are checkpointed.
[u, v )
One can construct a left-branching binary checkpoint tree
[p , v )
[u, p )
over the same root execution interval [u, v ) with the split (a)
[u, p )
[p , p )
points p0 , p1 , and p2 (Fig. 4a). This can also be viewed as
[u, p ) [p , p )
constructing an n-ary checkpoint tree where all children but
[u, v )
the rightmost are checkpointed (Fig. 4b). This leads to nested
checkpoints for the execution intervals [u, p0 ), [u, p1 ), and (b) [u, p ) [p , p ) [p , p ) [p , v)
[u, p2 ) as shown in Fig. 2(c). Since the starting execution
Figure 4.
point u is the same for these intervals, a single snapshot (gold) with longer lifetime
suffices. These checkpoints divide the forward (red) and reverse (blue) sweeps into four
stages. This allows the storage needed for the tape to be reduced arbitrarily (i.e., the
red and blue segments can be made arbitrarily short), by rerunning successively shorter
prefixes of the primal computation (green), without taping, running only short segments
(red) with taping. This requires an O(t) increase in time for (re)computation of the
[u, v )
primal (green).
Alternatively, one can construct a right-branching binary
[ p0 , v )
checkpoint tree over the same root execution interval [u, v ) [u, p0)
with the same split points p0 , p1 , and p2 (Fig. 5). This also di[ p0 , p1 )
[ p1 , v )
vides the forward (red) and reverse (blue) sweeps into four
[ p1 , p2 ) [ p2 , v )
stages. With this, the requisite tape length (the maximal
Figure 5.
length of the red and blue segments) can be reduced arbitrarily while running the primal
(green) just once, by saving more snapshots (gold and pink), as shown in Fig. 2(d), This
requires an O(t) increase in space for storage of the live snapshots (gold and pink).
Thus we see that divide-and-conquer checkpointing can make the requisite tape arbitrarily small with either left- or right-branching binary checkpoint trees. This involves
a space-time tradeoff. The left-branching binary checkpoint trees require a single snapshot but an O(t) increase in time for (re)computation of the primal (green). The right2

2

1

0

0

4

1

0

0

1

2

1

1

2

2

branching binary checkpoint trees require an O(t) increase in space for storage of the
[u, v )
live snapshots (gold and pink) but (re)run the primal only once.
One can also construct a complete binary checkpoint tree
[ p1 , v )
[u, p1)
over the same root execution interval [u, v ) with the same split
points p0 , p1 , and p2 (Fig. 6). This constitutes application of [u, p0) [p0, p1) [p1, p2) [p2, v)
Figure 6.
the approach from Fig. 2(b) in a divide-and-conquer fashion
as shown in Fig. 2(e). This also divides the forward (red) and reverse (blue) sweeps into
four stages. One can continue this divide-and-conquer process further, with more split
points, more snapshots, and more but shorter stages, as shown in Fig. 2(f). This leads
to an O(log t) increase in space for storage of the live snapshots (gold and pink) and an
O(log t) increase in time for (re)computation of the primal (green). Variations of this
technique can tradeoff between different improvements in space and/or time complexity,
leading to overhead in a variety of sublinear asymptotic complexity classes in one or
both. In order to apply this technique, we must be able to construct a checkpoint tree
of the desired shape with appropriate split points. This in turn requires the ability to
interrupt the primal computation at appropriate execution points, save the interrupted
execution state as a capsule, and restart the computation from the capsules, sometimes
repeatedly.2
[u, v )
Any given divide-and-conquer decomposition of the same
root execution interval with the same split points can be (a)
[u, p )
[p , v )
[u, p ) [p , p ) [p , p ) [p , v )
viewed as either a binary checkpoint tree or an n-ary check[u, v )
point tree. Thus Fig. 2(e) can be viewed as either Fig. 7(a)
or Fig. 7(b). Similarly, Fig. 2(f) can be viewed as either (b)
[u, p )
[p , p )
[p , v )
Fig. 7(c) or Fig. 7(d). Thus we distinguish between two al[p , p ) [p , v )
gorithms to perform divide-and-conquer checkpointing.
1

0

1

0

0

1

0

1

2

1

2

1

2

1

2

[u, v )

binary An algorithm that constructs a binary checkpoint
tree.
treeverse The algorithm from [9, Figs. 2 and 3] that constructs an n-ary checkpoint tree.

[u, p3 )

(c)
[u, p1 )

[u, p0 )

[ p0 , p1 )

[ p3 , v )
[ p1 , p3 )

[ p1 , p2 )

[ p2 , p3 )

[ p3 , p5 )

[ p3 , p4 )

[ p4 , p5 )

[ p5 , v )

[ p5 , p6 )

[ p6 , v )

[u, v )

(d)

[u, p0 )

[ p0 , p1 )

[ p1 , p3 )

[ p1 , p2 )

[ p2 , p3 )
[ p3 , p4 )

[ p3 , v )

There is, however, a simple correspondence between associated binary and n-ary checkpoint trees. The n-ary checkFigure 7.
point tree is derived from the binary checkpoint tree by coalescing each maximal sequence
of left branches into a single node. Thus we will see, in Section 5, that these two algorithms exhibit the same properties.
Note that (divide-and-conquer) checkpointing does not incur any space or time overhead in the forward or reverse sweeps themselves (i.e., the number of violet lines and the
total length of red and blue lines). Any space overhead results from the snapshots (gold
and pink) and any time overhead results from (re)computation of the primal (green).
Several design choices arise in the application of divide-and-conquer checkpointing in
addition to the choice of binary vs. n-ary checkpoint trees.
[ p4 , p5 )

[ p5 , v )

[ p5 , p6 )

[ p6 , v )

● What root execution interval(s) should be subject to divide-and-conquer checkpointing?
● Which execution points are candidate split points? The divide-and-conquer process of
constructing the checkpoint tree will select actual split points from these candidates.
● What is the shape or depth of the checkpoint tree, i.e., what is the termination criterion
for the divide-and-conquer process?
2

The correspondence between capsules and snapshots will be discussed in Section 3.10.

5

Since the leaf nodes of the checkpoint tree correspond to stages, the termination criterion and the number of evaluation steps in the stage at each leaf node (the length of a
pair of red and blue lines) are mutually constrained. The number of live snapshots at
a leaf (how many gold and pink lines are crossed by a horizontal line drawn leftward
from that stage, the pair of red and blue lines, to the root) depends on the depth of
the leaf and its position in the checkpoint tree. Different checkpoint trees, with different
shapes resulting from different termination criteria and split points, can lead to a different maximal number of live snapshots, resulting in different storage requirements. The
amount of (re)computation of the primal (the total length of the green lines) can also
depend on the shape of the checkpoint tree, thus different checkpoint trees, with different
shapes resulting from different termination criteria and split points, can lead to different compute-time requirements. Thus different strategies for specifying the termination
criterion and the split points can influence the space-time tradeoff.
We make a distinction between several different approaches to selecting root execution
intervals subject to divide-and-conquer checkpointing.
loop Execution intervals resulting from invocations of specified DO loops are subject to
divide-and-conquer checkpointing.
entire derivative calculation The execution interval for an entire specified derivative
calculation is subject to divide-and-conquer checkpointing.
We further make a distinction between several different approaches to selecting candidate
split points.
iteration boundary Iteration boundaries of the DO loop specified as the root execution
interval are taken as candidate split points.
arbitrary Any execution point inside the root execution interval can be taken as a
candidate split point.
We further make a distinction between several different approaches to specifying the
termination criterion and deciding which candidate split points to select as actual split
points.
bisection Split points are selected so as to divide the computation dominated by a
node in half as one progresses successively from right to left among children [9,
equation (12)]. One can employ a variety of termination criteria, including that
from [9, p. 46]. If the termination criterion is such that the total number of leaves
is a power of two, one obtains a complete binary checkpoint tree. A termination
criterion that bounds the number of evaluation steps in a leaf limits the size of
the tape and achieves logarithmic overhead in both asymptotic space and time
complexity compared with the primal.
binomial Split points are selected using the criterion from [9, equation (16)]. The termination criterion from [9, p. 46] is usually adopted to achieve the desired properties
discussed in [9]. Different termination criteria can be selected to control spacetime tradeoffs.
fixed space overhead One can bound the size of the tape and the number of
snapshots to obtain sublinear but superlogarithmic overhead in asymptotic
time complexity compared with the primal.
fixed time overhead One can bound the size of the tape and the
(re)computation of the primal to obtain sublinear but superlogarithmic overhead in asymptotic space complexity compared with the primal.
logarithmic space and time overhead One can bound the size of the tape

6

and obtain logarithmic overhead in both asymptotic space and time complexity compared with the primal. The constant factor is less than that of
bisection checkpointing.
We elaborate on the strategies for selecting actual split points from candidate split points
and the associated termination criteria in Section 5.
Divide-and-conquer checkpointing has only been provided to date in AD systems in
special cases. For example, Tapenade allows the user to select invocations of a specified
DO loop as the root execution interval for divide-and-conquer checkpointing with the
c$ad binomial-ckp pragma, taking iteration boundaries of that loop as candidate split
points. Tapenade employs binomial selection of split points and a fixed space overhead termination criterion. Note, however, that Tapenade only guarantees this fixed
space overhead property for DO loop bodies that take constant time. Similarly adolc [10] contains a nested taping mechanism for time-integration processes [17] that also
performs divide-and-conquer checkpointing. This only applies to code formulated as a
time-integration process.
Here, we present a framework for applying divide-and-conquer checkpointing to arbitrary code with no special annotation or refactoring required. An entire specified derivative calculation is taken as the root execution interval, rather than invocations of a
specified DO loop. Arbitrary execution points are taken as candidate split points, rather
than iteration boundaries. As discussed below in Section 5, both binary and n-ary (treeverse) checkpoint trees are supported. Furthermore, as discussed below in Section 5,
both bisection and binomial checkpointing are supported. Additionally, all of the above
termination criteria are supported: fixed space overhead, fixed time overhead, and logarithmic space and time overhead. Any combination of the above checkpoint-tree generation algorithms, split-point selection methods, and termination criteria are supported.
In order to apply this framework, we must be able to interrupt the primal computation
at appropriate execution points, save the interrupted execution state as a capsule, and
restart the computation from the capsules, sometimes repeatedly. This is accomplished
by building divide-and-conquer checkpointing on top of a general-purpose mechanism
for interrupting and resuming computation. This mechanism is similar to engines [13]
and is orthogonal to AD. We present several implementations of our framework which
we call checkpointVLAD. In Section 6, we compare the space and time usage of our
framework with that of Tapenade on an example.
Note that one cannot generally achieve the space and time guarantees of divide-andconquer checkpointing with program-interval, subroutine-call-site, or subroutine-body
checkpointing unless the call tree has the same shape as the requisite checkpoint tree.
Furthermore, one cannot generally achieve the space and time guarantees of divideand-conquer checkpointing for DO loops by specifying the loop body as a programinterval checkpoint, as that would lead to a right-branching checkpoint tree and behavior
analogous to Fig. 2(d). Moreover, if one allows split points at arbitrary execution points,
the resulting checkpoint execution intervals may not correspond to program intervals.
Some form of divide-and-conquer checkpointing is necessary. One may wish to take the
gradient of a long-running computation, even if it has low asymptotic time complexity.
The length of the tape required by reverse mode without divide-and-conquer checkpointing increases with increasing run time. Modern computers can execute several billion
floating point operations per second, even without GPUs and multiple cores, which only
exacerbate the problem. If each such operation required storage of a single eight-byte
double precision number, modern terabyte RAM sizes would fill up after a few seconds
of computation. Thus without some form of divide-and-conquer checkpointing, it would

7

not be possible to efficiently take the gradient of a computation that takes more than a
few seconds.
Machine learning methods in general, and deep learning methods in particular, require
taking gradients of long-running high-dimension computations, particularly when training deep neural networks in general or recurrent neural networks over long time series.
Thus variants of divide-and-conquer checkpointing have been rediscovered and deployed
by the machine learning community in this context [6, 11]. These implementations are far
from automatic, and depend on compile-time analysis of the static primal flow graphs.
The general strategy of divide-and-conquer checkpointing, the n-ary treeverse algorithm, the bisection and binomial strategies for selecting split points, and the termination
criteria that provide fixed space overhead, fixed time overhead, and logarithmic space and
time overhead were all presented in [9]. Furthermore, Tapenade has implemented divideand-conquer checkpointing with the n-ary treeverse algorithm, the binomial strategy for
selecting split points, and the termination criterion that provides fixed space overhead,
but only for root execution intervals corresponding to invocations of specified DO loops
that meet certain criteria with split points restricted to iteration boundaries of those
loops. To our knowledge, the binary checkpoint-tree algorithm presented here and the
framework for allowing it to achieve all of the same guarantees as the n-ary treeverse algorithm is new. However, our central novel contribution here is providing a framework for
supporting either the binary checkpoint-tree algorithm or the n-ary treeverse algorithm,
either bisection or binomial split point selection, and any of the termination criteria of
fixed space overhead, fixed time overhead, or logarithmic space and time overhead in a
way that supports taking the entire derivative calculation as the root execution interval and taking arbitrary execution points as candidate split points, by integrating the
framework into the language implementation.
Some earlier work [14, 15, 23] prophetically presaged the work here. This work seems
to have received far less exposure and attention than deserved. Perhaps because the ideas
therein were so advanced and intricate that it was difficult to communicate those ideas
clearly. Moreover, the authors report difficulties in getting their implementations to be
fully functional. Our work here formulates the requisite ideas and mechanisms carefully
and precisely, using methods from the programming-language community, like formulation of divide-and-conquer checkpointing of a function as divide-and-conquer application
of reverse mode to two functions whose composition is the original function, formulation
of the requisite decomposition as a precise and abstract interruption and resumption
interface, formulation of semantics precisely through specification of evaluators, use of
CPS evaluators to specify an implementation of the interruption and resumption interface, and systematic derivation of a compiler from that evaluator via CPS conversion, to
allow complete, correct, comprehensible, and fully general implementation.

2.

The Limitations of Divide-and-Conquer Checkpointing with Split Points
at Fixed Syntactic Program Points like Loop Iteration Boundaries

Consider the example in Fig. 8. This example, y = f (x; l, φ), while contrived, is a simple
caricature of a situation that arises commonly in practice: modeling a physical system
with an adaptive grid. An initial state vector x ∶ Rn is repeatedly transformed by a
state update process Rn → Rn and, upon termination, an aggregate property y of the
final state is computed by a function Rn → R. We wish to compute the gradient of that
property y relative to the initial state x. Here, the state update process first rotates the
value pairs at adjacent odd-even coordinates of the state x by an angle θ and then rotates
8

f u n c t i o n ilog2 ( l )
ilog2 = dlog ( real (l , 8))/ dlog (2.0 d0 )
end
s u b r o u t i n e rotate ( theta , x1 , x2 , x1p , x2p )
double p r e c i s i o n theta , x1 , x2 , x1p , x2p , c , s
c = cos ( theta )
s = sin ( theta )
x1p = c * x1 - s * x2
x2p = s * x1 + c * x2
end
s u b r o u t i n e rot1 ( theta , n , x )
double p r e c i s i o n theta , x , x1p , x2p
dimension x(n)
do k = 1 , n -1 , 2
call rotate ( theta , x ( k ) , x ( k +1) , x1p , x2p )
x ( k ) = x1p
x ( k +1) = x2p
end do
end

s u b r o u t i n e f (n , x , l , phi , y )
double p r e c i s i o n phi , x , y , x1
d i m e n s i o n x ( n ) , x1 (1000)
do k = 1 , n
x1 ( k ) = x ( k )
end do
c $ad binomial - ckp l +1 40 1
do i = 1 , l
m = 2**( ilog2 ( l ) +
ilog2 (1+ int ( mod ( 1 0 1 3 . 0 d0 *3.0 d0 ** phi * real (i , 8) ,
real (l , 8)))))
+
do j = 1 , m
call magsqr (n , x1 , y )
y = sqrt ( y )
call rot1 (1.2 d0 *y , n , x1 )
call rot2 (1.4 d0 *y , n , x1 )
end do
end do
call magsqr (n , x1 , y )
y = y /2.0 d0
end

s u b r o u t i n e rot2 ( theta , n , x )
double p r e c i s i o n theta , x , x1p , x2p
dimension x(n)
do k = 2 , n -2 , 2
call rotate ( theta , x ( k ) , x ( k +1) , x1p , x2p )
x ( k ) = x1p
x ( k +1) = x2p
end do
end

p r o g r a m main
double p r e c i s i o n phi , x , xb , y , yb
d i m e n s i o n x (1000) , xb (1000)
read * , n
read * , l
read * , phi
do k = 1 , n
x ( k ) = n +1 - k
xb ( k ) = 0.0 d0
end do
yb = 1.0 d0
call f (n , x , l , phi , y )
print * , y
call f_b (n , x , xb , l , phi , y , yb )
do k = 1 , n
print * , xb ( k )
end do
end

s u b r o u t i n e magsqr (n , x , y )
double p r e c i s i o n x , y
dimension x(n)
y = 0.0 d0
do k = 1 , n
y = y + x ( k )* x ( k )
end do
end

Figure 8. Fortran example. This example is rendered in checkpointVLAD in Fig. 28. Space and time overhead
of two variants of this example when run under Tapenade are presented in Fig. 30. The pragma used for the
variant with divide-and-conquer checkpointing of the outer DO loop is shown. This pragma is removed for the
variant with no checkpointing.

those at adjacent even-odd coordinates. The rotation θ is taken to be proportional to
the magnitude of x. The adaptive grid manifests in two nested update loops. The outer
loop has duration l, specified as an input hyperparameter. The duration m of the inner
loop varies wildly as some function of another
input hyperparameter φ and the outer
⌊lg(l)⌋−⌊lg(1+(1013⌊3x ⌋i mod l))⌋
loop index i, perhaps 2
, that is small on most iterations of
the outer loop but O(l) on a few iterations. If the split points were limited to iteration
boundaries of the outer loop, as would be common in existing implementations, the
increase in space or time requirements would grow larger than sublinearly. The issue is
that for the desired sublinear growth properties to hold, it must be possible to select
arbitrary execution points as split points. In other words, the granularity of the divideand-conquer decomposition must be primitive atomic computations, not loop iterations.
The distribution of run time across the program is not modularly reflected in the static
syntactic structure of the source code, in this case the loop structure. Often the user is
unaware of, or even unconcerned with, the micro-level structure of atomic computations,
and does not wish to break the modularity of the source code to expose it. Yet the user
may still wish to reap the sublinear space or time overhead benefits of divide-and-conquer
checkpointing. Moreover, the relative duration of different paths through a program
may vary from loop iteration to loop iteration in a fashion that is data dependent,
as shown by the above example, and not even statically determinable. We will now
proceed to discuss an implementation strategy for divide-and-conquer checkpointing that
does not constrain split points to loop iteration boundaries or other syntactic program
constructs and does not constrain checkpoints to program intervals or other syntactic
program constructs. Instead, it can take any arbitrary execution point as a split point
and introduce checkpoints at any resulting execution interval.

9

3.

Technical Details of our Method

Implementing divide-and-conquer checkpointing requires the capacity to
(1)
(2)
(3)
(4)

measure the length of the primal computation,
interrupt the primal computation at a portion of the measured length,
save the state of the interrupted computation as a capsule, and
resume an interrupted computation from a capsule.

For our purposes, the second and third operations are always performed together and can
be fused into a single operation. These can be difficult to implement efficiently as library
routines in an existing language implementation (see Section 7.1). Thus we design a new
language implementation, checkpointVLAD, with efficient support for these low-level
operations.
3.1

Core Language

checkpointVLAD adds builtin AD operators to a functional pre-AD core language.
In the actual implementation, this core language is provided with a Scheme-like surface
syntax.3 But nothing turns on this; the core language can be exposed with any surface
syntax. For expository purposes, we present the core language here in a simple, more
traditional, math-like notation.
checkpointVLAD employs the same functional core language as our earlier vlad
system [21]. Support for AD in general, and divide-and-conquer checkpointing in particular, is simplified in a functional programming language (see Section 7.2). Except for
this simplification, which can be eliminated with well-known techniques (e.g., monads
[26] and uniqueness types [1]) for supporting mutation in functional languages, nothing
turns on our choice of core language. We intend our core language as a simple expository
vehicle for the ideas presented here; they could be implemented in other core languages
(see Section 7.1).
Our core language contains the following constructs:
e ∶∶= c ∣ x ∣ λx.e ∣ e1 e2 ∣ if e1 then e2 else e3 ∣ ◇e ∣ e1 ● e2

(1)

Here, e denotes expressions, c denotes constants, x denotes variables, e1 e2 denotes function application, ◇ denotes builtin unary operators, and ● denotes builtin binary operators. For expository simplicity, the discussion of the core language here omits many
vagaries such as support for recursion and functions of multiple arguments; the actual
implementation supports these using standard mechanisms that are well known within
the programming-language community (e.g., tupling or Currying).
3.2

Direct-Style Evaluator for the Core Language

We start by formulating a simple evaluator for this core language (Fig. 9) and extend
it to perform AD and ultimately divide-and-conquer checkpointing. This evaluator is
written in what is known in the programming-language community as direct style, where
functions (in this case E , denoting ‘eval’ and A, denoting ‘apply’) take inputs as functioncall arguments and yield outputs as function-call return values [18]. While this evaluator
can be viewed as an interpreter, it is intended more as a description of the evaluation
3

The surface syntax employed differs slightly from Scheme in ways that are irrelevant to the issue at hand.

10

A ⟨(λx.e), ρ⟩ v = E ρ[x ↦ v ] e
E ρc=c
E ρx=ρx
E ρ (λx.e) = ⟨(λx.e), ρ⟩
E ρ (e1 e2 ) = A (E ρ e1 ) (E ρ e2 )
E ρ (if e1 then e2 else e3 ) = if (E ρ e1 ) then (E ρ e2 ) else (E ρ e3 )
E ρ (◇e) = ◇(E ρ e)
E ρ (e1 ● e2 ) = (E ρ e1 ) ● (E ρ e2 )

(2a)
(2b)
(2c)
(2d)
(2e)
(2f)
(2g)
(2h)

Figure 9. Direct-style evaluator for the core checkpointVLAD language.

mechanism; this mechanism could be the underlying hardware as exposed via a compiler.
Indeed, as described below in Section 3.14, we have written three implementations, one
an interpreter, one a hybrid compiler/interpreter, and one a compiler.
With any evaluator, one distinguishes between two language evaluation strata: the
target, the language being implemented and the process of evaluating programs in that
language, and the host, the language in which the evaluator is written and the process of
evaluating the evaluator itself. In our case, the target is checkpointVLAD, while the
host varies among our three implementations; for the first two it is Scheme while for
the third it is the underlying hardware, achieved by compilation to machine code via c.
In the evaluator in Fig. 9, ρ denotes an environment, a mapping from variables to
their values, ρ0 denotes the empty environment that does not map any variables, ρ x
denotes looking up the variable x in the environment ρ to obtain its value, ρ[x ↦ v ]
denotes augmenting an environment ρ to map the variable x to the value v, and E ρ e
denotes evaluating the expression e in the context of the environment ρ. There is a clause
for E in Fig. 9, (2b–h), for each construct in (1). Clause (2b) says that one evaluates a
constant by returning that constant. Clause (2c) says that one evaluates a variable by
returning its value in the environment. The notation ⟨e, ρ⟩ denotes a closure, a lambda
expression e together with an environment ρ containing values for the free variables in e.
Clause (2d) says that one evaluates a lambda expression by returning a closure with the
environment in the context that the lambda expression was evaluated in. Clause (2e)
says that one evaluates an application by evaluating the callee expression to obtain a
closure, evaluating the argument expression to obtain a value, and then applying the
closure to the value with A. A, as described in (2a), evaluates the body of the lambda
expression in the callee closure in the environment of that closure augmented with the
formal parameter of that lambda expression bound to the argument value. The remaining
clauses are all analogous to clause (2h), which says that one evaluates an expression e1 ● e2
in the target by evaluating e1 and e2 to obtain values and then applying ● in the host to
these values.
3.3

Adding AD Operators to the Core Language

Unlike many AD systems implemented as libraries, we provide support for AD by augmenting the core language to include builtin AD operators for both forward and reverse
mode [21]. This allows seamless integration of AD into the language in a completely
general fashion with no unimplemented or erroneous corner cases. In particular, it allows
Ð
→
←
Ð
nesting [19]. In checkpointVLAD, we adopt slight variants of the J and J operators
11

previously incorporated into vlad. (Nothing turns on this. The variants adopted here
are simpler, better suit our expository purposes, and allow us to focus on the issue at
hand.) In checkpointVLAD, these operators have the following signatures.

Ð
→
J ∶ f x x́ ↦ (y, ý)

←
Ð
J ∶ f x ỳ ↦ (y, x̀)

(3)

We use the notation x́ and x̀ to denote tangent or cotangent values associated with the
primal value x respectively, and the notation (x, y ) to denote a pair of values. Since
in checkpointVLAD, functions can take multiple arguments but only return a single
result, which can be an aggregate like a pair, the AD operators take the primal and the
associated (co)tangent as distinct arguments but return the primal and the associated
(co)tangent as a pair of values.
Ð
→
The J operator provides the portal to forward mode and calls a function f on a
←
Ð
primal x with a tangent x́ to yield a primal y and a tangent ý. The J operator provides
the portal to reverse mode and calls a function f on a primal x with a cotangent ỳ to
yield a primal y and a cotangent x̀.4 Unlike the original vlad, here, we restrict ourselves
to the case where (co)tangents are ground data values, i.e., reals and (arbitrary) data
structures containing reals and other scalar values, but not functions (i.e., closures).
Nothing turns on this; it allows us to focus on the issue at hand.
The implementations of vlad and checkpointVLAD are disjoint and use completely
different technology. The Stalin∇ [21] implementation of vlad is based on source-code
transformation, conceptually applied reflectively at run time but migrated to compile
time through partial evaluation. The implementation of checkpointVLAD uses something more akin to operator overloading. Again, nothing turns on this; this simplification
is for expository purposes and allows us to focus on the issue at hand (see Section 7.1).
In checkpointVLAD, AD is performed by overloading the arithmetic operations in
the host, in a fashion similar to fadbad++ [5]. The actual method used is that employed
by r6rs-ad5 and DiffSharp6 . The key difference is that fadbad++ uses c++ templates
to encode a hierarchy of distinct forward-mode types (e.g., F<double>, F<F<double> >,
. . . ), distinct reverse-mode types (e.g., B<double>, B<B<double> >, . . . ), and mixtures
thereof (e.g., F<B<double> >, B<F<double> >, . . . ) while here, we use a dynamic, runtime approach where numeric values are tagged with the nesting level [19, 24]. Template
instantiation at compile-time specializes code to different nesting levels. The dynamic
approach allows a single interpreter (host), formulated around unspecialized code, to
interpret different target programs with different nesting levels.
3.4

Augmenting the Direct-Style Evaluator to Support the AD Operators

We add AD into the target language as new constructs.

Ð
→
←
Ð
e ∶∶= J e1 e2 e3 ∣ J e1 e2 e3

(4)

We implement this functionality by augmenting the direct-style evaluator with new
Ð
→
←
Ð
clauses for E (Fig. 10), clause (2k) for J and clause (2l) for J . These clauses are all
Ð
→
←
Ð
analogous to clause (2h), formulated around J and J operators in the host. These are
Ð
→
←
Ð
In the implementation, J and J are named j* and *j respectively.
https://github.com/qobi/R6RS-AD and https://engineering.purdue.edu/~ qobi/stalingrad-examples2009/
6 http://diffsharp.github.io/DiffSharp/
4

5

12

Ð
→
J v1 v2 v́3 = let (v4 ⊳ v́5 ) = (A v1 (v2 ⊳ v́3 )) in (v4 , v́5 )
←
Ð
J v1 v2 v̀3 = let (v4 ⊲ v̀5 ) = ((A v1 v2 ) ⊲ v̀3 ) in (v4 , v̀5 )
Ð
→
Ð
→
E ρ ( J e1 e2 e3 ) = J (E ρ e1 ) (E ρ e2 ) (E ρ e3 )
←
Ð
←
Ð
E ρ ( J e1 e2 e3 ) = J (E ρ e1 ) (E ρ e2 ) (E ρ e3 )

(2i)
(2j)
(2k)
(2l)

Figure 10. Additions to the direct-style evaluator for checkpointVLAD to support AD.

Ð
→
←
Ð
defined in (2i, 2j). The J and J operators in the host behave like A except that they
level shift to perform AD. Just like (A f x) applies a target function f (closure) to a
Ð
→
target value x, ( J f x x́) performs forward mode by applying a target function f (clo←
Ð
sure) to a target primal value x and a target tangent value x́, while ( J f x ý ) performs
reverse mode by applying a target function f (closure) to a target primal value x and a
target cotangent value ỳ.
Ð
→
As described in (2i), J operates by recursively walking v2 , a data structure containing primals, in tandem with v́3 , a data structure containing tangents, to yield a single
data structure where each numeric leaf value is a dual number, a numeric primal value
associated with a numeric tangent value. This recursive walk is denoted as v2 ⊳ v́3 . A is
then used to apply the function (closure) v1 to the data structure produced by v2 ⊳ v́3 .
Since the input argument is level shifted and contains dual numbers instead of ordinary
reals, the underlying arithmetic operators invoked during the application perform forward mode by dispatching on the tags at run time. The call to A yields a result data
structure where each numeric leaf value is a dual number. This is then recursively walked
to separate out two data structures, one, v4 , containing the numeric primal result values,
and the other, v́5 , containing the numeric tangent result values, which are returned as a
pair (v4 , v́5 ) This recursive walk is denoted as let (v4 ⊳ v́5 ) = . . . in . . ..
←
Ð
As described in (2j), J operates by recursively walking v2 , a data structure containing
primals, to replace each numeric value with a tape node. A is then used to apply the
function (closure) v1 to this modified v2 . Since the input argument is level shifted and
contains tape nodes instead of ordinary reals, the underlying arithmetic operators invoked
during the application perform the forward sweep of reverse mode by dispatching on the
tags at run time. The call to A yields a result data structure where each numeric leaf
value is a tape node. A recursive walk is performed on this result data structure, in
tandem with a data structure v̀3 of associated cotangent values, to initiate the reverse
sweep of reverse mode. This combined operation is denoted as ((A v1 v2 ) ⊲ v̀3 ). The
result of the forward sweep is then recursively walked to replace each tape node with its
numeric primal value and the input value is recursively walked to replace each tape node
with the cotangent computed by the reverse sweep. These are returned as a pair (v4 , v̀5 ).
This combined operation is denoted as let (v4 ⊲ v̀5 ) = . . . in . . ..
3.5

An Operator to Perform Divide-and-Conquer Checkpointing in
Reverse-Mode AD
✓

We introduce a new AD operator J to perform divide-and-conquer checkpointing.7 (For
expository simplicity, we focus for now on binary bisection checkpointing. In Section 5
7

✓

In the implementation, J is named checkpoint-*j.

13

✓

To compute (y, x̀) = J f x ỳ:
←
Ð
base case (f x fast): (y, x̀) = J f x ỳ
inductive case:

h○g =f
z=gx

✓

(y, z̀) = J h z ỳ
✓

(z, x̀) = J g x z̀

(step 0)
(step 1)
(step 2)
(step 3)
(step 4)

Figure 11. Algorithm for binary checkpointing.

✓

below, we provide alternate implementations of J that perform treeverse and/or binomial
✓

checkpointing.) The crucial aspect of the design is that the signature (and semantics) of J
←
Ð
is identical to J ; they are completely interchangeable, differing only in the space/time
complexity tradeoffs. This means that code need not be modified to switch back and forth
between ordinary reverse mode and various forms of divide-and-conquer checkpointing,
←
Ð
✓
save interchanging calls to J and J .
✓

Conceptually, the behavior of J is shown in Fig. 11. In this inductive definition, a
function f is split into the composition of two functions g and h in step (1), the capsule z
is computed by applying g to the input x in step (2), and the cotangent is computed by
✓

recursively applying J to h and g in steps 3 and 4. This divide-and-conquer behavior is
terminated in a base case, when the function f is small, at which point the cotangent
←
Ð
is computed with J , in step (0). If step (1) splits a function f into two functions g
and h that take the same number of evaluation steps, and we terminate the recursion
when f takes a bounded number of steps, the recursive divide-and-conquer process yields
logarithmic asymptotic space/time overhead complexity.
The central difficulty in implementing the above is performing step (1), namely splitting
a function f into two functions g and h, such that f = h ○ g, ideally where we can specify
the split point, the number of evaluation steps through f where g transitions into h. A
sophisticated user can manually rewrite a subprogram f into two subprograms g and h.
A sufficiently powerful compiler or source transformation tool might also be able to do
so, with access to nonlocal program text. But an overloading system, with access only
to local information, would not be able to.
3.6

General-Purpose Interruption and Resumption Mechanism

We solve this problem by providing an interface to a general-purpose interruption and
resumption mechanism that is orthogonal to AD (Fig. 12). This interface allows (a) determining the number of evaluation steps of a computation, (b) interrupting a computation
after a specified number of steps, usually half the number of steps determined by the
mechanism in (a), and (c) resuming an interrupted computation to completion. A variety of implementation strategies for this interface are possible. We present two in detail
below, in Sections 3.8 and 3.12, and briefly discuss another in Section 7.1.
Irrespective of how one implements the general-purpose interruption and resumption
✓

interface, one can use it to implement the binary bisection variant of J in the host, as
14

primops f x ↦ l
interrupt f x l ↦ z
resume z ↦ y

Return the number l of evaluation steps needed to compute y = f (x).

Run the first l steps of the computation of f (x) and return
a capsule z.
If z = (interrupt f x l), return y = f (x).

Figure 12. General-purpose interruption and resumption interface.

✓

To compute (y, x̀) = J f x ỳ:
←
Ð
base case:
(y, x̀) = J f x ỳ
inductive case:

l = primops f x
z = interrupt f x
✓

(step 0)

(y, z̀) = J (λz.resume z ) z ỳ
✓

(step 1)

⌊ 2l ⌋

(z, x̀) = J (λx.interrupt f x

(step 2)
(step 3)

⌊ 2l ⌋)

x z̀

(step 4)

Figure 13. Binary bisection checkpointing via the general-purpose interruption and resumption interface. Step (1)
need only be performed once at the beginning of the recursion, with steps (2) and (4) taking l at the next
recursion level to be ⌊ 2l ⌋ and step (3) taking l at the next recursion level to be ⌈ 2l ⌉. As discussed in the text, this
implementation is not quite correct, because (λz.resume z) in step (3) and (λx.interrupt f x ⌊ 2l ⌋) in step (4)
are host closures but need to be target closures. A proper implementation is given in Fig. 18.

shown in Fig. 13. The function f is split into the composition of two functions g and h
by taking g as (λx.interrupt f x l), where l is half the number of steps determined by
(primops f x), and h as (λz.resume z ).
3.7

Continuation-Passing-Style Evaluator

One way of implementing the general-purpose interruption and resumption interface is
to convert the evaluator from direct style to what is known in the programming-language
community as continuation-passing style (CPS) [18], where functions (in this case E , A,
Ð
→
←
Ð
J , and J in the host) take an additional continuation input k and instead of yielding
outputs via function-call return, do so by calling the continuation with said output as
arguments (Figs. 14 and 15). In CPS, functions never return: they just call their continuation. With tail-call merging, this corresponds to a computed go to and does not
incur stack growth. This crucially allows an interruption to actually return a capsule
containing the saved state of the evaluator, including its continuation, allowing the evaluation to be resumed by calling the evaluator with this saved state. This ‘level shift’ of
return to calling a continuation, allowing an actual return to constitute interruption, is
analogous to the way backtracking is classically implemented in Prolog, with success
implemented as calling a continuation and failure implemented as actual return. In our
case, we further instrument the evaluator to thread two values as inputs and outputs:
the count n of the number of evaluation steps, which is incremented at each call to E ,
and the limit l of the number of steps, after which an interrupt is triggered.
Fig. 14 contains the portion of the CPS evaluator for the core language corresponding
to Fig. 9, while Fig. 15 contains the portion of the CPS evaluator for the AD constructs
corresponding to Fig. 10. Except for (5b), the equations in Figs. 9 and 10 are in one-toone correspondence to those in Figs. 14 and 15, in order. Clauses (5c–e) are analogous
15

A k n l ⟨(λx.e), ρ⟩ v = E k n l ρ[x ↦ v ] e
E k l l ρ e = ⟦k, ⟨(λ .e), ρ⟩⟧
E k n l ρ c = k (n + 1) l c
E k n l ρ x = k (n + 1) l (ρ x)
E k n l ρ (λx.e) = k (n + 1) l ⟨(λx.e), ρ⟩
E k n l ρ (e1 e2 ) = E (λn l v1 .
(E (λn l v2 .
(A k n l v1 v2 ))
n l ρ e2 ))
(n + 1) l ρ e1

E k n l ρ (if e1 then e2 else e3 ) = E (λn l v1 .
(if v1
then (E k n l ρ e2 )
else (E k n l ρ e3 )))
(n + 1) l ρ e1
E k n l ρ (◇e) = E (λn l v.
(k n l (◇v )))
(n + 1) l ρ e

E k n l ρ (e1 ● e2 ) = E (λn l v1 .
(E (λn l v2 .
(k n l (v1 ● v2 )))
n l ρ e2 ))
(n + 1) l ρ e1

(5a)
(5b)
(5c)
(5d)
(5e)
(5f)

(5g)

(5h)

(5i)

Figure 14. CPS evaluator for the core checkpointVLAD language.

to the corresponding clauses (2b–d) except that they call the continuation k with the
result, instead of returning that result. The remaining clauses for E in the CPS evaluator
are all variants of

E (λn l v1 .(E (λn l v2 .(k n l . . .)) n l ρ e2 )) (n + 1) l ρ e1

(6)

for one-, two-, or three-argument constructs. This evaluates the first argument e1 and
calls the continuation (λn l v1 . . . .) with its value v1 . This continuation then evaluates
the second argument e2 and calls the continuation (λn l v2 . . . .) with its value v2 . This
continuation computes something, denoted by . . ., and calls the continuation k with the
resulting value.
The CPS evaluator threads a step count n and a step limit l through the evaluation
process. Each clause of E increments the step count exactly once to provide a coherent
fine-grained measurement of the execution time. Clause (5b) of E implements interruption. When the step count reaches the step limit, a capsule containing the saved state of
the evaluator, denoted ⟦k, f ⟧, is returned. Here, f is a closure ⟨(λ .e), ρ⟩ containing the
environment ρ and the expression e at the time of interruption. This closure takes an
argument that is not used. The step count n must equal the step limit l at the time of
interruption. As will be discussed below, in Section 3.8, neither the step count nor the
step limit need to be saved in the capsule, as the computation is always resumed with
different step count and limit values.

16

Ð
→
J v1 v2 v́3 = A (λn l (v4 ⊳ v́5 ).
(v4, v́5 ))
0 ∞ v1 (v2 ⊳ v́3 )
←
Ð
J v1 v2 v̀3 = A (λn l v.
let (v4 ⊲ v̀5 ) = v ⊲ v̀3
in (v4 , v̀5 ))
0 ∞ v1 v2
Ð
→
E k n l ρ ( J e1 e2 e3 ) = E (λn l v1 .
(E (λn l v2 .
Ð
→
(E (λn l v3 .(k n l ( J v1 v2 v3 )))
n l ρ e3 ))
n l ρ e2 ))
(n + 1) l ρ e1
←
Ð
E k n l ρ ( J e1 e2 e3 ) = E (λn l v1 .
(E (λn l v2 .
←
Ð
(E (λn l v3 .(k n l ( J v1 v2 v3 )))
n l ρ e3 ))
n l ρ e2 ))
(n + 1) l ρ e1

(5j)

(5k)

(5l)

(5m)

Figure 15. Additions to the CPS evaluator for checkpointVLAD to support AD.

Several things about this CPS evaluator are of note. First, all builtin unary and binary
operators are assumed to take unit time. This follows from the fact that all clauses for E ,
as typified by (6), increment the step count by one. Second, the builtin unary and binary
operators in the host are implemented in direct style and are not passed a continuation.
This means that clauses (5h, 5i), as typified by (6), must call the continuation k on the
result of the unary and binary operators. Third, like all builtin operators, invocations
Ð
→
←
Ð
of the J and J operators, including the application of v1 , are assumed to take unit
time. This follows from the fact that clauses (5l, 5m), again as typified by (6), increment
Ð
→
←
Ð
the step count by one. Fourth, like all builtin operators, J and J in the host, in (5j,
5k), are implemented in direct style and are not passed a continuation. This means that
Ð
→
clauses (5l, 5m), as typified by (6), must call the continuation k on the result of J
←
Ð
Ð
→
←
Ð
and J . Finally, since J and J receive target functions (closures) for v1 , they must
apply these to their arguments with A. Since A is written in CPS in the CPS evaluator,
these calls to A in (5j, 5k) must be provided with a continuation k, a step count n, and
a step limit l as arguments. The continuation argument simply returns the result. The
step count, however, is restarted at zero, and the step limit is set to ∞. This means that
Ð
→
←
Ð
invocations of J and J are atomic and cannot be interrupted internally. We discuss this
further below in Section 7.3.
3.8

Implementing the General-Purpose Interruption and Resumption
Interface with the CPS Evaluator

With this CPS evaluator, it is possible to implement the general-purpose interruption and
resumption interface (Fig. 16). The implementation of primops (7a) calls the evaluator
17

primops f x = A (λn l v.n) 0 ∞ f x
interrupt f x l = A (λn l v.v ) 0 l f x
resume ⟦k, f ⟧ = A k 0 ∞ f 

(7a)
(7b)
(7c)

Figure 16. Implementation of the general-purpose interruption and resumption interface using the CPS evaluator.

with no step limit and simply counts the number of steps to completion. The implementation of interrupt (7b) calls the evaluator with a step limit that must be smaller
than that needed to complete so an interrupt is forced and the capsule ⟦k, ⟨(λ .e), ρ⟩⟧ is
returned. The implementation of resume (7c) calls the evaluator with arguments from
the saved capsule. Since the closure in the capsule does not use its argument, an arbitrary
value  is passed as that argument.
←
Ð
Ð
→
Note that the calls to A in J (5j), J (5k), primops (7a), interrupt (7b), and
resume (7c) are the only portals into the CPS evaluator. The only additional call to A
is in the evaluator itself, clause (5f) of E . All of the portals restart the step count at
zero. Except for the call in interrupt (7b), none of the portals call the evaluator with
a step limit. In particular, resume (7c) does not provide a step limit; other mechanisms
detailed below provide for interrupting a resumed capsule.
This implementation of the general-purpose interruption and resumption interface can✓

not be used to fully implement J in the host as depicted in Fig. 13. The reason is that
←
Ð
the calls to J in the base case, step (0), and interrupt in step (2), must take a target
←
Ð
function (closure) for f , because that is what is invoked by the calls to A in J (5k)
✓

and interrupt (7b). As written in Fig. 13, the recursive calls to J , namely steps (3)
and (4), pass (λz.resume z ) and (λx.interrupt f x ⌊ 2l ⌋) for f . There are two problems with this. First, these are host closures produced by host lambda expressions, not
target closures. Second, these call the host functions resume and interrupt that are
not available in the target. Thus it is not possible to formulate these as target closures
without additional machinery.
Examination of Fig. 13 reveals that the general-purpose interruption and resump✓

tion interface is invoked four times in the implementation of J . primops is invoked in
step (1), interrupt is invoked in steps (2) and (4), and resume is invoked in step (3).
Of these, primops is invoked only in the host, resume is invoked only in the target,
and interrupt is invoked in both the host and the target. Thus we need to expose
interrupt and resume to the target. We do not need to expose primops to the target; the implementation in Fig. 13 only uses it in the host. For interrupt, the call in
step (2) can use the host implementation (7b) in Fig. 16 but the call in step (4) must
use a new variant exposed to the target. For resume, the call in step (3) must also use
a new variant exposed to the target. The host implementation (7c) in Fig. 16 is never
used since resume is never invoked in the host.
We expose interrupt and resume to the target by adding them to the target language
as new constructs.
e ∶∶= interrupt e1 e2 e3 ∣ resume e

(8)

We implement this functionality by augmenting the CPS evaluator with new clauses
for E (Fig. 17), clause (6p) for interrupt and clause (6q) for resume. We discuss the
implementation of these below. But we first address several other issues.
18

I f l = ⟨(λx.(interrupt f x l)), ρ0[f ↦ f ][l ↦ l]⟩
R = ⟨(λz.(resume z )), ρ0 ⟩
E k n l ρ (interrupt e1 e2 e3 ) = E (λn l v1 .
(E (λn l v2 .
(E (λn l v3 .
if l = ∞
then (A k 0 v3 v1 v2 )
else let ⟦k, f ⟧ = (A k 0 l v1 v2 )
in ⟦k, (I f (v3 − l))⟧)
n l ρ e3 ))
n l ρ e2 ))
(n + 1) l ρ e1
E k n l ρ (resume e) = E (λn l ⟦k′ , f ⟧.(A k′ 0 l f )) (n + 1) l ρ e

(6n)
(6o)
(6p)

(6q)

Figure 17. Additions to the CPS evaluator for checkpointVLAD to expose the general-purpose interruption and
resumption interface to the target.

✓

To compute (y, x̀) = J f x ỳ:
←
Ð
base case:
(y, x̀) = J f x ỳ
inductive case:

l = primops f x
z = interrupt f x
✓

(y, z̀) = J R z ỳ
✓

(z, x̀) = J (I f

(step 0)

⌊ 2l ⌋

⌊ 2l ⌋)

(step 1)
(step 2)
(step 3)

x z̀

(step 4)

Figure 18. Binary bisection checkpointing in the CPS evaluator. This is a proper implementation of the algorithm
in Fig. 13 where the host closure (λz.resume z) in step (3) is replaced with the target closure R and the host
closure (λx.interrupt f x ⌊ 2l ⌋) in step (4) is replaced with the target closure (I f ⌊ 2l ⌋).

With appropriate implementations of interrupt and resume expressions in the target language, one can create target closures for the expressions (λz.resume z ) and
✓

(λx.interrupt f x ⌊ 2l ⌋), and use these to formulate a proper implementation of J in
the host. We formulate a target closure to correspond to (λz.resume z ) and denote this
as R. The definition is given in (6o) in Fig. 17. Note that since (λz.resume z ) does not
contain any free variables, the closure created by R is constructed from the empty environment ρ0 . Thus there is a single constant R. We similarly formulate a target closure to
correspond to (λx.interrupt f x l) and denote this as I . The definition is given in (6n)
in Fig. 17. Here, however, (λx.interrupt f x l) contains two free variables: f and l.
Thus the closure created by I contains a nonempty environment with values for these
two variables. To provide these values, I is formulated as a function that takes these
values as arguments.
✓

With (I f l) and R, it is now possible to reformulate the definition of J in the
host from Fig. 13, replacing the host closure (λz.resume z ) in step (3) with the target
closure R and the host closure (λx.interrupt f x ⌊ 2l ⌋) in step (4) with the target
✓

closure (I f ⌊ 2l ⌋). This new, proper, definition of J in the host is given in Fig. 18.
✓

In this proper implementation of J in the host, the interrupt and resume operations
19

✓

need to be able to nest, even without nesting of calls to J in the target. The recursive
✓

calls to J in the inductive case of Fig. 18 imply that it must be possible to interrupt
a resumed capsule. This happens when passing R for f in step (3) and then passing
(I f . . .) for f in step (4), i.e., the left branch of a right branch in the checkpoint tree.
The resulting function f = (I R . . .) will interrupt when applied to some capsule. It also
happens when passing (I f . . .) for f twice in succession in step (4), i.e., the left branch
of a left branch in the checkpoint tree. The resulting function f = (I (I f . . .) . . .) will
interrupt and the capsule produced will interrupt when resumed.
Consider all the ways that evaluations of interrupt and resume expressions can nest.
User code will never contain interrupt and resume expressions; they are created only
✓

by invocations of I and R. R is only invoked by step (3) of J in Fig. 18. I is invoked two
✓

ways: step (4) of J in Fig. 18 and a way that we have not yet encountered, evaluation
of nested interrupt expressions in the else branch of clause (6p) in Fig. 17.

✓

✓

Consider all the ways that evaluations of I and R can be invoked in J in Fig. 18.

J is invoked with some user code for f , i.e., code that does not contain interrupt
✓

and resume expressions. The inductive cases for J create a binary checkpoint tree of
invocations. The leaf nodes of this binary checkpoint tree correspond to the base case in
←
Ð
step (0) where the host J is invoked. At internal nodes, the host interrupt is invoked
←
Ð
in step (2). The target closure values that can be passed to the host J and interrupt
are constructed from f , I , and R in steps (3) and (4). What is the space of all possible
constructed target closures? The constructed target closures invoked along the left spine
of the binary checkpoint tree look like the following

(I (I . . . (I (I f l0 ) l1 ) . . . li−1 ) li )

(9)

with zero or more nested calls to I . In this case li < li−1 < ⋯ < l1 < l0 , because the recursive
✓

calls to J in step (4) of Fig. 18 always reduce l. The constructed target closures invoked
in any other node in the binary checkpoint tree look like the following

(I (I . . . (I (I R l0 ) l1 ) . . . li−1 ) li )

(10)

with zero or more nested calls to I . In this case, again, li < li−1 < ⋯ < l1 < l0 , for the same
←
Ð
reason. These are the possible target closures f passed to J in step (0) or interrupt
✓

in step (2) of J in Fig. 18. (We assume that the call to primops in step (1) is hoisted
out of the recursion.)
A string of calls to I as in (9) will result in a nested closure structure whose invocation
will lead to nested invocations of interrupt expressions.

⟨(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨. . .
(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨(λx.(interrupt f x l)),
ρ0 [f ↦ f ][l ↦ l0 ]⟩][l ↦ l1 ] . . .⟩][l ↦ li−1 ]⟩][l ↦ li ]⟩
20

(11)

A string of calls to I as in (10) will also result in a nested closure structure whose
invocation will lead to nested invocations of interrupt expressions.

⟨(λx.(interrupt f x l)),

ρ0 [f ↦ ⟨(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨. . .
(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨(λz.(resume z)), ρ0 ⟩]
[l ↦ l0 ]⟩][l ↦ l1 ] . . .⟩][l ↦ li−1 ]⟩][l ↦ li ]⟩

(12)

In both of these, li < li−1 < ⋯ < l1 < l0 , so the outermost interrupt expression will
interrupt first. Since the CPS evaluator only maintains a single step limit, li will be that
step limit during the execution of the innermost content of these nested closures, namely f
in (11) and ⟨(λz.(resume z )), ρ0 ⟩ in (12). None of the other intervening interrupt
expressions will enforce their step limits during this execution. Thus we need to arrange
for the capsule created when the step limit li is reached during the execution of f or
⟨(λz.(resume z )), ρ0 ⟩ to itself interrupt with the remaining step limits li−1 , . . . , l1 , l0 .
This is done by rewrapping the closure in a capsule with interrupt expressions. The
interruption of f or ⟨(λz.(resume z )), ρ0 ⟩ will produce a capsule that looks like the
following

⟦k, f ⟧

(13)

where the closure f contains only user code, i.e., no interrupt or resume expressions.
The f in (13) is wrapped with calls to I to reintroduce the step limits li−1 , . . . , l1 , l0 .

⟦k, (I . . . (I (I f l0 ) l1 ) . . . li−1 )⟧

(14)

This will yield a capsule that looks like the following

⟦k, ⟨(λx.(interrupt f x l)),
ρ0 [ f ↦ ⟨ . . .
(λx.(interrupt f x l)),
ρ0 [f ↦ ⟨(λx.(interrupt f x l)),
ρ0 [f ↦ f ][l ↦ l0 ]⟩][l ↦ l1 ] . . .⟩][l ↦ li−1 ]⟩⟧

(15)

which will interrupt upon resumption. Each such interruption will peel off one interrupt
expression. Note that since the closure f in a capsule (13) contains only user code, it
will not contain a resume expression. Further, since the wrapping process (15) only
introduces interrupt expressions via calls to I (14), and never introduces resume expressions, the closures in capsules, whether wrapped or not, will never contain resume
expressions.
When there is no contextual step limit, i.e., when l = ∞, the interrupt expression
must introduce v3 , the step limit specified as the argument to the interrupt expression,
as the step limit. This is handled by the then branch of clause (6p) in Fig. 17. When
there is a contextual step limit, i.e., when l =/ ∞, the interrupt expression must wrap
the returned capsule. This wrapping is handled by the else branch of clause (6p) in
Fig. 17. Since capsule resumption restarts the step count at zero, the wrapping that
handles nested step limits is relativized to this restart by the v3 − l in the else branch in
clause (6p).
21

Capsule resumption happens in one place, the call to A in clause (6q) in Fig. 17 for
a resume expression. Except for the contextual step limit l, this is the same as the call
to A in the implementation of resume in (7c) in Fig. 16. Said resumption is performed by
applying the capsule closure f , a target closure, to , since the lambda expression in the
capsule closure ignores its argument. This call to A is passed the capsule continuation k′
as its continuation. Unlike the implementation of resume in (7c), the step limit l is that
which is in effect for the execution of the resume expression. This is to allow capsule
resumption to itself interrupt. Because capsules are resumed with a step count of zero
and the step limit at the time of resumption, the step count and limit at the time of the
interruption need not be saved in the capsule.
As a result of this, all interrupt expressions will appear in one of two places. The first
is a preamble (11) or (12) wrapped around either a user function f by (9) or a resume
←
Ð
expression in R by (10), respectively. This will always be invoked either by J in the base
case, step (0), or by interrupt in step (2), of Fig. 18. The second is a preamble (15)
wrapped around the closure of a capsule by the else branch in clause (6p) of Fig. 17,
i.e., (14). This will always be invoked during capsule resumption, i.e., clause (6q) of
Fig. 17. We assume that the step limits are such that an interruption never occurs during
either of these preambles. This is enforced by ensuring that the termination criterion that
←
Ð
triggers the base case, step (0), of Fig. 18 is sufficiently long so that the calls to A in J in
step (0) and interrupt in step (2) won’t interrupt before completion of the preamble.
There is one further requirement to allow the CPS evaluator to support divide-and←
Ð
conquer checkpointing. The base case use of J in step (0) of Fig. 18 needs to be able to
produce cotangents z̀ of capsules z in step (3) and consume them in step (4). A capsule
⟦k, f ⟧ is the saved state of the evaluator. The value f is a target closure ⟨(λx.e), ρ⟩
←
Ð
which contains an environment with saved state. This state is visible to J . But the
continuation k is a host continuation, which is opaque. Any evaluator variables that
←
Ð
it closes over are not visible to J . Thus the implementation of host continuations in
the CPS evaluator must employ a mechanism to expose them. When we replace the
CPS evaluator with a direct-style evaluator applied to CPS-converted target code, in
Sections 3.11 and 3.12 below, this will no-longer be necessary since continuations will be
←
Ð
represented as target closures which are visible to J .
3.9

Augmenting the CPS Evaluator to Support Divide-and-Conquer
Checkpointing
✓

We can now add the J operator to the target language as a new construct.
✓

e ∶∶= J e1 e2 e3

(16)

We implement this functionality by augmenting the CPS evaluator with a new clause (6r)
for E (Fig. 19).
With this addition, target programs can perform divide-and-conquer checkpointing
←
Ð
✓
✓
simply by calling J instead of J . Note that it is not possible to add the J operator to
the direct-style evaluator because the implementation of binary bisection checkpointing
is built on the general-purpose interruption and resumption interface which is, in turn,
built on the CPS evaluator. We remove this limitation below in Section 3.12. Also note
that since the implementation of binary bisection checkpointing is built on the general22

✓

E k n l ρ (J e1 e2 e3 ) = E (λn l v1 .
(E (λn l v2 .
(E (λn l v3 .

(6r)

✓

(k n l (J v1 v2 v3 )))
n l ρ e3 ))
n l ρ e2 ))
(n + 1) l ρ e1

Figure 19. Addition to the CPS evaluator for checkpointVLAD to support divide-and-conquer checkpointing.

purpose interruption and resumption interface which is, in turn, built on an evaluator,
it is only available for programs that are evaluated, i.e., for programs in the target, but
not for programs in the host. We remove this limitation below as well, in Section 3.13.
3.10

Some Intuition

The algorithm in Fig. 18 corresponds to Fig. 2(b). The start of the computation of f
in Fig. 18 corresponds to u in Fig. 2(b). The computation state at u is x in Fig. 18.
Collectively, the combination of f and x in Fig. 18 comprises a snapshot, the gold line in
Fig. 2(b). The end of the computation of f in Fig. 18 corresponds to v in Fig. 2(b). The
computation state at v is y in Fig. 18. Step (1) computes ⌊ 2l ⌋ which corresponds to the
split point p in Fig. 2(b). Step (2) corresponds to the green line in Fig. 2(b), i.e., running
the primal without taping from the snapshot f and x at u until the split point p which
is ⌊ 2l ⌋. The capsule z in Fig. 18 corresponds to the computation state at p in Fig. 2(b).
Brown and pink lines in Fig. 2 denote capsules. If step (3) would incur the base case,
step (0), in the recursive call, it would correspond to the right stage (pair of red and
blue lines) in Fig. 2(b). If step (4) would incur the base case, step (0), in the recursive
call, it would correspond to the left stage (pair of red and blue lines) in Fig. 2(b). Note
that f and x is used both in steps (2) and (4). Referring to this as a snapshot is meant
to convey that the information must be saved across the execution of step (3). And it
must be possible to apply f to x twice, once in step (2) and once in step (4). In some
implementations, such a snapshot involves saving mutable state that must be restored.
In our formulation in a functional framework (Section 7.2), we need not explicitly save
and restore state; we simply apply a function twice. Nonetheless, the storage required
for the snapshot is implicit in the extended lifetime of the values f and x which extends
✓

from the entry into J , over step (3), until step (4).
✓

Note that recursive calls to J in step (4) extend the lifetime of a snapshot. These are
denoted as the black tick marks on the left of the gold and pink lines. In the treeverse
algorithm from [9, Figs. 2 and 3], the lifetime of one snapshot ends at a tick mark by
a call to retrieve in one recursive call to treeverse in the while loop of the parent
and the lifetime of a new snapshot begins by a call to snapshot in the next recursive
call to treeverse in the while loop of the parent. But since the state retrieved and
then immediately saved again as a new snapshot is the same, these adjacent snapshot
execution intervals can conceptually be merged.
✓

Also note that recursive calls to J in step (3) pass R and a capsule z as the f and x
of the recursive call. Thus capsules from one level of the recursion become snapshots at
the next level, for all but the base case step (0). Pink lines in Fig. 2 denote values that
23

2x∣k7 ↝ k x
2(λx.e)∣k7 ↝ k (λk′ x.2e∣k′ 7)
2(e1 e2 )∣k7 ↝ 2e1 ∣(λx1 .2e2 ∣(λx2 .(x1 k x2 ))7)7
e0 ↝ 2e0 ∣(λx.x)7

(17a)
(17b)
(17c)
(17d)

Figure 20. CPS conversion for the untyped lambda calculus.

are capsules at one level but snapshots at lower levels. Some, but not all, capsules are
snapshots. Some, but not all, snapshots are capsules. Gold lines in Fig. 2 denote snapshots
that are not capsules. Brown lines in Fig. 2 denote capsules that are not snapshots. Pink
z = (interrupt f x ⌊ ⌋)
lines in Fig. 2 denote values that are both snapshots and capsules.
It is now easy to see that the recursive call tree of the algorithm in Fig. 18 is isomorphic to a binary checkpoint tree. (z, x̀) = (✓J (I f ⌊ ⌋) x z̀) (y, z̀) = (✓J R z ỳ)
Figure 21.
The binary checkpoint tree in Fig. 3 corresponds to the call
tree in Fig. 21 produced by the algorithm in Fig. 18. This depicts just one level of the
recursion. If one unrolls this call tree, one obtains a binary checkpoint tree.
l
2

l
2

3.11

CPS Conversion

So far, we have formulated divide-and-conquer checkpointing via a CPS evaluator. This
can be—and has been—used to construct an interpreter. A compiler can be—and has
been—constructed by generating target code in CPS that is instrumented with step
counting, step limits, and limit checks that lead to interrupts. Code in direct style can be
automatically converted to CPS using a program transformation known in the programming language community as CPS conversion. Many existing compilers, such as sml/nj
for sml, perform CPS conversion as part of the compilation process [4].
We illustrate CPS conversion for the untyped lambda calculus (Fig. 20).
e ∶∶= x ∣ λx.e ∣ e1 e2

(18)

The notation 2e∣k7 denotes the transformation of the expression e to CPS so that it calls
the continuation k with the result. There is a clause for 2e∣k7 in Fig. 20, (17a–c), for
each construct in (18). Clause (17a) says that one converts a variable x by calling the
continuation k with the value of that variable. Clause (17b) says that one converts a
lambda expression (λx.e) by adding a continuation variable k′ to the lambda binder,
converting the body relative to that variable, and then calling the continuation k with
that lambda expression. Clause (17c) says that one converts an application (e1 e2 ) by
converting e1 with a continuation that receives the value x1 of e1 , then converts e2 with
a continuation that receives the value x2 of e2 , and then calls x1 with the continuation k
and x2 . Clause (17d) says that the top level expression e0 can be converted with the
identity function as the continuation.
This technique can be extended to thread a step count n and a step limit l through
the computation along with the continuation k, and to arrange for the step count to be
incremented appropriately. Further, this technique can be applied to the entire target
language (Fig. 22). Clauses (20a–j) correspond one-to-one to the checkpointVLAD
constructs in (1), (4), and (16). Since CPS conversion is only applied once at the beginning
of compilation, to the user program, and the user program does not contain interrupt
and resume expressions, since these only appear internally in the target closures created
by I and R, CPS conversion need not handle these constructs. Finally, ⟪e⟫k,n,l denotes
24

a limit check that interrupts and returns a capsule when the step count n reaches the
step limit l. The implementation of this limit check is given in (20k). Each of the clauses
(20a–j) is wrapped in a limit check.
3.12

Augmenting the Direct-Style Evaluator to Support CPS-Converted
Code and Divide-and-Conquer Checkpointing

The direct-style evaluator must be modified in several ways to support CPS-converted
code and divide-and-conquer checkpointing (Fig. 23). First, CPS conversion introduced
lambda expressions with multiple arguments and their corresponding applications. Continuations have three arguments and converted lambda expressions have four. Thus we
add several new constructs into the target language to replace the single argument lambda
expressions and applications from (1).
e ∶∶= λ3 n l x.e ∣ λ4 k n l x.e ∣ e1 e2 e3 e4 ∣ e1 e2 e3 e4 e5

(19)

Second, we need to modify E to support these new constructs. We replace clause (2a)
with clauses (21a, 21b) to update A and clauses (2d, 2e) with clauses (21c–f) to update E .
Third, we need to add support for interrupt and resume expressions, as is done with
clauses (21g, 21h). These are direct-style variants of clauses (6p, 6q) from the CPS evaluator and are needed to add support for the general-purpose interruption and resumption
interface to the direct-style evaluator when evaluating CPS code. Note that the calls
to A from (6p, 6q) are modified to use the converted form A4 of A (21b) in (21g, 21h).
Similarly, the calls to continuations from (6p, 6q) are modified to use the continuation
form A3 of A (21a) in (21g, 21h). Fourth, the calls to A4 must be modified in the host
Ð
→
←
Ð
implementations of the AD operators J and J , as is done with (21i, 21j). Note that
unlike the corresponding (2i, 2j), the calls to A4 here take target closures instead of host
closures. Fifth, the general-purpose interruption and resumption interface, (7a, 7b, 6n,
6o), must be migrated from the CPS evaluator to the direct-style evaluator as (21k–n).
In doing so, the calls to A4 in primops and interrupt are changed to use (21b), the
host continuations are modified to be target continuations in (21k, 21l), and the lambda
expressions in (21m, 21n) are CPS converted.
3.13

Compiling Direct-Style Code to C

One can compile target checkpointVLAD code, after CPS conversion, to c (Figs. 24
and 25). Modern implementations of c, like GCC, together with modern memory management technology, like the Boehm-Demers-Weiser garbage collector, allow the compilation process to be a straightforward mapping of each construct to a small fragment of c
code. In particular, garbage collection, GC malloc, eases the implementation of closures
and statement expressions, ({...}), together with nested functions, ease the implementation of lambda expressions. Furthermore, the flow analysis, inlining, and tail-call merging performed by GCC generates reasonably efficient code. In Figs. 24 and 25, S denotes
such a mapping from checkpointVLAD expressions e to c code fragments. Instead
of environments ρ, S takes π, a mapping from variables to indices in environment, the
run-time environment data structure. Here, π x denotes the index of x, πi denotes the
variable for index i, φ e denotes a mapping for the free variables in e, and N denotes
a mapping from a checkpointVLAD operator to the name of the c function that implements that operator. This, together with a library containing the typedef for thing,
25

2c∣k, n, l7 ↝ ⟪k (n + 1) l c⟫k,n,l
2x∣k, n, l7 ↝ ⟪k (n + 1) l x⟫k,n,l
2(λx.e)∣k, n, l7 ↝ ⟪k (n + 1) l (λ4 k n l x.2e∣k, n, l7)⟫k,n,l
2(e1 e2 )∣k, n, l7 ↝ ⟪2e1 ∣(λ3 n l x1 .
2e2 ∣(λ3 n l x2 .
(x1 k n l x2 )),
n, l7),
(n + 1), l7⟫k,n,l

2(if e1 then e2 else e3 )∣k, n, l7 ↝ ⟪2e1 ∣(λ3 n l x1 .
(if x1
then 2e2 ∣k, n, l7
else 2e3 ∣k, n, l7)),
(n + 1), l7⟫k,n,l
2(◇e)∣k, n, l7 ↝ ⟪2e∣(λ3 n l x.
(k n l (◇x))),
(n + 1), l7⟫k,n,l

(20e)

(20f)

2(e1 ● e2 )∣k, n, l7 ↝ ⟪2e1 ∣(λ3 n l x1 .
2e2 ∣(λ3 n l x2 .
(k n l (x1 ● x2 ))),
n, l7),
(n + 1), l7⟫k,n,l
Ð
→
2( J e1 e2 e3 )∣k, n, l7 ↝ ⟪2e1 ∣(λ3 n l x1 .
2e2 ∣(λ3 n l x2 .
2e3 ∣(λ3 n l x3 .
Ð
→
(k n l ( J x1 x2 x3 ))),
n, l7),
n, l7),
(n + 1), l7⟫k,n,l
←
Ð
2( J e1 e2 e3 )∣k, n, l7 ↝ ⟪2e1 ∣(λ3 n l x1 .
2e2 ∣(λ3 n l x2 .
2e3 ∣(λ3 n l x3 .
←
Ð
(k n l ( J x1 x2 x3 ))),
n, l7),
n, l7),
(n + 1), l7⟫k,n,l
✓

2(J e1 e2 e3 )∣k, n, l7 ↝ ⟪2e1 ∣(λ3 n l x1 .
2e2 ∣(λ3 n l x2 .
2e3 ∣(λ3 n l x3 .

(20a)
(20b)
(20c)
(20d)

(20g)

(20h)

(20i)

(20j)
✓

(k n l (J x1 x2 x3 ))),
n, l7),

n, l7),
(n + 1), l7⟫k,n,l
⟪e⟫k,n,l ↝ if n = l then ⟦k, λ4 k n l .e⟧ else e

(20k)

Figure 22. CPS conversion for the checkpointVLAD language that threads step counts and limits.

26

A3 ⟨(λ3 n l x.e), ρ⟩ n′ l′ v = E ρ[n ↦ n′ ][l ↦ l′ ][x ↦ v ] e
A4 ⟨(λ4 k n l x.e), ρ⟩ k′ n′ l′ v = E ρ[k ↦ k′ ][n ↦ n′ ][l ↦ l′ ][x ↦ v ] e
E ρ (λ3 n l x.e) = ⟨(λ3n l x.e), ρ⟩
E ρ (λ4 k n l x.e) = ⟨(λ4k n l x.e), ρ⟩
E ρ (e1 e2 e3 e4 ) = A3 (E ρ e1 ) (E ρ e2 ) (E ρ e3 ) (E ρ e4 )
E ρ (e1 e2 e3 e4 e5 ) = A4 (E ρ e1 ) (E ρ e2 ) (E ρ e3 ) (E ρ e4 ) (E ρ e5 )
E ρ (interrupt e1 e2 e3 ) = let v1 = (E ρ e1 )
v2 = (E ρ e2 )
v3 = (E ρ e3 )
k = ρ ‘k’
l = ρ ‘l’
in if l = ∞
then (A4 v1 k 0 v3 v2 )
else let ⟦k, f ⟧ = (A4 v1 k 0 l v2 )
in ⟦k, (I f (v3 − l))⟧

(21a)
(21b)
(21c)
(21d)
(21e)
(21f)
(21g)

E ρ (resume e) = let ⟦k′ , f ⟧ = (E ρ e)
(21h)
l = ρ ‘l’
in (A4 f k′ 0 l )
Ð
→
J v1 v2 v́3 = let (v4 ⊳ v́5 ) = (A4 v1 ⟨(λ3 n l v.v ), ρ0 ⟩ 0 ∞ (v2 ⊳ v́3 )) (21i)
in (v4 , v́5 )
←
Ð
J v1 v2 v̀3 = let (v4 ⊲ v̀5 ) = ((A4 v1 ⟨(λ3n l v.v ), ρ0 ⟩ 0 ∞ v2 ) ⊲ v̀3 ) (21j)
in (v4 , v̀5 )

primops f x = A4 f ⟨(λ3 n l v.n), ρ0 ⟩ 0 ∞ x
interrupt f l n = A4 f ⟨(λ3 n l v.v ), ρ0 ⟩ 0 l x
I f l = ⟨(λ4k n l x.(interrupt f x l)), ρ0[f ↦ f ][l ↦ l]⟩
R = ⟨(λ4k n l z.(resume z )), ρ0⟩

(21k)
(21l)
(21m)
(21n)

Figure 23. Extensions to the direct-style evaluator and the implementation of the general-purpose interruption
and resumption interface to support divide-and-conquer checkpointing on target code that has been converted to
CPS.

the enum for tag, definitions for null constant, true constant, false constant, cons,
as closure, set closure, continuation apply, converted apply, is false, and all
of the functions named by N (essentially a translation of r6rs-ad, the general-purpose
interruption and resumption interface from Fig. 23, and the implementation of binary
bisection checkpointing from Fig. 18 into c), allows arbitrary checkpointVLAD code
to be compiled to machine code, via c, with complete support for AD, including forward
mode, reverse mode, and binary bisection checkpointing.
3.14

Implementations

We have written three complete implementations of checkpointVLAD.8 All three accept exactly the same source language in its entirety and are able to run the example discussed in Section 6. The first implementation is an interpreter based on the CPS evaluator
(Figs. 14, 15, 17, and 19), where the evaluator, the operator overloading implementation
of AD, the general-purpose interruption and resumption mechanism (Fig. 16), and the
8 The

code for all three implementations will be released on github upon acceptance of this manuscript.

27

S π () = null constant
S π true = true constant
S π false = false constant
S π (c1 , c2 ) = cons((S π c1 ), (S π c1 ))
S πn=n
S π ‘k’ = continuation
S π ‘n’ = count
S π ‘l’ = limit
S π ‘x’ = argument
S π x = as closure(target)->environment[π x]
S π (λ3 n l x.e) = ({
thing function(thing target,
thing count,
thing limit,
thing argument) {
return (S (φ e) e);
}
thing lambda = (thing)GC malloc(sizeof(struct {
enum tag tag;
struct {
thing (*function)();
unsigned n;
thing environment[∣φ e∣];
}))
}
set closure(lambda);
as closure(lambda)->function = &function;
as closure(lambda)->n = ∣φ e∣;
as closure(lambda)->environment[0] = S π (φ e)0
⋮
as closure(lambda)->environment[∣φ e∣ − 1] = S π (φ e)∣φ e∣−1
lambda;
})

(22a)
(22b)
(22c)
(22d)
(22e)
(22f)
(22g)
(22h)
(22i)
(22j)
(22k)

Figure 24. Compiler for the checkpointVLAD language when in CPS. Part I.

binary bisection checkpointing driver (Fig. 18) are implemented in Scheme. The second
implementation is a hybrid compiler/interpreter that translates the checkpointVLAD
source program into CPS using CPS conversion (Fig. 22) and then interprets this with an
interpreter based on the direct-style evaluator (Figs. 9, 10, and 23), where the compiler,
the evaluator, the operator overloading implementation of AD, the general-purpose interruption and resumption mechanism (Fig. 23), and the binary bisection checkpointing
driver (Fig. 18) are implemented in Scheme. The third implementation is a compiler
that translates the checkpointVLAD source program into CPS using CPS conversion (Fig. 22) and then compiles this to machine code via c using GCC, where the
compiler (Figs. 24 and 25) is implemented in Scheme, the evaluator is the underlying
hardware, and the operator overloading implementation of AD, the general-purpose interruption and resumption mechanism (Fig. 23), and the binary bisection checkpointing
driver (Fig. 13) are implemented in c. The first implementation was used to generate the
results reported in [20] and presented at AD (2016). The techniques of Figs. 20 and 22
were presented at AD (2016). The third implementation was used to generate the results
reported here.

28

S π (λ4 k n l x.e) = ({

(22l)
thing function(thing target,
thing continuation,
thing count,
thing limit,
thing argument) {
return (S (φ e) e);
}
thing lambda = (thing)GC malloc(sizeof(struct {
enum tag tag;
struct {
thing (*function)();
unsigned n;
thing environment[∣φ e∣];
}))
}
set closure(lambda);
as closure(lambda)->function = &function;
as closure(lambda)->n = ∣φ e∣;
as closure(lambda)->environment[0] = S π (φ e)0
⋮
as closure(lambda)->environment[∣φ e∣ − 1] = S π (φ e)∣φ e∣−1
lambda;

})

S π (e1 e2 e3 e4 ) = continuation apply((S π e1 ),
(S π e2 ),
(S π e3 ),
(S π e4 ))

S π (e1 e2 e3 e4 e5 ) = converted apply((S π e1 ),
(S π e2 ),
(S π e3 ),
(S π e4 ),
(S π e5 ))
S π (if e1 then e2 else e3 ) = (!is false((S π e1 ))?(S π e2 ):(S π e3 ))
S π (◇e) = (N ◇)((S π e))
S π (e1 ● e2 ) = (N ●)((S π e1 ), (S π e2 ))
Ð
→
Ð
→
S π ( J e1 e2 e3 ) = (N J )((S π e1 ), (S π e2 ), (S π e3 ))
←
Ð
←
Ð
S π ( J e1 e2 e3 ) = (N J )((S π e1 ), (S π e2 ), (S π e3 ))
✓

✓

S π (J e1 e2 e3 ) = (N J )((S π e1 ), (S π e2 ), (S π e3 ))

(22m)

(22n)

(22o)
(22p)
(22q)
(22r)
(22s)
(22t)

Figure 25. Compiler for the checkpointVLAD language when in CPS. Part II.

4.

Complexity

The internal nodes of a binary checkpoint tree correspond to invocations of interrupt
in step (2). The right branches of each node correspond to step (3). The left branches
←
Ð
of each node correspond to step (4). The leaf nodes correspond to invocations of J in
the base case, step (0). Each leaf node corresponds to a stage, the red, blue, and violet
lines in Fig. 2(g). The checkpoint tree is traversed in depth-first right-to-left preorder.
In our implementation, we terminate the recursion when the step limit l is below a fixed
constant. Consider a general primal computation f that uses maximal live storage w
and that runs for t steps. This results in the following space and time complexities for
←
Ð
reverse mode without checkpointing, including our implementation of J , and for binary
29

✓

bisection checkpointing, including our implementation of J .
Complexity
space
Computation

without checkpointing
binary bisection checkpointing

time

primal

snapshots

tape

total

overhead

recomputation

O (w )
O (w )

O(w log t)

O ( t)
O(1)

O ( w + t)
O(w log t)

O( ww+t )
O(log t)

O(t log t)

forward
sweep

O ( t)
O ( t)

reverse
sweep

O ( t)
O ( t)

total

overhead

O ( t)
O(t log t)

O(1)
O(log t)

If we assume that O(t) ≥ O(w), i.e., that the computation uses all storage, the total
space requirement without checkpointing becomes O(t) and the overhead becomes O(t).
5.

Extensions to Support Treeverse and Binomial Checkpointing
✓

The general-purpose interruption and resumption interface allows implementation of J
using the treeverse algorithm from [9, Fig. 4] as shown in Fig. 26. This supports the full
functionality of that algorithm with the ability to select arbitrary execution points as
split points. By selecting the choice of mid as either [9, equation (12)] or [9, equation
(16)], one can select between bisection and binomial checkpointing. By selecting which
of d, t, and α the user specifies, computing the others from the ones specified, together
with n, using the methods described in [9] one can select the termination criterion to be
either fixed space overhead, fixed time overhead, or logarithmic space and time overhead.9
All of this functionality has been implemented in all three of our implementations: the
interpreter, the hybrid compiler/interpreter, and the compiler.
But it turns out that the binary checkpointing algorithm from Fig. 18 can be easily modified to support all of the functionality of the treeverse algorithm, including the
ability to select either bisection or binomial checkpointing and the ability to select any
of the termination criteria, including either fixed space overhead, fixed time overhead,
or logarithmic space and time overhead, with exactly the same guarantees as treeverse.
The idea is simple and follows from the observation that the right branch introduces a
snapshot and the left branch introduces (re)computation of the primal. One maintains
two counts, a right-branch count δ and a left-branch count τ , decrementing them as one
descends into a right or left branch respectively, to limit the number of snapshots or
the amount of (re)computation introduced. The base case is triggered when either gets
to zero or a specified constant bound on the number of steps to be taped is reached.
The binary checkpoint tree so produced corresponds to the associated n-ary checkpoint
tree produced by treeverse, as discussed in Section 1. Again, this supports the full functionality of treeverse with the ability to select arbitrary execution points as split points.
By selecting the choice of mid as either [9, equation (12)] or [9, equation (16)], one can
select between bisection and binomial checkpointing. By selecting which of d, t, and α
the user specifies, computing the others from the ones specified, together with n, using
the methods described in [9] one can select the termination criterion to be either fixed
space overhead, fixed time overhead, or logarithmic space and time overhead. All of this
functionality has been implemented in all three of our implementations: the interpreter,
the hybrid compiler/interpreter, and the compiler.
As per [9], with a binomial strategy for selecting split points, the termination criteria
can be implemented as follows. Given a measured number n of evaluation steps, d, t,
9

In Section 5 and Figs. 26 and 27 we use notation similar to that in [9] to facilitate understanding. Thus n and t
here means something different then elsewhere in this manuscript.

30

treeverse f x ỳ α δ τ β σ φ = if σ > β
then let z = interrupt f x (σ − β )
in first R z ỳ α (δ − 1) τ β σ φ
else first f x ỳ α δ τ β σ φ

first f x ỳ α δ τ β σ φ = if φ − σ > α ∧ δ =/ 0 ∧ τ =/ 0
then let κ = mid δ τ σ φ
(y, z̀) = treeverse f x ỳ α δ τ σ κ φ
in rest f x z̀ α y δ (τ − 1) β σ κ
←
Ð
else J (I f (φ − σ )) x ỳ

rest f x ỳ α y δ τ β σ φ = if φ − σ > α ∧ δ =/ 0 ∧ τ =/ 0
then let κ = mid δ τ σ φ
( , z̀ ) = treeverse f x ỳ α δ τ σ κ φ
in rest f x z̀ α y δ (τ − 1) β σ κ
←
Ð
else let ( , x̀) = J (I f (φ − σ )) x ỳ
in (y, x̀)
✓

J f x ỳ = let n = primops f x
pick α d t
in treeverse f x ỳ α d t 0 0 n

Figure 26. Implementation of treeverse from [9, Fig. 4] using the general-purpose interruption and resumption
interface, written in a functional style with no mutation. The variables δ, τ , β, σ, φ, n, d, and t have the same
meaning as in [9]. The variables f , x, x̀, y, ỳ, z, and z̀ have the same meaning as earlier in this manuscript.
The variable α denotes an upper bound on the number of evaluation steps for a leaf node. Different termination
criteria allow the user to specify some of α, d, and t and compute the remainder as a function of the ones specified,
together with n.

binary f x ỳ α δ τ φ = if φ ≤ α ∨ δ = 0 ∨ τ = 0
←
Ð
then J f x ỳ
else let κ = mid δ τ 0 φ
z = interrupt f x κ
(y, z̀ ) = binary R z ỳ (δ − 1) τ (φ − κ)
(z, x̀) = binary (I f κ) x z̀ δ (τ − 1) κ
in (y, x̀)
✓

J f x ỳ = let n = primops f x
pick α d t
in binary f x ỳ α d t n

Figure 27. Implementation of binary checkpointing using the general-purpose interruption and resumption interface in a fashion that supports all of the functionality of treeverse from [9, Fig. 4]. The variables δ, τ , φ, n, d,
and t have the same meaning as in [9]. The variables f , x, x̀, y, ỳ, z, and z̀ have the same meaning as earlier in this
manuscript. The variable α denotes an upper bound on the number of evaluation steps for a leaf node. Different
termination criteria allow the user to specify some of α, d, and t and compute the remainder as a function of the
ones specified, together with n.

and α are mutually constrained by a single constraint.
n = primops f x

η (d, t) = (

31

d+t
)
t

α=⌈

n
⌉
η (d, t)

(23)

( define ( car ( cons car cdr )) car )
( define ( cdr ( cons car cdr )) cdr )
( define ( first x ) ( car x ))
( define ( rest x ) ( cdr x ))
( define ( second x ) ( first ( rest x )))
( define ( iota n ) ( if ( zero ? n ) ’() ( cons n ( iota (- n 1)))))
( define ( ilog2 l ) ( floor (/ ( log l ) ( log 2))))
( define ( rotate theta x1 x2 )
( let (( c ( cos theta )) ( s ( sin theta )))
( cons (- (* c x1 ) (* s x2 )) (+ (* s x1 ) (* c x2 )))))
( define ( rot1 theta x )
( if ( or ( null ? x ) ( null ? ( rest x )))
x
( let (( x12 ( rotate theta ( first x ) ( second x ))))
( cons ( car x12 )
( cons ( cdr x12 ) ( rot1 theta ( rest ( rest x ) ) ) ) ) ) ) )
( define ( rot2 theta x )
( if ( null ? x ) x ( cons ( first x ) ( rot1 theta ( rest x )))))
( define ( magsqr x )
( if ( null ? x ) 0 (+ (* ( first x ) ( first x )) ( magsqr ( rest x )))))

( define ( w r i t e - v e c t o r v )
( if ( null ? v ) ’() ( cons ( w r i t e - r e a l ( first v )) ( w r i t e - v e c t o r ( rest v )))))
( define ( f x l phi )
( let outer (( i 1) ( x1 x ))
( if ( > i l )
(/ ( magsqr x1 ) 2)
( let (( m ( expt
2
(- ( ilog2 l )
( ilog2
(+ 1 ( modulo (* (* 1013 ( floor ( expt 3 phi ))) i ) l ) ) ) ) ) ) )
( let inner (( j 1) ( x1 x1 ))
( if ( > j m )
( outer (+ i 1) x1 )
( inner (+ j 1)
( let (( y ( sqrt ( magsqr x1 ))))
( rot2 (* 1.4 y ) ( rot1 (* 1.2 y ) x1 ) ) ) ) ) ) ) ) ) )
( let * (( n ( r e a d - r e a l ))
( l ( r e a d - r e a l ))
( phi ( r e a d - r e a l ))
( x ( iota n ))
( result ( c h e c k p o i n t - * j ( lambda ( x ) ( f x l phi )) x 1)))
( cons ( w r i t e - r e a l ( car result )) ( w r i t e - v e c t o r ( cdr result ))))

Figure 28. A rendering of the example from Fig. 8 in checkpointVLAD. Space and time overhead of two variants
of this example when run under checkpointVLAD are presented in Fig. 30. The variant for divide-and-conquer
checkpointing is shown. The variant with no checkpointing replaces checkpoint-*j with *j.

One can√ select any two and determine the third. Selecting d and α to determine
t = O ( d n) yields the fixed space overhead termination criterion. Selecting t and α
√
to determine d = O ( t n) yields the fixed time overhead termination criterion. Alternatively, one can further constrain d = t. With this, selecting α to determine d and t yields
the logarithmic space and time overhead termination criterion.

6.

An Example

As discussed in Section 2, existing implementations of divide-and-conquer checkpointing,
such as Tapenade, are limited to placing split points at execution points corresponding
to particular syntactic program points in the source code, i.e., loop iteration boundaries.
Our approach can place split points at arbitrary execution points. The example in Fig. 8
illustrates a situation where placing split points only at loop iteration boundaries can fail
to yield the sublinear space overhead of divide-and-conquer checkpointing while placement of split points at arbitrary execution points will yield the sublinear space overhead
of divide-and-conquer checkpointing. To illustrate this, we run the Fortran variant of
this example two different ways with Tapenade:
(1) without checkpointing, by removing all pragmas and
(2) with divide-and-conquer checkpointing, particularly the treeverse algorithm applied
to a root execution interval corresponding to the invocations of the outer DO loop,
split points selected with the binomial criterion from execution points corresponding
to iteration boundaries of the outer DO loop, and a fixed space overhead termination
criterion, by placing the c$ad binomial-ckp pragma as shown in Fig. 8.
For comparison, we reformulate this Fortran example in checkpointVLAD (Fig. 28)
and run it two different ways:
←
Ð
(1) without checkpointing, by calling J , written here as *j and
(2) with divide-and-conquer checkpointing, particularly the binary checkpointing algorithm applied to a root execution interval corresponding to the entire derivative
calculation, split points selected with the bisection criterion from arbitrary execution points, and a logarithmic space and time overhead termination criterion, by
✓

calling J , written here as checkpoint-*j.
For this example, n is the input dimension and l is the number of iterations of the
outer loop. Using the notation from Section 4, the maximal space usage of the primal
32

for this example should be w = O(n). The time required for the primal for this example
should be t = O(l), since there are l iterations of the outer loop and the inner loop has
average case O(1) iterations per iteration of the outer loop. The analysis in Sections 4
and 5 predicts the following asymptotic space and time complexity of the Tapenade
and checkpointVLAD variants that compute gradients on this particular example:
Complexity
Computation

space

time

O (n )

primal
Tapenade no checkpointing
Tapenade divide-and-conquer checkpointing
checkpointVLAD no checkpointing
checkpointVLAD divide-and-conquer checkpointing

O (n + l )
O(nl)10

O (n + l )
O(n log l)

O (l )

O(l)√
d
O (l l)

O (l )
O(l log l)

The efficacy of our method can be seen in the plots (Fig. 30) of the observed space and
time usage of the above two Fortran variants and the above two checkpointVLAD
variants with varying l and n = 1000. We observe that Tapenade space and time usage
←
Ð
grows with l for all cases. checkpointVLAD space and time usage grows with l with J .
✓

checkpointVLAD space usage is sublinear with J . checkpointVLAD time usage
✓

7.

Discussion

Our current implementations, as evaluated in Section 6, are expository prototypes and not intended as practical
artifacts. Nonetheless, technology exists
that can support construction of largescale practical and efficient implementations based on the conceptual ideas presented here.

Maximum resident set size (kbytes)

grows with l with J .
The crucial aspect of this example is that we observe sublinear space usage overhead
with divide-and-conquer checkpointing in checkpointVLAD but not in Tapenade.11
The reason that we fail to observe sublinear space usage overhead with divide-andconquer checkpointing in Tapenade is that the space overhead guarantees only hold
when the asymptotic time complexity of the loop body is constant. Since the asymptotic time complexity of the loop body is O(l), the requisite tape size grows with O(l)
even though the number of snapshots, and the size of those snapshots, is bounded by a
space
constant.
40000
Tapenade treeverse binomial checkpointing
checkpointVLAD binary bisection checkpointing

35000
30000
25000
20000
15000
10000
5000
0
0

500

1000

1500

2000

2500

3000

3500

l

Figure 29. Comparison of space usage of the the example in Figs. 8 and 28 when run with divide-andconquer checkpointing. checkpointVLAD achieves sublinear growth while Tapenade does not, thus for long
enough run times, the space usage of Tapenade exceeds
that of checkpointVLAD.
10
The space complexity of Tapenade with divide-and-conquer checkpointing would be O(n) if the inner DO loop
would have a constant number of iterations. The fact that it has average case O(1) iterations but O(l) worst case,
foils checkpointing and causes the space complexity to increase.
11
Technically, the space and time usage overhead of the checkpointVLAD variant of this example with divideand-conquer checkpointing should be logarithmic. It is difficult to see that precise overhead in the plots. Linear
regression does indeed fit logarithmic growth to the time usage better than linear growth. But the observed space
usage appears to be grow in steps. This is likely due to the coarse granularity of the measurement techniques that
are based on kernel memory page allocation and thus fail to measure the actual live fraction of the heap data
managed by the garbage collector.

33

space

time
9000

Tapenade no
Tapenade treeverse binomial
checkpointVLAD no
checkpointVLAD binary bisection

7x107
6x107

checkpointing
checkpointing
checkpointing
checkpointing

Tapenade no
Tapenade treeverse binomial
checkpointVLAD no
checkpointVLAD binary bisection

8000
User time (seconds)

Maximum resident set size (kbytes)

8x107

5x107
4x107
3x107
2x107

7000
6000
5000
4000
3000
2000

1x107

1000

0

0
0

500

1000

1500

2000

2500

3000

3500

0

500

1000

1500

2000

l
space

3000

3500

time
7

Tapenade no checkpointing
Tapenade treeverse binomial checkpointing
O(l)
O(l)

350000
300000
250000
200000
150000
100000
50000

5
4
3
2
1

0

0

-50000

-1
0

500

1000

1500

2000

2500

3000

Tapenade no checkpointing
Tapenade treeverse binomial checkpointing
O(l)
O(l**(1+1/40))

6
User time (seconds)

Maximum resident set size (kbytes)

2500

l

400000

3500

0

500

1000

1500

l

2000

2500

3000

3500

l

space

time

8x107

9000
checkpointVLAD no checkpointing
checkpointVLAD binary bisection checkpointing
O(l)
O(log l)

7x107
6x107

checkpointVLAD no checkpointing
checkpointVLAD binary bisection checkpointing
O(l)
O(l log l)

8000
7000
User time (seconds)

Maximum resident set size (kbytes)

checkpointing
checkpointing
checkpointing
checkpointing

5x107
4x107
3x107
2x107
1x107

6000
5000
4000
3000
2000
1000

0

0

-1x107

-1000
0

500

1000

1500

2000

2500

3000

3500

l

0

500

1000

1500

2000

2500

3000

3500

l

Figure 30. Space and time usage of reverse-mode AD with and without divide-and-conquer checkpointing for the
example in Figs. 8 and 28. Space and time usage was measured with /usr/bin/time --verbose The center and
bottom rows repeat the information from the top row just for Tapenade and checkpointVLAD, overlaid with
the theoretical asymptotic complexity fit to the actual data by linear regression.

7.1

Implementation Technologies

One can use posix fork() to implement the general-purpose interruption and resumption
interface, allowing it to apply in the host, rather than the target, and thus it could be
used to provide an overloaded implementation of divide-and-conquer checkpointing in a
fashion that was largely transparent to the user [9]. The last paragraph of [9] states:
For the sake of user convenience and computational efficiency, it would be ideal if reverse
automatic differentiation were implemented at the compiler level.

We have exhibited such a compiler. Our implementation, however, is not as efficient as
Tapenade. There can be a number of reasons for this. First, Fortran unboxes double precision numbers whereas checkpointVLAD boxes them. This introduces storage
allocation, reclamation, and access overhead for arithmetic operations. Second, array access and update in Fortran take constant time whereas access and update in checkpointVLAD take linear time. Array update in checkpointVLAD further involves
storage allocation and reclamation overhead. Third, Tapenade implements the base
34

case reverse mode of divide-and-conquer checkpointing using source-code transformation
whereas checkpointVLAD implements it with operator overloading. In particular, the
dynamic method for supporting nesting involves tag dispatch for every arithmetic operation.
We have exhibited a very aggressive compiler (Stalin∇) for vlad that ameliorates
some of these issues. It unboxes double precision numbers and implements AD via sourcecode transformation instead of operator overloading as is done by checkpointVLAD.
This allows it to have numerical performance rivaling Fortran. While it does not support constant-time array access and update, methods that are well-known in the programming languages community (e.g., monads and uniqueness types) can be used for that.
But it does not include support for divide-and-conquer checkpointing. However, there is
no barrier, in principle, to supporting divide-and-conquer checkpointing, of the sort described above, in an aggressive optimizing compiler that implements AD via source-code
transformation. One would simply need to reformulate the source-code transformations
that implement AD, along with the aggressive compiler optimizations, in CPS instead
of direct style. Moreover, the techniques presented here could be integrated into other
compilers for other languages that generate target code in CPS by instrumenting the
CPS conversion with step counting, step limits, and limit-check interruptions.12 A driver
✓

can be wrapped around such code to implement J . For example, existing compilers, like
sml/nj [4], for functional languages like sml, already generate target code in CPS, so it
seems feasible to adapt their techniques to the purpose of AD with divide-and-conquer
checkpointing. In fact, the overhead of the requisite instrumentation for step counting,
step limits, and limit-check interruptions need not be onerous because the step counting,
step limits, and limit-check interruptions for basic blocks can be factored, and those for
loops can be hoisted, much as is done for the instrumentation needed to support storage
allocation and garbage collection in implementations like MLton [27], for languages like
sml, that achieve very low overhead for automatic storage management.
7.2

Advantages of Functional Languages for Interruption and Resumption

Functional languages simplify interruption and resumption, allowing these to be much
more efficient. Two different capsules taken at two different execution points can share
common substructure, by way of pointers, without needing to copy that substructure.
Indeed, CPS in the checkpointVLAD implementation renders all program state, including the stack and variables in the environment, as closures, possibly nested. Creating a capsule simply involves saving a pointer to a closure. Resuming a capsule simply
involves invoking the saved closure, a simple function call that is passed the closure environment as its argument. The garbage collector can traverse the pointer structure of the
program state to determine the lifetime of a capsule. The interruption and resumption
framework need not do so itself. This simplicity and efficiency would be disrupted by
structure mutation or assignment.
7.3

Nesting of AD Operators

It has been argued that the ability to nest AD operators is important in many practical
domains [2, 3, 7]. Supporting nested use of AD operators involves many subtle issues
12 We

note that many optimizing compilers, for example GCC, use an intermediate program representation called
Single Static Assignment, or SSA, which is formally equivalent to CPS [16].

35

[19]. checkpointVLAD addresses these issues and fully supports nested use of AD
operators. One can write programs of the form
α (λx.(. . . (β f . . .))) . . .

(24)

Ð
→ ←
Ð
✓
where each of α and β can be any of J , J , and J . I.e., one can apply AD to one
function that, in turn, applies AD to another function. This allows one not only to do
forward-over-reverse, reverse-over-forward, and reverse-over-reverse, it also allows one to
do things like divide-and-conquer-reverse-over-forward, reverse-over-divide-and-conquerreverse, and even divide-and-conquer-reverse-over-divide-and-conquer-reverse.
There is one catch however. When one applies an AD operator, that application is
considered to be atomic by an application of a surrounding AD operator. This is evident
✓

by the n + 1 in (5l, 5m, 6r), What this means is that if α in (24) were J , a split point
for α could not occur inside f . While this does not affect the correctness of the result, it
could affect the space and time complexity.
The reason for this is that while the semantics of the AD operators are functional,
Ð
→
their internal implementation involves mutation. In particular, to support nesting, J ,
←
Ð
✓
J , and J internally maintain and update an ǫ tag as described in [19]. The ǫ tag
is incremented upon entry to an AD operator and decremented upon exit from that
invocation to keep track of the nesting level. All computation within a level must be
performed with the same ǫ tag. This requires that the entries to and exits from AD
✓

operator invocations obey last-in-first-out sequencing. If α were J , and the computation
of f were interrupted, then situations could arise where the last-in-first-out sequencing
of β was violated. Moreover, since divide-and-conquer checkpointing executes different
portions of the forward sweep different numbers of times, the number of entries into a
nested AD operator could exceed the number of exits from that operator.
Reverse mode involves a further kind of mutation. The forward sweep creates a tape
represented as a directed acyclic graph. The nodes in this tape contain slots for the
cotangent values associated with the corresponding primal values. The reverse sweep
operates by traversing this graph to accumulate the cotangents in these slots. Such
accumulation is done by mutation.
The above issues arise because of mutation in the implementation of AD operators.
Conceivably, these could be addressed using methods that are well-known in the programming languages community for supporting mutation in functional languages (e.g.,
✓

monads and uniqueness types). Issues arise beyond this, however. If both α and β were J ,
✓

and J was not atomic, situations could arise where f could interrupt for α instead of β.
Currently, interruption is indicated by returning instead of calling a continuation. If f
were to return instead of calling a continuation, there is no way to indicate that that
interruption was due to α instead of β. It is unclear whether this issue could be resolved.

8.

Conclusion

Reverse-mode AD with divide-and-conquer checkpointing is an enabling technology, allowing gradients to be efficiently calculated even where classical reverse mode imposes an
impractical storage overhead. We have shown that it is possible to provide an operator
that implements reverse-mode AD with divide-and-conquer checkpointing, implemented
36

as an interpreter, a hybrid compiler/interpreter, and a compiler, which

● has an identical API to the classical reverse-mode AD operator,
● requires no user annotation,
● takes the entire derivative calculation as the root execution interval, not just the execution intervals corresponding to the invocations of particular constructs such as DO
loops,
● takes arbitrary execution points as candidate split points, not just the execution points
corresponding to the program points at the boundaries of particular constructs like
the iteration boundaries of DO loops,
● supports both an algorithm that constructs binary checkpoint trees and the treeverse
algorithm that constructs n-ary checkpoint trees,
● supports selection of actual split points from candidate split points using both a bisection and a binomial criterion, and
● supports any of the termination criteria of fixed space overhead, fixed time overhead,
or logarithmic space and time overhead,
yet still provides the favorable storage requirements of reverse mode with divide-andconquer checkpointing, guaranteeing sublinear space and time overhead.

Acknowledgments
This manuscript is an expanded version of the extended abstract of a presentation at
AD (2016), which is available as a technical report [20]. Aspects of this work may be
subject to the intellectual property policies of Purdue University and/or Maynooth University. Code for the implementations described above and related materials are available
at https://github.com/qobi/checkpoint-VLAD/. This work was supported, in part,
by the US National Science Foundation under Grants 1522954-IIS and 1734938-IIS, by
the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/Interior Business Center (DOI/IBC) contract number D17PC00341, and by Science
Foundation Ireland under Grant SFI 09/IN.1/I2637. 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.

References
[1] P. Achten, J. Van Groningen, and R. Plasmeijer, High level specification of I/O in functional languages, in Functional Programming, Glasgow 1992, Springer, 1993, pp. 1–17.
[2] N. Agarwal, B. Bullins, and E. Hazan, Second order stochastic optimization in linear time (2016),
Available at https://arxiv.org/abs/1602.03943.
[3] M. Andrychowicz, M. Denil, S.G. Colmenarejo, M.W. Hoffman, D. Pfau, T. Schaul, and N. de
Freitas, Learning to learn by gradient descent by gradient descent, in Neural Information Processing
Systems, Available at https://arxiv.org/abs/1606.04474, NIPS, 2016.
[4] A.W. Appel, Compiling with continuations, Cambridge University Press, 2006.
[5] C. Bendtsen and O. Stauning, FADBAD, a flexible C++ package for automatic differentiation, Technical Report IMM–REP–1996–17, Department of Mathematical Modelling, Technical University of
Denmark, Lyngby, Denmark, 1996.
[6] T. Chen, B. Xu, Z. Zhang, and C. Guestrin, Training deep nets with sublinear memory cost (2016),
Available at https://arxiv.org/abs/1604.06174.

37

[7] B. Christianson, Reverse accumulation of functions containing gradients, Tech. Rep.
278, University of Hertfordshire Numerical Optimisation Centre, 1993, Available at
http://hdl.handle.net/2299/4337, presented at the Theory Institute Argonne National
Laboratory Illinois, Procs. of the Theory Institute on Combinatorial Challenges in Automatic
Differentiation.
[8] B. Dauvergne and L. Hascoët, The Data-Flow Equations of Checkpointing in Reverse Automatic
Differentiation, in Computational Science – ICCS 2006, Lecture Notes in Computer Science, Vol.
3994, Springer, Heidelberg, 2006, pp. 566–573.
[9] A. Griewank, Achieving logarithmic growth of temporal and spatial complexity in reverse automatic
differentiation, Optimization Methods and Software 1 (1992), pp. 35–54.
[10] A. Griewank, D. Juedes, and J. Utke, A package for the automatic differentiation of algorithms written in C/C++. User manual, Tech. Rep., Institute of Scientific Computing, Technical University of Dresden, Dresden, Germany, 1996, this version of the manual is superceded by
http://www.math.tu-dresden.de/~ adol-c/adolc110.ps.
[11] A. Gruslys, R. Munos, I. Danihelka, M. Lanctot, and A. Graves, Memory-Efficient
Backpropagation Through Time, in Neural Information Processing Systems, Available at
https://arxiv.org/abs/1606.03401, NIPS, 2016.
[12] L. Hascoët and V. Pascual, TAPENADE 2.1 user’s guide, Rapport technique 300, INRIA, 2004.
[13] C.T. Haynes and D.P. Friedman, Engines build process abstractions, in Proceedings of the 1984 ACM
Symposium on LISP and functional programming, ACM, 1984, pp. 18–24.
[14] R. Heller, Checkpointing without operating system intervention: Implementing griewank’s algorithm,
Master’s thesis, Ohio University, 1998.
[15] Y. Kang, Implementation of forward and reverse mode automatic differentiation for GNU Octave
applications, Master’s thesis, Ohio University, 2003.
[16] R.A. Kelsey, A correspondence between continuation passing style and static single assignment form,
ACM SIGPLAN Notices, Papers from the 1995 ACM SIGPLAN Workshop on Intermediate Representations 30 (1995), pp. 13–22.
[17] A. Kowarz and A. Walther, Optimal Checkpointing for Time-Stepping Procedures in ADOL-C, in
Computational Science – ICCS 2006, Lecture Notes in Computer Science, Vol. 3994, Springer,
Heidelberg, 2006, pp. 566–573.
[18] J.C. Reynolds, The discoveries of continuations, Lisp and Symbolic Computation 6 (1993), pp.
233–247.
[19] J.M. Siskind and B.A. Pearlmutter, Nesting forward-mode AD in a functional framework, HigherOrder and Symbolic Computation 21 (2008), pp. 361–376.
[20] J.M. Siskind and B.A. Pearlmutter, Binomial checkpointing for arbitrary programs with no user
annotation (2016), Available at https://arxiv.org/abs/1611.03410.
[21] J.M. Siskind and B.A. Pearlmutter, Efficient implementation of a higher-order language with built-in
AD (2016), Available at https://arxiv.org/abs/1611.03416.
[22] B. Speelpenning, Compiling fast partial derivatives of functions given by algorithms, Ph.D. diss.,
Department of Computer Science, University of Illinois at Urbana-Champaign, 1980.
[23] A. Stovboun, A tool for creating high-speed, memory efficient derivative codes for large scale applications, Master’s thesis, Ohio University, 2000.
[24] G.J. Sussman, J. Wisdom, and M.E. Mayer, Structure and interpretation of classical mechanics,
MIT Press, 2001.
[25] Y.M. Volin and G.M. Ostrovskii, Automatic computation of derivatives with the use of the multilevel
differentiating technique — I: Algorithmic basis, Computers and Mathematics with Applications 11
(1985), pp. 1099–1114.
[26] P.L. Wadler, Comprehending monads, in Proceedings of the 1990 ACM Conference on Lisp and
Functional Programming, ACM, 1990, pp. 61–78.
[27] S.
Weeks,
Whole-program
compilation
in
MLton
(2006),
Available
at
http://www.mlton.org/References.attachments/060916-mlton.pdf, workshop on ML.

38

