A spectral-in-time Newton–Krylov method for nonlinear PDE-constrained optimization

Abstract

We devise a method for nonlinear time-dependent partial-differential-equation-constrained optimization problems that uses a spectral-in-time representation of the residual, combined with a Newton–Krylov method to drive the residual to zero. We also propose a preconditioner to accelerate this scheme. Numerical results indicate that this method can achieve fast and accurate solution of nonlinear problems for a range of mesh sizes and problem parameters. The numbers of outer Newton and inner Krylov iterations required to reach the attainable accuracy of a spatial discretization are robust with respect to the number of collocation points in time and also do not change substantially when other problem parameters are varied.

PDE-constrained optimization, time-dependent PDE, residual-based method, Newton–Krylov method, coupled systems of PDEs, preconditioner

1. Introduction

Partial-differential-equation (PDE)-constrained optimization is an important class of mathematical problems (Lions, 1971; Ito & Kunisch, 2008; Tröltzsch, 2010), with a wide range of applications across science and engineering (see Arridge, 1999; Cheney et al., 1999; Griesse & Volkwein, 2005; Gunzberger, 2010, for instance). The fast and accurate solution of the first-order optimality conditions resulting from such problems is a significant challenge for researchers in these communities. For example, when an ‘all-at-once’ approach is applied to solve such conditions, one is faced with coupled linear systems of huge scale, particularly when standard finite difference or finite element schemes are used for the discretization procedure in the spatial coordinates. For time-dependent problems there is an additional challenge of how to appropriately discretize in the time variable. When fast solvers for the resulting linear system are required, in particular preconditioned iterative methods, it has been popular to utilize a backward Euler scheme for this purpose (see for example Pearson & Stoll, 2013; Pearson et al., 2012; Stoll & Breiten, 2015); however, this requires a very small step size in order to obtain an accurate numerical solution, and thus very large linear systems to solve.

Motivated by this, in Güttel & Pearson (2018) the authors devised deferred correction methods for linear PDE-constrained optimization problems, where equations for the errors and residuals at each deferred correction step were constructed in order to successively increase the order of the time-stepping scheme. As a result one is required to solve a sequence of much smaller linear systems in order to achieve the same accuracy of the numerical solution, compared to a method without a deferred correction acceleration. This approach was found to be highly effective for linear PDE-constrained optimization problems; however, a significant question remained about how to devise a related strategy for nonlinear PDE-constrained optimization. These problems possess a vastly increased level of difficulty, compared to linear problems, as the matrices describing the spatial behavior of the physical system are different at every time step.

In this paper we devise a residual-based approach for nonlinear PDE-constrained optimization problems; in particular this is based on using a spectral-in-time representation of the residual that is then linearized and solved by a Newton method. The use of a spectral method in the time variable means that high accuracy solutions can be obtained with only a small number of time steps, keeping the linearized Newton system of relatively small dimension. We implement our approach using a Newton–Krylov method, of the form described in Kelley (2003, Chapter 3) and Knoll & Keyes (2004). To make such a method numerically viable we suggest a general preconditioning strategy, which is found to substantially accelerate the Newton–Krylov scheme for the problems examined. Care has to be taken with the implementation of the preconditioner as the spectral time integration matrix is dense. We hence transform the Jacobian arising in the Newton system into a form that can be easily approximated by a Kronecker product, allowing for the application of a preconditioner whose cost is approximately proportional to the number of spatial degrees of freedom. Our approach mitigates a key difficulty encountered with alternative discretization techniques for nonlinear PDE-constrained optimization problems, specifically being required to solve linear systems arising from very large numbers of grid points in time to obtain even modest discretization error properties, and our new method is found to be an effective strategy for a number of examples, mesh sizes and problem parameters.

We highlight that there has been a substantial amount of research undertaken with the goal of achieving numerical solutions to time-dependent PDE-constrained optimization problems, with low discretization error and within reasonable computation time. Among many valuable references, we refer to Dolgov & Stoll (2017); Stoll & Breiten (2015) for low-rank solution methods, to Pearson & Stoll (2013), Pearson et al. (2012) for preconditioned iterative solvers, to Borzì & Schulz (2009), Hinze et al. (2012) for a discussion of multigrid approaches, to Heinkenschloss (2005) for a multiple shooting strategy, to Maday & Turinici (2002), Mathew et al. (2010) for parareal approaches and to Götschel & Minion (2019) for a recently developed time-parallel method with the computation of the adjoint gradient information performed using the PFASST framework. We also point to the body of work on Krylov deferred correction methods for initial value problems developed in Huang et al. (2007), Jia & Huang (2008). Indeed, for linear initial value problems, deferred correction can be interpreted as a preconditioned Newton-like method for solving the time collocation system (Huang et al., 2006; Weiser, 2015).

While the use of a spectral-in-time residual function in this work is inspired by the spectral deferred correction approach in Dutt et al. (2000), we found in extensive numerical tests that the direct Newton-based minimization of this function without the outer deferred correction loop is both conceptually simpler and indeed more efficient for nonlinear problems. This is the approach we would like to explore in this paper. Compared to our previous work (Güttel & Pearson, 2018) on linear PDE-constrained optimization, the method described here is not a deferred correction scheme. Our aim is to provide a general approach for nonlinear time-dependent problems, which, as we will demonstrate, can lead to very accurate solutions while being applicable with any discretization in the spatial variables of the user’s choice and which has the potential to be combined with a number of the approaches listed above. Our focus is on problems with a squared |$L^2$|-norm regularization term applied to the control variable in the objective function, along with first-order derivatives in time and a linear term involving the control within the PDE constraint, as such a formulation allows us to eliminate the control variable a priori and solve a coupled state–adjoint system. We believe that, with suitable modifications that we will describe, our method is more generalizable still.

This paper is structured as follows. In Section 2 we present the residual-based method for nonlinear PDE-constrained optimization. Specifically, in Section 2.1 we state the coupled systems of PDEs arising from the PDE-constrained optimization problems considered in this paper, in Section 2.2 we derive the residual-in-time representation and introduce the Newton–Krylov method we apply, then in Section 2.3 we present the preconditioner that is embedded within the Newton–Krylov scheme. Sections 2.4 and 2.5 further discuss the implementation and the stopping condition used. Section 3 presents the results of a range of numerical experiments, using test problems for which analytic solutions are given in Appendix  A. In Section 4 we provide our concluding remarks.

2. Residual-based Newton–Krylov method

In this section we present our residual-based method for nonlinear time-dependent PDE-constrained optimization problems. To focus our discussion we will present the methodology using two specific examples, although this approach is much more general. Firstly, motivated by research in literature such as Buchholz et al. (2013),Gunzburger et al. (1991), we consider the optimal control of the Schlögl equation given by

$$\begin{align} \ \min_{y,c}~~\frac{1}{2}\int_{0}^{T}\int_{\varOmega}\big(y-\widehat{y}\big)^{2}~\textrm{d}\varOmega&\,\textrm{d}t+\frac{\beta}{2}\int_{0}^{T}\int_{\varOmega} c^{2}~ \textrm{d}\varOmega\, \textrm{d}t \\[-3pt] \ \nonumber \text{s.}\,\text{t.}\quad~~\frac{\partial{}y}{\partial{}t}-\nabla^{2}y-y+y^{3}&{}=c+f \textrm{in}\ \varOmega\times(0,T),& \\[-3pt] \ \nonumber y(\vec{x},t)&{}=h(\vec{x},t) \textrm{on}\ \partial\varOmega\times(0,T),& \\[-3pt] \ \nonumber y(\vec{x},0)&{}=y_{0}(\vec{x}) \textrm{at}\ t=0.& \end{align}$$

(2.1)

Here, |$y$|⁠, |$\widehat{y}$| and |$c$| denote the state, desired state and control variables, respectively, |$f$| and |$h$| are given functions in space and time, |$y_{0}$| is a given initial condition, |$\beta>0$| denotes a regularization parameter (or Tikhonov parameter) and |$T>0$| is the final time. The PDE constraint is posed on a space–time domain, with |$(\vec{x},t)\in \varOmega \times (0,T)$|⁠, and with |$\partial \varOmega $| denoting the boundary of |$\varOmega $|⁠. We highlight that formulation (2.1) is similar in structure to the optimal control problem constrained by the Nagumo equation considered in Götschel & Minion (2019, Section 5.2).

Secondly, we consider a reaction–diffusion control problem, based on the formulation in Barthel et al. (2010):

