An Efficient Parallel-in-Time Method for Optimization with Parabolic PDEs

TL;DR

This paper presents a parallel-in-time method using PFASST to efficiently solve PDE-constrained optimization problems with parabolic equations, reducing computational costs and enabling faster solutions.

Contribution

It introduces a fully time-parallel algorithm applying PFASST to both state and adjoint equations, improving efficiency in PDE-constrained optimization.

Findings

Achieves significant parallel speedup in linear and nonlinear problems

Demonstrates efficiency gains from reusing previous iteration information

Validates approach with numerical experiments on reaction-diffusion problems

Abstract

To solve optimization problems with parabolic PDE constraints, often methods working on the reduced objective functional are used. They are computationally expensive due to the necessity of solving both the state equation and a backward-in-time adjoint equation to evaluate the reduced gradient in each iteration of the optimization method. In this study, we investigate the use of the parallel-in-time method PFASST in the setting of PDE constrained optimization. In order to develop an efficient fully time-parallel algorithm we discuss different options for applying PFASST to adjoint gradient computation, including the possibility of doing PFASST iterations on both the state and adjoint equations simultaneously. We also explore the additional gains in efficiency from reusing information from previous optimization iterations when solving each equation. Numerical results for both a linear and a non-linear reaction-diffusion optimal control problem demonstrate the parallel speedup and efficiency of different approaches.

Figures

Tables

Table 1. Table 1 : Strong scaling results for the heat example: first state, then adjoint approach with cold start. Speedup and parallel efficiency are compared to the sequential run with one processor.

#CPUs	time [s]		speedup	efficiency [%]
$1$	$7,808$	$.9$	–	$0.0$
$2$	$6,063$	$.5$	$1.3$	$65.0$
$5$	$3,797$	$.3$	$2.1$	$42.0$
$10$	$2,677$	$.8$	$2.9$	$29.0$
$20$	$1,679$	$.9$	$4.6$	$23.0$

Table 2. Table 2 : Strong scaling results for the heat example: mixed approach with cold start. Speedup computed with respect to the first state, then adjoint approach with one CPU.

#CPUs	time [s]		speedup	efficiency [%]
$1$	$8,888$	$.8$	–	$0.0$
$2$	$6,116$	$.0$	$1.3$	$65.0$
$5$	$3,518$	$.1$	$2.2$	$44.0$
$10$	$2,508$	$.1$	$3.1$	$31.0$
$20$	$1,606$	$.0$	$4.9$	$24.5$

Table 3. Table 3 : Strong scaling results for the heat example: first state, then adjoint with warm start. Total speedup computed with respect to first state, then adjoint approach with one processor and cold start. Speedup to cold denotes the speedup only due to warm start, i.e., compared to the time for cold start with the same number of CPUs. Efficiency is computed using the total speedup.

#CPUs	time [s]		total speedup	speedup to cold	efficiency [%]
$1$	$6,474$	$.1$	$1.2$	$1.2$	–
$2$	$5,087$	$.1$	$1.5$	$1.2$	$75.0$
$5$	$3,172$	$.1$	$2.5$	$1.2$	$50.0$
$10$	$2,034$	$.9$	$3.8$	$1.3$	$38.0$
$20$	$1,347$	$.5$	$5.8$	$1.2$	$29.0$

Table 4. Table 4 : Strong scaling results for the heat example: mixed approach with warm start. Total speedup computed with respect to the first state, then adjoint approach with one CPU and cold start. Speedup to cold denotes the speedup only due to warm start, i.e., compared to the time for the mixed approach using cold start with the same number of CPUs. Efficiency is computed using the total speedup.

#CPUs	time [s]		total speedup	speedup to cold	efficiency [%]
$1$	$7,187$	$.6$	$1.1$	$1.2$	–
$2$	$5,537$	$.4$	$1.4$	$1.1$	$70.0$
$5$	$3,424$	$.8$	$2.3$	$1$	$46.0$
$10$	$2,276$	$.3$	$3.4$	$1.1$	$34.0$
$20$	$1,494$	$.6$	$5.2$	$1.1$	$26.0$

Table 5. Table 5 : Influence of residual tolerance on parallel speedup for the first state then adjoint approach with cold start. Note that the sequential reference is computed using the same tolerances, otherwise total speedups are computed as in Tables  1 – 4 .

	speedup (20 CPUs)
residual tolerance	first state, then adjoint		mixed
	cold	warm	cold	warm
${10}^{-10}$	$4.9$	$5.8$	$4.9$	$5.2$
${10}^{-8}$	$4.9$	$5.9$	$4.6$	$5.3$
${10}^{-6}$	$4.5$	$5.5$	$4.3$	$4.9$
${10}^{-4}$	$4.1$	$5.0$	$3.9$	$4.4$

Table 6. Table 6 : IMEX vs. MISDC for the Nagumo example: required number of PFASST iterations (average per time step, rounded). For the nonlinear solver, Newton’s method is stopped when the discrete ℓ 2 subscript ℓ 2 \ell_{2} -norm of the correction is smaller than 10 − 12 superscript 10 12 10^{-12} .

$\gamma$	$T/\Delta t$	#CPUs	IMEX		MISDC lagged		MISDC nonlinear
			state	adjoint	state	adjoint	state	adjoint
1	32	1	-	-	14	9	14	9
		32	-	-	53	23	55	23
	64	1	7	6	7	6	7	6
		32	-	-	31	19	32	19
	128	1	5	4	5	4	5	4
		32	20	18	20	18	20	18
3	32	1	-	-	14	9	14	9
		32	-	-	53	24	54	24
	64	1	-	-	7	6	7	6
		32	-	-	31	19	32	19
	128	1	5	4	5	4	5	4
		32	20	18	21	18	21	18
5	32	1	-	-	14	9	14	9
		32	-	-	53	24	56	24
	64	1	-	-	7	6	7	6
		32	-	-	31	19	32	19
	128	1	5	4	5	4	5	4
		32	-	-	21	18	21	18
	256	1	4	4	4	4	4	4
		32	18	16	18	18	18	18

Table 7. Table 7 : Nagumo example: Strong scaling results for IMEX (first state, then adjoint; cold start). Speedup and parallel efficiency are compared to the sequential IMEX version with one CPU.

#CPUs	time [s]		speedup	efficiency [%]	rel. err. to $u_{\text{exact}}$
$1$	$169$	$.8$	–	$0.0$	$1.1\cdot {10}^{-1}$
$2$	$109$	$.0$	$1.6$	$80.0$	$1.3\cdot {10}^{-1}$
$4$	$73$	$.4$	$2.3$	$57.5$	$1.3\cdot {10}^{-1}$
$8$	$51$	$.0$	$3.3$	$41.3$	$1.1\cdot {10}^{-1}$
$16$	$37$	$.7$	$4.5$	$28.1$	$1.1\cdot {10}^{-1}$
$32$	$32$	$.0$	$5.3$	$16.6$	$1.2\cdot {10}^{-1}$

