Solving Parameter Estimation Problems with Discrete Adjoint Exponential Integrators

TL;DR

This paper develops discrete adjoint formulas for exponential integrators used in inverse problems, enabling efficient gradient computation and optimization for complex PDE models.

Contribution

It derives the discrete adjoint formulas for EPIRK-W exponential integrators, facilitating their use in inverse problems with efficient adjoint code generation.

Findings

Validated adjoint formulas with Lorenz-96 model

Demonstrated efficiency in inverse problem solving

Enabled Jacobian approximation independence from model state

Abstract

The solution of inverse problems in a variational setting finds best estimates of the model parameters by minimizing a cost function that penalizes the mismatch between model outputs and observations. The gradients required by the numerical optimization process are computed using adjoint models. Exponential integrators are a promising family of time discretizations for evolutionary partial differential equations. In order to allow the use of these discretizations in the context of inverse problems adjoints of exponential integrators are required. This work derives the discrete adjoint formulae for a W-type exponential propagation iterative methods of Runge-Kutta type (EPIRK-W). These methods allow arbitrary approximations of the Jacobian while maintaining the overall accuracy of the forward integration. The use of Jacobian approximation matrices that do not depend on the model state avoids the complex calculation of Hessians in the discrete adjoint formulae, and allows efficient adjoint code generation via algorithmic differentiation. We use the discrete EPIRK-W adjoints to solve inverse problems with the Lorenz-96 model and a computational magnetics benchmark test. Numerical results validate our theoretical derivations.

Tables

Table 1. Table 1: Measurements compared to background simulation and simulation after data assimilation in Wbs − 1 superscript Wbs 1 \mathrm{Wb}\mathrm{s}^{-1} .

	Measurement	Background	Optimum
$z_1$	7.656	7.450	7.660
$z_2$	-0.556	-0.377	-0.554
$z_3$	-9.565	-9.414	-9.571
$z_4$	-6.092	-4.004	-6.092
$z_5$	4.472	4.098	4.470

Equations

$$\frac{d\mathbf{y}}{dt}=\mathbf{f}(\mathbf{y}(t),\mathbf{\theta }),\mathbf{y}(t_0)={\mathbf{y}}_{\mathrm{ini}}(\mathbf{\theta }),t_0\le t\le t_F,$$

$${\mathbf{\theta }}_{\text{textsc}\mathbf{map}}=\arg \underset{\mathbf{\theta }}{\min }$$

$${\mathbf{y}}^{\prime }=\mathbf{f}(\mathbf{y}(t),\mathbf{\theta }),\mathbf{y}(t_0)={\mathbf{y}}_{\mathrm{ini}}(\mathbf{\theta }),t_0\le t\le t_F,$$

$$\begin{array}{rl}{\mathbf{Y}}_{n,i}& ={\mathbf{y}}_n+a_{i,1}{\mathbf{\psi }}_{i,1}(g_{i,1}h{\mathbf{J}}_n)h\mathbf{f}({\mathbf{y}}_n)\\ & +\sum _{j=2}^i{a}_{i,j}{\mathbf{\psi }}_{i,j}(g_{i,j}h{\mathbf{J}}_n)h{\Delta }^{(j-1)}\mathbf{r}({\mathbf{y}}_n),i=1\dots s-1,\end{array}$$

$${\mathbf{y}}_{n+1}$$

$${\mathbf{\psi }}_{i,j}(\mathbf{Z})=\sum _{k=1}^s{p}_{i,j,k}{\mathbf{\phi }}_k(\mathbf{Z}),$$

$${\phi }_k(z)={\int }_0^1{e}^{z(1-\theta )}\frac{{\theta }^{k-1}}{(k-1)!}d\theta =\sum _{i=0}^{\infty }\frac{z^i}{(i+k)!},k=1,2,\dots ,$$

$${\phi }_0(z)=e^z;{\phi }_{k+1}(z)=\frac{{\phi }_k(z)-1/k!}{z},k=1,2,\dots ;{\phi }_k(0)=\frac{1}{k!}.$$

$$\begin{array}{rl}\mathbf{r}(\mathbf{y})& :=\mathbf{f}(\mathbf{y})-\mathbf{f}({\mathbf{y}}_n)-{\mathbf{J}}_n(\mathbf{y}-{\mathbf{y}}_n),\\ {\Delta }^{(0)}\mathbf{r}(\mathbf{y})& :=\mathbf{r}(\mathbf{y}),{\mathbf{Y}}_{n,0}:={\mathbf{y}}_n,\\ {\Delta }^{(j)}\mathbf{r}({\mathbf{Y}}_{n,i})& ={\Delta }^{(j-1)}\mathbf{r}({\mathbf{Y}}_{n,i+1})-{\Delta }^{(j-1)}\mathbf{r}({\mathbf{Y}}_{n,i}),j\ge 1,\\ {\Delta }^{(j)}\mathbf{r}({\mathbf{y}}_n)& ={\Delta }^{(j)}\mathbf{r}({\mathbf{Y}}_{n,0})={\Delta }^{(j-1)}\mathbf{r}({\mathbf{Y}}_{n,1})-{\Delta }^{(j-1)}\mathbf{r}({\mathbf{Y}}_{n,0}),j\ge 1.\end{array}$$

$$\begin{array}{rl}{\mathbf{Y}}_{n,i}& ={\mathbf{y}}_n+a_{i,1}{\mathbf{\psi }}_{i,1}(g_{i,1}h{\mathbf{T}}_n)h\mathbf{f}({\mathbf{y}}_n)\\ & +\sum _{j=2}^i{a}_{i,j}{\mathbf{\psi }}_{i,j}(g_{i,j}h{\mathbf{T}}_n)h{\Delta }^{(j-1)}\mathbf{r}({\mathbf{y}}_n),i=1\dots s-1,\end{array}$$

$${\mathbf{y}}_{n+1}$$

$${\mathbf{A}}_{n,i,j}:=a_{i,j}{\mathbf{\psi }}_{i,j}(g_{i,j}h{\mathbf{T}}_n),{\mathbf{B}}_{n,j}:=b_j{\mathbf{\psi }}_{s,j}(g_{s,j}h{\mathbf{T}}_n).$$

$${\mathbf{Y}}_{n,i}$$

$${\mathbf{y}}_{n+1}$$

$${\mathbf{\theta }}_{\text{textsc}\mathbf{map}}=\arg \underset{\mathbf{\theta }}{\min }$$