$$\begin{align} \ \min_{y,z,c}~~\frac{\beta_{y}}{2}\int_{0}^{T}\int_{\varOmega}\big(y-\widehat{y}\big)^{2}~&\textrm{d}\varOmega\, \textrm{d}t+\frac{\beta_{z}}{2}\int_{0}^{T}\int_{\varOmega}\big(z-\widehat{z}\big)^{2}~\textrm{d}\varOmega\,\textrm{d}t+\frac{\beta_{c}}{2}\int_{0}^{T}\int_{\partial\varOmega} c^{2}~\textrm{d}s\, \textrm{d}t \\[-3pt] \ \nonumber \text{s.}\,\text{t.}\quad\frac{\partial{}y}{\partial{}t}-D_{1}\nabla^{2}y+k_{1}y&{}=-\gamma_{1}yz+f \textrm{in}\ \varOmega\times(0,T),& \\[-3pt] \ \nonumber \frac{\partial{}z}{\partial{}t}-D_{2}\nabla^{2}z+k_{2}z&{}=-\gamma_{2}yz+g \textrm{in}\ \varOmega\times(0,T),& \\[-3pt] \ \nonumber D_{1}\frac{\partial{}y}{\partial{}n}&{}=c \textrm{on}\ \partial\varOmega\times(0,T),& \\[-3pt] \ \nonumber D_{2}\frac{\partial{}z}{\partial{}n}&{}=0 \textrm{on}\ \partial\varOmega\times(0,T),& \end{align}$$

(2.2)

$$\begin{align}\nonumber\\[-40pt] \ \nonumber y(\vec{x},0)&{}=y_{0}(\vec{x}) \textrm{at}\ t=0,& \\ \ \nonumber z(\vec{x},0)&{}=z_{0}(\vec{x}) \textrm{at}\ t=0.&\end{align}$$

(2.2)

In this formulation |$y$| and |$z$| denote state variables with corresponding desired states |$\widehat{y}$| and |$\widehat{z}$|⁠, |$c$| is again a control variable, |$f$| and |$g$| are given functions in space and time, |$y_{0}$| and |$z_{0}$| are given initial conditions and |$T>0$| is again a final time. The parameters |$\beta _{c}$|⁠, |$D_{1}$| and |$D_{2}$| are positive, with |$\beta _{y}$|⁠, |$\beta _{z}$|⁠, |$k_{1}$|⁠, |$k_{2}$|⁠, |$\gamma _{1}$| and |$\gamma _{2}$| non-negative (and generally positive). In particular at least one of |$\beta _{y}$| and |$\beta _{z}$| must be positive. We highlight that a key structural difference between the two examples is that (2.1) is a distributed control problem, whereas (2.2) is a boundary control problem. In both examples the source terms |$f$| and |$g$| would often be zero in practice, but we allow nonzero functions so that test problems with analytic solutions are more readily available.

2.1 Derivation of coupled systems of PDEs

In order to apply our residual-based approach, the first step is to describe the solution of our time-dependent PDE-constrained optimization problems using a coupled system of PDEs, which define the first-order optimality conditions. A general approach, an outline of which is provided in Tröltzsch (2010, Chapter 3) for time-dependent (parabolic) PDE-constrained optimization problems, is to apply an optimize-then-discretize strategy. In that strategy one seeks the stationary points of a continuous Lagrangian involving the cost functional being minimized, as well as the PDE constraints and boundary conditions enforced by an adjoint variable (or Lagrange multiplier). One may then take Fréchet derivatives in the direction of the adjoint, control and state variables, and test the conditions that the derivatives must be equal to zero using functions with appropriate continuity and differentiability properties and which satisfy suitable boundary conditions. Taking the derivatives in this order gives rise to state equations, gradient equations and adjoint equations, respectively.

For instance, to briefly outline how this may be applied to the optimal control problem involving the Schlögl equation (2.1), the continuous Lagrangian reads

$$\begin{align*} \ \mathcal{L}(y,c,p)={}&\frac{1}{2}\int_{0}^{T}\int_{\varOmega}\big(y-\widehat{y}\big)^{2}~\textrm{d}\varOmega\, \textrm{d}t+\frac{\beta}{2}\int_{0}^{T}\int_{\varOmega} c^{2}~ \textrm{d}\varOmega\, \textrm{d}t \\ \ &+\int_{0}^{T}\int_{\varOmega}p_{\varOmega}\left(\frac{\partial{}y}{\partial{}t}-\nabla^{2}y-y+y^{3}-c-f\right)~\textrm{d}\varOmega\, \textrm{d}t+\int_{0}^{T}\int_{\partial\varOmega}p_{\partial\varOmega}\left(y-h\right)~\textrm{d}s\, \textrm{d}t, \end{align*}$$

where the adjoint variable is split into its components in the interior of |$\varOmega $| and boundary |$\partial \varOmega $|⁠, denoted |$p_{\varOmega }$| and |$p_{\partial \varOmega }$|⁠. Here, we have omitted the initial condition from the definition of the Lagrangian for ease of presentation. Using integration by parts and applying the divergence theorem allows us to rearrange the Lagrangian to obtain

$$\begin{align*} \ \mathcal{L}(y,c,p)={}&\frac{1}{2}\int_{0}^{T}\int_{\varOmega}\big(y-\widehat{y}\big)^{2}~\textrm{d}\varOmega\, \textrm{d}t+\frac{\beta}{2}\int_{0}^{T}\int_{\varOmega} c^{2}~ \textrm{d}\varOmega\, \textrm{d}t \\ \ &+\int_{0}^{T}\int_{\varOmega}y\left(-\frac{\partial{}p_{\varOmega}}{\partial{}t}-\nabla^{2}p_{\varOmega}-p_{\varOmega}+y^{2}p_{\varOmega}\right)~\textrm{d}\varOmega\, \textrm{d}t-\int_{0}^{T}\int_{\varOmega}p_{\varOmega}\left(c+f\right)~\textrm{d}\varOmega\, \textrm{d}t \\ \ &+\int_{\varOmega}\left[p_{\varOmega}y\right]_{t=0}^{t=T}~ \textrm{d}\varOmega+\int_{0}^{T}\int_{\partial\varOmega}\left(y\frac{\partial{}p_{\varOmega}}{\partial{}n}-p_{\varOmega}\frac{\partial{}y}{\partial{}n}\right)~\textrm{d}s\, \textrm{d}t+\int_{0}^{T}\int_{\partial\varOmega}p_{\partial\varOmega}\left(y-h\right)~\textrm{d}s\, \textrm{d}t. \end{align*}$$

Applying the Fréchet derivatives in the directions of |$p_{\varOmega }$|⁠, |$c$| and |$y$|⁠, testing with appropriate functions and relabeling the adjoint variable |$p_{\varOmega }$| as |$p$| gives the state, gradient and adjoint equations, respectively. Here,

$$\begin{align} \ \left.\begin{array}{rl} \frac{\partial{}y}{\partial{}t}-\nabla^{2}y-y+y^{3}=c+f & \textrm{in}\ \varOmega\times(0,T), \\ y(\vec{x},t)=h(\vec{x},t) & \textrm{on}\ \partial\varOmega\times(0,T), \\ y(\vec{x},0)=y_{0}(\vec{x}) & \textrm{at}\ t=0, \\ \end{array}\right\}&~~\text{state equation,} \end{align}$$

(2.3)

$$\begin{align} \ \left.\begin{array}{rl} \beta{}c-p=0 & \textrm{in}\ \varOmega\times(0,T), \\ \end{array}\right\}&~~\text{gradient equation,} \end{align}$$

(2.4)

$$\begin{align} \ \left.\begin{array}{rl} -\frac{\partial{}p}{\partial{}t}-\nabla^{2}p-p+3y^{2}p=\widehat{y}-y & \textrm{in}\ \varOmega\times(0,T), \\ p(\vec{x},t)=0 & \textrm{on}\ \partial\varOmega\times(0,T), \\ p(\vec{x},T)=0 & \textrm{at}\ t=T \\ \end{array}\right\}&~~\text{adjoint equation.} \end{align}$$

(2.5)

Applying a similar approach to the reaction–diffusion control problem (2.2), the first-order optimality conditions are as follows (see also Barthel et al., 2010):