Table 8. Table 8 : Nagumo example: Strong scaling results for IMEX (first state, then adjoint; warm start). Total speedup is compared to IMEX with one CPU and cold start, speedup to cold denotes the speedup only due to warm start, i.e., compared to the time for cold start with the same number of CPUs. Efficiency is computed using the total speedup.

#CPUs	time [s]		total speedup	speedup to cold	efficiency [%]	rel. err. to $u_{\text{exact}}$
$1$	$132$	$.7$	$1.3$	$1.3$	$130.0$	$1.3\cdot {10}^{-1}$
$2$	$84$	$.6$	$2$	$1.3$	$100.0$	$1.3\cdot {10}^{-1}$
$4$	$58$	$.0$	$2.9$	$1.3$	$72.5$	$1.2\cdot {10}^{-1}$
$8$	$39$	$.6$	$4.3$	$1.3$	$53.8$	$1.2\cdot {10}^{-1}$
$16$	$30$	$.4$	$5.6$	$1.2$	$35.0$	$1.2\cdot {10}^{-1}$
$32$	$22$	$.1$	$7.7$	$1.4$	$24.1$	$1.2\cdot {10}^{-1}$

Table 9. Table 9 : Runtimes of different sequential methods yielding similar optimization results.

method	$\Delta t$	time [s]	rel. err. to $u_{\text{exact}}$
IMEX-Euler	${10}^{-3}$	$102.5$	$1.2\cdot {10}^{-1}$
ARK-4	$2.5\cdot {10}^{-3}$	$183.0$	$1.4\cdot {10}^{-1}$
MLSDC cold	$1.5625\cdot {10}^{-1}$	$169.8$	$1.1\cdot {10}^{-1}$
MLSDC warm	$1.5625\cdot {10}^{-1}$	$132.7$	$1.3\cdot {10}^{-1}$

Equations

$$\underset{y\in Y,u\in U}{\min }J(y,u)\text{subject to}c(y,u)=0,$$

$$c_y(y,u)y^{\prime }(u)+c_u(y,u)=0.$$

$$\underset{u\in U}{\min }j(u):=J(y(u),u).$$

$$⟨j^{\prime }(u^{\star }),u-u^{\star }{⟩}_{U^{\star },U}\ge 0\forall u\in U_{\text{ad}}.$$

$$\mathcal{L}(y,u,p)=J(y,u)+⟨p,c(y,u){⟩}_{Z,Z^{\star }},$$

$$⟨j^{\prime }(u),\delta u{⟩}_{U^{\star },U}=⟨{\mathcal{L}}_y(y(u),u,p),y^{\prime }(u)\delta u{⟩}_{Y^{\star },Y}+⟨{\mathcal{L}}_u(y(u),u,p),\delta u{⟩}_{U^{\star },U}.$$

$$c_y(y(u),u{)}^{\star }p(u)=-J_y(y(u),u)$$

$${\mathcal{L}}_y(y(u),u,p(u))=J_y(y(u),u)+c_y(y(u),u{)}^{\star }p(u)=0.$$

$$⟨j^{\prime }(u),\delta u{⟩}_{U^{\star },U}=⟨{\mathcal{L}}_u(y(u),u,p),\delta u{⟩}_{U^{\star },U},$$

$$j^{\prime }(u)={\mathcal{L}}_u(y(u),u,p(u))=J_u(y(u),u)+c_u(y(u),u{)}^{\star }p(u).$$

$$\bigl{(}\delta u,\nabla j(u)\bigr{)}_{U}=j^{\prime}(u)\delta u\ \ \forall\delta u\in U.$$

$$J(y,u)=\underset{=J^{\Omega }(y,u)}{\underbrace{\frac{1}{2}{\int }_0^T{\Vert y-y_d\Vert }_{L^2(\Omega )}^{2}dt}}+\underset{=J^u(y,u)}{\underbrace{\frac{\lambda }{2}{\int }_0^T{\Vert u\Vert }_{L^2(\Omega )}dt}}$$

$$y_t-\kappa \Delta y+f(y)-u$$

$$y(0)-y_0$$

$$-p_t-\kappa \Delta p+f_y(\bar{y})p$$

$$p(T)$$

$$J^T(y,u)=\frac{\sigma }{2}{\Vert y(T)-y_d^T\Vert }_{L^2(\Omega )}^2$$

$$p(T)=-J_y^T(\bar{y},u)=-(y(T)-y_d^T).$$

$$u_{k+1}$$

$$d_{k+1}$$

$${\beta }_k^{\text{FR}}$$

$${\beta }_k^{\text{PRP}}$$

$${\beta }_k^{\text{DY}}$$

$$j(u_k+{\alpha }_k{d}_k)$$

$$|⟨j^{\prime }(u_k+{\alpha }_k{d}_k),d_k⟩|$$

$$y^{\prime }(t)=F(t,y),$$

$$y(t_m)=y_j+{\int }_{t_j}^{t_m}F(\tau ,y(\tau ))d\tau .$$

$$y_m\approx y_j+\Delta t\sum _{i=0}^M{w}_{m,i}F(t_i,y_i),$$

$$w_{m,i}=\frac{1}{\Delta t}{\int }_{t_j}^{t_m}{\ell }_i(s)ds,m=0,\dots ,M,i=0,\dots ,M,$$

$$y_m^{[k+1]}=y_j+\Delta t\sum _{i=0}^m{\tilde{w}}_{m,i}F(t_i,y_i^{[k+1]})+\Delta t\sum _{i=0}^M(w_{m,i}-{\tilde{w}}_{m,i})F(t_i,y_i^{[k]}).$$

$$(v,w)={\int }_0^T(v(t),w(t){)}_{H}dt.$$

$$(v,w)=\sum _{i=0}^N{\int }_{t_i}^{t_{i+1}}(v(t),w(t){)}_H.$$

$$-{\tilde{p}}_t-\kappa \Delta \tilde{p}+f_y(\bar{y})\tilde{p}$$

Full text

Keywords: PDE-constrained optimization, parallel-in-time methods, PFASST

AMS subject classifications: 65K10, 65M55, 65M70, 65Y05

1 Introduction

