Mathematics Physics Computer Science
Computer Science-2021.gz:143154
Computer Science
60399585a9e6eca70f2237274276a55bae3cda7b
A parallel-in-time approach for accelerating direct-adjoint studies

Parallel-in-time methods are developed to accelerate the direct-adjoint looping procedure. Particularly, we utilize the Paraexp algorithm, previously developed to integrate equations forward in time, to accelerate the direct-adjoint looping that arises from gradient-based optimization. We consider both linear and non-linear governing equations and exploit the linear, time-varying nature of the adjoint equations. Gains in efficiency are seen across all cases, showing that a parallel-in-time approach is feasible for the acceleration of direct-adjoint studies. This signifies a possible approach to further increase the run-time performance for optimization studies that either cannot be parallelized in space or are at their limit of efficiency gains for a parallel-in-space approach.


Introduction
The use of direct-adjoint looping techniques to compute the extremum of an objective functional, subject to a state variable satisfying complementary governing (direct) equations, is the subject of a wide range of optimization studies. Direct-adjoint looping is a technique that originated in control theory and was pioneered in fluid dynamics by Jameson [13] to find the optimal aerofoil shape that minimizes drag. These Lagrangian optimization problems are formulated with the governing equations included as constraints. The direct (governing) equations are integrated forward in time yielding a solution for a set of control parameters. Afterwards, the adjoint equations are integrated backwards in time, providing the gradient of the objective functional with respect to the parameters. This information can then be used as part of a gradient-based optimization routine to find the extremum of the objective functional via repeated applications of this direct-adjoint loop. Recent examples include, among others, the optimization of mixing in binary fluids [8,19,4,6,5], finding the minimal seed that triggers turbulence in pipe flow [21] and determining the optimal place to ignite a diffusion flame [22].
While the inclusion of the adjoint procedure significantly speeds up the calculation of the gradient information, several complications can arise due to the use of the adjoint. Firstly, the inherently iterative nature of direct-adjoint looping can cause optimization studies performed in this manner to become prohibitively slow, despite the speedup achieved through the efficient computation of the gradient. Another complication arises when the adjoint equations depend on the direct solution. This necessarily implies that the full direct solution must be saved and injected into the adjoint equation at specific points in time to obtain accurate solutions for the adjoint variables, leading to high memory costs. To circumvent these memory costs the checkpointing library revolve [11] can be used. In this way, the direct solution is saved at various checkpoints in time, and, at the expense of recomputing the solution, the direct-adjoint loop can then be solved between checkpoints leading to a reduced memory cost. These two complications, among many others, can result in run-times that may be unfeasible for some optimization studies. Traditionally, to alleviate these problems and retain viable run-times, parallel-in-space approaches are used to speed up the calculation, examples of which include work by Pekurovsky [20], Laizet and Vassilicos [17], Kallala et al. [14]. However, when encountering problems that cannot be spatially parallelized or where the parallelization in space has reached its maximum efficiency, gains in run-time efficiency must be achieved through other means. In our study, we consider one such alternative method to accelerate scenarios where more computational power is available: parallelization in time.
Parallel-in-time approaches have previously been used to accelerate direct studies, with examples including the work of Farhat and Chandesris [7] or Lions et al. [18]. For an in-depth review, we refer the reader to the review of J. Gander [12] which describes a large variety of parallel-in-time algorithms. A particular algorithm which will be utilized in this study is the Paraexp algorithm developed by Gander and Güttel [9]. This algorithm achieves parallelization-in-time for a linear equation by isolating the homogeneous and inhomogeneous parts. The time domain is then partitioned and the inhomogeneous equations are solved on each time partition with a zero initial condition. Solving the homogeneous equations allows for propagation of the correct initial conditions throughout the domain, with the full solution obtained as a superposition of all inhomogeneous and homogeneous parts. Although the homogeneous solutions can span more than one time partition, and hence more equations are solved than the original problem, the fact that the homogeneous equations can be efficiently solved via exponential time integration, examples of which can be seen in Schulze et al. [24] and Bergamaschi et al. [1], leads to an overall speedup. Furthermore, this algorithm can be extended to non-linear equations via an iterative approach [16,10].
This present work aims to accelerate direct-adjoint looping by extending the Paraexp algorithm to also include the adjoint component. We consider both linear and non-linear equations, as well as the possible inclusion of a checkpointing scheme. For the case of non-linear governing equations two algorithms are developed; one in which the direct equation is solved iteratively in parallel, and another in which the direct equation is solved in series. In both cases, the adjoint equation is solved in parallel. When the direct equation is solved in series we obtain a speedup by overlapping the direct and adjoint solutions using a particularly efficient time partition. For all algorithms we derive theoretical expressions for the speedups. Furthermore, we describe in a step-by-step fashion the developed parallel-in-time algorithms and apply them to the 2D advectiondiffusion equation and viscous Burgers' equation for the linear and nonlinear cases, respectively. In all cases we observe a marked improvement in run-time efficiency in line with the theoretical predictions.
The rest of this article is organized in the following manner; §2 gives a brief overview of adjoint-looping techniques before moving on to §3 where parallelin-time algorithms are presented for the direct and adjoint system. Next, §4 covers the numerical implementation, and the resulting numerical experiments are compared to theoretical scalings. Finally, conclusions are offered in section §5.