$$\begin{align} \ \left.\begin{array}{rl} \frac{\partial{}y}{\partial{}t}-D_{1}\nabla^{2}y+k_{1}y=-\gamma_{1}yz+f & \textrm{in}\ \varOmega\times(0,T), \\[4pt] \frac{\partial{}z}{\partial{}t}-D_{2}\nabla^{2}z+k_{2}z=-\gamma_{2}yz+g & \textrm{in}\ \varOmega\times(0,T), \\[4pt] D_{1}\frac{\partial{}y}{\partial{}n}=c & \textrm{on}\ \partial\varOmega\times(0,T), \\[4pt] D_{2}\frac{\partial{}z}{\partial{}n}=0 & \textrm{on}\ \partial\varOmega\times(0,T), \\[4pt] y(\vec{x},0)=y_{0}(\vec{x}) & \textrm{at}\ t=0, \\ z(\vec{x},0)=z_{0}(\vec{x}) & \textrm{at}\ t=0 \\ \end{array}\right\}&~~\text{state equations,} \end{align}$$

(2.6)

$$\begin{align} \ \left.\begin{array}{rl} \beta_{c}c-p=0 & \textrm{on}\ \partial\varOmega\times(0,T) \\ \end{array}\right\}&~~\text{gradient equation,} \end{align}$$

(2.7)

$$\begin{align} \ \left.\begin{array}{rl} -\frac{\partial{}p}{\partial{}t}-D_{1}\nabla^{2}p+k_{1}p+\gamma_{1}zp+\gamma_{2}zq=\beta_{y}\big(\widehat{y}-y\big) & \textrm{in}\ \varOmega\times(0,T), \\[4pt] -\frac{\partial{}q}{\partial{}t}-D_{2}\nabla^{2}q+k_{2}q+\gamma_{1}yp+\gamma_{2}yq=\beta_{z}\big(\widehat{z}-z\big) & \textrm{in}\ \varOmega\times(0,T), \\[4pt] D_{1}\frac{\partial{}p}{\partial{}n}=0 & \textrm{on}\ \partial\varOmega\times(0,T), \\[4pt] D_{2}\frac{\partial{}q}{\partial{}n}=0 & \textrm{on}\ \partial\varOmega\times(0,T), \\[2pt] p(\vec{x},T)=0 & \textrm{at}\ t=T, \\ q(\vec{x},T)=0 & \textrm{at}\ t=T \\ \end{array}\right\}&~~\text{adjoint equations.} \end{align}$$

(2.8)

Note how in both cases the first-order optimality conditions take the form of coupled time-dependent PDEs, with the state equations equipped with initial conditions, and the adjoint equations possessing final-time conditions.

2.2 Derivation of the Newton system

The formulations derived in the previous section, specifically (2.3)–(2.5) for the Schlögl problem and (2.6)–(2.8) for the reaction–diffusion control problem, may be discretized in space and result in coupled initial/final value problems in time, of the form

$$\begin{align*}\begin{cases} M_u \textbf{u}^{\prime}(t) = \textbf{f}(t,\textbf{u},\textbf{v}), \qquad M_u \textbf{u}(0) = M_u \textbf{u}_0\in \mathbb{R}^N\ \textrm{given},\qquad\qquad\qquad\qquad{(2.9)} \\ M_v \textbf{v}^{\prime}(t) = \textbf{g}(t,\textbf{u},\textbf{v}), \qquad M_v \textbf{v}(T) = M_v \textbf{v}_T\in \mathbb{R}^N\ \textrm{given},\qquad\qquad\qquad\quad\ \ {(2.10)}\end{cases} \end{align*}$$

(2.9)

for two vector-valued functions |$\textbf{u},\textbf{v}: [0,T] \mapsto \mathbb{R}^N$| and with matrices |$M_u,M_v\in \mathbb{R}^{N\times N}$|⁠. This system may be thought of as a coupled system of ordinary differential equations in the time variable. We highlight the key structures that allow us to rewrite the optimization problem in the form (2.9)–(2.10): (i) the spatial derivatives are separate from the first-order time derivatives within the PDE constraints; (ii) a squared |$L^2$|-norm term measures the control within the objective function and a linear term incorporates the control within the constraints, allowing us to derive a linear relation between control and adjoint variables (as above). We believe the second structure could be relaxed in order to allow a nonlinear relation between control and adjoint, which may be incorporated into (2.9)–(2.10) through an additional algebraic constraint (i.e., not involving a time derivative). In certain cases, higher-order PDEs in time could also be written as (2.9)–(2.10) by rewriting these as a system of first-order equations.

In the case of the Schlögl example, substituting the gradient equation (2.4) into (2.3) leads to a coupled system of the form (2.9)–(2.10) with

$$\begin{align*} \ &\textbf{u}(t):=u(t)\leftarrow{}y(t),\quad{}M_{u}\leftarrow\textrm{Id},\quad\textbf{f}(t,\textbf{u},\textbf{v})\leftarrow\nabla^{2}u+u-u^{3}+\frac{1}{\beta}v+f, \\ \ &\textbf{v}(t):=v(t)\leftarrow{}p(t),\quad{}M_{v}\leftarrow\textrm{Id},\quad\textbf{g}(t,\textbf{u},\textbf{v})\leftarrow-\nabla^{2}v-v+3u^{2}v+u-\widehat{y}. \end{align*}$$

We have used the notation ‘|$\leftarrow $|’, with the quantities on the left-hand side being the result of spatial discretization. The symbol |$\textrm{Id}$| refers to the discretized identity operator.

For the reaction–diffusion example we again have a system of the form (2.9)–(2.10) with

$$\begin{align*} \ &\textbf{u}(t):=\left[\begin{array}{c} u_1(t) \\ u_2(t) \\ \end{array}\right]\leftarrow\left[\begin{array}{c} y(t) \\ z(t) \\ \end{array}\right],\quad{}M_{u}\leftarrow\left[\begin{array}{cc} \textrm{Id} & 0 \\ 0 & \textrm{Id} \\ \end{array}\right],\quad\textbf{f}(t,\textbf{u},\textbf{v})\leftarrow\left[\begin{array}{c} D_{1}\nabla^{2}u_1-k_{1}u_1-\gamma_{1}u_1 u_2+f \\ D_{2}\nabla^{2}u_2-k_{2}u_2-\gamma_{2}u_1 u_2+g \\ \end{array}\right], \\ \ &\textbf{v}(t):=\left[\begin{array}{c} v_1(t) \\ v_2(t) \\ \end{array}\right]\leftarrow\left[\begin{array}{c} p(t) \\ q(t) \\ \end{array}\right],\quad{}M_{v}\leftarrow\left[\begin{array}{cc} \textrm{Id} & 0 \\ 0 & \textrm{Id} \\ \end{array}\right],\quad{} \\ \ &\textbf{g}(t,\textbf{u},\textbf{v})\leftarrow\left[\begin{array}{c} -D_{1}\nabla^{2}v_1+k_{1}v_1+\gamma_{1}u_2 v_1+\gamma_{2}u_2 v_2+\beta_{y}\big(u_1-\widehat{y}\big) \\ -D_{2}\nabla^{2}v_2+k_{2}v_2+\gamma_{1}u_1 v_1+\gamma_{2}u_1 v_2+\beta_{z}\big(u_2-\widehat{z}\big) \\ \end{array}\right]. \end{align*}$$

Having derived systems of the above forms we need to equip the relevant operators with boundary conditions. For example, in the reaction–diffusion example and with the notation above, using the gradient equation (2.7) gives that the first Neumann boundary condition reads |$D_{1}\frac{\partial{}u_1}{\partial{}n}=\frac{1}{\beta _{c}}v_1$|⁠.

Given a system of the form (2.9)–(2.10) let us assume that approximations |$\widetilde{\textbf{u}}_j$| and |$\widetilde{\textbf{v}}_j$| at time points |$0=t_0<t_1<\cdots < t_n=T$| are available and consider the interpolants

$$\begin{align}& \widetilde{\textbf{u}}(t) = \sum_{j=0}^n \ell_j(t) \widetilde{\textbf{u}}_j \quad \textrm{and} \quad \widetilde{\textbf{v}}(t) = \sum_{j=0}^n \ell_j(t) \widetilde{\textbf{v}}_j, \end{align}$$

(2.11)

where |$\ell _j(t)$| are differentiable Lagrange functions satisfying |$\ell _j(t_i) = \delta _{ij}$|⁠. We consider the associated residual functions

$$\begin{align*} \begin{cases}\displaystyle \textbf{r}_u(t) := \int_0^t \textbf{f}(\tau,\widetilde{\textbf{u}}(\tau),\widetilde{\textbf{v}}(\tau))\, \mathrm{d} \tau - M_u \widetilde{\textbf{u}}(t) + M_u \widetilde{\textbf{u}}(0),\qquad\qquad\qquad\qquad{(2.12)}\\ \displaystyle \ \textbf{r}_v(t) := \int_0^t \textbf{g}(\tau,\widetilde{\textbf{u}}(\tau),\widetilde{\textbf{v}}(\tau)) \,\mathrm{d} \tau - M_v \widetilde{\textbf{v}}(t) + M_v \widetilde{\textbf{v}}(0).\qquad\qquad\qquad\qquad{(2.13)} \end{cases}\end{align*}$$