$${\mathbf{Y}}_{n,i}={\mathbf{y}}_n+{\mathbf{A}}_{n,i,1}h\mathbf{f}({\mathbf{y}}_n,\mathbf{\theta })$$

$$+\sum _{j=2}^i{\mathbf{A}}_{n,i,j}h{\Delta }^{(j-1)}\mathbf{r}({\mathbf{y}}_n,\mathbf{\theta }),i=1,\dots ,s-1.$$

$${\mathbf{y}}_{n+1}={\mathbf{y}}_n+{\mathbf{B}}_{n,1}h\mathbf{f}({\mathbf{y}}_n,\mathbf{\theta })$$

$$+\sum _{j=2}^s{\mathbf{B}}_{n,j}h{\Delta }^{(j-1)}\mathbf{r}({\mathbf{y}}_n,\mathbf{\theta }),0\le n\le N-1,$$

$${\mathbf{y}}_0={\mathbf{y}}_{\mathrm{ini}}(\mathbf{\theta }).$$

$$\tilde{\Phi }(\mathbf{\theta })$$

$${\nabla }_{\mathbf{\theta }}\tilde{\Phi }(\mathbf{\theta })$$

$${\tilde{\mathbf{y}}}^{\prime }={\left[\begin{array}{c}\mathbf{y}\\ \mathbf{\tau }\end{array}\right]}^{\prime }=\left[\begin{array}{c}\mathbf{f}(\mathbf{y},\mathbf{\tau })\\ 0\end{array}\right]=\tilde{\mathbf{f}}(\tilde{\mathbf{y}}),\tilde{\mathbf{y}}(t_0)=\left[\begin{array}{c}\mathbf{y}(t_0)\\ \mathbf{\tau }(t_0)\end{array}\right]=\left[\begin{array}{c}{\mathbf{y}}_{\mathrm{ini}}\\ \mathbf{\theta }\end{array}\right]={\tilde{\mathbf{y}}}_{\mathrm{ini}}.$$

$${\mathbf{Y}}_n={\left({\mathbf{Y}}_{n,1}^T,\dots ,{\mathbf{Y}}_{n,s-1}^T\right)}^T,\hat{\mathbf{Y}}={\left({\mathbf{Y}}_0^T,\dots ,{\mathbf{Y}}_{N-1}^T\right)}^T,\hat{\mathbf{y}}={\left({\mathbf{y}}_{\mathrm{ini}}^T,\dots ,{\mathbf{y}}_N^T\right)}^T.$$

$$\mathcal{L}(\hat{\mathbf{y}},\hat{\mathbf{Y}},\hat{\boldsymbol{\lambda}},\hat{\boldsymbol{\Lambda}},\boldsymbol{\theta})=\sum_{k=0}^{N}q_{k}\big{(}\mathbf{y}_{k},\boldsymbol{\theta}\big{)}-\mathcal{L}_{\textnormal{ext}}\big{(}\hat{\mathbf{y}},\hat{\mathbf{Y}},\hat{\boldsymbol{\lambda}},\boldsymbol{\theta}\big{)}-\mathcal{L}_{\rm int}\big{(}\hat{\mathbf{y}},\hat{\mathbf{Y}},\hat{\boldsymbol{\Lambda}},\boldsymbol{\theta}\big{)},$$

$${\mathcal{L}}_{\text{ext}}$$

$$+\sum _{k=0}^{N-1}{\mathbf{\lambda }}_{k+1}^T\left({\mathbf{y}}_{k+1}-{\mathbf{y}}_k-{\mathbf{B}}_{k,1}h\mathbf{f}({\mathbf{y}}_k,\mathbf{\theta })-\sum _{j=2}^s{\mathbf{B}}_{k,j}h{\Delta }^{(j-1)}\mathbf{r}({\mathbf{y}}_k,\mathbf{\theta })\right),$$

$${\mathcal{L}}_{\mathrm{int}}$$

$${\Delta }^{(j-1)}\mathbf{r}({\mathbf{y}}_k,\mathbf{\theta })=\sum _{\ell =0}^{j-1}{C}_{\ell ,j}\mathbf{r}({\mathbf{Y}}_{k,j-1-\ell },\mathbf{\theta }),C_{\ell ,j}=(-1{)}^{\ell }(\genfrac{}{}{0ex}{}{j-1}{\ell }),$$

$$\mathbf{r}({\mathbf{Y}}_{k,j-1-\ell },\mathbf{\theta })=\mathbf{f}({\mathbf{Y}}_{k,j-1-\ell },\mathbf{\theta })-\mathbf{f}({\mathbf{y}}_k,\mathbf{\theta })-{\mathbf{T}}_k\cdot ({\mathbf{Y}}_{k,j-1-\ell }-{\mathbf{y}}_k).$$

$${\mathcal{L}}_{\text{ext}}$$

$$\displaystyle-\displaystyle\sum_{j=2}^{s}\mathbf{B}_{k,j}\,h\bigg{(}\sum_{\ell=0}^{j-1}C_{\ell,j}\big{(}\mathbf{f}(\mathbf{Y}_{k,j-1-\ell},\boldsymbol{\theta})-\mathbf{f}(\mathbf{y}_{k},\boldsymbol{\theta})-\mathbf{T}_{k}(\mathbf{Y}_{k,j-1-\ell}-\mathbf{y}_{k})\big{)}\bigg{)}\Bigg{)},$$

$${\mathcal{L}}_{\mathrm{int}}$$

$$\displaystyle-\displaystyle\sum_{j=2}^{m}\mathbf{A}_{k,m,j}\,h\bigg{(}\sum_{\ell=0}^{j-1}C_{\ell,j}\big{(}\mathbf{f}(\mathbf{Y}_{k,j-1-\ell},\boldsymbol{\theta})-\mathbf{f}(\mathbf{y}_{k},\boldsymbol{\theta})-\mathbf{T}_{k}(\mathbf{Y}_{k,j-1-\ell}-\mathbf{y}_{k})\big{)}\bigg{)}\Bigg{)}.$$

$$\displaystyle\frac{\partial\mathcal{L}_{\textnormal{ext}}}{\partial\mathbf{y}_{n}}=\boldsymbol{\lambda}_{n}^{T}\cdot\frac{\partial}{\partial\mathbf{y}_{n}}\Bigg{\{}\mathbf{y}_{n}-\mathbf{y}_{n-1}-\mathbf{B}_{n-1,1}\,h\mathbf{f}(\mathbf{y}_{n-1},\boldsymbol{\theta})$$

$$\displaystyle\quad-\sum_{j=2}^{s}\mathbf{B}_{n-1,j}\;h\,\sum_{\ell=0}^{j-1}C_{\ell,j}\bigg{[}\mathbf{f}(\mathbf{Y}_{n-1,j-1-\ell},\boldsymbol{\theta})-\mathbf{f}(\mathbf{y}_{n-1},\boldsymbol{\theta})-\mathbf{T}_{n-1}(\mathbf{Y}_{n-1,j-1-\ell}-\mathbf{y}_{n-1})\bigg{]}\Bigg{\}}$$