Large-scale PDE-constrained optimization problems occur in a multitude of applications, for example in solving inverse problems for non-destructive testing of materials and structures [23] or in individualized medicine [15]. More recently, the training of certain deep neural networks in machine learning, e.g, for image recognition or natural language processing, has been formulated as a dynamic optimal control problem [27, 41]. Algorithms for the solution of PDE-constrained optimization problems are computationally extremely demanding as they require to numerically solve multiple partial differential equations during the iterative optimization process. This is especially challenging for transient problems, where the solution of the associated optimality system requires information about the discretized variables on the whole space-time domain. For the solution of such optimization problems, methods working on the reduced objective functional are often employed to avoid a full spatio-temporal discretization of the problem. The evaluation of the reduced gradient then requires one solve of the state equation forward in time, and one backward-in-time solve of the adjoint equation. In order to tackle realistic applications, it is not only essential to devise efficient discretization schemes for optimization, but also to use advanced techniques to exploit computer architectures and decrease the time-to-solution, which otherwise is prohibitively long. One approach is to utilize the number of CPU cores of current and future many-core high performance computing systems by parallelizing the PDE solution method. In addition to well-established methods for parallelization in the spatial degrees of freedom, parallel-in-time methods have seen a growing interest in the last 15 years. Research into time parallelism dates back at least to the 1960s, and in 2001 the introduction of parareal by Lions et al. [34] sparked new research into time-parallel methods. The methods in this paper are based on the parallel full approximation scheme in space and time (PFASST) introduced by Emmett and Minion [13]. We will not attempt a thorough review of the field here, and the interested reader is encouraged to consult the survey article [17] for an overview of competing approaches.

More recently, the application of space-time parallel methods to the solution of optimization problems governed by PDEs has become an active research area, with approaches including multiple shooting (e.g., [29] and the references therein), Schwarz methods [5, 19], and the application of parareal preconditioners [35, 44]. A time-parallel gradient type method is presented in [10]. There the time interval of interest is subdivided into time steps, which are solved in parallel using quantities from the previous optimization iteration as input. This leads to jumps in the solutions of state and adjoint equation such that these equations are not satisfied during optimization. While they report excellent speedups and linear scaling up to 50 processors and show convergence if sufficiently small step sizes for updating the control are used, it is unclear how to automatically select such a step size. Alternatively, space-time parallel multigrid methods are applied to adjoint gradient computation and simultaneous optimization [24, 25] within the XBraid software library [4]. XBraid provides a non-intrusive framework adding time-parallelism to existing serial time stepping codes, and using simultaneous instead of reduced space optimization, a speedup of 19 using 256 time processors has been reported. The same method is also applied to perform layer-parallel training of neural networks [26].

In this paper we employ PFASST to provide a fully time-parallel reduced-space gradient- or nonlinear conjugate gradient method that allows using the usual line search criteria, e.g., the strong Wolfe conditions, for step size selection to guarantee convergence, and is thus non-intrusive with respect to the optimization algorithm. On the other hand this is an intrusive approach concerning the PDE solvers as it requires using multilevel spectral deferred correction methods as time steppers. The implementation effort is mitigated by the availability of libraries like LibPFASST [2] or dune-PFASST [1]. While the basic ideas of the approach are outlined in the short paper [20], here we provide more details, and develop additional approaches to increase speedup and efficiency. We demonstrate that using PFASST brings additional benefits, like flexibility in treating nonlinearities, and, more importantly, the possibility of warm starting the PDE-solutions required during the optimization. The remainder of the paper is organized as follows. In Sect. 2 we present the optimization problem and review optimality conditions as well as adjoint gradient computation. The PFASST method is introduced in Sect. 3; it is used to derive parallel-in-time methods for solving optimization problems with parabolic PDEs in Sect. 4. Finally, in Sect. 5 we present numerical examples, followed by a discussion of results and future improvements in Sect. 6.

2 Adjoint gradient computation for optimization with parabolic PDEs

Here we briefly summarize the mathematical approach to parabolic PDE-constrained optimization problems. For more details and generalizations, we refer to, e.g., [31].

We consider optimization problems of the form

[TABLE]

with the equality constraint $c:Y\times U\rightarrow Z^{\star}$ being a parabolic PDE on Hilbert spaces $Y,U,Z$. $Z^{\star},U^{\star}$ denote the dual spaces of $Z$ and $U$, respectively; the dual pairing between a space $X$ and its dual is denoted by $\langle\cdot,\cdot\rangle_{X^{\star},X}$, and $(\cdot,\cdot)_{X}$ denotes the scalar product on $X$. We drop the subscripts like $\cdot_{X}$ if the involved spaces are clear from the context. In the present setting, the constraint $c(y,u)=0$ means that the state $y$ satisfies a PDE where the control $u$ is a specific forcing term, occurring, e.g., as a source term, in the boundary conditions, or as some other parameter.

To derive optimality conditions, we assume that there exists a unique solution $y=y(u)\in Y$ of the state equation $c(y,u)=0$ for each control $u\in U$. We additionally assume that $c_{y}(y,u):Y\rightarrow Z^{\star}$ is continuously invertible. Then, by the implicit function theorem (see, e.g., [46, Section 4.7]), the control-to-state mapping is continuously differentiable, and the derivative $y^{\prime}(u)$ is given by the solution of

[TABLE]

Subscripts like $c_{u}()$ denote the partial derivatives with respect to the indicated variable. By inserting $y(u)$ into the optimization problem (1) we arrive at the reduced problem

[TABLE]

In this unconstrained setting, the following simple first-order necessary optimality condition holds. If $u^{\star}\in U$ is a local solution of the reduced problem (3) it is a zero of the reduced derivative, $j^{\prime}(u^{\star})=0.$ If the reduced functional $j$ is convex, this condition is also sufficient. If we allow control constraints, i.e., demand $u\in U_{\text{ad}}$ with $U_{\text{ad}}\subset U$ non-empty, convex and closed, the optimality condition changes to the variational inequality for the local minimizer $u^{\star}\in U_{\text{ad}}$

[TABLE]

To formally derive a representation for the reduced gradient, we define the Lagrange functional $\mathcal{L}:Y\times U\times Z\rightarrow\mathbb{R}$,

[TABLE]

where in the present context, the Lagrange multiplier $p\in Z$ is referred to as the adjoint. Clearly, inserting $y=y(u)$ into (5), we get $j(u)=\mathcal{L}(y(u),u,p)$ for arbitrary $p\in Z$. Differentiation in direction $\delta u\in U$ yields

[TABLE]

Choosing $p=p(u)$ such that the adjoint equation

[TABLE]

is fulfilled gives

[TABLE]

Inserting this into (6), the first term on the right hand side vanishes, and we get the reduced derivative $j^{\prime}(u)\in U^{\star}$ as

[TABLE]

i.e.,

[TABLE]

In the Hilbert space setting used here, for a given $u\in U$ the reduced gradient $\nabla j(u)\in U$ is then given as the Riesz representative of the reduced derivative $j^{\prime}(u)\in U^{\star}$, i.e., via,

[TABLE]

To be more concrete, we consider a distributed tracking-type objective functional $J(y,u)$ with an additional term penalizing the control cost for $\lambda\geq 0$,

[TABLE]

with $y_{d}$ representing a desired solution. The state $y$ is subject to a linear or semi-linear parabolic PDE with distributed control

[TABLE]