(2.12)

It is important to note that the residuals are zero for all solutions that satisfy the system (2.9) and (2.10), irrespective of the initial/final condition for |$\textbf{u}$|/|$\textbf{v}$|⁠. Hence, we still need to impose |$M_u \widetilde{\textbf{u}}(0) = M_u \textbf{u}_0$| and |$M_v \widetilde{\textbf{v}}(T) = M_v \textbf{v}_T$| explicitly.

Denoting by |$\textbf{r}_{u,j}$| and |$\textbf{r}_{v,j}$| the approximations to the residuals |$\textbf{r}_u(t_j)$| and |$\textbf{r}_v(t_j)$| at time points |$t_j$|⁠, respectively, we can write

$$\begin{align*} & [ \textbf{r}_{u,0},\textbf{r}_{u,1},\ldots,\textbf{r}_{u,n} ] = [ \textbf{f}_0,\textbf{f}_1,\ldots,\textbf{f}_n] Q + M_u [\widetilde{\textbf{u}}_0 - \widetilde{\textbf{u}}_0,\widetilde{\textbf{u}}_0 - \widetilde{\textbf{u}}_1,\ldots,\widetilde{\textbf{u}}_0 - \widetilde{\textbf{u}}_n] \end{align*}$$

and

$$\begin{align*} & [ \textbf{r}_{v,0},\textbf{r}_{v,1},\ldots,\textbf{r}_{v,n} ] = [ \textbf{g}_0,\textbf{g}_1,\ldots,\textbf{g}_n] Q + M_v [\widetilde{\textbf{v}}_0 - \widetilde{\textbf{v}}_0,\widetilde{\textbf{v}}_0 - \widetilde{\textbf{v}}_1,\ldots,\widetilde{\textbf{v}}_0 - \widetilde{\textbf{v}}_n]. \end{align*}$$

Here, the |$(n+1)$|-by-|$(n+1)$| collocation matrix |$Q$| corresponds to cumulative integration, i.e.,

$$\begin{align}& q_{ij} = \int_{0}^{t_{j}} \ell_i(\tau)\, \mathrm{d} \tau, \quad i,j=0,1,\ldots,n, \end{align}$$

(2.14)

with |$q_{00}$| being the top-left entry of |$Q$|⁠. (The name of that matrix can be remembered by thinking of ‘quadrature’.) The motivation for using forward integration for both (2.12) and (2.13), as opposed to using coupled forward–backward integration as in our previous work (Güttel & Pearson, 2018), is to avoid the introduction of separate time-collocation matrices |$Q_{\textbf{f}}$| and |$Q_{\textbf{g}}$| (in Güttel & Pearson, 2018, these matrices were called |$C_u$| and |$C_v$|⁠, respectively). We will see that using a common collocation matrix for both |$\textbf{f}$| and |$\textbf{g}$| results in a regular matrix structure that we can exploit for the construction of a computationally efficient preconditioner.

We define the global residual function|$\textbf{R}:\mathbb{R}^{2N(n+1)}\to \mathbb{R}^{2N(n+1)}$| as

$$\begin{align}& \textbf{R}: \underbrace{\left[\begin{array}{c} \widetilde{\textbf{u}}_0\\ \widetilde{\textbf{v}}_0\\ \widetilde{\textbf{u}}_1\\ \widetilde{\textbf{v}}_1\\ \widetilde{\textbf{u}}_2\\ \widetilde{\textbf{v}}_2\\ \vdots\\ \widetilde{\textbf{u}}_n\\ \widetilde{\textbf{v}}_n \end{array}\right]}_{=:\textbf{w}} \mapsto \left[\begin{array}{c} \textbf{0}\\ \textbf{0}\\ \sum_{i=0}^n q_{i1} \textbf{f}_i \\ \sum_{i=0}^n q_{i1} \textbf{g}_i \\ \sum_{i=0}^n q_{i2} \textbf{f}_i \\ \sum_{i=0}^n q_{i2} \textbf{g}_i\\ \vdots\\ \sum_{i=0}^n q_{in} \textbf{f}_i \\ \sum_{i=0}^n q_{in} \textbf{g}_i \end{array}\right] + \left[\begin{array}{c} \widetilde{\textbf{u}}_0 - \textbf{u}_0\\ \widetilde{\textbf{v}}_n - \textbf{v}_T\\ M_u (\widetilde{\textbf{u}}_0 - \widetilde{\textbf{u}}_1)\\ M_v (\widetilde{\textbf{v}}_0 - \widetilde{\textbf{v}}_1)\\ M_u (\widetilde{\textbf{u}}_0 - \widetilde{\textbf{u}}_2)\\ M_v (\widetilde{\textbf{v}}_0 - \widetilde{\textbf{v}}_2)\\ \vdots\\ M_u (\widetilde{\textbf{u}}_0 - \widetilde{\textbf{u}}_n)\\ M_v (\widetilde{\textbf{v}}_0 - \widetilde{\textbf{v}}_{n}) \end{array}\right]. \end{align}$$

(2.15)

We aim to apply the Newton method |$\textbf{w}^{(k+1)} = \textbf{w}^{(k)} - [\mathcal{D}\textbf{R}(\textbf{w}^{(k)})]^{-1}\textbf{R}(\textbf{w}^{(k)})$| in order to approximately solve |$\textbf{R}\textbf{w} = \textbf{0}$|⁠. The Jacobian |$\textbf{J} = \mathcal{D}\textbf{R}(\textbf{w})$| is readily computed as

$$\begin{equation} \textbf{J} = \left[\begin{array}{@{}cc|ccccc@{}} I_N & O & \cdots & \cdots &\cdots & O & O \\ O & O & \cdots & \cdots &\cdots & O & I_N \\\hline q_{01} \mathcal{D}_u\textbf{f}_0 + M_u & q_{01}\mathcal{D}_v\textbf{f}_0 & q_{11}\mathcal{D}_u\textbf{f}_1 - M_u & q_{11}\mathcal{D}_v\textbf{f}_1 & \cdots & q_{n1}\mathcal{D}_u\textbf{f}_n & q_{n1}\mathcal{D}_v\textbf{f}_n\\ q_{01} \mathcal{D}_u\textbf{g}_0 & q_{01}\mathcal{D}_v\textbf{g}_0 + M_v & q_{11}\mathcal{D}_u\textbf{g}_1 & q_{11}\mathcal{D}_v\textbf{g}_1 - M_v & \cdots & q_{n1}\mathcal{D}_u\textbf{g}_n & q_{n1}\mathcal{D}_v\textbf{g}_n\\ \vdots&\vdots & \vdots &\vdots & \ddots&\vdots&\vdots \\ q_{0n} \mathcal{D}_u\textbf{f}_0 + M_u & q_{0n}\mathcal{D}_v\textbf{f}_0 & q_{1n}\mathcal{D}_u\textbf{f}_1 & q_{1n}\mathcal{D}_v\textbf{f}_1 & \cdots & q_{nn}\mathcal{D}_u\textbf{f}_n - M_u & q_{nn}\mathcal{D}_v\textbf{f}_n\\ q_{0n} \mathcal{D}_u\textbf{g}_0 & q_{0n}\mathcal{D}_v\textbf{g}_0 + M_v & q_{1n}\mathcal{D}_u\textbf{g}_1 & q_{1n}\mathcal{D}_v\textbf{g}_1 & \cdots & q_{nn}\mathcal{D}_u\textbf{g}_n & q_{nn}\mathcal{D}_v\textbf{g}_n- M_v \end{array}\right], \end{equation}$$

(2.16)

with |$\mathcal{D}_u\textbf{f}_i\in \mathbb{R}^{N\times N}$| denoting the Jacobian of |$\textbf{f}(t,\textbf{u},\textbf{v})$| with respect to |$\textbf{u}$|⁠, evaluated at the linearization point |$(t_i,\widetilde{\textbf{u}}_i,\widetilde{\textbf{v}}_i)$|⁠. The other Jacobians |$\mathcal{D}_v\textbf{f}_i, \mathcal{D}_u\textbf{g}_i, \mathcal{D}_v\textbf{g}_i$| are denoted analogously, |$O$| denotes the |$N$|-by-|$N$| zero matrix.