Direct-adjoint looping
We begin by briefly outlining a general direct-adjoint problem [13]. Consider the non-linear governing (direct) equation where N (q, f ) represents the right-hand side of an equation dependent upon the state q and parameters f . We seek to minimize an objective functional J that also can depend on the state q and parameters f . To this end, we form the Lagrangian where the Lagrange multiplier q † is the adjoint variable. We take the inner products and the cost functional to be of the form where differing inner products and cost functionals are introduced for the discrete and continuous cases, respectively. In our formulation, Ω represents the spatial domain, and the cost functional is split into an integral part, J I , and a part that solely depends on the final time, J T . To minimize our objective functional subject to the constraint that our state satisfies the governing equation (1), we must ensure that all first variations of the Lagrangian are zero. Taking the variation with respect to the conjugate adjoint variable enforces the state equation (1). Before we take the variations with respect to the adjoint variable q † , we utilize integration by parts in time to rewrite the Lagrangian as follows Upon taking the variations of the Lagrangian in this form with respect to the state q, we arrive at the adjoint equation where (·) denotes the complex conjugate. The adjoint operator must be found such that the adjoint relation is satisfied for all vectors a and b. In the continuous case, this relation can be used to find the adjoint operator via integration by parts in space. For the discrete case, on the other hand, simple linear algebra manipulations can be used to find the adjoint operator As the adjoint equation (7) must be integrated backwards in time from the final time T , we require a final time condition for the adjoint variable, q † . This can found by ensuring the first variation with respect to q(T ) is zero, i.e., Finally, it remains to take the first variation with respect to f . In general, unless f is a local minimum of the Lagrangian, the first variation with respect to f will be non-zero. Instead, we obtain the gradient of the Lagrangian as Likewise, if the cost functional is to be optimized with respect to the initial condition, the associated gradient becomes In practice, achieving an analytic expressions for all required quantities is not feasible. Therefore, the governing equations and adjoint equations are solved explicitly, while some initial guess is employed for our control variable, f in the case of forcing or q 0 in the case of optimal initial condition. An iterative procedure, in conjunction with either equation (11) or (12) and a gradientdescent routine, is then utilized to drive the control to a value that achieves a local minimum for the chosen cost functional.

Linear governing equation
We now turn to the discussion of the parallel-in-time approach, particularly the Paraexp algorithm developed by Gander and Güttel [9]. The fundamental aspect of this approach is that a linear equation can be split into an inhomogeneous and homogeneous component. Consider the general linear equation where the inhomogeneity arises from the source term f . Following Kooij et al. [16], we make the substitution q = q − q 0 to obtain the equation which absorbs the initial condition into the governing equation. To solve this equation over the range t ∈ [0, T ] in parallel, we employ the Paraexp algorithm. For N processors the time domain is decomposed into the partition 0 = T 0 < T 1 < ... < T N = T . The p th processor, where available processors are labelled p ∈ {0, 1, ..., N − 1}, then solves the inhomogeneous equation on [T p , T p+1 ]. Thus the inhomogeneous part of the equation is solved on each time interval with a zero initial condition. Following this we solve the homogeneous problems (d) Finally, we sum all parts on each processor to find the full solution. in addition to the inhomogeneous evaluation covers the full domain. However, the speedup is achieved by recognizing that the homogeneous equations have the ability to be solved far more efficiently than the inhomogeneous equation, for example by using exponential time-stepping [23]. Hence, the slow-to-solve inhomogeneous parts are solved in parallel first and then the homogeneous parts are solved rapidly to propagate the correct initial conditions throughout the time domain. The full solution on the interval [T p , T p+1 ] can then be found by the sum In their implementation of the Paraexp algorithm, Gander and Güttel [9] suggest solving problems q I,p and q H,p+1 on the same processor so that communication between processors is only necessary when computing the sum (20) at the end. Here, however, we will keep the time domain completely separate, i.e., the p th processor will always solve equations on its time partition [T p , T p+1 ]. Therefore, whilst the inhomogeneous problems are separated by processors, an individual solution for the homogeneous part q H,j will span multiple processors, resulting in necessary communications between adjacent processors. The reason for this particular arrangement is due to the inclusion of the adjoint solver. The adjoint will require full state information on each time interval and therefore it is vital for efficiency that this is available without communicating a large amount of data across processors. Hence, we rewrite the homogeneous problems (16) as for p > 0 and j ∈ {0, ..., p − 1}, to be solved on [T p , T p+1 ] with the initial conditions Here the subscript p completely defines the processor that the relevant equation is solved on, and the p th processor solves p homogeneous problems. We note that the subscript j is not the same as in equation (16) but is instead used to denote the different homogeneous problems that must be solved for each processor. Each processor (except the zero th ) must solve the j = 0 homogeneous problem that stems from the final time inhomogeneous state from the previous processor. Then processor p must solve the additional homogeneous problems j = 1, ..., (p − 1) that arise from the homogeneous problems started on previous processors that must be integrated to the final time. By solving the equations this way the full solution on the interval [T p , T p+1 ] becomes Now that the algorithm for the direct equation has been presented we turn to the discussion on the inclusion of the adjoint solver. The adjoint equation to (13) for given cost functional (see section 2) is where we now have to integrate backwards in time, and the final time condition was absorbed into the equation using the substitution q † = q † − q † (T ). The linear nature of this equation means that we can use the Paraexp algorithm in exactly the same manner as previously described. The adjoint inhomogeneous problems that must be solved are on [T j+1 , 0] and for j ∈ {N − 2, ...0}. Once again, we keep the adjoint homogeneous problems local to each processor's time partition and instead solve the problems This time the p th processor solves (N − 1) − p homogeneous equations. Note that these equations can be solved forwards in time by making the substitution τ = T p + T p+1 − t on each time interval. The full adjoint solution on processor p is then given as the sum