and suitable boundary conditions. For simplicity, we consider a constant scalar diffusivity $\kappa$; $\Delta y$ denotes the Laplacian of $y$. Thus, for sufficiently smooth $f$, we choose $Y=U=Z=L^{2}(0,T;H)$ with $H=L^{2}(\Omega)$ in (1). Here, $\Omega\subset\mathbb{R}^{n}$, and $T>0$ denotes a given final time. The adjoint equation (7) for a given state solution $\bar{y}$ becomes

[TABLE]

with homogeneous boundary conditions of the same type as in the state equation. Generalizations to controls or observations on the boundary are straightforward. Having a terminal cost contribution in the objective, e.g.,

[TABLE]

leads to the modified terminal condition

[TABLE]

As the adjoint equation (7) is backward in time, due to the occurrence of $-J_{y}(y(u),u)$ as a source term, and—in the nonlinear case—the dependence of $c_{y}(y(u),u)$ on the state solution $y(u)$, adjoint gradient computation consists of three steps (see also Fig. 1):

Solve the PDE $c(y,u)=0$ for a give control $u$, and store the solution trajectory $y\in Y$. 2. 2.

Solve the adjoint PDE $c_{y}(y,u)^{\star}p=-J_{y}(y,u)$ for $p\in Z$. 3. 3.

Set $j^{\prime}(u)=J_{u}(y,u)+c_{u}(y,u)^{\star}p$ and compute the Riesz representation $\nabla j(u)$.

Thus, the computation of the reduced gradient requires the solution of two parabolic PDEs. For solving the optimization problem, in this work we consider nonlinear conjugate gradient (ncg)- and steepest descent (sd) methods, as they require only gradient information. Hence gradient-based methods for the optimization problems that iterate over these three steps require many forward and backward PDE solves. The goal of this paper is to exploit properties of the the PFASST parallel-in-time algorithm to reduce the computational complexity of the optimization procedure.

Remark 2.1**.**

State and adjoint equations can be solved as a coupled system, as is done, e.g., in [22]. In view of parallelizing in time this becomes more difficult, see the brief discussion in Sect. 4.2 and Fig. 4; our main focus thus is on gradient computation and using algorithms like steepest descent or nonlinear conjugate gradients to iteratively solve the optimization problem. This can easily be applied to, e.g., control constrained problems, and be extended to second order methods like semi-smooth Newton, or Newton-CG.

For the gradient-based methods considered here, the optimization iteration proceeds as

[TABLE]

where $d_{0}=\nabla j(u_{0})$, and the choice of $\beta_{k}$ defines the actual method. Abbreviating $\nabla j(u_{k})$ by $g_{k}$, specific variants of ncg include, e.g.,

[TABLE]

Using $\beta_{k}=0$ yields the usual steepest descent method.

To guarantee convergence under the usual assumptions, the step size $\alpha_{k}$ is chosen to satisfy the strong Wolfe conditions for ncg,

[TABLE]

$0<c_{1}<c_{2}<1$, and just the Armijo condition (18) for sd [38].

Remark 2.2**.**

In this work we follow the first optimize, then discretize approach. To implement the optimization methods based on the optimality conditions of the previous section, we need to discretize the arising parabolic PDEs in time and space. To facilitate using PFASST for time parallelism, this is done using the method of lines approach, so we discretize space first. In the examples in Sect. 5 we use a pseudo-spectral method, but other techniques like finite element, finite difference, or finite volume methods are possible as well. The resulting system of ODEs is then solved using PFASST.

3 SDC, MLSDC, and PFASST

In this section we give an overview of the parallel-in-time strategy used to solve the state and adjoint equations. The strategy is mainly based on the PFASST algorithm [13] with some modifications specific to PDE-constrained optimization problems. PFASST can be thought of as a time parallel variant of the multi-level spectral deferred correction method (MLSDC) [42], which in turn is constructed from the method of spectral deferred corrections (SDC) [12]. Since all three of these methods are well established, we give only a brief overview here skewed towards the numerical methods tested in Sect. 5.

3.1 SDC Methods

Consider the generic ODE over the time interval $[t_{j},t_{j+1}]$ representing one time step

[TABLE]

with initial condition $y(t_{j})=y_{j}$. Divide the interval $[t_{j},t_{j+1}]$ into $M$ smaller sub-intervals by choosing points $t_{m},\leavevmode\nobreak\ m=0,\dots,M$ corresponding to the Gauss-Lobatto quadrature nodes.

The exact solution of the ODE at each point $t_{m}$ is given by

[TABLE]

The collocation method is defined by the solution of the system of equations derived by applying quadrature rules to (21),

[TABLE]

where $\Delta t=t_{j+1}-t_{j}$ and $w_{m,i}$ are the quadrature weights,

[TABLE]

with Lagrange polynomials $\ell_{i}$ defined by the quadrature nodes $(t_{m})_{m=0,\dots,M}$ (here Lobatto IIIA). Equation (22) is a (typically nonlinear) system that is $M$ times larger than that of a single-step implicit method like backward Euler. It is also equivalent to a fully implicit Runge-Kutta method with the values $w_{m,i}$ corresponding to the matrix in the Butcher tableaux. Such methods have good stability properties and have formal order of accuracy $2M$ for $M+1$ Lobatto quadrature nodes. (See e.g. [28] for a more detailed discussion of collocation and implicit Runge-Kutta methods).

Spectral deferred correction methods (SDC) can be considered as a fixed point iteration to solve the collocation formulation (22). Each SDC iteration or sweep consists of stepping through the nodes $t_{m}$ and updating the solution at that node $y_{m}^{{[k]}}$. A typical version of SDC using implicit sub-stepping takes the form (with ${{[k]}}$ denoting iteration)

[TABLE]

Here the values ${\tilde{w}}_{m,i}=0$ for $i>m$, hence each of the substeps has the same computational complexity of a single backward Euler step.

One attractive aspect of SDC methods is that they can easily be extended to cases in which (20) can be split into stiff and nonstiff components. Semi-implicit or Implicit-Explicit (IMEX) methods that split the equation into stiff and nonstiff terms first appeared in [37]. So called multi-implicit or MISDC methods that treat two implicit terms in an operator splitting approach are introduced in [7]. Such splittings are explored in the numerical tests in Sect. 5.2 by using MISDC variants to treat the nonlinearity of the state equation. In the numerical examples, the values of ${\tilde{w}}_{m,i}$ for implicit terms are chosen following the $LU$ factorization method of Weiser [45], while those for explicit terms correspond to the usual forward Euler substepping.

3.2 MLSDC Methods

Higher-order SDC methods require a relatively large number of function evaluations (explicit and/or implicit) per time step. One method to reduce the cost of these iterations is to employ a multi-level formulation of iterations where SDC sweeps are done on a hierarchy of discretization levels. In [42], so called multi-level SDC (MLSDC) methods are studied where the levels are differentiated by the spatial and/or temporal order and resolution as well as the tolerance of implicit solvers. SDC sweeps are scheduled like V-cycles in multigrid, and coarse level problems are modified with a term analogous to a space-time full approximation scheme (FAS) correction term. In the current study, only the number of spatial degrees of freedom and the number of SDC quadrature nodes is varied on MLSDC or PFASST levels. MLSDC is the basic building block for the parallel-in-time method PFASST discussed next.

