Preconditioned pseudo-time continuation for parameterized inverse problems
Joseph Hart, Alen Alexanderian, and Bart van Bloemen Waanders
TL;DR
This paper introduces an efficient pseudo-time continuation method with a novel adaptive quasi-Newton preconditioner to solve parameterized inverse problems constrained by PDEs, significantly reducing computational costs for uncertainty quantification.
Contribution
We develop a new adaptive quasi-Newton preconditioner for pseudo-time continuation, improving the efficiency of solving parameterized inverse PDE problems.
Findings
Preconditioner accelerates convergence of inverse problem solutions.
Method effectively handles nonlinear inverse problems.
Framework enables efficient uncertainty quantification.
Abstract
We consider parameterized variational inverse problems that are constrained by partial differential equations (PDEs). We seek to efficiently compute the solution of the inverse problem when auxiliary model parameters, which appear in the governing PDE, are varied. Computing the solution of the inverse problem for different auxiliary parameter values is crucial for uncertainty quantification. This, however, is computationally challenging since it requires solving many optimization problems for different realizations of the auxiliary parameters. We leverage pseudo-time continuation and solve an initial value problem to evolve the optimal solution along an auxiliary parameter path. This article introduces the use of an adaptive quasi-Newton Hessian preconditioner to accelerate the computation. Our proposed preconditioner exploits properties of the pseudo-time continuation process to…
| Method | State | Gradient | Total Linear Solves | ||
|---|---|---|---|---|---|
| 7 | 7 | 6 | 567 | 1181 | |
| 7 | 7 | 6 | 1064 | 2175 | |
| 7 | 7 | 6 | 1511 | 3069 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
TopicsMatrix Theory and Algorithms · Numerical methods in inverse problems · Model Reduction and Neural Networks
\newsiamremark
remarkRemark \newsiamremarkhypothesisHypothesis
\newsiamthmclaimClaim \newsiamremarkfactFact
\headersPseudo-time ContinuationJ. Hart, A. Alexanderian, and B. van Bloemen Waanders
Preconditioned Pseudo-time Continuation for Parameterized Inverse Problems††thanks: Submitted to the editors DATE.
Joseph Hart Department of Scientific Machine Learning, Sandia National Laboratories, Albuquerque, NM (, https://www.sandia.gov/ccr/staff/joseph-lee-hart/).
Alen Alexanderian Department of Mathematics, North Carolina State University, Raleigh, NC (, https://aalexan3.math.ncsu.edu).
Bart van Bloemen Waanders Department of Scientific Machine Learning, Sandia National Laboratories, Albuquerque, NM (, https://www.sandia.gov/ccr/staff/bart-g-van-bloemen-waanders/)..)
Abstract
We consider parameterized variational inverse problems that are constrained by partial differential equations (PDEs). We seek to efficiently compute the solution of the inverse problem when auxiliary model parameters, which appear in the governing PDE, are varied. Computing the solution of the inverse problem for different auxiliary parameter values is crucial for uncertainty quantification. This, however, is computationally challenging since it requires solving many optimization problems for different realizations of the auxiliary parameters. We leverage pseudo-time continuation and solve an initial value problem to evolve the optimal solution along an auxiliary parameter path. This article introduces the use of an adaptive quasi-Newton Hessian preconditioner to accelerate the computation. Our proposed preconditioner exploits properties of the pseudo-time continuation process to achieve reliable and efficient computation. We elaborate our proposed framework and elucidate its properties for two nonlinear inverse problems.
keywords:
Optimization; post-optimal sensitivity analysis; inverse problems
{MSCcodes}
65B99, 65K10, 65L05, 65M32
1 Introduction
Inverse problems arise in a variety of applications where a quantity of interest cannot be measured directly. Rather, indirect observations are made and the quantity of interest is related to the observations via a mathematical model. Often, the model is a partial differential equation (PDE), the inversion parameter of interest is a model parameter, and the observable quantity is a function of the solution of the PDE. Various approaches exist to estimate the model parameters. All the approaches require solving the PDE for different model parameters and comparing the model output with the observed data by evaluating a data-misfit functional.
The PDE governing an inverse problem typically has a number of auxiliary parameters. These are model parameters that are needed to fully specify the governing PDE, but are not the inversion parameters. We focus on the challenge of understanding the impact of auxiliary parameters on the inverse problem. This is crucial for quantifying the uncertainty in the estimated inversion parameter. A first-principles approach is to jointly estimate the inversion parameter of interest and the auxiliary parameters [6]. However, this is frequently intractable due to the increased dimensionality induced by including the auxiliary parameters in the inversion. The auxiliary parameter dimension may be decreased via screening and sensitivity analysis techniques [24, 26]. Although such approaches may provide auxiliary parameter dimension reduction, the presence of influential auxiliary parameters is frequently inescapable and thus must be accounted for in the inverse problem.
In a Bayesian paradigm, auxiliary parameter uncertainty can be accounted for by the Bayesian approximation error (BAE) approach [22, 21]. This is achieved by pre-marginalizing the auxiliary parameter uncertainty, which results in a modified noise model. In this article, we focus on variational inverse problems. Such deterministic inverse problems require solving an optimization problem whose objective is the data-misfit—the difference between the model prediction and the observed data. Also, in ill-posed inverse problems, which are our primary focus, a regularization term is added to the data-misfit functional to define the objective function.
A pragmatic approach is to fix the auxiliary parameter to a nominal value and solve the optimization problem to determine the parameter of interest that results in the model agreeing with the observational data. However, to incorporate auxiliary parameter uncertainty, it is necessary to solve many such optimization problems for different auxiliary parameters. This is computationally prohibitive due to the many optimization solves required.
To avoid solving optimization problems repeatedly, hyper-differential sensitivity analysis (HDSA) was developed to compute the sensitivity of the nominal optimization problem’s solution with respect to perturbations of the auxiliary parameters [13, 14, 29, 30]. This provides a linear approximation of the auxiliary parameter to optimal solution mapping. The key benefit of HDSA is its computational scalability. Namely, HDSA provides an efficient and accurate framework for estimating the sensitivity of the estimated inversion parameter with respect to the auxiliary parameter in large-scale inverse problems with high-dimensional inversion and auxiliary parameters. However, this linear approximation may be insufficient for large perturbations, e.g., when there is significant auxiliary parameter uncertainty.
To overcome the limitations of linear approximation, continuation and path-following algorithms can be used to solve perturbed optimization problems. Early work focused on continuation methods to solve systems of nonlinear equations. See the books [2, 3], and the references therein, for foundational contributions in continuation methods, such as the predictor-corrector algorithm. These ideas also translate to optimization problems since the first-order optimality condition defines a system of nonlinear equations. Furthermore, optimization problems have additional structure, such as symmetry of the Hessian matrix, that can be exploited to develop efficient algorithms. See the book [12] for an exposition of continuation methods for parametric optimization problems and [17] for a treatment of dual degenerate problems. See also the recent work [18, 28], which demonstrates the utility of path-following—continuation along a path through the parameter space—to compute optimal parameters for statistical models. These works focus on regression problems where the primary source of complexity is in the regularization term arising from a Bayesian prior. On the other hand, the present work focuses on inverse problems where the greatest challenge arises from the PDE-constrained data-misfit term.
We build on the recent work [1] that introduced the idea of pseudo-time continuation, which is synonymous with path-following, in the context of parameterized PDE-constrained inverse problems. The present article differs from the breadth of continuation literature as a result of the computational challenges and mathematical structures that arise from PDE-constrained optimization. These include the necessity of iterative solvers for Hessian systems and the computation of Hessian-vector products via the adjoint method. Our contributions in this article are: (i) the design of adaptive Hessian preconditioners that improve computational efficiency, and (ii) demonstration on challenging PDE-constrained inverse problems that elucidate the key features of the proposed approach. We focus on efficiently computing the minimizer for one specific auxiliary parameter perturbation by tracking minimizers along a path from the nominal parameter value to the perturbed parameter value. There are natural extensions of this work to computing minimizers for many different parameter perturbations. This is discussed in our conclusions.
The article is organized as follows. Section 2 provides the mathematical framing of the problem and recalls relevant literature. The computational kernel of pseudo-time continuation is iterative linear solves to invert the Hessian. Our approach to preconditioning these solves is introduced in Section 3. Numerical results are presented in Section 4 to demonstrate the benefits of our proposed approach and give insight into the effect of the algorithm parameters (step sizes, tolerances, etc.). In Section 5, we provide concluding remarks and discuss potential extensions of the proposed approach.
2 Preliminaries
Consider parameterized variational inverse problems of the form
[TABLE]
Here, is the vector of inversion parameters, is a parameter-to-observable map, {\mathchoice{\mbox{\displaystyle{y}}}{\mbox{\textstyle{y}}}{\mbox{\scriptstyle{y}}}{\mbox{\scriptscriptstyle{y}}}}\in\mathbb{R}^{d} is measurement data, is a vector of auxiliary parameters, and {\mathchoice{\mbox{\displaystyle{\mathrm{R}}}}{\mbox{\textstyle{\mathrm{R}}}}{\mbox{\scriptstyle{\mathrm{R}}}}{\mbox{\scriptscriptstyle{\mathrm{R}}}}}\in\mathbb{R}^{n\times n} is a symmetric positive definite matrix that is used for regularization. We consider inverse problems governed by PDEs in which the auxiliary parameters correspond to model parameters that are not being estimated. This is why the parameter-to-observable map is written as a function of both and . The governing PDE must be solved each time {\mathchoice{\mbox{\displaystyle{\mathrm{F}}}}{\mbox{\textstyle{\mathrm{F}}}}{\mbox{\scriptstyle{\mathrm{F}}}}{\mbox{\scriptscriptstyle{\mathrm{F}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}};\boldsymbol{\theta}) is evaluated. Consequently, solving (1) is challenging due to the computational cost of solving the PDE repeatedly as {\mathchoice{\mbox{\displaystyle{\mathrm{F}}}}{\mbox{\textstyle{\mathrm{F}}}}{\mbox{\scriptstyle{\mathrm{F}}}}{\mbox{\scriptscriptstyle{\mathrm{F}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}};\boldsymbol{\theta}) is evaluated for different inputs.
We assume throughout that J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) is twice continuously differentiable with respect to ({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) and use \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J and \nabla_{{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}}J to denote the gradient and Hessian of with respect to ; similar notation is used to denote the derivatives. We consider problems for which matrix-vector products with the Hessian \nabla_{{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) may be computed and use Newton methods to solve (1). The Hessian information facilitates rapid convergence of the optimization iterates. Given a fixed and current optimization iterate {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k}, the Newton search direction is obtained by solving the system of linear equations
[TABLE]
where {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta})=\nabla_{{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}) is the Hessian of the objective function. The Newton step is {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k+1}={\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k}+{\mathchoice{\mbox{\displaystyle{s}}}{\mbox{\textstyle{s}}}{\mbox{\scriptstyle{s}}}{\mbox{\scriptscriptstyle{s}}}}_{k}. A globalization strategy, such as line search or trust region, is necessary to ensure convergence of the optimization algorithm. However, when sufficiently close to the minimizer, the Newton steps result in rapid (quadratic) convergence of the optimization iterations toward the minimizer.
When solving the inverse problem, traditionally is fixed at a vector of nominal auxiliary parameters. When the auxiliary parameters are subject to considerable uncertainty, it is important to quantify the uncertainty in the solution of the inverse problem. However, doing this by repeatedly solving (1) for different values of is computationally expensive. In what follows, we let {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}) denote a minimizer of . Note that for non-convex optimization problems, uniqueness of {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}) is ensured by restricting to a neighborhood of the local minimizer computed when ; see [1] for details. In this section, we recall the pseudo-time continuation approach and the predictor-corrector algorithm to compute {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}) with {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\overline{\theta}}) used as a starting point. Section 3 details this article’s contribution to enhance the algorithms through judicious preconditioning to compute {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}) more efficiently.
2.1 Pseudo-time continuation
Consider the Jacobian of {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*} with respect to , which we denote by \nabla_{\boldsymbol{\theta}}{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}. Applying the Implicit Function Theorem to the first-order optimality condition for (1) yields
[TABLE]
Here {\mathchoice{\mbox{\displaystyle{\mathrm{B}}}}{\mbox{\textstyle{\mathrm{B}}}}{\mbox{\scriptstyle{\mathrm{B}}}}{\mbox{\scriptscriptstyle{\mathrm{B}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}),\boldsymbol{\theta})=\nabla_{{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}),\boldsymbol{\theta}) is the mixed second derivative of the objective function and {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}),\boldsymbol{\theta}) is the Hessian previously introduced in the discussion on Newton’s method; see (2). We can consider a linear approximation model
[TABLE]
This approximation is only valid in a neighborhood of and may therefore be inaccurate if is large. To overcome this limitation, we consider a pseudo-time continuation approach. Let be a parameter vector for which we seek to compute {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\tilde{\theta}}) and parameterize a path from to via , . With an abuse of notation, we consider {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(t)\coloneqq{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}(t)). It follows that
[TABLE]
Let
[TABLE]
where . We may compute {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\tilde{\theta}}) by solving the initial value problem
[TABLE]
To ensure that {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}(t)),\boldsymbol{\theta}(t)) is invertible, and that {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}(t)) is tracking a strict minimizer, we make the following assumption:
Assumption 2.1**.**
{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}(t)),\boldsymbol{\theta}(t))* is positive definite for all . *
In general, tracking a stationary point {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}(t)) that satisfies \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}(t)),\boldsymbol{\theta}(t))={\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}} does not guarantee positive definiteness of the Hessian. However, this assumption is reasonable in practice for many inverse problems arising from physical systems, particularly since regularization is typically used in ill-posed problems. Lastly, we note that a common approach to addressing near-zero eigenvalues in the Hessian is to project onto a likelihood informed subspace [7]. If indefiniteness is expected at the nominal minimizer {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(0), we can consider optimization in only the likelihood informed subspace, which will ensure satisfaction of the assumption in almost all practical scenarios.
Observe that the linear approximation (3) corresponds to discretizing (5) with forward Euler and integrating with a single time step of size 1. However, by viewing the evolution of the minimizer via an ODE system (5), we may consider higher-fidelity time stepping algorithms to ensure more accurate solutions. To demonstrate the benefits of pseudo-time continuation relative to an optimization algorithm, we consider a simple example below.
Example 2.1**.**
Consider the parameterized optimization problem
[TABLE]
where is an uncertain auxiliary parameter. This illustrative example captures the common inverse problem scenario that the objective function has small curvature in some parameter directions [11]. We assume that (6) has been solved when to produce the minimizer . To compute the minimizer when , we take as the initial iterate and use Newton’s method to solve (6) with . For comparison, we use pseudo-time continuation taking a forward Euler discretization of (5) with time steps. Figure 1 displays the iterations produced by Newton’s method and the time steps produced by pseudo-time continuation.
We note that the pseudo-time continuation approach converges faster than Newton’s method. We conjecture that this improved convergence occurs for the following interrelated reasons:
In Newton’s method, we consider a sequence of iterations that minimize successive local quadratic approximations of the objective function. On the other hand, in pseudo-time continuation we consider (5), which describes the evolution of the minimizer as a function of . Subsequently, we approximate the numerical solution of that initial value problem with a time-stepping method. Hence, the pseudo-time continuation approach permutes the actions of approximation and minimization. 2. 2.
By taking forward Euler time steps, we perturb by at each step whereas Newton’s method takes a size perturbation () of initially. 3. 3.
Since our pseudo-time continuation algorithms changes gradually, we stay within the Newton basin such that a time step of the form (7) tracks the minimizer well. Whereas initializing at is outside the Newton basin and hence does not enjoy the rapid quadratic convergence that characterizes Newton’s method when it is “close enough” to the minimizer.
These observations motivate the study of pseudo-time continuation algorithms. Also, as demonstrated shortly, posing the perturbed optimization problem in terms of the initial value problem (5) provides other opportunities to accelerate computation.
2.2 Predictor-corrector algorithm
To obtain the minimizer corresponding to a perturbed auxiliary parameter vector , we can, in principle, solve the initial value problem (5) using a time-stepping method. However, the ODE system (5) has additional structure that is not exploited by time-stepping methods. Specifically, for any , {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(t) is a local minimizer of J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}(t)) and hence \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}(t),\boldsymbol{\theta}(t))={\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}}. If we integrate (5) via a traditional time-stepping method, the accumulation of error results in \|\nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}(t),\boldsymbol{\theta}(t))\| increasing as a function of time. Rather, we can use a predictor-corrector procedure [2, 3, 12] where a traditional time-stepping method predicts the solution and a Newton optimization iteration corrects the predicted solution to ensure that \|\nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}(t),\boldsymbol{\theta}(t))\| is less than a specified tolerance. Note that despite using the terminology, and having many conceptual similarities, predictor-corrector algorithms in the present context are different from the classical predictor-corrector ODE solvers. In particular, predictor-corrector ODE solvers use two different time stepping schemes, whereas the algorithms we consider use a time stepping scheme for prediction followed by a Newton step for correction. This corrector is unique to path-following of optimization solutions.
To facilitate the discussion, we first consider the setup of a time-stepping approach for (5). We focus on explicit methods and discretize (5) using time steps. Let denote the step size and let , for . We let {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k} denote the estimated value of {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(t_{k}) and for .
Consider the numerical solution of (5) approximated using forward Euler,
[TABLE]
The Euler-Newton predictor-corrector scheme consists of the steps:
[TABLE]
In the above procedure, the forward Euler step (8a) serves as a predictor. Subsequently, one step of Newton’s method is used as a corrector. In practice, the predictor step results in an approximation to {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k+1} that falls close to the basin of attraction for the {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(t_{k+1}). The corrector step then uses this as an initial guess for a Newton step. As noted later in the article, additional corrector steps may be needed to ensure optimality within a prescribed tolerance. However, the optimality typically holds with one or two corrector steps thanks to the rapid convergence of Newton’s method.
The Euler-Newton predictor-corrector scheme (8) is commonly used [3, 12]. To provide further intuition, we recall classical observations of similarity between the forward Euler predictor step (8a) and an optimization iterate using Newton’s method. Specifically, a first-order approximation of \Delta t{\mathchoice{\mbox{\displaystyle{\mathrm{B}}}}{\mbox{\textstyle{\mathrm{B}}}}{\mbox{\scriptstyle{\mathrm{B}}}}{\mbox{\scriptscriptstyle{\mathrm{B}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})\Delta\boldsymbol{\theta} yields
[TABLE]
where the latter approximation follows since \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}(t_{k}),\boldsymbol{\theta}_{k})={\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}} and {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k}\approx{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(t_{k}). Consequently, the forward Euler step (8a) is approximately given by
[TABLE]
To minimize J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}_{k+1}) with Newton’s method using an initial iterate {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k}, the iteration takes the form
[TABLE]
Hence, the predictor step may be viewed as a type of Newton step that is informed by the variability.
There is a lot to say about the procedure outlined in (8). In the first place, one may consider higher-order time-stepping methods for the predictor. In Section 2.3, we present and discuss the advantages of a second-order accurate method in the predictor step. Another important consideration is the cost of the Hessian solves in the predictor and corrector steps. In the case of optimization problems arising from discretized inverse problems constrained by PDEs, these Hessian solves are expensive and are tackled using Krylov iterative methods. Specifically, we rely on the Conjugate Gradient (CG) method. The cost of CG solves for these Hessian systems is dominated by expensive matrix-vector products with the Hessian—each Hessian-vector product requires two PDE solves. In this context, preconditioning is an essential tool to reduce the number of CG iterations. We outline our proposed preconditioners in Section 3. Our approach leverages properties of the pseudo-time continuation and predictor-corrector algorithms to evolve the preconditioner with time.
2.3 The predictor step
The role of the predictor step is to provide a sufficiently accurate initial guess for the corrector step to ensure rapid convergence. Much of the work on optimal solution path-following predictor-corrector methods rely on forward Euler time-stepping [3]. However, if the auxiliary parameter perturbation is large, we may need a small time step when using forward Euler.
Forward Euler is first-order accurate, that is, the truncation error of time discretization is . A second-order method has a truncation error that is . This has the potential to reduce computational cost by taking larger time steps. There are multiple second-order methods that we can consider. In this article, we focus on a modified Euler method defined by the steps:
[TABLE]
This modified Euler method corresponds to taking a step of size using forward Euler (9a) to predict an intermediate {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k+\frac{1}{2}}, followed by a step of size from time to using the slope \boldsymbol{f}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k+\frac{1}{2}},\boldsymbol{\theta}_{k+\frac{1}{2}}) evaluated at the midpoint predicted by the forward Euler step (9b). There are many potential second-order numerical integration schemes, see [10] for an overview. This modified Euler method is preferred to other second-order schemes because it increments time by in both (9a) and (9b).
Consider forward Euler time-stepping with time steps and modified Euler time-stepping with time steps (assuming is even). The time nodes in the modified Euler time discretization correspond exactly to the time nodes in the forward Euler time discretization. Hence, both methods require evaluations of and thus have approximately the same computational cost. However, the truncation error of forward Euler is whereas the truncation error of modified Euler is . If and the constants hidden in the are equal, then the modified Euler method has a smaller truncation error with the same computational cost. This simple analysis is instructive to see when and why higher-order methods may be superior. Our numerical results presented in Section 4 demonstrate the improvements that can be attained using modified Euler when the auxiliary parameter perturbation is large enough to mandate small forward Euler time steps. Although we omit discussion of higher-order (greater than second-order) methods now, we highlight their potential benefits in our concluding remarks; see Section 5.
3 Preconditioners for the Hessian solves
The pseudo-time continuation algorithm requires many evaluations of \boldsymbol{f}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}), which was defined in (4). The key challenge when evaluating is computing the action of {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})^{-1} on vectors. These Hessian solves are performed using the Conjugate Gradient (CG) method. Each CG iteration requires a Hessian apply, which involves solving the incremental state and incremental adjoint equations; see Appendix A for details. Hence, to achieve computational efficiency, it is crucial to precondition the Hessian solves so that we reduce the number of CG iterations.
In this section, we present novel quasi-Newton methods to build effective Hessian preconditioners within the pseudo-time continuation algorithm. We note at the onset that this is made possible because our approach starts at a local minimizer and tracks its changes as is perturbed. Furthermore, as we previously assumed (and argued is reasonable in practice), the Hessian matrices remain positive definite throughout the pseudo-time steps. In a general optimization setting, effective preconditioning of the Hessian system is more difficult. We note that the idea to design Hessian preconditioners to accelerate CG solves within a path-following algorithm was also considered in [28]. In that work, the optimization problem arose from Bayesian hierarchical statistical modeling where the complexity stemmed from the prior model. Specialized preconditioners were developed in [28] by exploiting structure associated with the statistical modeling problem. In our context, the complexity comes from the underlying PDE constraint. This motivates a different approach to design preconditioners that adapt using information from the PDE (and associated adjoint) solves.
In preparation for our contributions, we first recall the Preconditioned Conjugate Gradient (PCG) method in Algorithm 1, and refer the reader to [16] for a more detailed discussion. Note that we present PCG using Euclidean inner products. In general, if the linear system arises from the discretization of an linear operator defined on function spaces, then the PCG solver should use a weighted inner product that arises from discretizing the function space inner product [19]. However, in the context of the Hessian systems in Newton optimization algorithms, the inner product weighting is already embedded in the discretized Hessian so that the algorithms may use Euclidean inner products.
Let denote the preconditioner, which must be symmetric positive definite and should approximate the inverse Hessian {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})^{-1}. Ensuring that the condition number of \mathbf{E}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) is small improves the convergence of the PCG iteration. An alternative perspective on preconditioning is that we desire to have the eigenvalues of \mathbf{E}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})\mathbf{E}^{\frac{1}{2}} clustered. In exact arithmetic, if \mathbf{E}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})\mathbf{E}^{\frac{1}{2}} has distinct eigenvalues then PCG will converge in at most iterations. From the perspective of clustering eigenvalues, an attractive preconditioner is , where is the regularization operator in (1) [11]. To see this, let {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}_{\text{M}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) be the Hessian of the data-misfit term \frac{1}{2}\|{\mathchoice{\mbox{\displaystyle{F}}}{\mbox{\textstyle{F}}}{\mbox{\scriptstyle{F}}}{\mbox{\scriptscriptstyle{F}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})-{\mathchoice{\mbox{\displaystyle{y}}}{\mbox{\textstyle{y}}}{\mbox{\scriptstyle{y}}}{\mbox{\scriptscriptstyle{y}}}}\|^{2} and observe that
[TABLE]
where denotes the identity matrix. In many applications, \mathbf{R}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}_{\text{M}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})\mathbf{R}^{\frac{1}{2}} is approximately low-rank as a result of data sparsity and smoothing properties of the parameter-to-observable map [11]. Hence the regularization-preconditioned data-misfit Hessian \mathbf{R}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})\mathbf{R}^{\frac{1}{2}}+\mathbf{I} will have many of its eigenvalues clustered near . While is a commonly used preconditioner in practice, we demonstrate that superior convergence can be achieved by adapting the preconditioner within the pseudo-time continuation algorithm. To this end, we introduce two novel quasi-Newton updates that exploit the pseudo-time continuation structure.
Quasi-Newton methods are commonly used for numerical optimization in problems where gradients are available, but Hessians are not. The idea of a quasi-Newton method is to begin the optimization algorithm with an initial estimate of the Hessian (frequently a poor one) and update it at each optimization iteration using gradient evaluations. We build on this idea and propose an approach wherein we initialize an estimate of the Hessian and update it after each time step within the predictor corrector iterations (8).
In the subsections that follow, we present two quasi-Newton preconditioner updates. The first one is based on introducing a secant equation that captures the variations. We refer to this as the parametric quasi-Newton update. The second one is an adaptation of the block BFGS method [9], which leverages the PCG iteration history to enhance the preconditioner in future time steps. We refer to this second update as the block quasi-Newton update. Central to both updates is that only existing gradients and Hessian-vector products are used to improve the approximation. Hence, additional PDE solves are not required. However, storage of additional vectors is needed that would otherwise be removed from memory.
3.1 Parametric Quasi-Newton Update
When solving (5), we require a sequence of {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})^{-1} applies, where ({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}) changes gradually. Traditional quasi-Newton methods, which rely on gradient differences, can be used to precondition these linear solves. However, this ignores the dependence of {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k} on . In the traditional setting, the Hessian only depends on . Here, we derive a quasi-Newton update that accounts for the parameter variations.
To initialize the preconditioned pseudo-time continuation algorithm, we require that a preconditioner for {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{0},\boldsymbol{\theta}_{0}) is provided. In the context of inverse problems, there are two natural options to define an initial preconditioner: (i) use the inverse of the matrix in (1), or (ii) use a low-rank approximation of {\mathchoice{\mbox{\displaystyle{R}}}{\mbox{\textstyle{R}}}{\mbox{\scriptstyle{R}}}{\mbox{\scriptscriptstyle{R}}}}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}_{\text{M}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{0},\boldsymbol{\theta}_{0}){\mathchoice{\mbox{\displaystyle{R}}}{\mbox{\textstyle{R}}}{\mbox{\scriptstyle{R}}}{\mbox{\scriptscriptstyle{R}}}}^{\frac{1}{2}} to approximate {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{0},\boldsymbol{\theta}_{0})^{-1}=({\mathchoice{\mbox{\displaystyle{R}}}{\mbox{\textstyle{R}}}{\mbox{\scriptstyle{R}}}{\mbox{\scriptscriptstyle{R}}}}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}_{\text{M}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{0},\boldsymbol{\theta}_{0}){\mathchoice{\mbox{\displaystyle{R}}}{\mbox{\textstyle{R}}}{\mbox{\scriptstyle{R}}}{\mbox{\scriptscriptstyle{R}}}}^{\frac{1}{2}}+\mathbf{I})^{-1}. We provide a detailed derivation of the latter in Appendix B. For the analysis that follows, we let {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k} denote our preconditioner that approximates {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})^{-1}, for , where {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{0} is provided via one of the options discussed above.
We assume that a preconditioner {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k-1} is available and seek to update it to produce a preconditioner {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k} for the next time step. To derive a quasi-Newton update, we introduce a secant equation in the context of pseudo-time continuation. Let
[TABLE]
be a local quadratic approximation of the objective function J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}_{k}) centered at {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k}. Here, \widehat{{\mathchoice{\mbox{\displaystyle{H}}}{\mbox{\textstyle{H}}}{\mbox{\scriptstyle{H}}}{\mbox{\scriptscriptstyle{H}}}}}_{k} denotes an approximate Hessian. We seek to determine \widehat{{\mathchoice{\mbox{\displaystyle{H}}}{\mbox{\textstyle{H}}}{\mbox{\scriptstyle{H}}}{\mbox{\scriptscriptstyle{H}}}}}_{k} and define {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k}=\widehat{{\mathchoice{\mbox{\displaystyle{H}}}{\mbox{\textstyle{H}}}{\mbox{\scriptstyle{H}}}{\mbox{\scriptscriptstyle{H}}}}}_{k}^{-1}.
In traditional optimization, is fixed and only {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k} changes with each iteration, whereas in our context both {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k} and change with each time step. We enforce the secant conditions
[TABLE]
The former is satisfied by construction. The latter gives
[TABLE]
The gradient \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}) is computed as a by-product of the pseudo-time continuation algorithm; see Appendix A. However, computing \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k-1},\boldsymbol{\theta}_{k}) would require an additional forward and adjoint solve. To avoid this, we consider a first-order Taylor expansion in ,
[TABLE]
This implies the approximate secant condition
[TABLE]
where {\mathchoice{\mbox{\displaystyle{z}}}{\mbox{\textstyle{z}}}{\mbox{\scriptstyle{z}}}{\mbox{\scriptscriptstyle{z}}}}_{k}={\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k}-{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k-1} and
[TABLE]
Note that {\mathchoice{\mbox{\displaystyle{y}}}{\mbox{\textstyle{y}}}{\mbox{\scriptstyle{y}}}{\mbox{\scriptscriptstyle{y}}}}_{k} is defined in terms of vectors that are already computed.
Mimicking the BFGS update [23, Chapter 6], rather than requiring (10), we multiply (10) by {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k}=\widehat{{\mathchoice{\mbox{\displaystyle{H}}}{\mbox{\textstyle{H}}}{\mbox{\scriptstyle{H}}}{\mbox{\scriptscriptstyle{H}}}}}_{k}^{-1} and enforce the secant condition
[TABLE]
We seek to determine a symmetric update matrix \Delta{\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k-1} such that {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k}={\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k-1}+\Delta{\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k-1} and {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k} satisfies (11). To determine an optimal \Delta{\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k-1}, we consider
[TABLE]
where {\mathchoice{\mbox{\displaystyle{x}}}{\mbox{\textstyle{x}}}{\mbox{\scriptstyle{x}}}{\mbox{\scriptscriptstyle{x}}}}_{k}={\mathchoice{\mbox{\displaystyle{z}}}{\mbox{\textstyle{z}}}{\mbox{\scriptstyle{z}}}{\mbox{\scriptscriptstyle{z}}}}_{k}-{\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k-1}{\mathchoice{\mbox{\displaystyle{y}}}{\mbox{\textstyle{y}}}{\mbox{\scriptstyle{y}}}{\mbox{\scriptscriptstyle{y}}}}_{k} and is suitably defined matrix norm. As detailed in [23, Chapter 6], (12) admits an analytic solution, which leads to the update
[TABLE]
where \rho_{k}\coloneq 1/({\mathchoice{\mbox{\displaystyle{y}}}{\mbox{\textstyle{y}}}{\mbox{\scriptstyle{y}}}{\mbox{\scriptscriptstyle{y}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{z}}}{\mbox{\textstyle{z}}}{\mbox{\scriptstyle{z}}}{\mbox{\scriptscriptstyle{z}}}}_{k}). Here, we require that the curvature condition holds. If it is not satisfied, the update (13) is omitted. However, empirical observations suggest that is satisfied as long as the pseudo-time continuation procedure is faithfully tracking a minimizer about which \nabla_{{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}}J is positive definite. Note that the form of (13) is identical to the classical BFGS update. The difference is in the definition of {\mathchoice{\mbox{\displaystyle{y}}}{\mbox{\textstyle{y}}}{\mbox{\scriptstyle{y}}}{\mbox{\scriptscriptstyle{y}}}}_{k}.
3.2 Block Quasi-Newton Update
The parametric quasi-Newton update introduced in the previous subsection incorporates the time evolution of directly into the preconditioner. This alone is beneficial for preconditioning the CG iterations. However, we can further enhance the quasi-Newton update by using the Hessian-vector products computed within PCG iterations. Specifically, after executing the PCG solve at time step , we have a set of matrix-vector products {\mathchoice{\mbox{\displaystyle{w}}}{\mbox{\textstyle{w}}}{\mbox{\scriptstyle{w}}}{\mbox{\scriptscriptstyle{w}}}}_{i}={\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}){\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}, (see Line 11 of Algorithm 1), where we assume PCG has converged in iterations. The data \{({\mathchoice{\mbox{\displaystyle{w}}}{\mbox{\textstyle{w}}}{\mbox{\scriptstyle{w}}}{\mbox{\scriptscriptstyle{w}}}}_{i},{\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i})\}_{i=1}^{\ell} can be used to enhance {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k} before applying the parametric quasi-Newton update (13) to obtain {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k+1}. We discuss the details of this process next.
Traditionally, quasi-Newton methods consider rank 1 or rank 2 updates that satisfy a single secant equation. To incorporate the data
[TABLE]
we seek to enforce secant equations as stated in the following block form,
[TABLE]
This amounts to requiring that the quasi-Newton Hessian approximation \widehat{{\mathchoice{\mbox{\displaystyle{H}}}{\mbox{\textstyle{H}}}{\mbox{\scriptstyle{H}}}{\mbox{\scriptscriptstyle{H}}}}}_{k} agrees with the full Hessian on the Krylov space explored in the PCG iterations. As shown in [27], a symmetric matrix satisfying (14) exists if and only if {\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k} is symmetric. In the present setting, this condition is satisfied since {\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}){\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}. We leverage this observation and define a block BFGS update following the approach introduced in [9]. Specifically, given a preconditioner {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k} and PCG data ({\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k},{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}), the symmetric matrix {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k}^{\text{up}} that is closest to {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k} and satisfies {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k}^{\text{up}}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k} is
[TABLE]
where {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}\in\mathbb{R}^{\ell\times\ell}.
Note that the PCG search directions \{{\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}\}_{i=1}^{\ell} satisfy the orthogonality property {\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}^{\top}{\mathchoice{\mbox{\displaystyle{w}}}{\mbox{\textstyle{w}}}{\mbox{\scriptstyle{w}}}{\mbox{\scriptscriptstyle{w}}}}_{j}=0 for . Therefore, in exact arithmetic, {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k} is a diagonal matrix. Since the block BFGS update (15) depends on {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k}^{-1}, care must be taken if a diagonal entry {\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}^{\top}{\mathchoice{\mbox{\displaystyle{w}}}{\mbox{\textstyle{w}}}{\mbox{\scriptstyle{w}}}{\mbox{\scriptscriptstyle{w}}}}_{i} is small. Following the filtering idea presented in [9], we discard vectors for which {\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}^{\top}{\mathchoice{\mbox{\displaystyle{w}}}{\mbox{\textstyle{w}}}{\mbox{\scriptstyle{w}}}{\mbox{\scriptscriptstyle{w}}}}_{i} is small. This is done by specifying a tolerance and discarding any vector for which {\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}^{\top}{\mathchoice{\mbox{\displaystyle{w}}}{\mbox{\textstyle{w}}}{\mbox{\scriptstyle{w}}}{\mbox{\scriptscriptstyle{w}}}}_{i}<\tau\|{\mathchoice{\mbox{\displaystyle{p}}}{\mbox{\textstyle{p}}}{\mbox{\scriptstyle{p}}}{\mbox{\scriptscriptstyle{p}}}}_{i}\|^{2}. If , we retain all vectors and as we discard all vectors and hence do not perform the update. We used as a default value in our numerical results.
However, in practice the aforementioned orthogonality of PCG search directions is not preserved numerically. This loss of orthogonality, which is a well-known phenomenon, causes the matrix {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k} to be non-diagonal. Nonetheless, {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k} is symmetric up to numerical precision. Therefore, in our implementation, we compute the eigenvalue decomposition {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k}={\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}{\mathchoice{\mbox{\displaystyle{\Lambda}}}{\mbox{\textstyle{\Lambda}}}{\mbox{\scriptstyle{\Lambda}}}{\mbox{\scriptscriptstyle{\Lambda}}}}_{k}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}^{\top}, and modify (15) by replacing {\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}, {\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}, and {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k} with {\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\prime}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}, {\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}^{\prime}={\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}, and {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k}^{\prime}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\prime\top}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}^{\prime}, respectively. Since the matrix {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k} is typically low-dimensional, the eigenvalue decomposition may be computed with standard dense linear algebra routines. Rotating {\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k} and {\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k} with the eigenvectors {\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k} is equivalent to enforcing the secant condition {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{k}^{\text{up}}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}^{\prime}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\prime} and leads to the numerically diagonal matrix {\mathchoice{\mbox{\displaystyle{D}}}{\mbox{\textstyle{D}}}{\mbox{\scriptstyle{D}}}{\mbox{\scriptscriptstyle{D}}}}_{k}^{\prime}={\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\prime\top}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}^{\prime}={\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{P}}}{\mbox{\textstyle{P}}}{\mbox{\scriptstyle{P}}}{\mbox{\scriptscriptstyle{P}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{W}}}{\mbox{\textstyle{W}}}{\mbox{\scriptstyle{W}}}{\mbox{\scriptscriptstyle{W}}}}_{k}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}={\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}{\mathchoice{\mbox{\displaystyle{\Lambda}}}{\mbox{\textstyle{\Lambda}}}{\mbox{\scriptstyle{\Lambda}}}{\mbox{\scriptscriptstyle{\Lambda}}}}_{k}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}^{\top}{\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}_{k}={\mathchoice{\mbox{\displaystyle{\Lambda}}}{\mbox{\textstyle{\Lambda}}}{\mbox{\scriptstyle{\Lambda}}}{\mbox{\scriptscriptstyle{\Lambda}}}}_{k}.
3.3 Algorithm Summary
Algorithm 2 summarizes preconditioned forward Euler pseudo-time continuation and highlights the preconditioner adaptation. In the loop over time steps , we first compute a Forward Euler time step to produce {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k+1}^{\text{pred}}\approx{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}_{k+1}). Next, we use a block quasi-Newton update to improve the time step preconditioner using {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}) matrix-vector products, and subsequently produce a time step preconditioner via a parametric quasi-Newton update. The corrector step executes a Newton optimization iteration at time step , and subsequently updates the time step preconditioner via a block quasi-Newton update using {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k+1}^{\text{pred}},\boldsymbol{\theta}_{k+1}) matrix-vector products. The tolerance satisfaction step may not be required as the predictor/corrector scheme typically predicts {\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}^{*}(\boldsymbol{\theta}_{k+1}) effectively. However, if the time step size is large, the predictor/corrector scheme may be inaccurate. The subsequent tolerance satisfaction step executes Newton optimization iterations until convergence. The PCG tolerance may be reduced in the tolerance satisfaction step to ensure Newton convergence properties. However, in our numerical experiments, tightening this tolerance yielded no benefits.
The initial preconditioner rank specifies {\mathchoice{\mbox{\displaystyle{\mathrm{E}}}}{\mbox{\textstyle{\mathrm{E}}}}{\mbox{\scriptstyle{\mathrm{E}}}}{\mbox{\scriptscriptstyle{\mathrm{E}}}}}_{0} as: (i) the inverse regularization Hessian if , or (ii) an approximate inverse based on a rank approximation of {\mathchoice{\mbox{\displaystyle{R}}}{\mbox{\textstyle{R}}}{\mbox{\scriptstyle{R}}}{\mbox{\scriptscriptstyle{R}}}}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}_{\text{M}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{0},\boldsymbol{\theta}_{0}){\mathchoice{\mbox{\displaystyle{R}}}{\mbox{\textstyle{R}}}{\mbox{\scriptstyle{R}}}{\mbox{\scriptscriptstyle{R}}}}^{\frac{1}{2}} if . Thanks to the block quasi-Newton update, our proposed approach has a self-correcting property. If the preconditioner is a poor approximation of the inverse Hessian, then more PCG iterations will be required. These iterations provide data to enhance the preconditioner so that it better approximates the inverse Hessian and hence enhances the future preconditioners. The maximum update rank , along with the tolerance , controls the rank of the block Quasi-Newton update. In practice, is specified based on storage limitations.
We can also consider a version of Algorithm 2 that uses modified Euler in place of forward Euler in the predictor step. This involves two Hessian inversions and correspondingly two preconditioner updates within the predictor step. However, this additional cost may be offset by the higher-order predictor method requiring fewer time steps, PCG iterations, and/or tolerance satisfaction iterations.
As a final remark, note that we check the norm of the gradient before executing the corrector step. If the optimality tolerance is satisfied, we forego the corrector step. This can lead to considerable computational savings. We revisit this issue in our numerical experiments in Section 4.2.
4 Numerical results
In this section, we present numerical results for two nonlinear inverse problems. Section 4.1 presents a study on the effect of algorithm parameters in our preconditioned pseudo-time continuation algorithm. These results focus on providing intuition about the algorithm. Section 4.2 demonstrates the use of preconditioned pseudo-time continuation for a challenging ice sheet bedrock inversion problem. The results in this section highlight the utility of our proposed approach for nonlinear and computationally intensive inverse problems.
4.1 Diffusion coefficient estimation
We consider the Poisson equation
[TABLE]
defined on the domain whose boundary consists of , , , and . Here, denotes the outward-facing normal vector to the boundary. In this problem is a prescribed (parameterized) source term and and are predefined functions. We use and . Consider the inverse problem of estimating the spatially distributed log-diffusion coefficient .
In this problem, the source term is parameterized by an auxiliary parameter vector as follows:
[TABLE]
To solve the inverse problem, we discretize (16) on a rectangular mesh. Hence, the discretized state variable and inversion parameter have dimension . In the present study, we use the ground-truth parameter depicted in the left panel of Figure 2 and use a ground-truth auxiliary parameter of , which will be specified below. Subsequently, we generate synthetic data by solving (16) with and extracting solution values on a grid of measurement points depicted in the right panel of Figure 2. To avoid the inverse crime, the data is generated by solving the forward problem on a rectangular mesh and contaminated with Gaussian noise.
We pose an inverse problem in the form of (1) where corresponds to the composition of the (discretized) solution operator of (16) and the observation operator maps the state solution to sparse observations. Furthermore, corresponds to the square of the discretized elliptic operator , with constants and , where and correspond to the Laplacian and identity operators, respectively. This choice is common in practice as it corresponds to using the squared inverse elliptic operator to define the prior covariance in the Bayesian formulation of the inverse problem; see [5].
When solving the inverse problem, we use a nominal source function parameterization rather than the ground-truth auxiliary parameter that was used to generate the data. Subsequently, we use pseudo-time continuation to compute the solution of the inverse problem when , where is an auxiliary parameter vector that approximates .
Varying the number of time steps
In our first numerical experiment, we compare the use of forward Euler and modified Euler time-stepping algorithms for the predictor. We fix , the canonical basis vector whose first entry is and all others are [math], and consider three perturbations of varying magnitudes. Specifically, we consider , where is the vector whose entries are all , and . Note that . For each scenario, we compare the performance of five potential approaches. Specifically, we establish a performance baseline by solving the perturbed optimization problem (1) with and the initial condition given by the optimal solution when . We refer to this as “re-optimization”. In addition, we execute the pseudo-time continuation procedure in four different configurations: using either forward Euler or modified Euler for the predictor steps, and using either the inverse regularization Hessian as the preconditioner for all solves or using our proposed adaptive preconditioner. In the four pseudo-time continuation configurations, we fix the PCG tolerance to and vary the numbers of time-steps; specifically . For our proposed adaptive preconditioner, we use the inverse regularization Hessian as the initial preconditioner and set the maximum rank of each block Quasi-Newton update to be . The results are shown in Figure 3, where we measure the computational cost via the total number of PDE solves required to achieve a gradient norm tolerance of .
We make several observations from Figure 3.
- •
A modified Euler predictor is generally superior to a forward Euler predictor; however, the difference in performance is less pronounced when the perturbation magnitude is small and the number of time steps is large.
- •
Using a fixed preconditioner results in considerable variability in computational cost as the number of time steps varies. In contrast, using the adaptive preconditioner results in small changes in computational cost as the number of time steps varies.
- •
Taking too few time steps, particularly with a forward Euler predictor step, may result in increased computational cost as a result of requiring multiple iterations to satisfy the gradient tolerance.
- •
The computational cost of the pseudo-time continuation algorithms grows as the magnitude of the perturbation increases.
- •
The computational cost of re-optimization is always larger than pseudo-time continuation with a modified Euler predictor and an adaptive preconditioner.
- •
The computational cost of re-optimization is always less than pseudo-time continuation with a forward Euler predictor and a fixed preconditioner. This is because re-optimization takes less costly steps (fewer Hessian-vector products per iteration) when it is initially far from the minimizer.
These observations highlight the benefit of using a modified Euler predictor with the adaptive preconditioner. However, the results above correspond to a fixed PCG tolerance and preconditioner update rank. In what follows, we vary these algorithm parameters while fixing the number of time steps () and the perturbation magnitude ().
Varying the PCG tolerance
In our second numerical experiment, we consider the effect of changing the PCG tolerance. We set the initial preconditioner as the inverse regularization Hessian and the maximum rank of the block Quasi-Newton update to be . Figure 4 displays the computational cost as a function of the PCG tolerance. As before, the cost is measured in terms of the total number of PDE solves required to achieve a gradient norm tolerance of .
We observe that the optimal PCG tolerance in this example is around . In this case, modified Euler with the adaptive preconditioner is slightly superior to the other three pseudo-time continuation approaches, which all have comparable computational costs. We observed that the tolerance could be effective, but the time stepping procedure may not converge if the loose tolerance leads to ineffective predictor steps that fail to track the Newton basin of attraction. Because of this lack of robustness, we omitted from the plot. The success of our pseudo-time continuation with such a loose PCG tolerance is explained through a combination of several factors. First, it is unsurprising that a loose tolerance is beneficial in the predictor step as the goal is to account for variations in and get sufficiently close to the minimizer that the corrector step will be successful. However, surprisingly, a loose tolerance on the corrector step is superior in this example. Although additional tolerance satisfaction steps were required, the total cost was lowest when using a loose tolerance.
It is noteworthy that using a fixed preconditioner results in a significant increase in computational cost as the the tolerance is tightened. This is because the number of Hessian-vector products required in the PCG solve increases. The adaptive preconditioner helps to mitigate this increased cost. However, in our later numerical example, we will observe that the computational cost using the adaptive preconditioner can vary nontrivially with changes in the PCG tolerance.
Varying the preconditioner rank parameters
and
In our final numerical experiment, we focus on the performance of pseudo-time continuation with the modified Euler predictor, as and are varied. In particular, we initialize the preconditioner using a rank approximation of the regularization-preconditioned data-misfit Hessian and set the maximum block Quasi-Newton update to be of rank . We consider all pairs . For each , we execute the pseudo-time continuation procedure with a modified Euler predictor, time steps, and a PCG tolerance of . Figure 4 displays the computational cost (left panel) and the number of vectors stored (right panel) for various rank parameters.
We make several observations from Figure 5.
- •
There is a general trend of decreasing cost with increasing the block Quasi-Newton update rank , although diminishing returns are observed after .
- •
There is a general trend of increasing storage with increasing the block Quasi-Newton update rank ; however, the rate of storage growth slows for larger .
- •
Setting (i.e., not using the block Quasi-Newton update) results in a considerable increase in computational cost.
- •
Investing in the initial preconditioner (i.e., taking a larger ) is advantageous when , but has little benefit for larger .
Scalability and mesh-independence
The proposed computational framework exploits specific problem structures that are independent of the discretization. Namely, the spectral properties of the Hessian operator depend primarily on the regularity of the parameter-to-observable map as well as smoothing properties of the regularization operator. The Newton–Krylov methods exploit these structures to achieve a nearly mesh-independent behavior, as observed in various studies; see, e.g., [11] for discussions and references.
To assess the mesh-independence properties of our proposed framework, we fixed all algorithm parameters and varied the mesh resolution. We observed negligible variation in the computational cost as the mesh resolution varied, thus indicating that the algorithm’s convergence behavior depends on the inherent spectral properties of the second-order derivative operators rather than the dimension of the discretized problem. The results are omitted for brevity.
4.2 Ice sheet bedrock inversion
We consider estimating the bedrock topography on the bottom of the Greenland Ice Sheet. For simplicity, we focus on a km region of the ice sheet (denoted as ) and use the shallow ice approximation [20, 15] with the isothermal assumption [4]. In what follows, we summarize the inverse problem and use it as a testbed to analyze computational aspects of our pseudo-time continuation approach. We refer the reader to [25], where the inverse problem and its sensitivities were studied, for a more detailed discussion.
The forward model is
[TABLE]
where is the surface height of the ice sheet (relative to sea level), is a velocity field, is the log basal friction on the interface of the ice sheet and the bedrock beneath it, and is a forcing term modeling the accumulation/ablation of ice. Scalar parameters include: the density of ice , acceleration due to gravity , and a flow rate factor . The forward model (17) is transient. We integrate over the time interval from [math] to years.
The inverse problem is to estimate the bedrock topography using observations of the surface ice velocity, which is defined as a function of the PDE solution via
[TABLE]
The log basal friction and forcing are spatially heterogeneous uncertain auxiliary parameters that we model as
[TABLE]
where
[TABLE]
are parameterized perturbations; are linear finite element basis functions defined on a rectangular mesh of the domain . This implies that the auxiliary parameter dimension is . The nominal auxiliary parameters correspond to taking \overline{\boldsymbol{\theta}}={\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}}\in\mathbb{R}^{1922} so that \beta_{{\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}}}=\overline{\beta} and h_{{\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}}}=\overline{h}.
After discretization, the inverse problem takes the form of (1) where , the parameter-to-observable map, is defined via composing the discretized solution operator for (17) and velocity model (18). In practice, observational data comes from satellite measurements whose spatial resolution is greater than the computational mesh resolution. Hence, observations are assumed to be available at all spatial locations. The regularization is defined as the negative log of the prior covariance in a Bayesian formulation of the inverse problem, so that corresponds to the square of the discretized elliptic operator , with constants and .
Synthetic data is generated from the model using the “true” bedrock topography, which is subsequently considered unknown as we solve the inverse problem to reconstruct it. The bedrock, ice thickness, accumulation/ablation, and log basal friction are taken from [31]. To avoid an “inverse crime,” we generate data by solving (17) on a mesh with time steps and interpolate the data onto a mesh with time steps to solve the inverse problem. Additionally, Gaussian noise is added to emulate the noise in satellite observations. Consequently, the discretized state dimension, which corresponds to the observational data dimension, is . The inversion parameter has dimension . Note, however, that this is still an ill-posed inverse problem due to smoothing properties of the forward operator. Furthermore, this is a challenging inverse problem due to the complexity of the nonlinear forward model and the auxiliary parameter uncertainty.
Figure 6 displays the terminal time ( years) surface height, the “true” bedrock topography used to synthesize data, and the estimated bedrock topography that results from solving the inverse problem with . The surface height is small near as the ice sheet is tapering down to the coastline on the western side of Greenland. The large bedrock topography in the top left corner corresponds to a mountain along the coast. We observe that the inversion correctly identifies the mountain, although it cannot resolve all of its local topography.
Preconditioned pseudo-time continuation
The inverse problem is computationally intensive as a result of the nonlinearities of the forward model and the large data dimension. The optimization algorithm required 26 iterations using a truncated CG Newton trust region method. To better understand uncertainty due to the auxiliary parameters, we seek to solve the inverse problem for perturbations . Since , many perturbations must be computed; consequently, solving the perturbed inverse problem solution quickly is crucial. For this article, we focus on one that perturbs the basal friction with the spatial field and the forcing term with the field .
Given our observations from the previous example, we use the modified Euler method as the predictor and take time steps. The preconditioner is initialized using the inverse of the regularization Hessian (i.e., ) and we set . This larger update rank is determined based on our observations from the optimization iteration history that implies the Hessian is not amenable for fast CG convergence (many optimization iterations required CG iterations). However, we found that only a total of vectors are stored over the three time steps as a result of filtering (with tolerance ) that reduces the number of vectors retained for the preconditioner.
We conduct three experiments using PCG tolerances . Figure 7 displays the number of PCG iterations required for the Hessian inversions. In each of the 3 time steps, we execute two Hessian inversions in the Predictor Step. As noted in Section 3.3, we can check for optimality after the predictor step to determine whether the corrector step is necessary. In this case, we found that no corrector steps were required.
For a proper computational cost comparison, note that nonlinear state solves are more expensive than the linear adjoint and incremental solves used for gradient and Hessian-vector product computation. Each state solve was executed using Newton’s method with an average of linear solves required per nonlinear solve. Hence, we will count cost in units of linear solves with a state solve counting as linear solves, a gradient evaluation as linear solve, and matrix-vector products with and as linear solves. Note that all linear solves have a comparable computational complexity since they involve similar coefficient matrices; see Appendix A for details regarding the coefficient matrices that arise in the adjoint method. The computational cost is summarized in Table 1. In this example, we observe that the computational cost is proportional to the PCG tolerance, with a loose tolerance being superior.
The comparison presented in Table 1 highlights that the computational cost does vary nontrivially with the selection of the PCG tolerance. However, it is noteworthy that the preconditioned pseudo-time continuation method converged with a coarse time step resolution and loose PCG tolerance. No algorithm parameter tuning was performed other than running with three different PCG tolerances for comparison. By contrast, we compared with re-optimization to solve the perturbed optimization problem. In this instance, the trust region radius and iteration limit for the perturbed optimization problem required tuning. We used the optimization algorithm parameters that were successful in the optimization problem. However, this resulted in a large optimization step which led to a lack of convergence from the forward solve. We reduced the trust region radius and increased the maximum number of optimization iterates to ensure convergence. Such challenges in algorithm parameter specification are common in practice for nonlinear optimization problems. Since preconditioned pseudo-time continuation takes the parameter perturbation gradually, it is less prone to convergence issues that require additional algorithm parameter tuning. This robustness is particularly important in uncertainty quantification contexts when minimizers are sought for many different perturbations .
5 Conclusion
The preconditioned pseudo-time continuation approach is superior to executing an optimization algorithm on the perturbed problem. The core concept is to gradually perturb the auxiliary parameter and adapt the optimization solution and Hessian preconditioner. This ensures desirable mathematical properties such as the gradient equaling zero (up to a tolerance) and the Hessian being positive definite at each time step. When combined with a predictor-corrector scheme, this produces a computationally efficient and reliable algorithm that requires little effort tuning algorithm parameters.
We observed that looser PCG tolerances generally yield lower computational costs. This algorithm parameter, denoted as , appears to be the most important consideration in our approach. Taking too small generally increases the computational cost. Taking too large may result in a poor step and thus require multiple Newton iterations for correction. This can increase the total computational cost. However, the predictor-corrector method converges as long as the time step is sufficiently small. Adaptive time stepping is another avenue of future research to enhance the algorithms. In particular, the adaptive algorithms of [8, chap. 5] are a promising avenue to leverage computations within the steps to estimate Lipschitz constants needed to automatically adapt the time step size.
One may be tempted to use the quasi-Newton inverse Hessian approximation in place of running PCG for the predictor step, corrector step, or both. However, error in the Hessian approximation may lead to running many corrector steps that exhibit slower convergence. Note also that in the context of optimization problems constrained by PDEs, the linear solves required to compute Hessian-vector products share common structures. Hence, the cost of computing multiple Hessian-vector products can be amortized by reusing matrix factorizations that arise from solving the incremental state and adjoint equations. Consequently, it is advantageous to take fewer iterations with more Hessian-vector products per iteration rather than using an approximate Hessian that requires many corrector iterations to ensure gradient norm tolerance satisfaction.
Our preconditioning approach reduces the computational cost measured in terms of the number of linearized PDE solves. However, it does introduce an additional storage requirement that is generally not required by a PCG implementation. For storage-limited applications, our approach may not be as effective. Nonetheless, the filtering tolerance, denoted as in our article, provides a natural mechanism to strike a balance between computational cost and storage. Furthermore, as we demonstrated in our numerical results, a modest storage budget is often sufficient to realize considerable benefits in computation cost reduction.
We demonstrated the benefit of using a second-order time integration method for the predictor step. The benefit is marginal in some instances, but we observed a significant benefit in scenarios with large auxiliary parameter perturbations and a small number of time steps. Higher-order time integration methods may also be considered. Based on the insight gained from our numerical results, we conjecture that higher-order time integrators may provide further improvements when the auxiliary parameter perturbation is very large. Exploration of such scenarios is an interesting line of future work.
Lastly, we note the potential benefit of our approach when many perturbed optimization problems must be solved repeatedly. For instance, in applications where auxiliary parameter samples are drawn and corresponding optimization problems must be solved. There is potential to design continuation paths and reuse preconditioners across samples so that we may solve the perturbed optimization problems at a computational cost that is much less than solving the problems independently. This is a topic of ongoing research as we explore ways to find a minimum number of paths necessary to pass through a set of auxiliary parameter samples. Using nonlinear paths may offer a considerable reduction in the number of paths that must be followed. However, the path designs must also consider the effects of nonlinearity in the ODE dynamics and the preconditioner updating.
Appendix A Computations within a time step
Within each time step, we require computation of {\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{k}={\mathchoice{\mbox{\displaystyle{\mathrm{B}}}}{\mbox{\textstyle{\mathrm{B}}}}{\mbox{\scriptstyle{\mathrm{B}}}}{\mbox{\scriptscriptstyle{\mathrm{B}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})\Delta\boldsymbol{\theta}, \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}), and {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})^{-1}{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{k}, where the latter is computed via PCG. In this appendix we outline these computations using the adjoint method and discuss the associated costs.
Assume that the parameter-to-observable map {\mathchoice{\mbox{\displaystyle{\mathrm{F}}}}{\mbox{\textstyle{\mathrm{F}}}}{\mbox{\scriptstyle{\mathrm{F}}}}{\mbox{\scriptscriptstyle{\mathrm{F}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) is defined in terms of an observation operator composed with the solution operator of a PDE. In particular, let {\mathchoice{\mbox{\displaystyle{u}}}{\mbox{\textstyle{u}}}{\mbox{\scriptstyle{u}}}{\mbox{\scriptscriptstyle{u}}}}\in\mathbb{R}^{m} be the discretized state variable, the discretized solution operator, and the observation operator. This way, {\mathchoice{\mbox{\displaystyle{\mathrm{F}}}}{\mbox{\textstyle{\mathrm{F}}}}{\mbox{\scriptstyle{\mathrm{F}}}}{\mbox{\scriptscriptstyle{\mathrm{F}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})=\mathcal{Q}(\mathcal{S}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})). Furthermore, we define the full-space objective as
[TABLE]
so that the reduced objective J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}) in (1) satisfies J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})=\mathcal{J}(\mathcal{S}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}),{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}).
Let {\mathchoice{\mbox{\displaystyle{c}}}{\mbox{\textstyle{c}}}{\mbox{\scriptstyle{c}}}{\mbox{\scriptscriptstyle{c}}}}({\mathchoice{\mbox{\displaystyle{u}}}{\mbox{\textstyle{u}}}{\mbox{\scriptstyle{u}}}{\mbox{\scriptscriptstyle{u}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})={\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}}\in\mathbb{R}^{m} denote the discretized PDE residual such {\mathchoice{\mbox{\displaystyle{c}}}{\mbox{\textstyle{c}}}{\mbox{\scriptstyle{c}}}{\mbox{\scriptscriptstyle{c}}}}(\mathcal{S}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}),\boldsymbol{\theta}) ={\mathchoice{\mbox{\displaystyle{0}}}{\mbox{\textstyle{0}}}{\mbox{\scriptstyle{0}}}{\mbox{\scriptscriptstyle{0}}}} that for all ({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta}). Then we may express the computation of \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}) and {\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{k}={\mathchoice{\mbox{\displaystyle{\mathrm{B}}}}{\mbox{\textstyle{\mathrm{B}}}}{\mbox{\scriptstyle{\mathrm{B}}}}{\mbox{\scriptscriptstyle{\mathrm{B}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})\Delta\boldsymbol{\theta} in terms of solves involving the Jacobian of with respect to , \nabla_{{\mathchoice{\mbox{\displaystyle{u}}}{\mbox{\textstyle{u}}}{\mbox{\scriptstyle{u}}}{\mbox{\scriptscriptstyle{u}}}}}{\mathchoice{\mbox{\displaystyle{c}}}{\mbox{\textstyle{c}}}{\mbox{\scriptstyle{c}}}{\mbox{\scriptscriptstyle{c}}}}; see Algorithm 3. In what follows, we denote the action of the PDE residual on the adjoint variable by {\mathchoice{\mbox{\displaystyle{c}}}{\mbox{\textstyle{c}}}{\mbox{\scriptstyle{c}}}{\mbox{\scriptscriptstyle{c}}}}({\mathchoice{\mbox{\displaystyle{u}}}{\mbox{\textstyle{u}}}{\mbox{\scriptstyle{u}}}{\mbox{\scriptscriptstyle{u}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})[{\mathchoice{\mbox{\displaystyle{\lambda}}}{\mbox{\textstyle{\lambda}}}{\mbox{\scriptstyle{\lambda}}}{\mbox{\scriptscriptstyle{\lambda}}}}]\vcentcolon={\mathchoice{\mbox{\displaystyle{c}}}{\mbox{\textstyle{c}}}{\mbox{\scriptstyle{c}}}{\mbox{\scriptscriptstyle{c}}}}({\mathchoice{\mbox{\displaystyle{u}}}{\mbox{\textstyle{u}}}{\mbox{\scriptstyle{u}}}{\mbox{\scriptscriptstyle{u}}}},{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}},\boldsymbol{\theta})^{\top}{\mathchoice{\mbox{\displaystyle{\lambda}}}{\mbox{\textstyle{\lambda}}}{\mbox{\scriptstyle{\lambda}}}{\mbox{\scriptscriptstyle{\lambda}}}}.
The computational cost of Algorithm 3 is approximately 1 state solve and 3 linearized (incremental) PDE solves, since the cost of other computations is negligible.
After executing Algorithm 3, we use \nabla_{\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}J({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k}), along with vectors computed in the previous time step, to perform the parametric quasi-Newton update (13). The computation of {\mathchoice{\mbox{\displaystyle{\mathrm{H}}}}{\mbox{\textstyle{\mathrm{H}}}}{\mbox{\scriptstyle{\mathrm{H}}}}{\mbox{\scriptscriptstyle{\mathrm{H}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})^{-1}{\mathchoice{\mbox{\displaystyle{\mathrm{B}}}}{\mbox{\textstyle{\mathrm{B}}}}{\mbox{\scriptstyle{\mathrm{B}}}}{\mbox{\scriptscriptstyle{\mathrm{B}}}}}({\mathchoice{\mbox{\displaystyle{m}}}{\mbox{\textstyle{m}}}{\mbox{\scriptstyle{m}}}{\mbox{\scriptscriptstyle{m}}}}_{k},\boldsymbol{\theta}_{k})\Delta\boldsymbol{\theta} via PCG requires Hessian-vector products, where is the number of PCG iterations. Algorithm 4 details the computation of one Hessian-vector product, whose total cost is 2 linearized PDE solves.
Appendix B Preconditioner initialization using low-rank approximation
Consider the generalized eigenvalue problem
[TABLE]
In ill-posed inverse problems, typically the eigenvalues decay rapidly. Also, since we have only finite-dimensional data, the number of measurements also constrains the data-misfit Hessian rank. Let us define the weighted inner product \left\langle{\cdot},{\cdot}\right\rangle_{{\mathchoice{\mbox{\displaystyle{\mathrm{R}}}}{\mbox{\textstyle{\mathrm{R}}}}{\mbox{\scriptstyle{\mathrm{R}}}}{\mbox{\scriptscriptstyle{\mathrm{R}}}}}} as
[TABLE]
Assuming the eigenvectors are normalized, i.e., \left\langle{{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{i}},{{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{j}}\right\rangle_{{\mathchoice{\mbox{\displaystyle{\mathrm{R}}}}{\mbox{\textstyle{\mathrm{R}}}}{\mbox{\scriptstyle{\mathrm{R}}}}{\mbox{\scriptscriptstyle{\mathrm{R}}}}}}=\delta_{ij}, for , we have
[TABLE]
Observing that \{{\mathchoice{\mbox{\displaystyle{\mathrm{R}}}}{\mbox{\textstyle{\mathrm{R}}}}{\mbox{\scriptstyle{\mathrm{R}}}}{\mbox{\scriptscriptstyle{\mathrm{R}}}}}^{\frac{1}{2}}{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{j}\}_{j=1}^{n} is an orthonormal basis and using the Sherman–Morrison–Woodbury formula, we obtain
[TABLE]
where {\mathchoice{\mbox{\displaystyle{V}}}{\mbox{\textstyle{V}}}{\mbox{\scriptstyle{V}}}{\mbox{\scriptscriptstyle{V}}}}=\begin{pmatrix}{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{1}&{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{2}&\cdots&{\mathchoice{\mbox{\displaystyle{v}}}{\mbox{\textstyle{v}}}{\mbox{\scriptstyle{v}}}{\mbox{\scriptscriptstyle{v}}}}_{n}\end{pmatrix}\in\mathbb{R}^{n\times n} and {\mathchoice{\mbox{\displaystyle{\Gamma}}}{\mbox{\textstyle{\Gamma}}}{\mbox{\scriptstyle{\Gamma}}}{\mbox{\scriptscriptstyle{\Gamma}}}}\in\mathbb{R}^{n\times n} is the diagonal matrix whose entry is .
Given a rank truncation of (20), where is chosen such that , we use the initial preconditioner
[TABLE]
where
[TABLE]
Acknowledgments
A. Alexanderian was supported in part by US National Science Foundation grant DMS #2111044. This work was supported by the Laboratory Directed Research and Development program at Sandia National Laboratories, a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc. for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This written work is authored by an employee of NTESS. The employee, not NTESS, owns the right, title and interest in and to the written work and is responsible for its contents. Any subjective views or opinions that might be expressed in the written work do not necessarily represent the views of the U.S. Government. The publisher acknowledges that the U.S. Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this written work or allow others to do so, for U.S. Government purposes. The DOE will provide public access to results of federally sponsored research in accordance with the DOE Public Access Plan.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] A. Alexanderian, J. Hart, and M. Stevens , A new perspective on parameter study of optimization problems , Applied Mathematics Letters, 140 (2023), p. 108548.
- 2[2] E. L. Allgower and K. Georg , Numerical Continuation Methods: An Introduction , Springer Series in Computational Mathematics, Springer Berlin, Heidelberg, 1 ed., 1990.
- 3[3] E. L. Allgower and K. Georg , Introduction to Numerical Continuation Methods , Classics in Applied Mathematics, SIAM, 2003.
- 4[4] E. Bueler, C. S. Lingle, J. A. Kallen-Brown, D. N. Covey, and L. N. Bowman , Exact solutions and verification of numerical models for isothermal ice sheets , Journal of Glaciology, 51 (2005), pp. 291–306.
- 5[5] T. Bui-Thanh, O. Ghattas, J. Martin, and G. Stadler , A computational framework for infinite-dimensional Bayesian inverse problems Part I: The linearized case, with application to global seismic inversion , SIAM Journal on Scientific Computing, 35 (2013), pp. A 2494–A 2523.
- 6[6] B. Crestel, G. Stadler, and O. Ghattas , A comparative study of structural similarity and regularization for joint inverse problems governed by PD Es , Inverse Problems, 35 (2018).
- 7[7] T. Cui, J. Martin, Y. M. Marzouk, A. Solonen, and A. Spantini , Likelihood-informed dimension reduction for nonlinear inverse problems , Inverse Problems, 30 (2014), pp. 1–28.
- 8[8] P. Deuflhard , Newton Methods for Nonlinear Problems , Springer Series in Computational Mathematics, Springer, 2006.