Non-linear governing equation solved in parallel
We now extend the parallel-in-time techniques to a direct-adjoint loop for which the governing equation is non-linear. As the adjoint equation is always a linear (and possibly linear, time-varying) equation, the Paraexp algorithm as presented in section 3.1 can always be used to solve the adjoint equation in parallel. To obtain speedup for the direct equation we present two approaches. The first uses a non-linear extension of the Paraexp algorithm, namely the method of Kooij et al. [16], to accelerate the direct equation. Note that other non-linear extensions of the Paraexp algorithm are available, see for example the proceedings of Gander et al. [10], but we do not consider them here. Specifically, the approaches of Gander et al. [10] and Kooij et al. [16] are both iterative approaches and therefore should yield fairly similar scalings.
For the non-linear governing equation we take the general equation (1) from section 2. To parallelize the time-integration of this equation we proceed using the method of Kooij et al. [16] by rewriting it in the form where we have split the operator N into its linear part A and genuine non-linear part N(q). We then introduce the following recurrence relation and note that if we achieve convergence for our series of q i , i.e. q i+1 ≈ q i , then (28) is equivalent to equation (27). The benefit of this altered equation, is that this equation is linear in q i+1 , as the non-linearity only acts on the previous iterate and becomes a source term only. To achieve the convergence necessary for our recurrence relation to equate to the full equation, Kooij et al. [16] augment this equation with where the average Jacobian has been added. Once again the substitution q = q − q 0 is used to absorb the initial condition into the equation. The non-linear equations have now been reduced to a set of equivalent linear ones and can be solved using the Paraexp method described in the previous section, giving the inhomogeneous problems as on the intervals [T p , T p+1 ], where the Jacobians are now averaged on each time interval. Similarly, the homogeneous equations become with the initial conditions Once the homogeneous equations (33) are solved the next iterate of the direct solution can be found via the sum where is a user-specified tolerance, we obtain the full non-linear solution. Once the non-linear solution has been determined, the adjoint equation is then solved using the algorithm presented in section 3.1.