3.3 PFASST

The parallel full approximation scheme in space and time (PFASST) [13] is, as the name suggests, a method for exploiting both spatial and temporal parallelism in a manner similar to FAS multigrid methods for nonlinear problems (see [6] for a multigrid perspective of PFASST). As mentioned above, PFASST can be considered to be a pipelined version of the MLSDC method, with each time step being assigned to a separate processor (or groups of processors), such that MLSDC sweeps are performed on multiple time steps in parallel. “Pipelined” here refers to the idea that each processor begins a step of the algorithm as soon as new initial conditions are passed forward in time from the previous processor [36].

As all processors handling time steps $[t_{j},t_{j+1}]$ need an initial condition $y_{j}$, which for $j>0$ is not known, PFASST starts with a predictor phase. Usually, this is done by integrating the equations at the coarsest level in serial (with or without the FAS correction term) with a low-order method. Alternatively, this step can be conducted without coarse level communication since after distributing the initial value $y_{0}$ to all processors, the $j$th processor can compute $j$ time steps sequentially rather than wait for the first $j-1$ time steps to be completed before beginning (burn-in). Note, however, in the context of PDE-constrained optimization, the latter is not possible, as each processor only has data (control, desired state) for the time steps it actually computes, and not for the other time intervals. Instead, except for the first iteration of optimization method, one can initialize all processors using the solution from the previous optimization iteration. We refer to this procedure as warm starting PFASST iterations. In this case, each processor will start doing SDC sweeps with communication at the coarse level propagating the solution forward in time to receive and updated initial value. The predictor step is represented at the bottom of Fig. 2, with warm starting discussed further in Sect. 4.3.

Having finished the predictor phase, each processor performs MLSDC sweeps, where after each SDC sweep on each level of the MLSDC hierarchy, the solution is communicated forward in time to the next processor. This communication overlaps with computation except at the coarsest level (see Fig. 2 and [14]). A more detailed description can be found in [13].

4 Time-parallel PDE-constrained Optimization

In this section, we discuss three components to produce an efficient fully time-parallel method for PDE-constrained optimization. First we briefly describe parallelizing the outer optimization loop, before coming to the more involved parallel-in-time computation of the reduced gradient. Finally, as the time integrators for solving state and adjoint equations are iterative, we discuss using previous solutions to warm start the time integration.

For parallelization in time, the time domain $[0,T]$ is subdivided into $N$ time steps $0=t_{0}<\dots<t_{N}=T$, which are distributed on $R$ processors. For examples with relatively fewer time steps, the shortest computation times may result from choosing $R=N$. On the other hand, when the number of processors used is smaller than the total number of time steps, (e.g., in the strong scaling studies included below), the PFASST algorithm is applied sequentially to blocks of time steps where $N$ is an integer multiple of $R$. Note that here we consider only temporal parallelism. Usually, this is used in combination with spatial parallelism, and multiplies the speedup gained from spatial parallelism. This means that the $R$ processors used for parallelizing in time can be considered as groups of processors with additional spatial parallelism. The best strategy for distributing processors between time and space parallelism is in general problem- and machine-dependent and not considered here.

4.1 Parallelization of the optimization loop

Time-parallelization of the optimization loop (16), (17) requires parallel-in-time computation of the reduced gradient, as well as evaluation of inner products. In the present optimal control setting for the prototype problem (10)–(14) with the usual $L^{2}$-regularity in time, $(\cdot,\cdot)$ denotes the inner product in $L^{2}(0,T;H)$, so

[TABLE]

The evaluation of these inner products, and thus computation of $\beta_{k}$ as well as evaluation of the objective function $j(u_{k}),j(u_{k}+\alpha_{k}d_{k})$ can easily be done time-parallel, as

[TABLE]

After discretization, each processor evaluates the discrete spatial inner product and integrates in time only for its own time steps. In the numerical examples below, evaluation of the time integral is done using a simple trapezoidal rule; of course using other quadrature rules, like re-using the spectral quadrature matrices provided by PFASST, is possible as well. In terms of communication, for each inner product only one scalar value per processor has to be transmitted to one master processor collecting the results. The essential ingredient for a fully parallel optimization method is the computation of the reduced gradient $\nabla j(u_{k})$, which is discussed next.

4.2 Time-parallel adjoint gradient computation

Here we discuss three different ways in which PFASST can be used to enable a time-parallel computation of the reduced gradient $\nabla j(u_{k})$. In the simplest case, PFASST is used to first solve the state equation, then to solve the adjoint equation afterwards (Sect. 4.2.1). Alternatively, if $R=N$, the adjoint can be solved simultaneously with the state equation (Sect. 4.2.2). A third variant is inspired by the paraexp method [18]. In this variant, we make use of the linearity of the adjoint equation to split the adjoint solve into an inhomogeneous equation with homogeneous terminal conditions on each time step (without communication), and a subsequent propagation of the correct terminal conditions backward-in-time (Sect. 4.2.3).

4.2.1 First state, then adjoint

In this straightforward approach, PFASST is used first to solve the state equation (11), (12). The state solution is recorded at the quadrature nodes, and then used in a subsequent PFASST run to solve the adjoint equation (13), (14), see Fig. 3 for a sketch. This has already been used in the preliminary study [20].

Compared to sequential time stepping, the storage overhead is larger here, as the discretized-in-space state solution has to be stored on quadrature nodes, not only on the time steps. This is somewhat mitigated by the fact that the time intervals are distributed across several processors, typically giving access to more memory. Using PFASST gives some flexibility in reducing memory requirements. Besides compressed storage, e.g., [21] for adaptive lossy compression of finite element solutions, options include storing the state on the coarse level only, and reducing the compression error by additional state sweeps without communication for the adjoint solve. A detailed analysis of storage requirements and storage reduction techniques will be reported elsewhere.

4.2.2 Simultaneous approach

Alternatively, if $R=N$, the adjoint can be solved simultaneously with the state equation (Fig. 4), requiring communication of updated initial values for the state equation forward-in-time as well as updated terminal conditions for the adjoint backward-in-time. This cross-communication makes an efficient implementation with overlapping communication and computation difficult (see Fig. 5). As our numerical experiments show that this induces severe wait times, thus rendering the method inefficient, we do not consider this approach further.

4.2.3 Mixed approach

Making use of the linearity of the adjoint equation, the adjoint $p$ described in equation (13), (15) can split into two terms $p=\tilde{p}+\delta$, where $\tilde{p}$ satisfies the same equation as $p$ but with homogeneous terminal conditions, and the defect $\delta$ which corrects the terminal conditions. Specifically, for a given state solution $\bar{y}$, let