At each Newton iteration we need to solve a linear system with the matrix |$\textbf{J}$| defined in (2.16). In order to simplify its structure and make it more amenable to preconditioning we apply a block-column transformation by right-multiplication with a matrix |$K$|⁠, adding trailing block-columns to the first and second block-columns of |$\textbf{J}$| and thereby eliminating all appearances of |$M_u$| and |$M_v$| in the first two block-columns and creating an invertible pivot in the leading 2-by-2 block. More precisely the matrix |$K$| is given by

$$\begin{align*} & K = \left[\begin{array}{@{}ccccc@{}} 1 & & \\ 1 & 1 & \\ 1 & & 1 \\ \vdots & & & \ddots \\ 1 & & & & 1 \end{array}\right]\otimes I_{2N}. \end{align*}$$

We denote the resulting matrix by |$\widetilde{\textbf{J}} = \textbf{J}\cdot K$|⁠,

$$\begin{equation*} \widetilde{\textbf{J}} = \left[\begin{array}{@{}cc|ccccccc@{}} I_N & O & \cdots & \cdots & \cdots& O & O \\ O & I_N & \cdots & \cdots & \cdots & O & I_N \\ \hline \mathcal{D}_u \widetilde{\textbf{f}}_1 & \mathcal{D}_v \widetilde{\textbf{f}}_1 & q_{11}\mathcal{D}_u\textbf{f}_1 - M_u & q_{11}\mathcal{D}_v\textbf{f}_1 & \cdots & q_{n1}\mathcal{D}_u\textbf{f}_n & q_{n1}\mathcal{D}_v\textbf{f}_n\\ \mathcal{D}_u \widetilde{\textbf{g}}_1 & \mathcal{D}_v \widetilde{\textbf{g}}_1 & q_{11}\mathcal{D}_u\textbf{g}_1 & q_{11}\mathcal{D}_v\textbf{g}_1 - M_v & \cdots & q_{n1}\mathcal{D}_u\textbf{g}_n & q_{n1}\mathcal{D}_v\textbf{g}_n\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \mathcal{D}_u \widetilde{\textbf{f}}_n & \mathcal{D}_v \widetilde{\textbf{f}}_n & q_{1n}\mathcal{D}_u\textbf{f}_1 & q_{1n}\mathcal{D}_v\textbf{f}_1& \cdots & q_{nn}\mathcal{D}_u\textbf{f}_n - M_u & q_{nn}\mathcal{D}_v\textbf{f}_n\\ \mathcal{D}_u \widetilde{\textbf{g}}_n & \mathcal{D}_v \widetilde{\textbf{g}}_n & q_{1n}\mathcal{D}_u\textbf{g}_1 & q_{1n}\mathcal{D}_v\textbf{g}_1 &\cdots & q_{nn}\mathcal{D}_u\textbf{g}_n & q_{nn}\mathcal{D}_v\textbf{g}_n- M_v \end{array}\right], \end{equation*}$$

with |$\mathcal{D}_u \widetilde{\textbf{f}}_i:= \sum _{j=0}^n q_{ji} \mathcal{D}_u\textbf{f}_j$|⁠, |$\mathcal{D}_u \widetilde{\textbf{g}}_i:= \sum _{j=0}^n q_{ji} \mathcal{D}_u\textbf{g}_j$|⁠, |$\mathcal{D}_v \widetilde{\textbf{f}}_i:= \sum _{j=0}^n q_{ji} \mathcal{D}_v\textbf{f}_j$| and |$\mathcal{D}_v \widetilde{\textbf{g}}_i:= \sum _{j=0}^n q_{ji} \mathcal{D}_v\textbf{g}_j$|⁠, for |$i=1,\ldots ,n$|⁠. Let us partition this matrix as

$$\begin{align*} & \widetilde{\textbf{J}} = \left[\begin{array}{@{}cc@{}} A & B \\ C & D \end{array}\right], \quad A = I_{2N}. \end{align*}$$

Within each Newton step we then solve

$$\begin{align*} & \textbf{J} \left[\begin{array}{@{}c@{}} \bar{x}\\ \bar{y} \end{array}\right] = \left[\begin{array}{@{}cc@{}} A & B \\ C & D \end{array}\right] \left[\begin{array}{@{}c@{}} \widetilde x\\ \widetilde y \end{array}\right] = \left[\begin{array}{@{}c@{}} a\\ b \end{array}\right], \end{align*}$$

or equivalently, using the Schur complement |$S = D - C A^{-1} B = D - C B$|⁠,

$$\begin{align}& \left[\begin{array}{@{}cc@{}} A & B \\ O & S \end{array}\right] \left[\begin{array}{c} \widetilde x\\ \widetilde y \end{array}\right] = \left[\begin{array}{@{}c@{}} a \\ b - Ca \end{array}\right]. \end{align}$$

(2.17)

A simplification is obtained by noting that, within the Newton iteration, we have as right-hand side vector the global residual evaluated for the previous iterate,

$ \left[\begin{array}{@{}c@{}} a\\ b \end{array} \right] = \textbf{R}(\textbf{w}^{(k)}). $

If |$\textbf{w}^{(k)}$| is such that |$\widetilde{\textbf{u}}_0 = \textbf{u}_0$| and |$\widetilde{\textbf{v}}_n = \textbf{v}_T$| then |$a=0$|⁠. As a consequence the right-hand side of (2.17) does not require the evaluation of |$Ca$|⁠, which involves the submatrix |$C$| of the Jacobian |$\widetilde{\textbf{J}}$| that may not be explicitly available and may need to be approximated by finite differencing.

Pseudocode for the resulting residual-based Newton method is given in Fig. 1. Note that the matrix computations in steps 6–7 can be performed cheaply: |$B$| is a very sparse matrix of size |$2N$|-by-|$2Nn$| whose only nonzeros correspond to an |$N$|-by-|$N$| identity matrix in its bottom right, while |$K$| can be applied efficiently using its Kronecker representation. The main computational cost of this Newton method is concentrated in step 5, the solution of the Schur complement system. In what follows we will introduce a preconditioner for that problem, allowing the development of an efficient Newton–Krylov solver.

Fig. 1.

Pseudocode for the residual-based Newton solver. The notation follows that of Section 2.2.

2.3 Preconditioner

As a preconditioner for the Schur complement |$S$| arising in (2.17) (and step 5 of Fig. 1) we consider the following matrix:

$$\begin{align*} & P = \widehat D - \widehat C B. \end{align*}$$

Here, |$\widehat D$| and |$\widehat C$| are approximations to |$D$| and |$C$|⁠, respectively, defined as

$$\begin{align*} & \widehat D = \left[\begin{array}{@{}cccc@{}} q_{11} & q_{21} & \cdots & q_{n1} \\ q_{12} & q_{22} & \cdots & q_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ q_{1n} & q_{2n} & \cdots & q_{nn} \\ \end{array}\right]\otimes G \quad \textrm{and} \quad \widehat C = \left[\begin{array}{@{}c@{}} \hat q_{1} \\ \hat q_{2} \\ \vdots \\ \hat q_{n} \\ \end{array}\right]\otimes G, \end{align*}$$

with |$\hat q_i = \sum _{j=0}^n q_{ji}$| and the matrix |$G$| corresponding to

$$\begin{align*} & G: = \left[\begin{array}{@{}cc@{}} \mathcal{D}_u \textbf{f} - M_u & \mathcal{D}_v \textbf{f} \\ \mathcal{D}_u \textbf{g} & \mathcal{D}_v \textbf{g} - M_v \end{array}\right] \in\mathbb{R}^{2N\times 2N}, \end{align*}$$

where the Jacobians of |$\textbf{f}$| and |$\textbf{g}$| are evaluated for the current Newton iterate at some time point in |$[0,T]$|⁠. (In our experiments we use the midpoint |$t=T/2$| as the time evaluation point.) In order to use this preconditioner efficiently we need to be able to apply |$P^{-1}$| cheaply. By the Sherman–Morrison formula we have

$$\begin{align*} & P^{-1} = \widehat D^{-1} + \widehat D^{-1} \widehat C ( I_{2N} - B \widehat D^{-1} \widehat C)^{-1} B \widehat D^{-1} = \underbrace{[ I_{2Nn} + \widehat D^{-1} \widehat C ( I_{2N} - B \widehat D^{-1} \widehat C)^{-1} B ]}_{=:F} \widehat D^{-1}. \end{align*}$$

Now using the facts that

$$\begin{align*} & \widehat D^{-1} = \left[\begin{array}{@{}cccc@{}} q_{11} & q_{21} & \cdots & q_{n1} \\ q_{12} & q_{22} & \cdots & q_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ q_{1n} & q_{2n} & \cdots & q_{nn} \\ \end{array}\right]^{-1}\otimes G^{-1} \end{align*}$$