Governing equation solved in series
Whilst the previous section shows that the governing equation can be solved in parallel, problems may nevertheless arise. As there is no guarantee that the iterations will converge, it is possible that the direct solution may never be obtained. Additionally, even if the iterations do converge, the number of iterations required might outweigh the benefit obtained from solving the equation in parallel. In these cases an alternative approach is needed; we propose an algorithm in which the non-linear equation is solved in series, whilst accommodating the time-parallelization of the adjoint equation. This gives rise to the hybrid serial-direct-parallel-adjoint method presented in this section.
For the direct equation we partition the time interval so that the p th processor only solves equations on the interval [T p , T p+1 ]. The direct equation is solved explicitly by integrating the following equations in sequence where the notation q p is used to denote the full solution q solved on processor p. As the p th processor needs an initial condition from the (p − 1) th processor, later processors must wait until previous processors have finished their direct solutions. However, unlike in the previous algorithms, the full direct solution is available on the interval [T p , T p+1 ] as each processor p finishes solving equation (36). Therefore, while later processors are computing the direct solutions, earlier processors can begin solving the adjoint inhomogeneous equation (22). The algorithm proceeds in this manner with the direct solution being solved across all the processors, with the adjoint inhomogeneous equations (22) beginning on each processor once the direct solution has been computed. After the adjoint inhomogeneous equations have been solved, all that remains is to compute the adjoint homogeneous equations. In contrast to the previous algorithms in which an equidistant time partition would be a sensible choice, the hybrid algorithm requires slightly more attention. Using the hybrid scheme means that the inhomogeneous adjoint equations start on earlier processors first. If an equidistant time partition were used, then these earlier processors would terminate first and consequently would idle for information from subsequent processors to begin solving the homogeneous adjoint equation. This time spent waiting limits the efficiency of the algorithm, as unused computing power is squandered. Therefore, it is imperative to choose the time partitioning strategically, thus avoiding the inefficiency of the equallyspaced time grid.
To prevent earlier processors waiting on later processors, all the inhomogeneous adjoint equations being solved need to terminate at the same time, i.e., the time for the full solution to be completed on one processor after having completed on the previous ones plus the subsequent adjoint inhomogeneous evaluation should remain constant on each processor. By introducing the time τ I taken per time unit for the direct solve and similarly the time τ † I taken per time unit for the adjoint inhomogeneous solve, we can write this relationship in terms of the time partition for processors n and n + 1 as This relationship can be rearranged to produce with k = τ † I /τ I . This expression necessarily implies that the simultaneous termination of the inhomogeneous evaluations is achieved by making subsequent time partitions finer as we approach the final time. Equation (39) is a secondorder difference equation, with boundary conditions T 0 = 0, and T N = T which can be solved to yield giving an analytic expression for the time partition in terms of the ratio k, which can be measured numerically.
One important difference of this hybrid serial-direct-parallel-adjoint algorithm to the previous approaches must be pointed out. In the previous cases, an equidistant time partition can be used. This means that, assuming the adjoint homogeneous equations (23) take the same amount of time to be solved on each processor, that they can be solved by instead solving the series of equations (24). This keeps the solved equations local to each processor's time partition, meaning that the direct solution does not need to be distributed to other processors. However, in this serial-direct-parallel-adjoint approach the time partition is inherently non-equidistant. Solving the adjoint homogeneous part using the series of equations (24) will mean that later processors, which have a smaller time partition, will need to wait to pass their states to lower processors which solve equations on increasingly longer time intervals. The thus introduced lag scales with the number of processors, negating any speedup obtained from overlapping the direct and adjoint solutions. In order to circumvent this problem, the adjoint homogeneous equations must instead be solved in their original form given by equation (23). Although this means that the direct solution must be distributed to all processors, this is a one-time cost that must be paid, as it enables the algorithm to achieve an overall speedup. We choose to solve the j th adjoint homogeneous equation (23) on processor p. The full adjoint solution on [T p , T p+1 ] can then be found by the sum  figure 2c shows that, when the direct integration finishes for p = 1, the direct integration continues on processor two, whilst the inhomogeneous integration begins on processor one and continues on processor zero. Finally, figure 2d shows that after the direct integration has finished, all inhomogeneous adjoint equations terminate at the same time. The adjoint solution can then be completed by solving the homogeneous equations (23) in parallel (not shown).

Step-by-step summary of the algorithms and scaling analysis
Now that the three algorithms have been presented, we provide a step-bystep summary for ease of implementation. A scaling analysis is also conducted to assess the theoretical speedups possible.

Linear governing equation
Algorithm 1 summarizes the linear case discussed in section 3.1. To quantify the speedup for the algorithm, we will introduce the time taken per time unit to solve the direct inhomogeneous and adjoint equations, denoted τ I and τ H , respectively. Similarly, we introduce adjoint counterparts τ † I , τ † H . For every processor, we assume that step 4 of algorithm 1 takes the same amount of time to solve and occurs perfectly in sync. This means that the time taken to solve the inhomogeneous direct equations is (T /N )τ I . After the direct inhomogeneous equations are solved, the direct homogeneous equations must be solved. This is achieved by performing the for-loop beginning in step 11. Although each processor solves a different amount of direct homogeneous integrations, it is processor N − 1 that solves the most. Therefore, the time for this loop is (T (N − 1)/N )τ H . Hence, the total time to solve the direct equation in parallel is Performing the same analysis on the adjoint part of the algorithm gives the time taken to solve the adjoint equation as (a) First, the direct equation is solved on the zeroth processor.
(b) Whilst the first processor takes over the direct solution, the adjoint equation begins to be solved on the zeroth processor.
(c) The direct solution is completed on the second processor. Simultaneously, the adjoint inhomogeneous solution is continued on processor zero, and begins on processor one.

22:
Solve the adjoint inhomogeneous equation (22)  Adding the two parts together gives the total time for one direct-adjoint loop In series, the time taken for the algorithm is simply T S L = T (τ I + τ † I ). Hence, the speedup s = T S L /T L (shown as the reciprocal speedup for convenience) is The speedup (45) shows us that there is a theoretical maximum to the possible speedup. Indeed, as N → ∞, s → (τ I + τ † I )/(τ H + τ † H ). This highlights the importance of solving the homogeneous equation faster than the inhomogeneous equation. The faster the homogeneous integration over the inhomogeneous evaluation, the more efficient the parallel-in-time approach becomes.