[TABLE]

and

[TABLE]

To allow for an efficient time-parallel solution, we have to further modify the equations on the time intervals. Denote by the superscript $j$, e.g., $\tilde{p}^{j},\delta^{j}$, the respective solution on the $j$th time interval $[t_{j},t_{j+1}]$. With this notation, let $\tilde{p}^{j}$ solve

[TABLE]

Note the homogeneous boundary condition on the interval. For the defect $\delta^{j}$, the respective equation then is

[TABLE]

with terminal conditions

[TABLE]

With this, equations (30), (31) can be solved in parallel together with the state equation on the interval without requiring communication backwards-in-time (first step in Fig. 6, the dash-dotted line denotes the modified adjoint equation). In the second step, the defect equation is solved (dashed line in Fig. 6). Note that on each time interval, the adjoint solvers only need to access the state solution computed on the same interval (even if each processor computes multiple time steps), thus requiring no additional communication of state values.

Remark 4.1**.**

In case of no terminal cost in the objective function, the splitting of the adjoint is still required, since the homogeneous terminal condition only holds at time $t=T$. After time discretization, the terminal condition of the timestep, i.e., $\tilde{p}^{j}(t_{j+1})$, is in general non-zero. To avoid backward-in-time communication these terminal conditions on the respective intervals are treated by the defect equation (32), (33).

For the solution of (32), (33) we note that formally, defining the operator $L=-\kappa\Delta+f_{y}(\bar{y}),$ the solution is given by

[TABLE]

For linear state equations we have $f_{y}=0$, and equation (34) can be solved efficiently using suitable approximations of the matrix exponential, see the discussion in [18]. Thus, for the numerical experiments we restrict ourselves to the linear case. Besides having a certain relevance on their own, such equations occur, e.g., as subproblems in an SQP method, or in the evaluation of Hessian-times-vector products for Newton-CG methods, where an additional linearized state equation and a corresponding second adjoint equation need to be solved. For nonlinear state equations, the differential operator and/or coefficients of the adjoint are time-dependent (as they depend on the state solution). The extension of the mixed approach to this case (e.g., facilitating the Magnus expansion along the lines of [33], where temporal parallelism of Magnus integrators is explored) is left for future work.

Remark 4.2**.**

The mixed approach is to some extent related to the paraexp algorithm [18] for linear initial value problems. While the treatment of the source term on time intervals $[t_{j},t_{j+1}]$ is the same, the paraexp method uses a near-optimal exponential integrator to compute in parallel solutions with the correct initial value on intervals $[t_{j},T],j=0,\dots,N-1$, i.e., up to the final time $T$, and then sums up the individual solutions, requiring communication between several processors for computing the superposition at the required time points. In contrast, we solve the defect equation to propagate the correct terminal values on the time interval $[t_{j},t_{j+1}],j=0,\dots,N-1$ only, such that only solution values at $t_{j},\,j=N,\dots,1$ have to be communicated to the adjacent processor.

4.3 Warm starts

An important benefit of using SDC sweeps is the option to re-use information from the previous optimization iteration to initialize the subsequent PDE solves. For this, the solution of state and/or adjoint at the quadrature nodes is stored in one optimization iteration, and used as an initial guess instead of the predictor step during the following state/adjoint solve with an updated control. As the optimization is converging, $\left\|u_{k+1}-u_{k}\right\|$ gets smaller, so that the previous solution of the state and adjoint equations are a suitable initialization, and the PFASST algorithm will take fewer iterations to converge. This incurs an overhead in storage that can be reduced by either using compression, or, as it is only used as initialization, values can be stored only on the coarsest level, and interpolated.

For the numerical examples below, to sweep on state or adjoint in optimization iteration $k+1$ we load the stored fine-level solution of the previous iteration $k$ to initialize at the fine quadrature nodes. Instead of invoking the usual PFASST predictor, we restrict the solution values and function evaluations to the coarser levels, and perform sweeps on the coarse level with communication as described in Sect. 3.3.

5 Numerical Results

In this section, we test the developed methods on two optimal control problems. First, we consider optimal control for a linear heat equation in Sect. 5.1 to compare the first state, then adjoint approach and the mixed approach for cold and warm starts. As a second example, an optimal control problem governed by the nonlinear Nagumo reaction-diffusion equation is considered (Sect. 5.2). There we compare different methods of treating the nonlinearity and again consider cold and warm starting of the time integration, using the first state, then adjoint approach. All examples are implemented using the LibPFASST library [2], and are timed on Intel Xeon E5-4640 CPUs clocked at 2.4GHz.

Before discussing the results in detail, two remarks are in order. First, we report speedup and parallel efficiency with respect to the sequential versions of the described methods, i.e., MLSDC for state and adjoint solution. In particular, the sequential integration method has the same temporal order as the parallel version. We do not perform a thorough comparison of our methods with other commonly used sequential temporal integrators, as it is not clear what the “best” sequential method is for the respective problems. However, in Sect. 5.2, we briefly report results using IMEX-Euler as well as a $4$th-order additive Runge-Kutta method as sequential references. Second, scaling experiments are performed for a relatively low number of processors (20 and 32 in the two examples). Usually, parallel-in-time methods are combined with parallelization in the spatial domain, and multiply the speedup (and the required number of processors). As a rule of thumb, time-parallel methods come into play when spatial parallelism saturates, i.e., most of the available processors will be used for spatial parallelism.

5.1 Linear problem: Heat equation

For the first numerical test, we consider the linear-quadratic optimal control problem to minimize $J(y,u)$ given by (10) subject to

[TABLE]

with periodic boundary conditions. The spatial domain $\Omega=[0,1]^{3}$, with initial conditions

[TABLE]

and target solution

[TABLE]

This implies the exact solution

[TABLE]

For strong scaling results we solve the above problem, with control cost parameter $\lambda=0.05$ in the objective (10). We use a three level PFASST scheme, and a pseudo-spectral discretization in space with $(16/32/64)^{3}$ degrees of freedom on the respective levels. The implicit linear solves are done via the FFT. In time, we use $T=2$ and $20$ time steps with $2/3/5$ Lobatto IIIA quadrature rules, yielding a temporal order $8$. To solve the optimization problem we apply gradient descent with the Armijo step size rule, and stop the optimization after $50$ iterations to compare the optimization progress, computed controls, and timings of the different algorithmic variants (first state, then adjoint vs. mixed, cold vs. warm start). PFASST iterations are stopped when the absolute or relative residual drops below $10^{-10}$. We note that this residual tolerance is too strict for the initial optimization iterations, where less exact gradients would be sufficient. Future research will be concerned with a thorough analysis of accuracy requirements, and the influence of inexactness on the convergence of the optimization methods along the lines of [21]. A brief comparison of different residual tolerances is shown in Table 5. We note that while naturally the wallclock time decreases, the influence of increased residual tolerances on parallel efficiency is small, as both sequential and parallel versions are accelerated similarly.