and hence

$$\begin{align*} & \widehat D^{-1} \widehat C = \left( \left[\begin{array}{@{}cccc@{}} q_{11} & q_{21} & \cdots & q_{n1} \\ q_{12} & q_{22} & \cdots & q_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ q_{1n} & q_{2n} & \cdots & q_{nn} \\ \end{array}\right]^{-1} \left[\begin{array}{@{}c@{}} \hat q_{1} \\ \hat q_{2} \\ \vdots \\ \hat q_{n} \\ \end{array}\right]\right)\otimes \left[\begin{array}{@{}cc@{}} I_N & O \\ O & I_N \end{array}\right] =: \left[\begin{array}{@{}c@{}} \widetilde q_{1} \\ \widetilde q_{2} \\ \vdots \\ \widetilde q_{n} \\ \end{array}\right] \otimes \left[\begin{array}{@{}cc@{}} I_N & O \\ O & I_N \end{array}\right], \end{align*}$$

a simple calculation shows that

$$\begin{align*} & F = I_{2Nn} + \widehat D^{-1} \widehat C (I_{2N} - B \widehat D^{-1} \widehat C)^{-1} B = \left[\begin{array}{cc|cc|cc|cc} I_N & & & & & & O & O\\ & I_N & & & & & O & \widetilde q_1 I_N /(1-\widetilde q_n) \\ \hline & & I_N & & & & O & O \\ & & & I_N & & & O & \widetilde q_2 I_N /(1-\widetilde q_n) \\\hline & & & & \ddots & & \vdots & \vdots \\ & & & & & \ddots & \vdots & \vdots \\\hline & & & & & & I_N & O \\ & & & & & & O & (1+\widetilde q_n/(1-\widetilde q_n)) I_N \\ \end{array}\right]. \end{align*}$$

This allows us to apply |$P^{-1} = F \widehat D^{-1}$| cheaply. In particular we have

$$\begin{align}& P^{-1} \textbf{v} = F \widehat D^{-1} \textbf{v} = F \, {\texttt{vec}} \left(G^{-1} {\texttt{mat}}(\textbf{v}) \left[\begin{array}{@{}cccc@{}} q_{11} & q_{21} & \cdots & q_{n1} \\ q_{12} & q_{22} & \cdots & q_{n2} \\ \vdots & \vdots & \ddots & \vdots \\ q_{1n} & q_{2n} & \cdots & q_{nn} \\ \end{array}\right]^{-\textrm{T}}\right), \end{align}$$

(2.18)

with the mat operation transforming a vector |$\textbf{v}\in \mathbb{R}^{2N n}$| into a |$2N$|-by-|$n$| matrix, and vec being the inverse operation. All computational costs involved in applying the preconditioner are linear with respect to |$N$|⁠, with the exception of creating a factorization of the |$2N$|-by-|$2N$| Jacobian matrix |$G$|⁠, which needs to be done once per Newton iteration.

2.4 Newton–Krylov implementation

We do not wish to form the Schur complement |$S = D - CB$| explicitly, but we will instead approximate its action within a Krylov method using finite differencing. The resulting Newton–Krylov method is standard and described in detail, e.g., in Kelley (2003, Chapter 3). We will therefore only outline the key parts required for its implementation.

Recall from (2.16) and the subsequent derivations that the matrix |$D$|⁠, which is the lower-right |$2Nn$|-by-|$2Nn$| submatrix of |$\widetilde{\textbf{J}}= \textbf{J}\cdot K$|⁠, is the same as the lower-right submatrix of |$\textbf{J} = \mathcal{D}\textbf{R}(\textbf{w})$|⁠. We can therefore obtain an approximation to the matrix–vector product |$Dx$| as a finite difference

$$\begin{align*} & \left[ \frac{\textbf{R}(\textbf{w}+ h [0_{2N\times 1} ; x]) - \textbf{R}(\textbf{w})}{h} \right]_{2N+1:2N(n+1)} \approx Dx, \end{align*}$$

where again ‘;’ corresponds to row-wise concatenation and the MATLAB operator ‘:’ was used to extract vector entries by their indices. For the appropriate choice of the parameter |$h$| in floating point arithmetic we refer to Kelley (2003, Section 3.2.1). Similarly, the action of the lower-left |$2Nn$|-by-|$2N$| submatrix |$C$| of |$\widetilde{\textbf{J}} = \textbf{J}\cdot K$| on a vector |$Bx$| can be approximated as

$$\begin{align*} & \left[ \frac{\textbf{R}(\textbf{w}+ h K [Bx ; 0_{2Nn\times 1}]) - \textbf{R}(\textbf{w})}{h} \right]_{2N+1:2N(n+1)} \approx C(Bx). \end{align*}$$

Overall, three evaluations of |$\textbf{R}(\cdot )$| are required to obtain an approximation of |$Sx = Dx - C(Bx)$|⁠. If multiple such approximations with different vectors |$x$| are required consecutively, e.g., within a Krylov iteration, then the evaluation of |$\textbf{R}(\textbf{w})$| can be reused for all of them.

The left-preconditioned system |$P^{-1} S x = P^{-1} b$| is solved using GMRES (Saad & Schultz, 1986). For these inner solves the standard Eisenstat–Walker residual stopping criterion based on the norm of the residual, |$\|\textbf{R}(\textbf{w}^{(k)})\|$|⁠, is used (Eisenstat & Walker, 1996); see also Kelley (2003, Section 3.2.3).

2.5 Stopping criterion

There are several options to stop our method, with the most obvious one being a norm criterion for the residual evaluated in step 3 of the algorithm in Fig. 1: stop at Newton iteration |$k$| if |$\|\textbf{R}(\textbf{w}^{(k)}) \| \leq {\texttt{tol}}\,{\|\textbf{R}(\textbf{w}^{(0)}) \|}$| for some user-specified tolerance |${\texttt{tol}}$|⁠.

An alternative option is to estimate the error in the iterate

$$\begin{align*} &\textbf{w}^{(k)} = \left[\begin{array}{c} \widetilde{\textbf{u}}^{(k)}\\ \widetilde{\textbf{v}}^{(k)} \end{array}\right] \end{align*}$$

from the difference |$\textbf{w}^{(k+s)} - \textbf{w}^{(k)}$| for some integer lag parameter |$s\geq 1$|⁠. This strategy gives good results if the method converges rapidly, as in this case |$\textbf{w}^{(k+s)}$| will be an ‘accurate’ (relative to |$\textbf{w}^{(k)}$|⁠) approximation for the vector to which the algorithm converges. As the state (⁠|$\widetilde{\textbf{u}}^{(k)}$|⁠) and adjoint (⁠|$\widetilde{\textbf{v}}^{(k)}$|⁠) components in |$\textbf{w}^{(k)}$| often vary significantly in magnitude it is advisable to treat them separately when estimating the error. We therefore suggest using the following estimate for relative error of |$\widetilde{\textbf{u}}^{(k)}$|⁠:

$$\begin{align*} & {\texttt{errest}}_{\textbf{u}}^{(k)} = \frac{\max_{j=0,1,\ldots,n} \| {\widetilde{\textbf{u}}_{j}^{(k+s)}} - {{\widetilde{\textbf{u}}}_j^{(k)}} \|_\infty}{\max_{j=0,1,\ldots,n} \| {\widetilde{\textbf{u}}_j^{(k+s)}} \|_\infty}, \end{align*}$$

and an analogous estimate for the relative error of |$\widetilde{\textbf{v}}^{(k)}$|⁠. Together with an estimate for the expected stagnation error level (dependent on the space and time discretization) these estimators provide a practical stopping criterion. We will illustrate this in Section 3.

3. Numerical experiments

Our numerical experiments are guided by two nonlinear model problems, the two-dimensional Schlögl problem and a two-dimensional reaction–diffusion problem, both posed on the spatial domain |$\varOmega =(-1,1)^2$| and with analytic solutions given in Appendix  A (with |$\alpha =0.1$| for the Schlögl problem). The experiments have been performed in Matlab 2019a on a Windows 10 laptop with 8 GB RAM and an Intel(R) Core(TM) i7-8650U CPU running at 2 GHz. A Matlab implementation of our method, including the preconditioner proposed in Section 2.3, as well as scripts reproducing these experiments, can be downloaded from https://github.com/nla-group/pdeoptim. We use Chebfun (Driscoll et al., 2014) to generate the spectral collocation matrices in space and time.

Fig. 2.