Non-linear governing equation solved in parallel
The case of a non-linear governing equation, see section (3.2), solved in parallel is summarized in algorithm 2. It has a strong resemblance to algorithm 1 except for the fact that the direct equation is now solved in an iterative fashion with the while-loop starting on line 9. For the purposes of assessing the speedup, we will include the cost of calculating the Jacobian terms on line 27 with the direct inhomogeneous equation. Furthermore, as the inhomogeneous direct equation in parallel is not the same as the direct equation in series, we introduce the time taken per time unit τ S D for the direct equation in series. This gives the time taken for a direct-adjoint loop in series as T S = T (τ S D + τ † I ). We introduce here the number of iterations K(N ) we require to achieve convergence for our non-linear scheme. Note the dependence of the iterations while err > tol do 10:

12:
Send the state q I,p,iter+1 (T p+1 ) to processor p + 1 MPI communication 13: end if 14: if p > 0 then 15: Receive the state q I,p−1,iter+1 (T p ) from processor p − 1 MPI communication 16: end if 17: for j ∈ {0, ..., p − 1} do 18: if j > 0 then 19: Receive the state q H,p−1,j−1,iter+1 (T p ) from processor p−1 MPI communication 20: end if 21: Solve equation (33)  end for 50: end procedure on the number of processors; the fineness of the time partition will affect the number of iterations. As each iteration involves the same steps as the linear case we obtain that the time taken for the direct equation is simply K(N ) times equation (42), i.e., The time taken to solve the adjoint equation is still given by equation (43), due to the required steps being identical to algorithm 1. Hence, the time taken to perform one direct-adjoint loop in parallel is giving the reciprocal speedup as Equivalently to the linear case, equation (48) shows the importance of the homogeneous solvers being faster than the inhomogeneous ones. However, unlike the linear case, the factor of K in the speedup equation implies that small N can lead to s < 1, signifying that no speedup is possible using a parallel approach. Requiring the speedup to be greater than one produces the condition that Applying condition (49) to the linear case, where K = 1 and τ S D = τ I simply gives the condition N > 1, implying that a speedup is achievable for any number of processors. In the non-linear case, this no longer holds true, and N must satisfy (49) in order to be faster than the serial case. Once again, the algorithm has a theoretical maximum speedup; as N → ∞, s → (τ S D + τ † I )/(Kτ H + τ † H ). This is similar to the linear case, except for the presence of a factor K in the denominator which implies that as the number of iterations increases, the maximum speedup decreases.

Governing equation solved in series
Lastly, we consider the hybrid direct-serial-parallel-adjoint case introduced in section 3.3. The procedure is summarized in algorithm 3. To obtain the speedup, we initially consider the time needed to solve the direct equation together with the adjoint inhomogeneous equations. Due to the overlapping nature of the algorithm, and the way in which we choose the non-equidistant partitioning, all inhomogeneous parts oof the adjoint equation should finish at the same time. First, the direct equation is integrated in series giving a time of τ I T . Following this, the adjoint inhomogeneous equations must all be integrated for Estimate the value of k in (39) by performing short direct and adjoint homogeneous calculations.
Dividing this expression by the time it takes to solve the system in series (T S L from the linear discussion in section (3.4.1)) gives the reciprocal speedup as . (51) Whilst this theoretical speedup is considerably more complicated than the other cases, we can still consider the maximum speedup achievable. We note that with the time partitioning based on equation (40), we obtain T N −1 → T and (T N − T N −1 ) → 0 as N → ∞. Therefore, we have s → (τ I + τ † I )/(τ I + τ † H ), showing that effectively we achieve a speedup by replacing an inhomogeneous adjoint equation with a homogeneous one.

Checkpointing
At this point it is important to note that, owing to the addition of the forcing term in the adjoint formulation, there is an explicit dependence of the adjoint equations on the direct variables. Therefore, the direct variables need to be stored at all time steps during the forward sweep and injected, at the appropriate time, into the adjoint equations. For high-resolution cases and large time horizons, we cannot afford to store all necessary direct variables, as the memory requirements may exceed the necessary amount of RAM of our computational system. In this case a checkpointing scheme is needed.
A graphical representation of the checkpointing process, where four checkpoints t 1 , ..., t 4 are used, can be seen in figure 3. Initially, the direct equation In view of our parallel-in-time algorithms we incorporate a checkpointing scheme as follows. We begin by defining the chosen location of our checkpoints. Then, the direct or direct-adjoint loops between each checkpoints are carried out in parallel. For example, in the four-checkpoint example we would first carry out the direct-Rough integration by integrating in parallel (or serially in the hybrid case) the equation in the interval [t 0 , t 1 ], followed by the interval [t 1 , t 2 ] and so on, until the direct-Rough integration is finished. When this is accomplished, the direct-adjoint loops are performed using the parallel algorithms developed so far on the checkpointed intervals, preceding to the first time interval. Although the checkpoint times can be chosen optimally [see 11], we proceed by using equispaced checkpoints for ease of implementation.