Results are obtained using the first state, then adjoint approach as well as the mixed approach. For each variant we compare the usual PFASST predictor (cold start) to initialize the MLSDC sweepers on each time step, and warm starts, using the previously computed solution on the fine level to initialize the sweeper. Tables 1–4 show speedup and efficiency for strong scaling using up to $20$ processors for the different approaches. In Tables 1 and 2, the usual PFASST predictor is used, spreading the given initial value on each parallel time step to the quadrature nodes. While speedup in Table 1 shows the gain in computation time due to using PFASST instead of MLSDC, the speedup in Table 2 is computed with respect to the sequential first state, then adjoint approach, i.e., contains speedup due to parallelism as well as splitting the adjoint solve. Tables 3 and 4 show results for warm starting the respective methods. Total speedup is computed with respect to the sequential first state, then adjoint base method. Speedup to cold denotes the speedup compared to the same method and same number of processors, but using the standard PFASST predictor (burn-in instead of warm start, see Sect. 3.3). In this example, warm starts increase the speedup by 10–20%. For the first state, then adjoint variant, the 10% gain due to warm start in the 20 processor setting corresponds to a reduction in MLSDC sweeps of 14% (state) and 34% (adjoint). Thus, for problems where the sweeps are more expensive, we expect warm starting to be more effective (compare also Sect. 5.2).

Fig. 7 shows the progress of the optimization, which is similar for all approaches. The individual wallclock times are compared again in Fig. 8. For cold starting, the mixed approach performs better than the first state, then adjoint approach for most settings, but the advantage is rather small. Here, the mixed approach mainly saves computation time due to reduced communication, as both the propagation of the terminal condition in the mixed approach and the SDC sweeps in the plain first state, then adjoint approach are computed using the FFT and have a similar computational burden. Thus for other space discretizations, e.g., finite differences or finite elements, we expect the mixed approach to offer larger gains in computation times. Warm starting offers a moderate gain in speedup and efficiency for both variants (first state then adjoint, and mixed). As the gain is slightly higher in the first state, then adjoint approach, this variant provides the best performance in this test, achieving a parallel efficiency of $29\%$ on $20$ processors.

For completeness, the third variant, solving state and adjoint simultaneously (Sect. 4.2.2), leads to the same optimization progress (shown in Fig. 7), but gives a speedup of merely $1.8$ on $20$ processors, which is significantly worse than other variants.

5.2 Nonlinear problem: Nagumo equation

Next we consider the following optimal control problem (see also [8, 3, 30]):

[TABLE]

i.e., as in the linear case, minimize (10), subject to

[TABLE]

The parameter $\gamma$ steers the influence of the nonlinear reaction terms, and thus the stiffness due to the reaction term, and will be used below to compare IMEX and MISDC (multi-implicit SDC) formulations for solving the state equation. The spatial domain $\Omega=(0,L)$, with $L=20$, and the equation is solved up to final time $T=5$. The initial condition is

[TABLE]

and target solution

[TABLE]

where $y_{\text{nat}}$ denotes the solution to the PDE (35) for $u\equiv 0$. Here, for $\lambda=0$, an exact optimal control is known,

[TABLE]

which will be used as a comparison for the computed optimal controls. For this example, the adjoint equation is

[TABLE]

and the reduced gradient is given as $\nabla j(u)=\lambda u-p$. For solving state and adjoint equation we do not consider the mixed method here, but always use the first state, then adjoint approach, since the adjoint equation contains time-dependent terms in the differential operator (see the discussion in Sect. 4.2.3).

State equation and adjoint solve: MISDC vs. IMEX

Before coming to the actual optimization, let us briefly compare MISDC and IMEX-MLSDC for solving the state and adjoint equations (35), (36) with zero control and varying $\gamma$. This demonstrates not only the flexibility of PFASST, but allows us to choose a suitable method in the scaling study for the optimization problem.

The IMEX and MISDC approach are explained by examining a single substep of an SDC sweep. With $k$ denoting the SDC iteration, $m$ the substep index, and $D^{2}$ the discretization of the second derivative term, the fully implicit SDC version of an SDC substep (23) for (35) takes the form

[TABLE]

where the term $S^{[k]}_{m}$ contains terms that either depend on the previous iteration $[k]$ or values at iteration $[k+1]$ already computed at substep $i<m$, including the control terms arising from the discretization of $u(x,t)$. The implicit equation couples nonlinear reaction and diffusion terms and hence requires a global nonlinear solver in each substep. For problems in which the reaction terms are non-stiff and can be treated explicitly, the reaction terms at node $m$ do not appear in the implicit equation, giving the form

[TABLE]

Each IMEX substep now requires only the solution of a linear implicit equation, and thus is computationally cheaper than the fully implicit approach, provided that the explicit treatment of the reaction term does not impose an additional time step restriction. When the reaction term is stiff, and hence it is advantageous to treat it implicitly, a standard MISDC approach applies an operator splitting between diffusion and reaction in the correction equation. For example,

[TABLE]

We use Newton’s method with damping and the natural monotonicity test according to [11] to solve the nonlinear equation (40) in each sweep. To reduce the computational effort, the MISDC approach is further modified so that the nonlinear solve for reaction in (40) is made linear by lagging terms in the splitting:

[TABLE]

This form creates an implicit solve with roughly the same cost as treating reaction explicitly but is more stable.

In all cases, a $3$-level PFASST scheme is applied to a method of lines discretization. We use a pseudo-spectral discretization with 64/128/256 degrees of freedom in space, with spatial derivatives computed spectrally with the FFT. In time, $32$ uniform time steps are used, with 3/5/9 Lobatto IIIA collocation rules. Thus, all variants have the same order, and converge to the same collocation solution, if they converge. Sweeps are stopped when the relative or absolute residual drop below $10^{-11}$. Table 6 shows the required number of sweeps for convergence for several values of $\gamma$, time steps, and number of processors. While both MISDC variants stably converge for all values of $\gamma$ and time steps in sequential and parallel, IMEX requires finer time discretizations to converge, especially in parallel. In the sequential runs, the stability of the IMEX scheme requires a time step small enough so that the IMEX SDC iterations converge despite the explicit treatment of the reaction term. On the other hand, in PFASST, the sequential coarse grid IMEX sweeps are low-order and have a more restrictive stability constraint than the IMEX serial schemes. When the solvers converge all variants require approximately the same number of sweeps on average, hence the IMEX scheme would be expected to need the least computational time. Fig. 9 shows the solution for $\gamma=5$ and $32$ parallel time steps, computed with MISDC with lagging, and the difference to the solution obtained using MISDC with a nonlinear solver. While in most settings nonlinear and lagged versions need the same number of iterations, here the nonlinear variant requires on average $3$ iterations more than the lagged version, but the solutions agree to the convergence tolerance.