$$\displaystyle\quad+\boldsymbol{\lambda}_{n+1}^{T}\cdot\frac{\partial}{\partial\mathbf{y}_{n}}\Bigg{\{}\mathbf{y}_{n+1}-\mathbf{y}_{n}-\mathbf{B}_{n,1}h\mathbf{f}(\mathbf{y}_{n},\boldsymbol{\theta})$$

$$\displaystyle\quad-\sum_{j=2}^{s}\mathbf{B}_{n,j}\;h\sum_{\ell=0}^{j-1}C_{\ell,j}\bigg{[}\mathbf{f}(\mathbf{Y}_{n,j-1-\ell},\boldsymbol{\theta})-\mathbf{f}(\mathbf{y}_{n},\boldsymbol{\theta})-\mathbf{T}_{n}(\mathbf{Y}_{n,j-1-\ell}-\mathbf{y}_{n})\bigg{]}\Bigg{\}}$$

$$\displaystyle=\boldsymbol{\lambda}_{n}^{T}\big{\{}\mathbf{I}-0\,-0-0\big{\}}+\boldsymbol{\lambda}_{n+1}^{T}\cdot\Bigg{\{}0-\mathbf{I}\,-\mathbf{B}_{n,1}\,h\mathbf{J}(\mathbf{y}_{n},\boldsymbol{\theta})$$

$$\displaystyle\quad-\sum_{j=2}^{s}\mathbf{B}_{n,j}\,h\sum_{\ell=0}^{j-2}C_{\ell,j}\big{[}0-\mathbf{J}(\mathbf{y}_{n},\boldsymbol{\theta})-\mathbf{T}_{n}(0-\mathbf{I})\big{]}\Bigg{\}}$$

Full text

{classcode}

34H05;34K29;34K35

1. 1 Introduction

2. 2 Formulation of the Discrete Parameter Estimation Problem

3. 2.1 EPIRK-W Time Discretization Methods

4. 2.2 The Discrete Formulation of the Parameter Estimation Problem

5. 2.3 4D-Var Data Assimilation

6. 3 First Order Optimality Conditions and the Discrete EPIRK-W Adjoint

7. 3.1 The Discrete Lagrangian

8. 3.2 Derivation of the Discrete Adjoint Equation

9. 3.2.1 Computing $\partial\mathcal{L}_{\textnormal{ext}}/\partial\mathbf{y}_{n}$

10. 3.2.2 Computing $\partial\mathcal{L}_{\textnormal{ext}}/\partial\mathbf{Y}_{n,i}$

11. 3.2.3 Computing $\partial\mathcal{L}_{\rm int}/\partial\mathbf{y}_{n}$

12. 3.2.4 Computing $\partial\mathcal{L}_{\rm int}/\partial\mathbf{Y}_{n,i}$

13. 3.3 First Order Optimality System

14. 3.4 Algorithmic Differentiation

15. 4 Applications

16. 4.1 Lorenz-96 Model

17. 4.2 Parameter Estimation in Magnetic Field Problems

18. 5 Conclusions

1 Introduction

Differential equations are widely used to model the dynamics of physical processes. Even if the form of the equations perfectly captures the physical effects under consideration, the predictive capability of the mathematical model depends on the availability of accurate parameter values and initial conditions. Data assimilation, roughly defined as the solution of inverse problems with models defined by differential equations, fuses information from model outputs and (noisy) physical measurements to produce better estimates of parameter values.

In this work we consider inverse problems where the dynamics is modeled by ordinary differential equations or time dependent partial differential equations, and where a maximum a posteriori (MAP) estimate of the parameter values is computed. Specifically, we consider four-dimensional variational (4D-Var) data assimilation problems where physical measurements are taken at different points in time, and the parameters are updated using all available measurements in the given time window and imposing the system dynamics as strong constraints [21, 37, 13, 30, 39, 26, 34, 38, 36, 15, 46]. The time dimension is discretized using a W-method [43] formulation of a class of exponential integrators that can achieve a high order of accuracy with low stage count and have been shown to be capable of dealing with many stiff problems [44, 45]. To maximize the posterior distribution gradient based optimization techniques are employed. The gradients required in the optimization process are given by the discrete adjoint approach, i.e., are obtained by algorithmic differentiation of the chosen exponential integration method.

Consider the following initial value problem in autonomous form:

[TABLE]

where $\mathbf{y}(t)\in\mathbb{R}^{K}$ is the model state, $\boldsymbol{\theta}\in\mathbb{R}^{P}$ is the vector of model parameters, and the right-hand side function $\mathbf{f}:\mathbb{R}^{K}\times\mathbb{R}^{P}\rightarrow\mathbb{R}^{K}$ is assumed to be smooth and continuously differentiable.

A general inverse problem to estimate the uncertain parameters $\boldsymbol{\theta}$ is formulated as follows:

[TABLE]

where the goal function $\Phi:\mathbb{R}^{P}\rightarrow\mathbb{R}$ is expressed as a sum of an integral involving a nonlinear function of the model state and parameters, $q:\mathbb{R}^{K}\times\mathbb{R}^{P}\rightarrow\mathbb{R}$ and nonlinear functions of the initial and final states and the model parameters, $u,w:\mathbb{R}^{K}\times\mathbb{R}^{P}\rightarrow\mathbb{R}$.

Two different approaches exist to solve problem (2). In the first optimize then discretize approach, the first-order-optimality conditions for the continuous problem (2) are derived first, and then a discretization scheme is applied to numerically solve the resulting equations [6]. This approach affords the freedom to choose different discretization schemes for the forward and adjoint problems. In the first discretize and then optimize approach, the original problem (2) is first discretized using an appropriate time stepping method to handle the constraint (2b) and approximating the integral in (2a) with a quadrature. The first-order-optimality conditions for the discrete optimization problem are then derived [6, 7, 31]. In this paper we study the latter approach, as it yields the exact gradient of the discrete system [34, 35] and can be partially automated with the help of algorithmic differentiation [17, 29]. Other reasons for not preferring the former approach may include the difficulty in obtaining the necessary conditions for problems that are non-trivial and the need to re-derive these for each new problem [8, 6], and the additional overhead involved in the derivation and implementation of the corresponding discrete system.