Implementation details
As revealed by the scaling analysis, it is critically important that the homogeneous equations are solved faster than the inhomogeneous ones. Indeed, we can see from the theoretical scalings (45), (48) and (51) that the maximum speedup possible for each algorithm crucially depends on solving the homogeneous equation more efficiently. This implies a natural restriction on the kind of equations that we are able to parallelize in time with this approach. Utilizing the approach of Gander and Güttel [9], we identify stiff equations as good candidates for time parallelization. The homogeneous components that arise from these stiff problems can be solved efficiently and quickly with exponential time-steppers, an example of which was derived and applied to stiff equations by Cox and Matthews [3]. As exponential time-steppers are based on the exact integration of the homogeneous equation, via the matrix exponential or approximations thereof, they alleviate the stiffness of the problem and consequently are able to take larger time-steps than traditional time-stepping methods based on linear multistep or multistage formula. This section begins with a brief discussion of exponential integrators. Next, we demonstrate the application of each of the three algorithms developed to a direct-adjoint looping procedure. Special attention is given to how well the scaling of the algorithm agrees with the ones derived theoretically, and how the maximum speedup depends on the stiffness of the problem.

Exponential integrators
Exponential integration is based on the analytic solution to a linear, autonomous, homogeneous problem in terms of the matrix exponential. For the linear problemẏ = Ay with initial condition y(0) = y 0 , the solution at time ∆t is y = exp(∆tA)y 0 . Exponential integrators rely on approximating the matrixvector product exp(∆tA)v, which advances the state v by ∆t time units. Even though this approach is based on the exact integration for a linear equation, this approach can still be used for non-linear homogeneous systems, albeit with smaller time-steps. In either case, exponential time-steppers can perform better than alternative methods for stiff problems [see, e.g., 15].
One way to approximate the matrix exponential B = exp(∆tA)v is via the orthogonalized, order-m Krylov subspace of B [23]. By computing the Krylov subspace of the matrix B we obtain a lower-dimensional representation of the matrix, namely B ≈ V m H m V T m where H m is an m × m matrix. If m is much smaller than the dimension of B, the matrix exponential exp(B) ≈ V m exp(H m )V T m can be computed significantly more quickly and efficiently. If the Krylov representation is accurate for small m, then this method provides an efficient way of calculating the matrix exponential. However, if m needs to be large (e.g., if B has a large spectral radius) then this method becomes timeconsuming and memory-intensive, as m Krylov vectors must be calculated and stored.
An alternative to the Krylov-based method is based on Newton interpolation of the exponential function φ(z) = (exp(z) − 1)/z [2]. Once the interpolation is found, the matrix exponential exp(B) can easily be retrieved from φ(B). When computing the Newton interpolation the choice of interpolation points {ζ j } N j=0 is critically important, as a clever choice of ζ j reduces the amount of points (N + 1) needed for an accurate interpolation. One way to achieve superlinear convergence of the matrix polynomial to the matrix exponential φ(Z) is to use Chebyshev nodes on the real focal interval [a, b], computed such that 'the "minimal" ellipse of the confocal family that contains the spectrum (or the field of values) of the matrix is not too "large" ' [1]. However, the locations of Chebyshev nodes are dependent on the interpolation degree N + 1. Therefore, if the interpolation degree needs to be increased in order to achieve a user-specified error, then the entire interpolation must be recomputed. This cost can be circumvented by instead using real Leja points, which still achieve a superlinear convergence. The advantage of Leja points over Chebyshev points is that they do not depend on the degree of interpolation. Hence, if a smaller error is required, the degree can be increased without recomputing the interpolation. Newton interpolation of the exponential function φ(z) is called the real Leja points method (ReLPM) and is summarized by Bergamaschi et al. [1]. In the case of large sparse matrices the ReLPM can become more advantageous over Krylov-based methods due to its lower memory cost. We choose to use the ReLPM method for our implementation so that the developed code is more easily portable to future large-scale problems. For a more detailed discussion of exponential integrators see Skene [25], which also contains a preliminary version of this work.