The conclusions of this comparison are two-fold. First, the benefits of using an MISDC scheme versus an IMEX scheme will depend on the relative stiffness of the two terms and the relative cost of the nonlinear versus linear implicit equations. Secondly, at least for the problems considered here, the MISDC with lagging method performs as good or better than an IMEX splitting in terms of iteration count with only a small computational overhead relative to the IMEX variant.

Optimization

To evaluate the performance of the parallel-in-time optimization method we use again $3$-level PFASST with an IMEX sweeper for solving state and adjoint equations up to a residual tolerance of $10^{-11}$. We set $\gamma=1$ in the state equation, and use 32/64/128 degrees of freedom in space as well as 32 time steps. In contrast to the experiments above, IMEX is converging in this setting due to the coarser space discretization, and has the smallest computation time. In the objective function the regularization parameter $\lambda=10^{-6}$ is used. Note that this leads to a difference between the computed optimal control and $u_{\text{exact}}$, as the latter is valid only for $\lambda=0$. The optimization is done using the nonlinear conjugate gradient method by Dai and Yuan [9] in combination with linesearch to satisfy the the strong Wolfe conditions. As reported in [8], the optimization progress with several ncg variants is slow, so we stop the optimization after 200 iterations. We again consider strong scaling with up to 32 CPUs. Due to the less severe nonlinearity in this test case and in view of the comparison above, only results using the IMEX variant are shown here.

Tables 7 and 8 show timings, speedup, and parallel efficiency for IMEX with cold and warm starts. For cold starts, a speedup of more than $5$ for $32$ processors is achieved, corresponding to a parallel efficiency of $16.5\%$.

For the case with 32 CPUs, warm starting reduces the number of required state sweeps by $65\%$, while the reduction in adjoint sweeps is less pronounced. Overall, this translates to a decrease in computation time by nearly $40\%$, leading to a parallel efficiency of around $24\%$. The gain here is significantly larger than in the linear heat example (Sect. 5.1). This is true despite the tests being in one spatial dimension: the PFASST method typically produces better parallel efficiencies for problems in higher dimensions since spatial coarsening then produces a greater relative reduction in computational cost on coarse levels.

Similar results (not reported here) are achieved using MISDC with lagging, albeit with overall higher computation times.

Figure 10 shows the objective function value over the optimization process with 1 and 32 processors as well as using cold and warm starts. We note that there are slight variances in the optimization progress and the computed controls, as also seen in the final error of the control in Table 7. The difference is that the sequential version does SDC sweeps on the later time steps with the correct initial value, while PFASST iteratively corrects the incorrect initial values. For warm starts, in addition the initial function evaluations as well as coarse grid corrections are inexact. Thus the iterative solution of the PDEs progresses differently, leading to slightly different final solutions and thus different gradients and step sizes, which accumulate during the optimization.

Comparison with classical time steppers

The results reported above use MLSDC on one processor as the sequential reference for the computation of speedups. This allows the use of the same time discretization for all setups in the strong scaling study. However, due to the relatively high computational cost, MLSDC may not be the best sequential method. For comparison we consider a standard IMEX-Euler method as well as a fourth-order additive Runge-Kutta (ARK-4) method of Kennedy and Carpenter [32]. The IMEX-Euler method treats the reaction part explicitly, the diffusion part implicitly; ARK-4 combines an explicit Runge-Kutta method for the reaction part with an ESDIRK method for the diffusion. To obtain a reference solution for state and adjoint equations, we fix the spatial discretization to 128 degrees of freedom as before and use SDC with 9 Lobatto IIIA collocation nodes on a very fine temporal mesh. Time step sizes $\Delta t$ for the different methods are chosen such that the time discretization error compared to the reference solution is approximately the same for each case, which would result in similar progress of the full optimization process. All other parameters are chosen as before.

One complication in this comparison is that computing the stage values in the Runge-Kutta method requires evaluating the control $u$ and desired state $y_{d}$ at intermediate times (as does the MLSDC method require this data at the collocation nodes). For MLSDC we chose symmetric collocation nodes, such that forward-in-time and backward-in-time use the same nodes. Thus, $u,y_{d}$ have distinct values at the collocation nodes. However, in the ARK-4 method, the stage values are not symmetric, and interpolation or dense output would be required to compute these values. For simplicity, we therefore (artificially) set $u$ and $y_{d}$ constant on the time steps for this example. Results are reported in table 9.

Despite the small time step size, IMEX-Euler is the fastest of the sequential methods. However, using PFASST with warm starting, speedup is achieved already by adding the second processor (speedup $1.2$). Compared to IMEX-Euler, on $32$ processors, PFASST with cold starting achieved a speedup/parallel efficiency of $3.2$/$10\%$, with warm starts this increases to $4.6$/$14\%$. As noted before, the benefit of MLSDC compared to SDC is greater for problems in more than one-dimension, hence these results should not be considered as a general comparison.

6 Conclusions and Outlook

In this paper, we introduce an efficient fully time-parallel strategy for gradient-based optimization with parabolic PDEs. In the most critical component, the computation of the reduced gradient, we discuss and compare three competing approaches for applying the PFASST algorithm to the solution of the state and adjoint equations. While simultaneously solving state and adjoint induces severe communication and wait times in the implementation, both the first state, then adjoint approach and the mixed approach yield comparatively good speedups and parallel efficiency of up to $29\%$. The warm starting capability of MLSDC/PFASST is an important feature, and makes the method competitive with standard time stepping methods even in a sequential run. It is also noteworthy that speedup is obtained already with two processors.

Although the presented results are promising, there are several additional avenues to further increase the parallel speedup or efficiency of the approach. For example, the test cases studied here were set in simple periodic domains using the FFT for implicit solves. In many cases, implicit solves are done iteratively and hence solver tolerances can be adjusted dynamically to do less work when less accuracy is needed [43]. The same is true of the tolerances used to terminate PFASST iterations. Likewise, the spatial or temporal order of the PDE solver could be changed dynamically during the optimization process. Hence a dynamic, adaptive strategy for controlling these tolerances and accuracy in the optimization loop is a promising direction for future study. Similarly, analyzing the impact of inexact storage of solution values for warm starting will be necessary for larger scale applications. Finding suitable extensions of the mixed approach to nonlinear state equations is also of interest. For this approach, it might also be beneficial to solve state and modified adjoint equation in the first step concurrently, making use of additional processors during SDC sweeps.

Acknowledgment

S.G. gratefully acknowledges partial funding by the Deutsche Forschungsgemeinschaft (DFG), Project WE 2937/6-1. The work of M.M. was supported by the Applied Mathematics Program of the DOE Office of Advanced Scientific Computing Research under the U.S. Department of Energy under contract DE-AC02-05CH11231, and the Alexander von Humboldt Foundation.