The paper is structured as follows. In Section 2 the exponential integrator and the discrete parameter estimation problem are introduced. Section 3 derives the first-order-optimality equations and the discrete adjoint of the exponential integrator. Numerical experiments are carried out in Section 4 for two different test problems. Conclusions are drawn in Section 5.

2 Formulation of the Discrete Parameter Estimation Problem

In this section we describe the discrete parameter estimation problem based on a time discretization that uses a W-formulation [43] of a class of exponential methods.

2.1 EPIRK-W Time Discretization Methods

Exponential Propagation Iterative Methods of Runge-Kutta type (EPIRK) [44, 45] use a very general ansatz among exponential time-stepping schemes that allows the construction of methods with low number of stages and high-order. A general $s$-stage EPIRK method reads:

[TABLE]

where $\mathbf{y}_{n}$ is the state at the current time $t_{n}$, $\mathbf{J}_{n}=\partial\mathbf{f}/\partial\mathbf{y}|_{t_{n}}$ is the Jacobian matrix of the right hand side function (1), $\mathbf{Y}_{n,1}\dots\mathbf{Y}_{n,s-1}$ are the intermediate stages, and $\mathbf{y}_{n+1}$ is the next-step solution. Each $\boldsymbol{\psi}_{i,j}(\cdot)$ matrix function is a linear combination of matrix functions $\boldsymbol{\varphi}_{k}(\cdot)$:

[TABLE]

where the scalar analytic function $\varphi_{k}(z)$ is defined as:

[TABLE]

with the matrix function counterparts defined using a series expansion. The functions $\varphi_{k}(z)$ satisfy the following recurrence relation:

[TABLE]

The major computational cost of EPIRK methods is the evaluation of matrix-function-times-vector products of the form \boldsymbol{\psi}_{i,j}\big{(}h\,\gamma\,\mathbf{J}_{n})\cdot\mathbf{v}, which are linear combinations of products \boldsymbol{\varphi}_{k}\big{(}h\,\gamma\,\mathbf{J}_{n})\cdot\mathbf{v}. For large systems Krylov-subspace methods are the preferred choice for the evaluation of these products [20]. It has been shown that for certain kinds of problems, Krylov-based approximation of matrix exponential products converge faster than the corresponding Krylov-based linear system solves [20]. Here we use a Krylov-subspace based algorithm to compute the matrix exponential-like vector products as explained in Section 3.4.

In (3) the $\boldsymbol{\psi}$ matrix functions multiply the right-hand side function evaluated at the current state, $\mathbf{f}(\mathbf{y}_{n})$, and the forward difference vector $\Delta^{(j-1)}\mathbf{r}(\mathbf{y}_{n})$. The forward difference is defined recursively with the help of the remainder function $\mathbf{r}(\mathbf{y})$ as follows:

[TABLE]

Although classical EPIRK methods perform well on stiff problems [45], using them in the first discretize and then optimize approach will require the computation of the Hessian of the right-hand side function. In order to circumvent this, we resort to a W-method formulation of the EPIRK method [27]. W-methods first introduced in [43] for Rosenbrock-type methods admit the use of arbitrary approximations in place of the exact Jacobian, while maintaining the derived order of convergence.

In [27] the authors have recently developed W-type EPIRK methods. An $s$-stage EPIRK-W method reads:

[TABLE]

where the exact Jacobian $\mathbf{J}_{n}$ in (3) has been replaced with an arbitrary approximation $\mathbf{T}_{n}$. The exact implication of this will become clear from the discussion that follows. To simplify notation, we define:

[TABLE]

With this notational change, one step of an EPIRK-W method reads:

[TABLE]

2.2 The Discrete Formulation of the Parameter Estimation Problem

In the first discretize and then optimize approach to solve (2) we partition the time interval into $N$ subintervals with nodes $t_{0},t_{1},\ldots,t_{N}={t_{F}}$. The constraint ODE (2b) is replaced by its EPIRK-W discretization (9) on each subinterval $[t_{n},t_{n+1}]$. The integral term of the goal function (2a) is also discretized by a quadrature rule evaluated on these time nodes. The discrete version of the parameter estimation problem (2) reads:

[TABLE]

The definition of the discrete cost function (10a) includes the initial term $u(\mathbf{y}_{\rm ini},\boldsymbol{\theta})$ in $q_{0}(\mathbf{y}_{\rm ini},\boldsymbol{\theta})$, and the final term $w(\mathbf{y}_{N},\boldsymbol{\theta})$ in $q_{N}(\mathbf{y}_{N},\boldsymbol{\theta})$.

We proceed by describing a special case of the optimization problem (10), that will also be used in the numerical examples of Section 4.

2.3 4D-Var Data Assimilation

The four-dimensional variational data assimilation problem is a special case of (10). Adopting the notation of [42], let $\boldsymbol{\theta_{\textsc{true}}}$ denote the true but unknown parameter vector and $\boldsymbol{\theta}_{\textsc{b}}$ its prior, which is also referred to as background. We assume that $\mathcal{E}_{\textsc{b}}=\boldsymbol{\theta}_{\textsc{b}}-\boldsymbol{\theta_{\textsc{true}}}$ is normally distributed with zero mean and covariance matrix $\mathbf{B}$. Let measurements $\mathbf{z}_{i}$ be given at time points $t_{j_{i}},\,i=1,\ldots,N_{\textsc{obs}}$, which are assumed to be a subset of the nodes of the time interval, for simplicity. Let $\mathcal{H}$ denote an observation operator that maps a state $\mathbf{y}_{j_{i}}$, at time $t_{j_{i}}$, to the space of observations. The measurement error $\mathcal{E}^{\mathrm{obs}}_{i}=\mathcal{H}(\mathbf{y}_{\textsc{true},j_{i}})-\mathbf{z}_{i}$, where $\mathbf{y}_{\textsc{true},j_{i}}$ represents the state associated to $\boldsymbol{\theta_{\textsc{true}}}$ at time $t_{j_{i}}$.We assume a normal distribution for $\mathcal{E}^{\mathrm{obs}}_{i}$, with mean zero and covariance matrix $\mathbf{R}_{i}$. Assuming no model errors, the 4D-Var cost function and gradient are defined as follows:

[TABLE]