Linear governing equation
To test the linear algorithm 1, we consider a forced advection-diffusion equation.
where the speed of advection is governed by a and the diffusion coefficient is D.
We take the domain to be Ω = [0, 2π] × [0, 2π], and apply periodic boundary conditions in both directions. By discretizing the derivatives using second-order central differences in space we can write this equation in the discretized form where the vectors q and f denote the state and forcing at each grid point, respectively. The state matrix A encapsulates the homogeneous component of the original system together with the boundary conditions. To formulate a direct-adjoint looping problem we consider choosing a force f = f true and evolving the system over the time interval [0, T ] starting from a zero initial condition. This results in the solution q true . We seek to recreate the forcing f true using only the observed solution, q true . To this end, let the cost functional be giving the adjoint equation where we have made the substitution τ = T − t to obtain an equation we integrate forward in the time variable τ . Once integrated, we can use the adjoint solution to evaluate the gradient of the cost functional with respect to the current forcing f as ∂L ∂f This gradient can then be used as part of a gradient-based optimization routine to minimize the cost functional and hence find f true . As mentioned previously, the algorithms developed perform better for stiffer systems. One way of introducing stiffness into our advection-diffusion problem is through the diffusion parameter D. This can be observed by studying the eigenvalues of A. Increasing D stretches the spectrum along the negative real axis and thus makes our system progressively more stiff. For our implementation we set a = 1, ω = 1, T = 10 and vary D ∈ {0.01, 0.1, 1, 10} and determine the average time to perform one direct-adjoint loop over 3 runs for a range of processor numbers. In this manner, we can observe to what extent the theoretical scalings agree with the scalings found from the numerical experiments.
The code is written in python, and a fourth-order Runge-Kutta scheme (with adaptive time-steps based on a fifth-order Runge-Kutta error approximation) is used for the inhomogeneous integrations. This time-stepper is implemented via the inbuilt RK45 method of the solve ivp function contained in the integrate submodule of the scipy library. For homogeneous integrations we use our implementation of the ReLPM as described by Bergamaschi et al. [1]. The relative error tolerance for each of these integrators is kept the same at rtol = 10 −3 .  Figure 4 shows the speedup obtained for the linear algorithm. As predicted, when D increases the speedup becomes larger, highlighting the observation that the homogeneous equation needs to be solved significantly faster than the inhomogeneous one. Furthermore, due to the extra cost in forming the matrix exponential which in turn causes the homogeneous integration to be slower than the inhomogeneous one, D = 0.01 and D = 0.1 report a worse performance for the parallel algorithm when compared to the serial approach. For D = 1 and D = 10 the stiffness introduced into the system is circumvented through the exponential time-stepping and thus, for the homogeneous equation, we are able to take far larger time steps than the inhomogeneous integrator, off-setting the cost of forming the matrix exponential.
The figure also shows good agreement between the numerically obtained scaling and the one predicted by the theoretical value. This shows that the effectiveness of a parallel-in-time approach for a specific problem can be assessed simply by performing a short simulation to determine the parameters τ I , τ † I , τ † H and , τ † H . Furthermore, by assessing the theoretical scaling the number of processors can best be chosen. In all cases, we can clearly see the scaling approaching its maximum value as the number of processors is increased. Initially, there are larger gains in speedup, followed by progressively smaller gains as more processors are included. For example, taking the D = 1 case for instance: using six processors we can expect a speedup of roughly 2.5, while increasing the number of processors only provides a speed up of approximately 3 -a relatively small speedup for more than double the number of processors. Hence, the user should initially consider the theoretical scaling and choose the number of processors that best balances speedup gains versus expended resources.

Non-linear governing equation
To test the two algorithms for non-linear direct equations we consider the harmonically forced viscous Burgers' equation for the state q = (U, V ) T on the same domain as the advection-diffusion equation of the previous section. Once again, we discretise using second-order centered finite differences and, together with periodic boundary conditions, arrive at The equations linearised about the direct solution q = (U, V) T are for the linearised state q = (u , v ) T . These equations are also discretized using second-order centered finite differences with periodic boundary conditions. This leads to the linearized equations in the form with the corresponding adjoint equation We have made explicit the dependence of the state matrix A on the direct solution and have, as before, introduced the time variable τ = T − t to integrate the adjoint equations forwards in the time variable τ.

Non-linear governing equation solved in parallel
Now we turn to timing the direct-adjoint loop, i.e., we solve both the nonlinear governing equation and the adjoint in parallel as described by algorithm 2. Figure 5 shows the comparison between the numerically obtained speedup and the speedup predicted by our theoretical scaling (48). As four iterations are used to converge the direct solution, we set K = 4 when computing the theoretical scaling. We note again that there is good agreement between the numerically obtained speedup and that predicted by our theoretical arguments. Initially the speedup is less than one due to the cost of performing four iterations negating any possible speedup. However, for N > 4 this cost is overcome, and we obtain a marked speedup.

Non-linear governing equation solved in series
Next, we time the same direct-adjoint loop, but using algorithm 3 in which the direct equation is solved in series. By performing short direct and adjoint simulations we obtain the value k = 2.11 for calculating the size of the time partitions using equation (40). The speedup obtained, compared to the theoretical value, is shown in figure 6.
The figure shows that the speedup converges to its theoretical maximum much faster than in the previous two algorithms. Indeed, about four processors are needed for the speedup to obtain its theoretical maximum. We also see a larger discrepancy in the obtained maximum efficiency gain and the one predicted by our scaling equation. This can mainly be attributed to the fact that the scaling argument assumes that, as the width of the time partition gets smaller, the inhomogeneous and homogeneous equations can be solved increas-