Eigenvalues of the Schur complement |$S$| and its preconditioned counterpart |$P^{-1}S$| for the two-dimensional Schlögl problem with finite difference discretization. Only the first three Newton iterations are shown since the subsequent ones look qualitatively similar. The condition numbers of the Schur complement |$S$| in the first seven iterations (until stagnation of the Newton method as determined in Section 3.2) are |$[2.5, 2.3, 2.4, 2.7, 3.0, 3.0, 3.0]\times 10^3$|⁠, and those of the preconditioned Schur complement |$P^{-1}S$| are |$[1.8,5.3,14.2,24.3,27.0,27.4,27.4]$|⁠.

Fig. 3.

Error in the computed state (left) and adjoint (right) variables for the Schlögl problem using a spectral space discretization, for varying regularization parameter |$\beta $|⁠.

3.1 Visual inspection of the preconditioner

We first explore visually the properties of our preconditioner from Section 2.3 with the help of the Schlögl problem, discretized with finite differences in space. The following parameters have been chosen: final time |$T=2$|⁠, regularization parameter |$\beta =0.01$|⁠, |$n_x=16$| equidistant interior grid points in each spatial dimension and degree |$n=5$| Chebyshev collocation on |$[0,T]$|⁠; that is, we use the time collocation points

$$\begin{align*} & t_j = T \, \frac{1-\cos(j\pi /n)}{2}, \quad j=0,1,\ldots,n. \end{align*}$$

In Fig. 2 we show the eigenvalues of both the unpreconditioned (left) and preconditioned (right) Schur complements, |$S$| and |$P^{-1}S$|⁠, at the first three Newton iterations. We observe that |$P^{-1}S$| has tightly clustered eigenvalues away from the origin, which serves as an indication (though not a guarantee) that fast GMRES convergence can be expected for this preconditioned system. We also show, at the bottom left of each eigenvalue plot, the condition numbers of |$S$| and |$P^{-1}S$|⁠, respectively. Again, the reduced condition number of the preconditioned system indicates that the preconditioner is effective. We will now turn to some concrete numerical tests to confirm this.

3.2 Schlögl problem, spectral in space

To test the achievable accuracy of our solver to the extreme, as well as the robustness with respect to the regularization parameter |$\beta $|⁠, we consider the two-dimensional Schlögl example with a spectral space discretization. We use |$n_x=12$| interior grid points in each spatial dimension and degree |$n=10$| Chebyshev collocation on |$[0,T\!=\!2]$| (i.e., |$n+1\!=\!11$| time collocation points). The errors of the computed state and adjoint solutions are obtained by taking the absolute deviation on the spectral collocation grid over all time steps, i.e.,

$$\begin{align*} & {\texttt{err}}_{\textbf{u}} = \frac{\max_{j=0,1,\ldots,n} \| \textbf{u}_j - {\widetilde{\textbf{u}}_j} \|_\infty}{\max_{j=0,1,\ldots,n} \| {\textbf{u}_j} \|_\infty}, \end{align*}$$

where |${\widetilde{\textbf{u}}_j}$| is the computed approximation to |${\textbf{u}}(t_j)$|⁠, the latter of which is known analytically. Likewise, the error in the adjoint, |$\texttt{err}_{\textbf{v}}$|⁠, is computed. The results are shown in Fig. 3, with the number of Newton iterations varying from |$0$| to |$10$|⁠. The zeroth Newton iteration corresponds to using the constant initial guess |$\widetilde{\textbf{u}}_j:= \textbf{u}_0$| and |$\widetilde{\textbf{v}}_j:= \textbf{v}_T$| for |$j=0,1,\ldots ,n$|⁠.

Figure 3 also indicates an error level of |$10^{-10}$|⁠, which we have used in Table 1 to evaluate the number of outer Newton iterations and inner GMRES iterations, as well as timings, required to achieve that accuracy. More precisely, we consider the algorithm to have converged to the target accuracy if

$$\begin{align*} &\texttt{err}:= \max \{ {\texttt{err}}_{\textbf{u}}, {\texttt{err}}_{\textbf{v}} \} \leq 10^{-10}.\end{align*}$$

Table 1

Outer Newton and inner GMRES iterations required to solve the Schlögl problem using a spectral space discretization, to an accuracy below |$10^{-10}$|

|$\beta $|	Outer iterations	Inner iterations	Total inner	Final accuracy	Total time (s)
|$1e\!-\!1$|	4	4, 14, 18, 30	66	|$6.92e\!-\!12$|	0.405
|$1e\!-\!2$|	7	1, 4, 6, 8, 12, 16, 25	72	|$6.27e\!-\!12$|	0.398
|$1e\!-\!3$|	8	1, 2, 4, 6, 9, 12, 20, 35	89	|$3.35e\!-\!13$|	0.480
|$1e\!-\!4$|	8	1, 2, 4, 6, 9, 7, 12, 30	71	|$8.49e\!-\!13$|	0.441
|$1e\!-\!5$|	8	1, 2, 4, 6, 9, 9, 10, 25	66	|$1.51e\!-\!12$|	0.409

As expected the total computation time is roughly proportional to the number of inner GMRES iterations performed, and this number increases steadily as the Newton iterates approach the sought minimum of the residual norm. Overall, the number of inner iterations and total time does not seem to depend much on the value of the regularization parameter |$\beta $|⁠, and we can solve complex optimal control problems in a fraction of a second using GMRES, up to spectral accuracy, for every |$\beta $| tested.

3.3 Schlögl problem, finite differences in space

We now consider the Schlögl example with a finite difference discretization in space, while still using a spectral method in time. This example intends to demonstrate that mixing these two discretizations is a viable option, as the spectral time discretization helps to keep the number of solution vectors |$\widetilde{\textbf{u}}_j,\widetilde{\textbf{v}}_j\in \mathbb{R}^N$||$(j=0,1,\ldots ,n)$| small, thereby reducing the overall memory consumption and arithmetic cost. This is in contrast to previous approaches that have used a backward Euler approximation in time and hence required a much larger number of time steps, at the cost of increased linear algebra complexity.

In this test we use degree |$n=5$| Chebyshev collocation in time and vary the number of equidistant spatial interior grid points |$n_x=4,8,16,32,64,128$| in both coordinate directions. The time interval is again |$[0,T=2]$|⁠, and the errors of the computed state and adjoint solutions are obtained as before.

The results are presented in Fig. 4 and Table 2. We first observe that as the number of spatial grid points is doubled the finally attained accuracy improves by a factor of about |$4$|⁠. This is consistent with the second-order centred finite difference scheme and also indicates that the time discretization is sufficiently fine. We further observe that, independently of |$n_x$|⁠, the error of the space discretization is reached after about six Newton iterations. By inspecting the number of inner GMRES iterations in Table 2 we find that the method is robust with respect to refinements in space, testament to an effective preconditioner. As a consequence the total computation time scales close to linearly in the number of overall spatial grid points, |$n_x^2$|⁠. We highlight that if |$n_x$| becomes very large the factorization of the matrix |$G$| in (2.18) at each Newton iteration could become a computational bottleneck, although our results demonstrate that high accuracy can be achieved for moderate |$n_x$|⁠, thus mitigating this issue.

Fig. 4.

Error in the computed state (left) and adjoint (right) variables for the Schlögl problem using a finite difference space discretization, with a varying number of grid points |$n_x$|⁠.

Table 2

Outer Newton and inner GMRES iterations required to solve the Schlögl problem using a finite difference space discretization

|$n_x$|	Outer iterations	Inner iterations	Total inner	Final accuracy	Total time (s)
4	6	1, 3, 4, 6, 9, 10	33	|$4.13e\!-\!03$|	0.0298
8	6	1, 3, 4, 6, 9, 9	32	|$1.22e\!-\!03$|	0.0466
16	6	1, 3, 4, 6, 9, 9	32	|$3.37e\!-\!04$|	0.0923
32	6	1, 3, 4, 6, 9, 9	32	|$8.92e\!-\!05$|	0.388
64	6	1, 3, 4, 6, 9, 9	32	|$2.35e\!-\!05$|	1.43
128	6	1, 3, 4, 6, 9, 9	32	|$7.67e\!-\!06$|	6.71

We also show in Fig. 5 the error at each Newton iteration when |$n_x=32$| interior grid points are used in space, as well as the error estimate discussed in Section 2.5, using the lag parameter |$s=1,2,3$|⁠. Note how the error estimates follow the actual error very closely until the method stagnates on the level of the spatial/temporal discretization error. Together with estimates for the expected discretization error this error estimate is a useful stopping criterion for our algorithm, and the user can specify a desired tolerance based on a reasonable expected discretization error given the |$n$| and |$n_x$| chosen.