The value $\boldsymbol{\theta_{\textsc{map}}}$ minimizing (10) with the goal function given by (11), represents the maximum-likelihood estimate of $\boldsymbol{\theta_{\textsc{true}}}$. In a data assimilation context, $\boldsymbol{\theta_{\textsc{map}}}$ is also referred to as analysis.

Remark 1*.*

In parameter estimation problem (2), the goal function and ODE constraint depend on a vector of unknown parameters $\boldsymbol{\theta}\in\mathbb{R}^{P}$ that need to be inferred from measurements. In the general case, computation of gradient of the cost function, as shown in ((C)), involves complicated derivatives with respect to $\boldsymbol{\theta}$. This can be circumvented by converting this into an initial state estimation problem. To this end, the ODE system can be extended with the addition of parameters $\boldsymbol{\tau}(t)=\boldsymbol{\theta}$ as formal variables to the state vector as follows:

[TABLE]

3 First Order Optimality Conditions and the Discrete EPIRK-W Adjoint

This section is devoted to the first-order-optimality system associated with the discrete inverse problem (10), which involves discrete adjoints for the EPIRK-W method. Various authors have addressed discrete adjoints in optimal control, constrained by ordinary differential equations. In [18] discrete adjoints for Runge-Kutta methods were derived together with control specific order conditions. A similar approach was taken in [22] for Runge-Kutta W-methods. In both works, the discrete adjoint method could be reformulated again as a Runge-Kutta method and this was exploited in the subsequent error analysis. A general discussion on the use of discrete time integration adjoints in the solution of inverse problems is presented in [36]. Discrete adjoints of Rosenbrock methods are formulated in [11, 12]. Theoretical properties of general discrete Runge-Kutta adjoints are proved in [34], and theoretical properties of discrete adjoints of linear multistep methods in [35]. Discrete adjoints of variable-step integrators are studied in [2]. Efficient implementations of Runge-Kutta adjoints are provided by software developed by the authors [26, 1, 41, 38, 46, 4].

This work is concerned with deriving the discrete adjoints of EPIRK-W methods. As we will see, the discrete adjoint equations do not readily show the structure of an exponential method. Discrete adjoints for another closely related class of exponential methods can be found in [33].

3.1 The Discrete Lagrangian

The first-order-optimality system is obtained by seeking a stationary point of the Lagrangian of the discrete optimal control problem (10). We introduce the compact notation:

[TABLE]

Let $\boldsymbol{\Lambda}_{n,i}$, $\boldsymbol{\lambda}_{n}$ denote the Lagrange multipliers associated to (10b) and (10c), respectively. The Lagrangian of the discrete optimal control problem (10) reads:

[TABLE]

where:

[TABLE]

In the above equations, the forward difference operator $\Delta^{(j-1)}\mathbf{r}(\mathbf{y}_{k},\boldsymbol{\theta})$ can be written in a closed form expression as

[TABLE]

where $\mathbf{r}(\mathbf{Y}_{k,j-1-\ell},\boldsymbol{\theta})$ is the remainder term of the first-order Taylor expansion of $\mathbf{f}$ around $\mathbf{y}_{k}$, evaluated at $\mathbf{Y}_{k,j-1-\ell}$ with the exact Jacobian $\mathbf{J}_{k}$ replaced by the approximation $\mathbf{T}_{k}$:

[TABLE]

It should be pointed out that in the above expressions $\mathbf{Y}_{k,0}=\mathbf{y}_{k}$ was introduced for brevity in later expressions.

Substituting (16) into (15) and (15b), respectively, we obtain:

[TABLE]

3.2 Derivation of the Discrete Adjoint Equation

The adjoint equations (20) are obtained by setting $\partial\mathcal{L}/\partial\mathbf{y}_{n}=0$ and $\partial\mathcal{L}/\partial\mathbf{Y}_{n,i}=0$ for all $n>0$. For a better readability, we compute various derivatives separately. We thereby adopt the convention that all terms for which subscripts turn out to be negative are automatically dropped or set to zero.

3.2.1 Computing $\partial\mathcal{L}_{\textnormal{ext}}/\partial\mathbf{y}_{n}$

For $0<n<N$, only the terms $k=n-1$ and $k=n$ in the sum $\sum_{k=0}^{N-1}$ give non-zero contributions, as follows:

[TABLE]

where $\mathbf{I}$ and $\mathbf{J}_{n}$ are the $\mathbb{R}^{K}\times\mathbb{R}^{K}$ identity matrix and the Jacobian matrix evaluated at $\mathbf{y}_{n}$, respectively. Also the definition $\mathbf{Y}_{k,0}=\mathbf{y}_{k}$ was used in the derivation above.

For $n=N$, only the term with index $k=N-1$ remains and we obtain:

[TABLE]

For $n=0$ we compute:

[TABLE]

3.2.2 Computing $\partial\mathcal{L}_{\textnormal{ext}}/\partial\mathbf{Y}_{n,i}$

The same nonzero terms in the sum remain when we differentiate with respect to $\mathbf{Y}_{n,i}$ as when we differentiate with respect to $\mathbf{y}_{n}$. For $0\leq n<N$ we obtain:

[TABLE]

For $n=N$ the derivative is zero.

3.2.3 Computing $\partial\mathcal{L}_{\rm int}/\partial\mathbf{y}_{n}$

As parts of the derivation were already given in Sec. 3.2.1 we omit intermediate steps. For $0\leq n<N$ we obtain:

[TABLE]

whereas there is no contribution for $n=N$.

3.2.4 Computing $\partial\mathcal{L}_{\rm int}/\partial\mathbf{Y}_{n,i}$

Again, omitting several intermediate results due to the similarity to Sec. 3.2.2 we obtain for $0\leq n<N$:

[TABLE]

Finally, we observe that at the stationary point $(\hat{\mathbf{y}}^{*},\hat{\mathbf{Y}}^{*},\hat{\boldsymbol{\lambda}}^{*},\hat{\boldsymbol{\Lambda}}^{*})$ of the Lagrangian, both $\mathcal{L}_{\textnormal{ext}}$ and $\mathcal{L}_{\textrm{int}}$ vanish and hence

[TABLE]

Differentiating this expression with respect to the parameter $\boldsymbol{\theta}$ yields

[TABLE]

from which we obtain the expression for the gradient, given in (C) below, by inserting the respective expressions for $\mathcal{L}_{\textnormal{ext}}$ and $\mathcal{L}_{\rm int}$.

3.3 First Order Optimality System

The first-order-optimality conditions of (10) are obtained from the results of Sec. 3.2.1–3.2.4 by enforcing the derivatives of $\mathcal{L}$ to zero:

1. (A)

The original EPIRK-W method, i.e., (10b), (10c) and (10d), is obtained by equating $\partial\mathcal{L}/\partial\boldsymbol{\lambda}_{n}=0$ and $\partial\mathcal{L}/\partial\boldsymbol{\Lambda}_{n,i}=0$, respectively. 2. (B)

The discrete adjoint system is obtained as follows.

[TABLE]

It should be emphasized that (20) is solved from ${t_{F}}$ to ${t_{0}}$ in reverse time direction. 3. (C)

The optimality condition reads:

[TABLE]

3.4 Algorithmic Differentiation

In the previous section, we derived the first-order-optimality conditions for discrete optimization problem in (10). To derive the adjoint we use algorithmic differentiation (AD) [29], a technique that generates code to compute the sensitivities of an output with respect to an input for a given program, to compute gradient of the Lagrangian. The site www.autodiff.org lists numerous tools that perform AD for programs written in two dozen languages. Excellent resources for the theory and implementation techniques of AD include [28, 17].

Two different modes of AD exist for first-order sensitivities - forward and reverse. The forward mode AD produces the tangent linear model of the input program, whereas the reverse mode produces the adjoint model. Mathematically, the forward mode computes Jacobian-vector products for the given program, and adjoint mode Jacobian-transposed-vector products, where the vector is some seed direction. Depending on the desire to compute a combination of the rows or columns of the Jacobian one uses the forward or adjoint mode, respectively. For a detailed discussion of cost in terms of operations and memory for each of the two modes we refer to [17]. It suffices to say that for the computation of the gradient of the Lagrangian, where the number of inputs far exceeds that of the outputs, the reverse mode is more economical in terms of operation count.

A drawback of the reverse mode is the necessity to store the entire forward program trajectory that will be required while “reversing” the program to compute the sensitivities. If the given program has a loop structure, then all the variables that affect the output and change from one iteration to the next have to be stored in a stack like data structure to be used during reversal of the loop in the program. As a result, storage requirements will grow with long running loop based programs. One technique to alleviate the situation is to use checkpointing [17], where periodically the stack is written to a file before being flushed for reuse. The checkpoint files are read in the reverse order while “reversing” the computations to calculate sensitivities. Technique are available to balance the costs of checkpointing and recomputation [16].

In this work we implemented a single-step EPIRK-W integrator in Fortran 90 making use of the following major operations:

Linear combinations of vectors (axpy) are carried out using the built-in array operation syntax.

• -

Matrix vector products (gemv) are implemented as an explicit subroutine call instead of using the Fortran intrinsic (MATMUL).

• -

The right-hand side $\mathbf{f}(\mathbf{y})$ is an external callback function provided by the caller.

• -

Computation of $\boldsymbol{\psi}_{i,j}$ function products is rather involved and is explained in detail below.

The $\boldsymbol{\psi}_{i,j}$ products can be written as a linear combination of individual $\boldsymbol{\varphi}_{k}$ products, where we approximate each $\boldsymbol{\varphi}_{k}$ product in the Krylov-subspace [44] as \boldsymbol{\varphi}_{k}\big{(}h\,\gamma\,\mathbf{T}_{n}\big{)}\mathbf{b}\approx\left\lVert\mathbf{b}\right\rVert\mathbf{V}\boldsymbol{\varphi}_{k}(h\,\gamma\,\mathbf{H})\mathbf{e}_{1}, where $\mathbf{V}$ is the matrix containing orthogonal basis vectors of the Krylov-subspace $\mathcal{K}_{m}=\;\textnormal{span}\,\{\mathbf{b},\mathbf{T}_{n}\mathbf{b},\mathbf{T}_{n}^{2}\mathbf{b},\ldots\}$ and $\mathbf{H}$ is the upper-Hessenberg matrix resulting from an Arnoldi process. Following [40, Theorem1], $\boldsymbol{\varphi}_{k}(h\,\gamma\,\mathbf{H})\mathbf{e}_{1}$ products are computed by constructing an augmented matrix and exponentiating it using Expokit [40].