Checkpointing
Lastly, we turn our attention to the inclusion of a checkpointing scheme in performing a direct-adjoint loop. We consider the same problem as in section 4.3.2, but we include five checkpoints using the method outlined in section 3.5. The speedup obtained is shown in figure 7.
Similarly to the hybrid example with no checkpoints, we obtain the maximum speedup at about four processors. The maximum speedup has decreased, but this is to be expected for the hybrid algorithm as no speedup is gained in solving the forward equation by itself. Only direct-adjoint loops are accelerated. The fact that a speedup is still obtained with the inclusion of checkpoints demonstrates the compatibility of our developed algorithms with checkpointing schemes. For the cases of a linear governing equation, and a non-linear direct equation solved in parallel, we can expect less of a penalty for using a checkpointing scheme as, in this case, the forward equations are also accelerated.

Conclusions
In this article we have presented three separate algorithms that aim to extend the Paraexp algorithm of Gander and Güttel [9] to adjoint-looping studies. Theoretical speed-ups have been derived for all cases and verified through a set of numerical experiments. The verifications of the theoretical predictions indicate that studies utilizing these approaches can a-priori assess the achievable speedups as well as the necessary processors required. All algorithms sped up the calculation of the relevant system in line with the predictions we presented, signifying that the adjoint Paraexp algorithm, and its non-linear extension, can be used to accelerate direct-adjoint optimization studies. This is particularly pertinent due to the linear, time-varying nature of the adjoint which lends itself naturally to this parallelization approach, and thus can always be parallelized in time.
We first considered the extension of the Paraexp algorithm to direct-adjoint loops that stem from a linear governing equation. As the adjoint equation is also a linear equation this can be achieved simply by integrating the adjoint equation using the Paraexp algorithm. By keeping solutions on each time partition local to each processor we were able to minimize the processor communications required for direct-adjoint looping. The theoretical scaling of the algorithm was derived, showing that this algorithm has a theoretical maximum that can be obtained as the number of processors is increased. The linear algorithm was demonstrated on a two-dimensional advection-diffusion equation, where good agreement with theoretical scalings was obtained. These numerical experiments, as well as the theoretical scalings, highlight an important consideration that must be made when using the Paraexp algorithm, namely that the stiffness of the governing equations should be taken into account. Only when the homogeneous component can be solved significantly faster than the inhomogeneous one, can speedups be obtained. Keeping this caveat in mind, we expect this approach to be particularly applicable to systems governing by a stiff (multi-scale) dynamics.
Further to the linear case, we also considered how a parallel-in-time approach can be used for non-linear direct equations. This led to the development of two algorithms; an iterative algorithm in which the direct equation is solved in parallel, and a hybrid algorithm in which the direct equation is solved in series. In both cases the adjoint equation is integrated in parallel using the Paraexp algorithm. The scalings for the iterative approach showed that due to iterations being carried out, care must be taken to ascertain whether a speedup is possible for a given number of processors as the cost of iterations can readily outweigh the cost of solving the equations in parallel. Also, similarly to the linear case, there is a maximum speedup possible. Therefore, an initial analysis of the governing equations may be necessary to ascertain if this approach is applicable and viable. As in the iterative non-linear algorithm the iterations may not converge, or too many processors may be needed in order to obtain a speedup, we also proposed a hybrid approach in which the direct equation is solved in series whilst the adjoint is solved in parallel. By using a non-equidistant time-partitioning and overlapping the direct and adjoint solves, a speedup can be obtained. The time-partitioning is chosen so that all the homogeneous adjoint solves begin at the same time. In this way, the theoretical scaling shows that, as the number of processors increases, we obtain a maximum speedup where the inhomogeneous adjoint integration is effectively replaced by a homogeneous one. The scaling and subsequent numerical verification demonstrated that in this approach the maximum speedup, although lower than the other algorithms, is obtained with a smaller number of processors. Again, we note that for the non-linear algorithms the speedups rely on the homogeneous equations being able to be solved more efficiently than their inhomogeneous counterparts, for example as in the case of a stiff system.
Lastly, we demonstrated that all developed algorithms are easily implementable with checkpointing regimes. Hence, we believe that parallel-in-time approaches can provide a viable option for accelerating direct-adjoint loops stemming from both linear and non-linear governing equations, with the theoretical scalings providing an easy 'a-priori check' of the speedups available. We emphasise that this approach is particularly aimed at studies that have exhausted spatial parallelization but still have additional computational power available. In this case, a significant impact can be made by parallelizing in time, adding valuable efficiency to these optimization studies. The optimal overall balance of spatial versus temporal parallelism for a given number of processors and a given computer architecture is an interesting extension of our study and will be addressed in a future effort.