Fig. 5.

Error and error estimates for the computed state (left) and adjoint (right) for the Schlögl problem using a finite difference space discretization, with |$n_x = 32$| interior grid points. The lag parameter for the error estimation is varied between |$s=1,2,3$|⁠.

Fig. 6.

Error in the computed state (left) and adjoint (right) variables for the reaction–diffusion problem using a finite difference space discretization, with a varying number of grid points |$n_x$|⁠.

3.4 Reaction–diffusion problem, finite differences in space

We now turn our attention to a nonlinear reaction–diffusion problem discretized using finite differences in space. We use degree |$n=10$| Chebyshev collocation on the time interval |$[0,T=1]$| and parameters |$\beta _y = 1$|⁠, |$\beta _z = 1$|⁠, |$\beta _c = 0.01$|⁠, |$D_1 = 0.5$|⁠, |$D_2 = 1$|⁠, |$k_1 = 1$|⁠, |$k_2 = 1$|⁠, |$\gamma _1 = 0.4$|⁠, |$\gamma _2 = 0.6$|⁠, noting that the method demonstrates robustness with respect to reasonable modifications of these parameters.

The results of the Newton–Krylov method are shown in Fig. 6 and Table 3. We can see that the error of the space discretization is reached after about five Newton iterations and that a slightly larger |$n_x$| is required to obtain a satisfactory numerical solution than for the Schögl problem. As shown in Fig. 6 the very low value of |$n_x=4$| (with |$n_x$| including boundary points for this problem, due to the imposition of Neumann boundary conditions) leads to worse error properties in the computed adjoint than an arbitrary initial guess; however, increasing the number of spatial grid points gives very high accuracy once again. We may see from Table 3 that, for sufficiently large |$n_x$|⁠, doubling this value leads to an improved accuracy by a factor of roughly |$4$| once more, as is expected. We also observe that the computation time scales close to linearly with respect to number of spatial grid points, with the only nonlinearity arising from the factorization of the matrix |$G$|⁠: this nonlinear scaling becomes visible for a smaller |$n_x$| than for the Schlögl problem, as |$G$| is a larger matrix for the reaction–diffusion problem; however, all computation times are low considering the high complexity of the problem being solved.

Table 3

Outer Newton and inner GMRES iterations required to solve the reaction–diffusion problem using a finite difference space discretization

|$n_x$|	Outer iterations	Inner iterations	Total inner	Final accuracy	Total time (s)
4	4	12, 23, 33, 60	128	|$3.73e\!+\!00$|	0.556
8	4	18, 28, 39, 68	153	|$1.83e\!-\!01$|	0.915
16	4	18, 30, 38, 60	146	|$1.84e\!-\!02$|	2.62
32	4	18, 31, 37, 55	141	|$3.23e\!-\!03$|	7.90
64	4	18, 31, 36, 50	135	|$6.63e\!-\!04$|	38.4
128	4	18, 32, 36, 51	137	|$1.61e\!-\!04$|	221

A visualization of the computed solution components |$y$| and |$q$| at time |$t=0.5$| using a finite difference space discretization with |$n_x=128$| spatial grid points is shown in Fig. 7. We conclude this section by highlighting that we can also solve this reaction–diffusion control example using a spectral method in space, and indeed our approach can in principle be applied to any space discretization, with all that is required being the matrices arising from the discretization of the linearized PDE operators at a given time step. For instance, one may apply a finite element discretization in space, and take advantage of the finite element method’s greater flexibility when it comes to grid and mesh structure.

Fig. 7.

Computed solution of the reaction–diffusion problem at time |$t=0.5$|⁠, using a finite difference space discretization with |$n_x=128$| spatial grid points.

4. Concluding remarks

In this paper we have derived a Newton–Krylov method to allow the spectral-in-time solution of nonlinear time-dependent PDE-constrained optimization problems, suggested a preconditioner to be applied within this method and validated our approach numerically on test problems relating to the optimal control of the Schlögl equation as well as a reaction–diffusion system. Thanks to the spectral time discretization, as opposed to more commonly used low-order time discretizations (such as backward Euler in Pearson & Stoll, 2013; Pearson et al., 2012; Stoll & Breiten, 2015), our method allows for a fast and accurate solution of such problems with relatively low memory requirements (as the solution is only collocated at a small number of time nodes). The method is also amenable to a range of strategies for the space discretization, and to apply the method the user only needs to provide matrices arising from their preferred discretization of the linearized PDE operators at a given time step. Possible extensions of this approach include the solution of more sophisticated systems of PDEs or integro-PDEs, problems with additional algebraic constraints on the state and control variables, and connecting with time-parallel methods for parabolic control problems, for instance the work presented in Götschel & Minion (2019). The derivation of more sophisticated preconditioning strategies and stopping criteria for the Newton–Krylov solver, taking account of the specific structures of the PDE operators involved, would also be of interest.

Acknowledgements

We are grateful to two anonymous referees for many insightful comments and suggestions that have improved this work.

Funding

Engineering and Physical Sciences Research Council (EP/S027785/1 to J.W.P.); Fellowships from the Alan Turing Institute (EP/N510129/1 to S.G. and J.W.P.).

A. Test problems with analytic solutions

Our first test problem involves the optimal control of the Schlögl equation, as given in (2.1). The continuous optimality conditions (2.3)–(2.5) are solved by the following functions:

$$\begin{align*} \ y={}&\alpha\left(\frac{4}{(d\pi^{2}-4)\beta}e^{T}-\frac{4}{d\pi^{2}\beta}e^{t}\right)~\prod_{k=1}^{d}\cos\left(\frac{\pi{}x_{k}}{2}\right), \\ \ p={}&\alpha(e^{T}-e^{t})~\prod_{k=1}^{d}\cos\left(\frac{\pi{}x_{k}}{2}\right),\quad\quad\left[c=\frac{1}{\beta}p\right] \\ \ \widehat{y}={}&\alpha\left(\left[\frac{4}{(d\pi^{2}-4)\beta}+\frac{d\pi^{2}}{4}-1\right]e^{T}+\left[2-\frac{4}{d\pi^{2}\beta}-\frac{d\pi^{2}}{4}\right]e^{t}\right)~\prod_{k=1}^{d}\cos\left(\frac{\pi{}x_{k}}{2}\right) \\ \ &+3\alpha^{3}\left(\frac{4}{(d\pi^{2}-4)\beta}e^{T}-\frac{4}{d\pi^{2}\beta}e^{t}\right)^{2}(e^{T}-e^{t})~\left[\prod_{k=1}^{d}\cos\left(\frac{\pi{}x_{k}}{2}\right)\right]^{3}, \\ \ f={}&\alpha^{3}\left(\frac{4}{(d\pi^{2}-4)\beta}e^{T}-\frac{4}{d\pi^{2}\beta}e^{t}\right)^{3}~\left[\prod_{k=1}^{d}\cos\left(\frac{\pi{}x_{k}}{2}\right)\right]^{3}, \\ \ y_{0}={}&\alpha\left(\frac{4}{(d\pi^{2}-4)\beta}e^{T}-\frac{4}{d\pi^{2}\beta}\right)~\prod_{k=1}^{d}\cos\left(\frac{\pi{}x_{k}}{2}\right), \\ \ h={}&0, \end{align*}$$

where |$d\in \{2,3\}$| is the dimension of the problem, solved on the space–time domain |$(-1,1)^{d}\times (0,T)$|⁠, and |$\alpha \ne 0$| denotes an arbitrary constant.

Our second test problem involves the reaction–diffusion control problem (2.2). The continuous optimality conditions (2.6)–(2.8) are solved by the following functions:

$$\begin{align*} \ y=&\,\frac{1}{8}(e^{T}-e^{t})~\sum_{k=1}^{d}x_{k}^{2}, \\ \ z=&\,\frac{1}{8}(e^{T}-e^{t}), \\ \ p=&\,\frac{\beta_{c}D_{1}}{4}(e^{T}-e^{t}),\quad\quad\left[c=\frac{1}{\beta_c}p~\ \textrm{on}\ \partial\varOmega\times(0,T)\right]\\q=&\,\frac{\beta_{c}D_{1}}{4}(e^{T}-e^{t})~\prod_{k=1}^{d}\left[1+\cos(\pi{}x_{k})\right], \end{align*}$$

with |$d\in \{2,3\}$| again the dimension of the problem, solved on the space–time domain |$(-1,1)^{d}\times (0,T)$|⁠.