The adjoint of a single-step EPIRK-W integrator is obtained with the help of Tapenade [19]. We differentiate through both axpy and gemv operations, but use black-box routines (http://www-sop.inria.fr/tropics/tapenade/faq.html#Libs1) for both $\mathbf{f}(\mathbf{y})$ and $\boldsymbol{\psi}_{i,j}(h\,\gamma\,\mathbf{T}_{n})\mathbf{v}$. The adjoint of $\mathbf{f}(\mathbf{y})$ is an external callback provided by the caller. The adjoint implementation of $\boldsymbol{\psi}_{i,j}(h\,\gamma\,\mathbf{T}_{n})\,\mathbf{v}$ is based on the following result.

Lemma 3.1** (Adjoint of matrix-function-times-vector operation).**

Let $\mathbf{y}_{n}$ be the independent and $\mathbf{y}_{n+1}$ the dependent of a single-step EPIRK-W method (8). From equations (4) and (5) we have:

[TABLE]

where $\mathbf{T}_{n}$ is an arbitrary approximate Jacobian that does not depend on $\mathbf{y}_{n}$. The adjoint of the operation:

[TABLE]

where $\mathbf{v}$ is a vector independent of $\mathbf{y}_{n}$, reads:

[TABLE]

where the adjoint variables are $\bar{\mathbf{v}}=\bar{\mathbf{y}}_{n+1}^{T}\frac{\partial\mathbf{y}_{n+1}}{\partial\mathbf{v}}$ and $\boldsymbol{\bar{\alpha}}=\bar{\mathbf{y}}_{n+1}^{T}\frac{\partial\mathbf{y}_{n+1}}{\partial\boldsymbol{\alpha}}$.

Proof.

Given an arbitrary matrix $\mathbf{X}$ and a vector $\mathbf{v}$ that is a function of the independent $\mathbf{y}_{n}$, the adjoint of the product $\mathbf{z}\leftarrow\mathbf{X}\mathbf{v}$ is the operation $\bar{\mathbf{v}}\leftarrow{\mathbf{X}^{T}}\bar{\mathbf{z}}$. Consequently, since $\mathbf{T}_{n}$ does not depend on $\mathbf{y}_{n}$, we have:

[TABLE]

It remains to show that

[TABLE]

We are given

[TABLE]

Taking the transpose then yields

[TABLE]

Since

[TABLE]

we have

[TABLE]

Taking the transpose then yields

[TABLE]

From equations (23) and (25), it is not too difficult to infer

[TABLE]

thus proving the result. ∎

Remark 2*.*

If the matrix argument is the exact Jacobian $\mathbf{J}_{n}$, or any matrix that depends on $\mathbf{y}_{n}$, then the adjoint of the operation

[TABLE]

involves the additional computation of the Hessian tensor $\partial\mathbf{J}_{n}/\partial\mathbf{y}_{n}$. This is the main motivation for using a W-method formulation of EPIRK in our work: the assumption that matrices $\mathbf{T}_{n}$ are independent of $\mathbf{y}_{n}$ allows the adjoint operations to be implemented in accordance with Lemma 3.1, while the W-property allows maintaining the overall accuracy of the forward simulation.

In the following section we consider two test problems where a 4D-Var problem is solved for parameter/initial state estimation. The discrete adjoint EPIRK-W integrator is used to compute the gradients needed in optimization.

4 Applications

In this section, two different numerical examples will be given. As a first test case, we consider the Lorenz-96 model, which is a discrete nonlinear example, where the parameter is the initial condition. As a second test case, an inductor will be considered. The original equation is obtained by applying the finite element method to the magnetoquasistatic partial differential equations. Then, data assimilation is carried out for parameters characterizing the nonlinear magnetic material in the iron core. For both examples a Python version of the limited memory variant of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm [10] is used for optimization. We use the third-order EPIRK-W method derived in [27, Figure 2] in our experiments. The integrator and its adjoint is made callable from Python using f2py and checkpoints are stored in-memory.

4.1 Lorenz-96 Model

The model considered here was originally proposed in [24] and used in a data assimilation context in [42], where a detailed description of the setup can be found. Its main points are summarized here to enhance readability. The deterministic model is given by the set of equations:

[TABLE]

where $K=40$ and $y=(y^{1},\ldots,y^{K})^{T}$. Let $t_{0}=0$, ${t_{F}}=0.3$ and the time interval be $\Delta t=0.015$ time units (time units).

We ran convergence test for the Lorenz-96 model of both the forward and the adjoint integrator with the Jacobian initialized to a random matrix at the beginning, and using a fixed random seed for the adjoint integrator. Reference solutions for both forward and adjoint problem were computed using the MATLODE [4] ERK forward and adjoint integrators, respectively, with both absolute and relative tolerance set to $1.010-12$.

As Figure 1 shows, both forward and adjoint EPIRK-W integrators achieve full third order convergence in the Euclidean norm, as desired. It should be noted, however, that either the exact Jacobian or a very good approximation of it may be needed for better stability of the integrators. In [27], the convergence behavior of the forward EPIRK-W integrator was treated in detail. A complete treatment of the convergence behavior of the adjoint EPIRK-W integrator will be considered in a future work. We continue to describe the example in which we apply these integrators.

In this example, $\boldsymbol{\theta}$ parametrizes the initial condition as $\mathbf{y}_{\rm ini}=\boldsymbol{\theta}$, whereas $\mathbf{f}$ is a function of $\mathbf{y}$, solely. Following [42], the reference $\boldsymbol{\theta}_{\textsc{true}}$ is obtained by integrating $y^{j}(-10\Delta t)=1+0.1\mathrm{mod}(j,5)$ in time until $t=0$ (time units). The observation operator $\mathcal{H}$ is defined as:

[TABLE]

We model the background and observation standard deviation as:

[TABLE]

where $\mathbf{y}_{j_{i}}^{\textsc{true}}$ refers to the solution at step $j_{i}=i100$, $i=1,\ldots,N_{\textsc{obs}}=10$, with initial condition $\boldsymbol{\theta_{\textsc{true}}}$. These standard deviations are used to define the background and observation covariance matrices as

[TABLE]

with $\alpha=0.1$, $L=4$, $(\mathbf{D})_{i,j}=\min(|i-j|,K-|i-j|)$. We emphasize that $\mathbf{R}$ is modeled to be independent of time in this case. Also, in this section, $\mathbf{I}$ refers to the $K\times K$ identity matrix. The background prior $\boldsymbol{\theta}_{\textsc{b}}$ is generated as a pseudo-random realization according to the normal distribution of $\mathcal{E}_{\textsc{b}}$ with mean value $\boldsymbol{\theta}_{\textsc{true}}$ and covariance $\mathbf{B}$. Correspondingly, the measurement $\mathbf{z}_{i}$ is generated as a pseudo-random realization around the observation reference according to the normal distribution of $\mathcal{E}_{\textsc{obs}}$, with mean value $\mathcal{H}(\mathbf{y}_{j_{i}}(\boldsymbol{\theta}_{\textsc{true}}))$ and covariance $\mathbf{R}$.

We perform data assimilation by solving (10) using the goal function and gradient given in (11) and (12), respectively. The discrete adjoints are computed using algorithmic differentiation as outlined in Section 3.4 and supplied to the optimization routine, with the approximation $\mathbf{T}_{n}=\mathbf{J}_{n}$ used in the forward and adjoint integrators. As can be seen from Figure 2 the algorithm converges quickly. More precisely, the cost function does not change significantly after two iterations, whereas a gradient of below $1.010-4$ is obtained within 8 iterations, for a step size of $\Delta t=0.0003$ time units. Figure 3 depicts the difference of the initial condition at the current iteration and the true initial condition in the Euclidean norm.

4.2 Parameter Estimation in Magnetic Field Problems

Before describing the data assimilation problem for the inductor, we give some background of the underlying partial differential equations. The magnetoquasistatic problem in a bounded computational domain $\Omega$ reads as

[TABLE]

and is widely used in the simulation of electrical machines, magnets and transformers. In (30), $\vec{E}$ denotes the electric field, $\kappa$ the conductivity, $\vec{J}$ the imposed source current density and $\vec{n}$ the outer unit normal. Also, the magnetic flux density $\vec{B}$ and the magnet field strength $\vec{H}$ are related through the nonlinear material relation $\vec{H}=\nu(|\vec{B}|)\vec{B}$. Here, we restrict ourselves to parametric nonlinearities of the type

[TABLE]

see [32], where $\nu(0)$ is obtained by taking the limit $|\vec{B}|\rightarrow 0$ in the above equation. The parameters $\boldsymbol{\theta}$ of the nonlinearity are often deduced from measurements. Equation (30) is discretized in space by applying the finite element method. To this end, as a first step, from (30) we derive the second order problem

[TABLE]

with the magnetic vector potential $\vec{A}$ such that $\vec{B}=\operatorname{\nabla\times}\vec{A}$ and the initial condition $\vec{A}(0)=\vec{A}_{0}$ in $\Omega$, see, e.g., [9]. We further simplify the problem by considering a two-dimensional domain. By inserting $\vec{A}=(0,0,u)^{T}$ in (32), where the third coordinate represents the direction perpendicular to the two-dimensional domain, we obtain

[TABLE]

Difficulties arise in the solution of problem (33), as $\kappa$ vanishes in non-conducting sub-regions of $\Omega$. A remedy consists in using $\kappa_{\text{reg}}=\max(\kappa,\kappa_{\text{air}})$, with artificial conductivity $\kappa_{\text{air}}>0$, in (32) instead. The error of this approximation as a function of $\kappa_{\text{air}}$ is well understood, see, e.g., [5]. We assume that $\kappa$ is replaced with $\kappa_{\text{reg}}$ from now on and omit the subscript for simplicity.

To state the weak formulation of (33) we introduce the usual Sobolev space

[TABLE]

Then, we seek $u:[t_{0},{t_{F}}]\rightarrow V$ such that

[TABLE]

The finite element space $V_{h}\subset V$ consists of piecewise linear polynomials on a quasi-uniform triangulation of $\Omega$. This results in an ordinary differential equation

[TABLE]

where $\mathbf{M},\mathbf{K},\mathbf{y}_{\text{FE}},\mathbf{f}_{\text{FE}}$ refer to the mass and stiffness matrix and to the vector of degrees of freedom and the source vector, respectively. Problem (36) can be recast as

[TABLE]

where $\mathbf{f}(\mathbf{y}_{\text{FE}},\boldsymbol{\theta})=\mathbf{M}^{-1}(\mathbf{f}_{\text{FE}}-\mathbf{K}(\mathbf{y}_{\text{FE}},\boldsymbol{\theta})\mathbf{y}_{\text{FE}})$. From (37) it is clearly visible that the parameter dependence is bound to $\mathbf{f}$ solely.

The example considered here is an inductor, introduced in [25], which is depicted in Figure 4. The observed (scalar) quantity is the flux linkage defined as

[TABLE]

where $N_{\text{turns}}$, $\Omega_{\text{coil}}$ and $|\Omega_{\text{coil}}|$ refer to the number of turns in the coil, the coil domain (red in Figure 4) and its area, respectively.

The parameter estimation problem is set up as follows. We assume that the true parameters of the nonlinearity (31) are given by:

[TABLE]

Here, the background and measurement covariance are modeled as:

[TABLE]

i.e., both the parameters and the measurements are assumed to be uncorrelated. We consider a sinusoidal current excitation $I(t)=I_{0}\sin(4\pi t)$ on the time interval $[0,0.5]$ s. Measurements are taken at time steps $j_{i}=10i$, for $i=1,\ldots,5=N_{\textsc{obs}}$. The measurement standard deviation is modeled to be $0.1$ percent of the effective measurement value:

[TABLE]

Then, the background prior $\boldsymbol{\theta}_{\textsc{b}}$ is generated as a pseudo-random realization according to a normal distribution with mean $\boldsymbol{\theta_{\textsc{true}}}$ and covariance $\mathbf{B}$. Accordingly, measurements $z_{i}$ are generated as pseudo-random realizations, according to a normal distribution with covariance $\mathbf{R}$ around the true trajectory.

With the goal function and gradient as in equations (11) and (12), we apply the EPRIK-W method to solve (37). This problem is very stiff and it warrants the use of the exact Jacobian, i.e. $\mathbf{T}_{n}=\mathbf{J}_{n}$, to enhance stability. We allow for a 20 percent variation of the coefficients during optimization. For the simulation a step size of $\Delta t=0.01$ s is employed. Using symmetry only the right half of the configuration is used in the computations. A current of $I_{0}=150$ A per turn, with $N_{\text{turns}}=66$ turns in total, is imposed to the coil. The associated flux distribution in the core is depicted in Figure 4 on the right. Meshing and finite element analysis are carried out with Gmsh [14] and FEniCS [3, 23], respectively. To reduce the stiffness of the problem we set $\kappa_{\text{air}}=$1.0106$\ \mathrm{S}\mathrm{m}^{-1}$, which has a negligible effect on the flux linkage.

Figure 5 shows the cost function and the gradient over the iteration steps of the optimization. Again, we observe a fast convergence of the cost function. In Figure 6 we plot the difference of parameter $\theta_{4}$ and its true value, during optimization. The error in $\theta_{4}$ decays quickly, while the other parameters remain unchanged. This can be attributed to the fact that the derivative of the objective function with respect to $\theta_{4}$ is orders of magnitudes larger than the derivatives with respect to $\theta_{1},\theta_{2}$ and $\theta_{3}$. The observation is further supported by the fact that $\theta_{4}$ models the saturated range of the nonlinearity [32], which contributes significantly to the shape of the waveform depicted in Figure 7. As this shape is modified during data assimilation to fit the measurements, $\theta_{4}$ plays a key role here. Figure 7 depicts the flux linkage over the time interval for the prior and the parameters after assimilation together with the measurements. The trajectory after data assimilation is in better accordance with the measurement data. Finally, Table 1 summarizes the measurement values as well as the numerical predictions before and after data assimilation.

5 Conclusions

This work derives the discrete adjoint formulae for a general class of exponential integrators, EPIRK-W. The choice of methods with the W-property allows to use arbitrary approximations of the Jacobian as arguments of the matrix functions while maintaining the overall accuracy of the forward integration. The use of matrices that do not depend on the model state avoids the complex calculation of Hessians (derivatives of the Jacobian with respect to state variables) in the discrete adjoint formulae. The simplified discrete adjoint can computed via algorithmic differentiation and then supplied to optimization routines for solving the inverse problem at hand. We were able to empirically verify that the convergence order of the discrete adjoint of the EPIRK-W method matches the order of convergence of the forward EPIRK-W integrator. The methodology was applied to estimate the initial condition of the Lorenz-96 model from synthetic measurement data. The methodology was also applied to estimate magnetic material parameters for computational magnetics problems, where a system of stiff ordinary differential equations is obtained by a method of lines approach. Exponential integrators are a promising alternative for this class of problems and this topic merits further investigation.

Acknowledgements

The work of M. Narayanamurthi and A. Sandu was supported in part by awards NSF DMS–1419003, NSF CCF–1613905, AFOSR DDDAS 15RT1037, and by the Computational Science Laboratory at Virginia Tech. Part of the work was performed during U. Römer’s visit at Virginia Tech.