# Hyper-differential sensitivity analysis with respect to model discrepancy:
optimal solution updating

Joseph Hart11footnotemark: 1<EMAIL_ADDRESS>Bart van Bloemen
Waanders22footnotemark: 2<EMAIL_ADDRESS>
###### Abstract

A common goal throughout science and engineering is to solve optimization
problems constrained by computational models. However, in many cases a high-
fidelity numerical emulation of systems cannot be optimized due to code
complexity and computational costs which prohibit the use of intrusive and
many query algorithms. Rather, lower-fidelity models are constructed to enable
intrusive algorithms for large-scale optimization. As a result of the
discrepancy between high and low-fidelity models, optimal solutions determined
using low-fidelity models are frequently far from true optimality. In this
article we introduce a novel approach that uses post-optimality sensitivities
with respect to model discrepancy to update the optimization solution. Limited
high-fidelity data is used to calibrate the model discrepancy in a Bayesian
framework which in turn is propagated through post-optimality sensitivities of
the low-fidelity optimization problem. Our formulation exploits structure in
the post-optimality sensitivity operator to achieve computational scalability.
Numerical results demonstrate how an optimal solution computed using a low-
fidelity model may be significantly improved with limited evaluations of a
high-fidelity model.

###### keywords:

Hyper-differential sensitivity analysis, post-optimality sensitivity analysis,
PDE-constrained optimization, model discrepancy

††journal: Computer Methods in Applied Mechanics and Engineering

[1]organization=Sandia National Laboratories, addressline=P.O. Box 5800,
Albuquerque, NM 87185 , country=United States

[2]organization=Sandia National Laboratories, addressline=P.O. Box 5800,
Albuquerque, NM 87185 , country=United States

##  1 Introduction

Optimization problems constrained by computational models are ubiquitous in
science and engineering and have been the topic of extensive research [44, 5,
22, 45, 8, 9, 31, 28, 26, 7, 11, 29, 6, 3]. Nonetheless, optimization problems
are only as useful as the models that constrain them. In many applications,
high-fidelity models cannot be instrumented for optimization due to their
software complexity and computational cost. This is common in production codes
which are developed over many years with the goal of achieving model fidelity
rather than enabling optimization or other intrusive analysis approaches. In
such cases, a logical option is to construct low-fidelity models and thereby
enable efficient optimization. Assuming that the optimization problem can be
solved using a low-fidelity model, we pose the question: “How can a small
number of evaluations from the high-fidelity model be used to improve the
optimal solution?”. This would benefit a variety of complex optimization
problems for applications with multi-scale, multi-physics, nonlinear, and
multi-component features (such as climate, additive manufacturing, and
experimental fusion).

In this article we introduce a novel extension of hyper-differential
sensitivity analysis (HDSA) to compute the sensitivity of optimization
problems with respect to model discrepancy, i.e. the difference in the high
and low-fidelity models. HDSA is built on the foundation of post-optimality
sensitivity analysis [10, 17] which was originally developed in the context of
operations research and then extended to optimization constrained by partial
differential equations (PDEs) targeting stability analysis [21, 20, 12]. HDSA
was specifically created to evaluate the influence of uncertainties by scaling
post-optimality sensitivities to high-dimensions through a coupling of tools
from PDE-constrained optimization and numerical linear algebra [24, 40, 43,
23, 42, 25].

Our extension of HDSA provides a tool to predict how the optimal solution will
change given a perturbation of the model discrepancy. However, since the space
of possible model discrepancies is large (the discrepancy is an operator
between function spaces), it is critical to constrain the sensitivity
calculation by high-fidelity data. Accordingly we formulate a Bayesian inverse
problem to estimate the model discrepancy and compute sensitivities to update
the optimal solution in the direction of the maximum a posterior probability
point. A systematic integration of prior knowledge, high-fidelity data, and
low-fidelity optimization is enabled in a pragmatic way for large-scale
applications.

Calibration of model discrepancy in a Bayesian framework was proposed in the
seminal work of Kennedy and O’Hagan [30] to account for uncertainty and
correlations between the discrepancy and model parameters being estimated.
This was followed by many works which analyzed and extended their approach,
[37, 13, 27, 4, 34, 32, 18]. In the calibration context, the model discrepancy
is accessed through observed data with an implicit dependence on the
calibration parameters which are not known a priori. Our work is inspired by
these concepts but is significantly different in our assumption that a high-
fidelity model may be queried for different values of the optimization
variables as specified by the user. We focus on optimal control and design
problems, although model calibration may still be considered from a different
perspective than Kennedy and O’Hagan.

Our motivation to incorporate limited high-fidelity model evaluations to
improve analysis using a low-fidelity model is not new. Multi-fidelity methods
are motivated by a need to enable outer loop analysis such as optimization and
uncertainty quantification by combining evaluations of high and low-fidelity
models to achieve accuracy commensurate to high-fidelity analysis using a
smaller number of high-fidelity model evaluations [38]. There has been
considerable interest in the context of accelerating Monte Carlo type analysis
[19] for general uncertainty quantification applications as well as in the
specific context of optimization under uncertainty [36, 35]. There has also
been work to combine high and low-fidelity model evaluations within each
iterate of a deterministic optimization algorithm [14, 2, 1, 33]. Analogously,
we assume access to both a high and low-fidelity model. However, unlike the
previous references we seek a framework where the optimization problem does
not require access to the high-fidelity model. We also assume that derivatives
of the high-fidelity model are not available. Furthermore, our focus is on
cases where the number of high-fidelity solves is significantly constrained
and must be performed offline. We cannot provide theoretical guarantees which
multi-fidelity methods ensure through recourse on the high-fidelity model. Our
framework is suited for applications where multi-fidelity optimization is
infeasible due to the complexity of the high-fidelity model.

Our contributions in this article consist of:

  1. $\bullet$

In section 3 we provide a formulation and derivation to compute post-
optimality sensitivities with respect to model discrepancy. Furthermore, we
enable computational scalability by exploiting structure in the post-
optimality sensitivities to constraint the form of the model discrepancy.

  2. $\bullet$

In section 4 we pose a Bayesian inverse problem to estimate the discrepancy
(which is an operator) from high-fidelity data, determine a computationally
efficient expression for the posterior discrepancy, and propagate it through
the post-optimality sensitivities.

  3. $\bullet$

In sections 5 and 6 we present an approach to specify the prior and noise
model.

  4. $\bullet$

In section 6 we demonstrate the proposed approach on an illustrative example
and a more complex fluid flow example where the Stokes equation is used as a
low-fidelity approximation of the Naiver-Stokes equation.

##  2 Optimization Formulation

Consider optimization problems of the form

$\displaystyle\min_{z\in{\mathcal{Z}}}J(S(z),z)$ (1)

where $z$ denotes optimization variables in the Hilbert space ${\mathcal{Z}}$,
$S:{\mathcal{Z}}\to{\mathcal{U}}$ denotes the solution operator for a PDE with
state variable $u$ in an infinite dimensional Hilbert space ${\mathcal{U}}$,
and $J:{\mathcal{U}}\times{\mathcal{Z}}\to{\mathbb{R}}$ is the objective
function. This formulation is applicable to design, control, and inverse
problems. For simplicity we will refer to $z$ as the controller, but emphasize
that our formulation is general. We focus on problems where ${\mathcal{Z}}$ is
infinite dimensional.

In many applications, evaluating $S(z)$, i.e. solving the PDE, is
computationally intensive, so in practice we solve

$\displaystyle\min_{z\in{\mathcal{Z}}}J(\tilde{S}(z),z)$ (2)

where $\tilde{S}:{\mathcal{Z}}\to{\mathcal{U}}$ is the solution operator for a
lower-fidelity PDE. Common practice is to solve (2) and use limited
evaluations of $S(z)$ to assess the quality of the solution. Our goal in this
article is to use the high-fidelity model evaluations to improve the solution
computed in (2).

##  3 Post-optimality sensitivity with respect to model discrepancy

Our derivations will leverage the PDE discretization to define the model
discrepancy. To this end, we briefly introduce the essential components of the
discretization. Let $\\{\phi_{i}\\}_{i=1}^{m}$ be a basis for a finite
dimensional subspace ${\mathcal{U}}_{h}\subset{\mathcal{U}}$ and
$\\{\psi_{j}\\}_{j=1}^{n}$ be a basis for a finite dimensional subspace of
${\mathcal{Z}}_{h}\subset{\mathcal{Z}}$. These may, for instance, be finite
element basis functions. Let ${\bm{u}}\in{\mathbb{R}}^{m}$ and
${\bm{z}}\in{\mathbb{R}}^{n}$ denote coordinates of $u$ and $z$, respectively,
in these bases. The solution operators are discretized by
${\bm{S}}:{\mathbb{R}}^{n}\to{\mathbb{R}}^{m}$ and
$\tilde{{\bm{S}}}:{\mathbb{R}}^{n}\to{\mathbb{R}}^{m}$ and the objective
function is discretized by
${\bm{J}}:{\mathbb{R}}^{m}\times{\mathbb{R}}^{n}\to{\mathbb{R}}$.

To propagate model discrepancy through the optimization problem, let
${\bm{\delta}}({\bm{z}},{\bm{\theta}})$ be a function of ${\bm{z}}$,
parameterized by ${\bm{\theta}}$ (which will be defined below), and consider

$\displaystyle\min_{{\bm{z}}\in{\mathbb{R}}^{n}}\hskip
2.84526pt\widehat{{\bm{J}}}({\bm{z}},{\bm{\theta}}):={\bm{J}}(\tilde{{\bm{S}}}({\bm{z}})+{\bm{\delta}}({\bm{z}},{\bm{\theta}}),{\bm{z}}).$
(3)

The parameterized optimization problem (3) coincides with the low-fidelity
optimization problem when ${\bm{\delta}}\equiv\bm{0}$ and the high-fidelity
optimization problem when ${\bm{\delta}}={\bm{S}}-\tilde{{\bm{S}}}$. Let
${\bm{\theta}}_{\text{nom}}=\bm{0}$ and assume that ${\bm{\delta}}$ is
parameterized such that
${\bm{\delta}}({\bm{z}},{\bm{\theta}}_{\text{nom}})=\bm{0}$ $\forall{\bm{z}}$.

To facilitate our subsequent analysis, we assume that $\widehat{{\bm{J}}}$ is
twice continuously differentiable with respect to $({\bm{z}},{\bm{\theta}})$.
Let $\tilde{{\bm{z}}}\in{\mathbb{R}}^{n}$ denote a local minimum of the low-
fidelity problem, i.e. a minimizer of the parameterized problem (3) when
${\bm{\theta}}={\bm{\theta}}_{\text{nom}}$. We assume that $\tilde{{\bm{z}}}$
satisfies the first and second order optimality conditions,

$\displaystyle\nabla_{\bm{z}}\widehat{{\bm{J}}}(\tilde{{\bm{z}}},{\bm{\theta}}_{\text{nom}})=\bm{0}$
(4)

and
$\nabla_{{\bm{z}},{\bm{z}}}\widehat{{\bm{J}}}(\tilde{{\bm{z}}},{\bm{\theta}}_{\text{nom}})$
is positive definite, respectively, where $\nabla_{z}$ and
$\nabla_{{\bm{z}},{\bm{z}}}$ denote the gradient and Hessian with respect to
${\bm{z}}$, respectively.

Applying the Implicit Function Theorem to the first order optimality condition
(4), we consider an operator
${\bm{F}}:\mathcal{N}({\bm{\theta}}_{\text{nom}})\to\mathcal{N}(\tilde{{\bm{z}}})$,
defined on neighborhoods of ${\bm{\theta}}_{\text{nom}}$ and
$\tilde{{\bm{z}}}$, such that ${\bm{F}}({\bm{\theta}})$ is a stationary point
of (3), i.e.

$\nabla_{\bm{z}}\widehat{{\bm{J}}}({\bm{F}}({\bm{\theta}}),{\bm{\theta}})=\bm{0}\qquad\forall{\bm{\theta}}\in\mathcal{N}({\bm{\theta}}_{\text{nom}}).$

Assuming that
$\nabla_{{\bm{z}},{\bm{z}}}\widehat{{\bm{J}}}({\bm{F}}({\bm{\theta}}),{\bm{\theta}})$
is positive definite for all
${\bm{\theta}}\in\mathcal{N}({\bm{\theta}}_{\text{nom}})$, ${\bm{F}}$ maps the
model discrepancy to the optimal solution. The Jacobian of ${\bm{F}}$ with
respect to ${\bm{\theta}}$, i.e. the sensitivity of the optimal solution to
the model discrepancy, evaluated at ${\bm{\theta}}=\bm{0}$, is given by

$\displaystyle\nabla_{\bm{\theta}}{\bm{F}}(\bm{0})=-{\bm{H}}^{-1}{\bm{B}}\in{\mathbb{R}}^{n\times
p}$ (5)

where
${\bm{H}}=\nabla_{{\bm{z}},{\bm{z}}}\widehat{{\bm{J}}}(\tilde{{\bm{z}}},\bm{0})\in{\mathbb{R}}^{n\times
n}$ is the Hessian of $\widehat{{\bm{J}}}$ and
${\bm{B}}=\nabla_{{\bm{z}},{\bm{\theta}}}\widehat{{\bm{J}}}(\tilde{{\bm{z}}},\bm{0})\in{\mathbb{R}}^{n\times
p}$ is the Jacobian of $\nabla_{\bm{z}}\widehat{{\bm{J}}}$ with respect to
${\bm{\theta}}$. Both are evaluated at the low-fidelity solution
$(\tilde{{\bm{z}}},\bm{0})$. We interpret the Jacobian matrix
$\nabla_{{\bm{\theta}}}{\bm{F}}(\bm{0})$ as a Newton update of the optimal
solution $\tilde{{\bm{z}}}$ when the model discrepancy is perturbed.

To determine an appropriate parameterization of
${\bm{\delta}}({\bm{z}},{\bm{\theta}})$, we consider (in the function spaces)
a general operator mapping from $\mathcal{Z}_{h}$ to $\mathcal{U}_{h}$. Such a
operator can be written as

$\displaystyle z\mapsto\sum\limits_{i=1}^{m}f_{i}(z)\phi_{i}$

where $f_{i}:{\mathcal{Z}}_{h}\to{\mathbb{R}}$, $i=1,2,\dots,m$, denotes
functionals, which will be parameterized by ${\bm{\theta}}$. We realize a
simplification by observing that the sensitivity operator (5) only depends on
$({\bm{z}},{\bm{z}})$ and $({\bm{z}},{\bm{\theta}})$ derivatives of the
objective $\widehat{{\bm{J}}}$ evaluated at $(\tilde{{\bm{z}}},\bm{0})$, so
without loss of generality we assume that $f_{i}(z)$, $i=1,2,\dots,m$, are
affine functions of $z$. Using the Riesz representation for bounded linear
functionals we have a general form

$\displaystyle
z\mapsto\sum\limits_{i=1}^{m}(\theta_{i,0}+(z,a_{i}({\bm{\theta}}))_{\mathcal{Z}})\phi_{i}$

where $\theta_{i,0}\in{\mathbb{R}}$, $i=1,2,\dots,m$ and
$a_{i}({\bm{\theta}})\in Z_{h}$, $i=1,2,\dots,m$, are the Riesz representation
elements, which we parameterize with ${\bm{\theta}}$. A general expression for
elements in $Z_{h}$ is given by writing them as a linear combination of basis
functions as

$\displaystyle
a_{i}({\bm{\theta}})=\sum\limits_{j=1}^{n}\theta_{i,j}\psi_{j},$

where $\theta_{i,j}\in{\mathbb{R}}$. This yields an operator of the form

$\displaystyle
z\mapsto\sum\limits_{i=1}^{m}\left(\theta_{i,0}+\sum\limits_{j=1}^{n}\theta_{i,j}(z,\psi_{j})_{\mathcal{Z}}\right)\phi_{i}$

where the vector of coefficients is defined as
${\bm{\theta}}=({\bm{\theta}}_{0}^{T},{\bm{\theta}}_{1}^{T},\dots,{\bm{\theta}}_{m}^{T})^{T}\in{\mathbb{R}}^{p}$,
$p=m(n+1)$, where
${\bm{\theta}}_{0}=(\theta_{0,0},\theta_{1,0},\dots,\theta_{m,0})^{T}\in{\mathbb{R}}^{m}$
corresponds to the $m$ intercept terms and
${\bm{\theta}}_{i}=(\theta_{i,1},\theta_{i,2},\dots,\theta_{i,n})^{T}\in{\mathbb{R}}^{n}$,
$i=1,2,\dots,m$, corresponds to the $m$ linear functionals.

Transforming into the coordinate spaces, the model discrepancy
${\bm{\delta}}:{\mathbb{R}}^{n}\times{\mathbb{R}}^{p}\to{\mathbb{R}}^{m}$ can
be written in a convenient Kronecker product representation

$\displaystyle{\bm{\delta}}({\bm{z}},{\bm{\theta}})=\left(\begin{array}[]{cc}{\bm{I}}_{m}&{\bm{I}}_{m}\otimes{\bm{z}}^{T}{\bm{M}}_{z}\end{array}\right){\bm{\theta}}$
(7)

where ${\bm{I}}_{m}\in{\mathbb{R}}^{m\times m}$ is the identity matrix and
${\bm{M}}_{z}\in{\mathbb{R}}^{n\times n}$ is the mass matrix whose $(i,j)$
entry is $(\psi_{i},\psi_{j})_{Z}$. The Kronecker product structure of (7)
will prove critical in the subsequent analysis to achieve computational
scalability as $p$, the dimension of ${\bm{\theta}}$, is extremely large.

The post-optimality sensitivity operator (5) depends on two matrices, the
inverse Hessian ${\bm{H}}^{-1}$ and the mixed second derivative ${\bm{B}}$.
The Hessian ${\bm{H}}$ does not depend on our representation of the
discrepancy since it is evaluated at ${\bm{\theta}}=\bm{0}$ and
${\bm{\delta}}({\bm{z}},\bm{0})\equiv\bm{0}$. To determine
${\bm{B}}=\nabla_{{\bm{z}},{\bm{\theta}}}\widehat{{\bm{J}}}(\tilde{{\bm{z}}},\bm{0})$
we apply the Chain and Product rules to differentiate
$\widehat{{\bm{J}}}({\bm{z}},{\bm{\theta}})={\bm{J}}(\tilde{{\bm{S}}}({\bm{z}})+{\bm{\delta}}({\bm{z}},{\bm{\theta}}),{\bm{z}})$.
This gives

$\displaystyle\nabla_{{\bm{z}},{\bm{\theta}}}\widehat{{\bm{J}}}=\nabla_{\bm{z}}\tilde{{\bm{S}}}^{T}\nabla_{{\bm{u}},{\bm{u}}}{\bm{J}}\nabla_{\bm{\theta}}{\bm{\delta}}+\nabla_{\bm{z}}{\bm{\delta}}^{T}\nabla_{{\bm{u}},{\bm{u}}}J\nabla_{\bm{\theta}}{\bm{\delta}}+\nabla_{\bm{u}}{\bm{J}}\nabla_{{\bm{z}},{\bm{\theta}}}{\bm{\delta}}+\nabla_{{\bm{z}},{\bm{u}}}{\bm{J}}\nabla_{\bm{\theta}}{\bm{\delta}},$
(8)

where we adopt the convention that gradients are row vectors.

To simplify the subsequent analysis we assume333This is common on a wide range
of problems when the objective function admits the structure
$J(u,z)=J_{mis}(u)+J_{reg}(z)$ where $J_{mis}$ is a state misfit or design
criteria and $J_{reg}$ is an optimization variable regularization. that
$\nabla_{{\bm{z}},{\bm{u}}}{\bm{J}}=0$. Our derivations can be done in the
more general case but such additional algebraic manipulations do not add value
to the exposition. Observe that
$\nabla_{\bm{z}}{\bm{\delta}}(\tilde{{\bm{z}}},\bm{0})=\bm{0}$ since
${\bm{\delta}}({\bm{z}},\bm{0})\equiv\bm{0}$. The remaining nonzero
derivatives of ${\bm{\delta}}$ in (8) can be computed thanks to the form of
the discrepancy (7). They are given by

$\displaystyle\nabla_{\bm{\theta}}{\bm{\delta}}({\bm{z}},{\bm{\theta}})=\left(\begin{array}[]{cc}{\bm{I}}_{m}&{\bm{I}}_{m}\otimes{\bm{z}}^{T}{\bm{M}}_{z}\end{array}\right)\qquad\text{and}\qquad\nabla_{\bm{u}}{\bm{J}}\nabla_{{\bm{z}},{\bm{\theta}}}{\bm{\delta}}({\bm{z}},{\bm{\theta}})=\left(\begin{array}[]{cc}\bm{0}&\nabla_{\bm{u}}{\bm{J}}\otimes{\bm{M}}_{z}\end{array}\right).$

It follows that

$\displaystyle{\bm{B}}=\nabla_{\bm{z}}\tilde{S}^{T}\nabla_{{\bm{u}},{\bm{u}}}{\bm{J}}\left(\begin{array}[]{cc}{\bm{I}}_{m}&{\bm{I}}_{m}\otimes\tilde{{\bm{z}}}^{T}{\bm{M}}_{z}\end{array}\right)+\left(\begin{array}[]{cc}\bm{0}&\nabla_{\bm{u}}{\bm{J}}\otimes{\bm{M}}_{z}\end{array}\right)\in{\mathbb{R}}^{n\times
p}$ (12)

where $\nabla_{\bm{z}}\tilde{\bm{S}}$, $\nabla_{\bm{u}}{\bm{J}}$, and
$\nabla_{{\bm{u}},{\bm{u}}}{\bm{J}}$ are evaluated at
$(\tilde{{\bm{S}}}(\tilde{{\bm{z}}}),\tilde{{\bm{z}}})$. We emphasize that the
Kronecker product structure of ${\bm{\delta}}$ is preserved in ${\bm{B}}$.

##  4 Updating the optimal solution

Assume that the low-fidelity optimization problem has been solved to determine
$\tilde{{\bm{z}}}$, and that we have access to $N$ forward solves of the high-
fidelity model ${\bm{S}}$ for different inputs
$\\{{\bm{z}}_{\ell}\\}_{\ell=1}^{N}$. Let
${\bm{y}}_{\ell}={\bm{S}}({\bm{z}}_{\ell})-\tilde{{\bm{S}}}({\bm{z}}_{\ell})$,
$\ell=1,2,\dots,N,$ denote the corresponding evaluations of the model
discrepancy at these inputs. The number of evaluations $N$ will be small since
evaluating ${\bm{S}}({\bm{z}}_{\ell})$ is computationally costly. Our goal is
to use these limited high-fidelity evaluations to improve the optimal
solution.

If $N<<n$, the number of high-fidelity solves is much less than the dimension
of ${\bm{z}}$, then estimating the model discrepancy is ill-posed. We
formulate a Bayesian inverse problem to incorporate domain expertise (prior
knowledge) and uncertainty quantification into our estimation of the model
discrepancy parameters ${\bm{\theta}}$. We use $\bm{\Theta}$ to denote a
random vector taking values in ${\mathbb{R}}^{p}$ which models uncertainty in
the discrepancy parameters ${\bm{\theta}}$ and
${\bm{\delta}}({\bm{z}},\bm{\Theta})$ to denote the stochastic operator which
corresponds to the push forward of $\bm{\Theta}$ through the deterministic
operator ${\bm{\delta}}({\bm{z}},{\bm{\theta}})$. We likewise consider the
push forward of $\bm{\Theta}$ through the post-optimality sensitivity
operator, shifted by the low-fidelity optimal solution,

$\displaystyle\bm{U}=\tilde{{\bm{z}}}+\nabla_{{\bm{\theta}}}{\bm{F}}(\tilde{{\bm{z}}},\bm{0})\bm{\Theta}.$

The random vector $\bm{U}$ takes values in ${\mathbb{R}}^{n}$. The mean of
$\bm{U}$ is our updated optimal solution given the high-fidelity data. This
article focuses on computing the posterior mean of $\bm{\Theta}$ and
propagating it through the post-optimality sensitivities to compute the
posterior mean of $\bm{U}$. A follow on article will explore posterior
sampling to quantify uncertainty.

Figure 1 depicts our proposed approach to combine Bayesian inversion and post-
optimality sensitivities to update the optimal solution. In this section, we
detail our contributions to enable this workflow, which include:

  1. 1.

defining a prior distribution for $\bm{\Theta}$,

  2. 2.

formulating a Bayesian inverse problem to estimate the discrepancy,

  3. 3.

determining an efficient closed form expression for the posterior mean of
$\bm{\Theta}$, and

  4. 4.

propagating the posterior mean through the post-optimality sensitivity
operator to compute the mean of $\bm{U}$.

Due to the high dimensionality of the model discrepancy parameterization,
${\bm{\theta}}\in{\mathbb{R}}^{p}$, it is critical that the Kronecker product
structure of the discrepancy is preserved in the posterior of $\bm{\Theta}$.
We make judicious choices in the problem formulation and manipulate the linear
algebra to achieve this. The approach never requires computation in
${\mathbb{R}}^{p}$, and hence is efficient for large-scale optimization
problems.

Figure 1: Diagram of the proposed analysis pipeline. We denote the discrepancy
prior covariance and posterior mean as ${\bm{M}}_{\bm{\theta}}^{-1}$ and
$\overline{\bm{\theta}}$, respectively. $\tilde{{\bm{z}}}$ and
${\overline{{\bm{z}}}}$ denotes the low-fidelity and updated optimal
solutions, respectively.

###  4.1 Defining the prior

The parameters $\bm{\theta}$ correspond to the discretization of an operator
mapping between infinite dimensional function spaces. Infinite dimensional
Bayesian inverse problems [15, 39] have seen great advances in recent years.
However, most existing research focuses on estimating functions rather than
operators (a mapping between functions). This subsection details our approach
to define a prior on the space of operators by exploiting the Kronecker
product structure of ${\bm{\delta}}$, which is crucial to ensure computational
scalability.

To facilitate computation and ease specification of prior parameters, we take
a Gaussian prior for $\bm{\Theta}$ with mean $\bm{0}$ and a covariance matrix
to incorporate application knowledge into the discrepancy. Such priors are
common in infinite dimensional Bayesian inverse problems. Recall that a mean
zero Gaussian random vector $\bm{X}$ with covariance $\bm{\Pi}$ has a
probability density function proportional to
$\exp\left(-\frac{1}{2}\bm{x}^{T}\bm{\Pi}^{-1}\bm{x}\right)$. The highest
probability regions correspond to the points $\bm{x}$ for which the weighted
norm $||\bm{x}||_{\bm{\Pi}^{-1}}=\sqrt{\bm{x}^{T}\bm{\Pi}^{-1}\bm{x}}$ is
small. Hence, we define a prior covariance for $\bm{\Theta}$ by determining a
symmetric positive definite matrix which defines a norm on the space of
operators.

Since ${\bm{\theta}}\in{\mathbb{R}}^{p}$ does not have a physical
interpretation, we start by defining a prior in ${\mathbb{R}}^{m}$ (or
equivalently ${\mathcal{U}}_{h}$). Let ${\bm{L}}\in{\mathbb{R}}^{m\times m}$
be a symmetric positive definite matrix which defines a weighted norm on the
state space for which small norm elements will conform to discrepancy
characteristics specified from domain expertise. This may embed assumptions
about the smoothness of ${\bm{\delta}}$, boundary conditions, or penalize
deviations from known conservation properties. A common approach is to define
${\bm{L}}$ as the square of an elliptic differential operator so that
${\bm{L}}^{-1}$ corresponds to the discretization of a well-defined covariance
operator in the infinite dimensional state space [15, 39, 41].

Computing the ${\bm{L}}$-weighted inner product of
${\bm{\delta}}({\bm{z}},{\bm{\theta}})$ with itself and exploiting Kronecker
product structure we have

$\displaystyle({\bm{\delta}}({\bm{z}},{\bm{\theta}}),{\bm{\delta}}({\bm{z}},{\bm{\theta}}))_{\bm{L}}={\bm{\theta}}^{T}\left(\begin{array}[]{cc}{\bm{L}}&{\bm{L}}\otimes{\bm{z}}^{T}{\bm{M}}_{z}\\\
{\bm{L}}\otimes{\bm{M}}_{z}{\bm{z}}&{\bm{L}}\otimes{\bm{M}}_{\bm{z}}{\bm{z}}{\bm{z}}^{T}{\bm{M}}_{\bm{z}}\end{array}\right){\bm{\theta}}.$

To remove dependence on ${\bm{z}}$, we integrate
$({\bm{\delta}}({\bm{z}},{\bm{\theta}}),{\bm{\delta}}({\bm{z}},{\bm{\theta}}))_{\bm{L}}$
with respect to a Gaussian measure on ${\mathbb{R}}^{n}$ with mean
$\tilde{{\bm{z}}}$ and covariance matrix $\zeta^{2}{\bm{M}}_{\bm{z}}^{-1}$.
This measures the average size of the model discrepancy giving the greatest
weight to the low-fidelity solution $\tilde{{\bm{z}}}$ and decreases the
weight proportionally (with constant $\zeta$) with the controllers’ distance
(in ${\mathcal{Z}}_{h}$) from $\tilde{{\bm{z}}}$. Recalling the property of
Gaussians that the expectation of ${\bm{z}}{\bm{z}}^{T}$ is
$\zeta^{2}{\bm{M}}_{\bm{z}}^{-1}+\tilde{{\bm{z}}}\tilde{{\bm{z}}}^{T}$, and
simplifying algebra, yields the symmetric positive definite matrix

$\displaystyle{\bm{M}}_{\bm{\theta}}=\left(\begin{array}[]{cc}{\bm{L}}&{\bm{L}}\otimes\tilde{{\bm{z}}}^{T}{\bm{M}}_{z}\\\
{\bm{L}}\otimes{\bm{M}}_{z}\tilde{{\bm{z}}}&{\bm{L}}\otimes(\zeta^{2}{\bm{M}}_{\bm{z}}+({\bm{M}}_{\bm{z}}\tilde{{\bm{z}}})({\bm{M}}_{\bm{z}}\tilde{{\bm{z}}})^{T})\end{array}\right)\in{\mathbb{R}}^{p\times
p},$ (15)

for which

$\displaystyle{\bm{\theta}}^{T}{\bm{M}}_{\bm{\theta}}{\bm{\theta}}=\int_{{\mathbb{R}}^{n}}||{\bm{\delta}}({\bm{z}},{\bm{\theta}})||_{{\bm{L}}}^{2}d\mu({\bm{z}}),$
(16)

where $\mu$ is the Gaussian measure on ${\mathbb{R}}^{n}$ with mean
$\tilde{{\bm{z}}}$ and covariance $\zeta^{2}{\bm{M}}_{\bm{z}}^{-1}$.

Hence ${\bm{M}}_{\bm{\theta}}$ defines an inner product for ${\bm{\theta}}$ to
measure the size of the model discrepancy
${\bm{\delta}}({\bm{z}},{\bm{\theta}})$ according to our prior knowledge
imposed by ${\bm{L}}$. Let ${\bm{M}}_{\bm{\theta}}^{-1}$ be the prior
covariance for $\bm{\Theta}$, which ensures that the prior is rooted in
established theory relating known physical properties to covariance matrices
which can be efficiently manipulated in the state space.

As a building block for our subsequent analysis, note that we can decompose
${\bm{M}}_{\bm{\theta}}=\bm{C}\bm{C}^{T}$ where

$\displaystyle\bm{C}=\left(\begin{array}[]{cc}{\bm{L}}^{\frac{1}{2}}&\bm{0}\\\
{\bm{L}}^{\frac{1}{2}}\otimes{\bm{M}}_{\bm{z}}\tilde{{\bm{z}}}&{\bm{L}}^{\frac{1}{2}}\otimes\zeta{\bm{M}}_{\bm{z}}^{\frac{1}{2}}\end{array}\right)$
(19)

and

$\displaystyle\bm{C}^{-1}=\left(\begin{array}[]{cc}{\bm{L}}^{-\frac{1}{2}}&\bm{0}\\\
{\bm{L}}^{-\frac{1}{2}}\otimes\left(-\zeta^{-1}{\bm{M}}_{\bm{z}}^{\frac{1}{2}}\tilde{{\bm{z}}}\right)&{\bm{L}}^{-\frac{1}{2}}\otimes\zeta^{-1}{\bm{M}}_{\bm{z}}^{-\frac{1}{2}}\end{array}\right).$
(22)

This factorization enables efficient manipulations in what follows.

###  4.2 Bayesian inverse problem

We formulate a Bayesian inverse problem to estimate ${\bm{\delta}}$ using the
data pairs $\\{{\bm{z}}_{\ell},{\bm{y}}_{\ell}\\}_{\ell=1}^{N}$. To facilitate
the analysis we assume:

###### Assumption 1

$\\{{\bm{z}}_{\ell}\\}_{\ell=1}^{N}$ is a linearly independent set of vectors
in ${\mathbb{R}}^{n}$.

This is usually satisfied in practice as one explores the model space in
different directions since the number of model evaluations $N$ is small
relative to the dimension of ${\bm{z}}$. For notational simplicity, we define

$\displaystyle{\bm{A}}_{\ell}=\begin{pmatrix}{\bm{I}}_{m}&{\bm{I}}_{m}\otimes{\bm{z}}_{\ell}^{T}{\bm{M}}_{z}\end{pmatrix}\in{\mathbb{R}}^{m\times
p},\qquad\ell=1,2,\dots,N,$ (23)

so that
${\bm{\delta}}({\bm{z}}_{\ell},{\bm{\theta}})={\bm{A}}_{\ell}{\bm{\theta}}$,
and the concatenation of these matrices

$\displaystyle{\bm{A}}=\left(\begin{array}[]{c}{\bm{A}}_{1}\\\ {\bm{A}}_{2}\\\
\vdots\\\ {\bm{A}}_{N}\end{array}\right)\in{\mathbb{R}}^{mN\times p}.$

${\bm{A}}{\bm{\theta}}\in{\mathbb{R}}^{mN}$ corresponds to the evaluation of
${\bm{\delta}}({\bm{z}},{\bm{\theta}})$ for all of the controller data
$\\{{\bm{z}}_{\ell}\\}_{\ell=1}^{N}$. In an analogous fashion, define
$\bm{b}\in{\mathbb{R}}^{mN}$ by stacking ${\bm{y}}_{\ell}$,
$\ell=1,2,\dots,N$, into a vector. Then we seek
${\bm{A}}{\bm{\theta}}\approx\bm{b}$.

To enable a closed form expression for the posterior we consider an additive
Gaussian noise model with mean $\bm{0}$ and covariance $\alpha\bm{I}$, where
$\alpha>0$ is specified by the user. Typically the noise covariance in a
Bayesian inverse problem is defined using knowledge of the data collection
process. In cases where data is collected from simulation we do not have a
clear definition of noise. However, ${\bm{\delta}}$ is defined as a linear
approximation of ${\bm{S}}({\bm{z}})-\tilde{{\bm{S}}}({\bm{z}})$, which is, in
general, a nonlinear operator. Hence the noise may be interpreted as the
approximation error due to the linearization of the discrepancy. Nonetheless,
one can show that the range space of ${\bm{A}}$ equals ${\mathbb{R}}^{mN}$, so
there exists infinitely many ${\bm{\theta}}\in{\mathbb{R}}^{p}$ such that
${\bm{A}}{\bm{\theta}}=\bm{b}$. The choice of $\alpha$ dictates the weight
given to the data misfit relative to the prior. Taking a small $\alpha$ will
drive the inverse problem toward interpolation with little regard for the
prior. We may experiment with different $\alpha$’s at a modest computational
cost.

Given our formulation with Gaussian prior and noise models, and the linearity
of ${\bm{\delta}}({\bm{z}},{\bm{\theta}})$, the posterior is Gaussian with a
negative log probability density (unnormalized) function

$\displaystyle\frac{1}{2\alpha}\left({\bm{A}}{\bm{\theta}}-\bm{b}\right)^{T}\left({\bm{A}}{\bm{\theta}}-\bm{b}\right)+\frac{1}{2}{\bm{\theta}}^{T}{\bm{M}}_{\bm{\theta}}{\bm{\theta}}.$
(25)

The posterior mean and covariance are given by

$\displaystyle\overline{{\bm{\theta}}}=\frac{1}{\alpha}\bm{\Sigma}{\bm{A}}^{T}\bm{b}\qquad\text{and}\qquad\bm{\Sigma}=\left({\bm{M}}_{\bm{\theta}}+\frac{1}{\alpha}{\bm{A}}^{T}{\bm{A}}\right)^{-1},$
(26)

respectively. Next we derive an explicit and computationally efficient
expression for $\overline{{\bm{\theta}}}$ and its image under the post-
optimality sensitivity operator.

###  4.3 Expression for the posterior mean

The posterior mean, given in (26), involves the posterior covariance matrix
which in turn involves the inverse of a sum of matrices in
${\mathbb{R}}^{p\times p}$. Since these matrices cannot be manipulated
directly, we seek to derive a computable expression for the posterior
covariance $\bm{\Sigma}$ and in turn the posterior mean
$\overline{{\bm{\theta}}}$. Emphasis is given to preserving the Kronecker
product structure throughout our analysis so that computations in
${\mathbb{R}}^{p}$ are avoided. The derivation is described in three steps in
the following subsections:

  1. 1.

factorize ${\bm{A}}$ in a ${\bm{M}}_{\bm{\theta}}$ orthogonal basis to rewrite
$\bm{\Sigma}^{-1}$,

  2. 2.

invert $\bm{\Sigma}^{-1}$ by exploiting orthogonality to invert a sum,

  3. 3.

compute
$\overline{{\bm{\theta}}}=\frac{1}{\alpha}\bm{\Sigma}{\bm{A}}^{T}\bm{b}$ using
Kronecker algebra.

#### Factorize ${\bm{A}}$

To rewrite the inverse posterior covariance,
${\bm{M}}_{\bm{\theta}}+\frac{1}{\alpha}{\bm{A}}^{T}{\bm{A}}$, in a form
amenable for inversion, we decompose ${\bm{A}}$ with respect to a chosen inner
product. In particular, the generalized singular value decomposition (GSVD) of
${\bm{A}}$ is used with the inner product on the column space weighted by
${\bm{M}}_{\bm{\theta}}$. The GSVD is given by

$\displaystyle{\bm{A}}=\bm{\Xi}\bm{\Phi}\bm{\Psi}^{T}{\bm{M}}_{\bm{\theta}}$

where $\bm{\Phi}$ is the diagonal matrix of singular values, and $\bm{\Xi}$
and $\bm{\Psi}$ are matrices containing the left and right singular vectors,
respectively, which satisfy $\bm{\Xi}^{T}\bm{\Xi}={\bm{I}}$ and
$\bm{\Psi}^{T}{\bm{M}}_{\bm{\theta}}\bm{\Psi}={\bm{I}}$. To determine the
singular vectors, note that
${\bm{A}}{\bm{M}}_{\bm{\theta}}^{-1}{\bm{A}}^{T}=\bm{\Xi}\bm{\Phi}^{2}\bm{\Xi}^{T}$.
Using the expressions (23) for ${\bm{A}}$ and (22) for
${\bm{M}}_{\bm{\theta}}^{-1}={\bm{C}}^{-T}{\bm{C}}^{-1}$, we can write

$\displaystyle{\bm{A}}{\bm{M}}_{\bm{\theta}}^{-1}{\bm{A}}^{T}=\bm{G}\otimes{\bm{L}}^{-1}$

where

$\displaystyle{\bm{G}}={\bm{e}}{\bm{e}}^{T}+\zeta^{-2}(\bm{Z}-\tilde{{\bm{z}}}{\bm{e}}^{T})^{T}{\bm{M}}_{\bm{z}}(\bm{Z}-\tilde{{\bm{z}}}{\bm{e}}^{T})\in{\mathbb{R}}^{N\times
N},$

$\bm{Z}=\begin{pmatrix}{\bm{z}}_{1}&{\bm{z}}_{2}&\dots&{\bm{z}}_{N}\end{pmatrix}$
is the matrix of controller data, and ${\bm{e}}\in{\mathbb{R}}^{N}$ is the
vector of ones. Hence the left singular vectors $\bm{\Xi}$ correspond to the
eigenvectors of $\bm{G}\otimes{\bm{L}}^{-1}$ and the squared singular values
$\bm{\Phi}^{2}$ correspond to the eigenvalues.

Denoting the eigenvectors and eigenvalues of $\bm{G}$ and ${\bm{L}}$ with
$(\bm{g}_{i},\lambda_{i})$ and $(\bm{l}_{j},\rho_{j})$, respectively, and
recalling properties of the eigenvalue decomposition of a Kronecker product,
we observe that the squared generalized singular values of $\bm{A}$ are given
by $\frac{\lambda_{i}}{\rho_{j}}$ and are associated with the left singular
vectors $\bm{\xi}_{i,j}=\bm{g}_{i}\otimes\bm{l}_{j}$. Rewriting the GSVD to
solve for $\bm{\Psi}$, the right singular vector associated with
$\bm{\xi}_{i,j}$ is given by

$\displaystyle\bm{\psi}_{i,j}=\frac{1}{\sqrt{\lambda_{i}\rho_{j}}}\left(\begin{array}[]{cc}s_{i}\bm{l}_{j}\\\
\bm{l}_{j}\otimes\zeta^{-2}\bm{w}_{i}\end{array}\right),\qquad i=1,2,\dots,N\
\ j=1,2,\dots,m,$ (29)

where

$\displaystyle\bm{w}_{i}=\bm{Z}\bm{g}_{i}-({\bm{e}}^{T}\bm{g}_{i})\tilde{{\bm{z}}}\qquad\text{and}\qquad
s_{i}=({\bm{e}}^{T}\bm{g}_{i})-\zeta^{-2}\bm{w}_{i}^{T}{\bm{M}}_{\bm{z}}\tilde{{\bm{z}}}.$
(30)

#### Invert $\bm{\Sigma}^{-1}$

Given this decomposition of ${\bm{A}}$ and our previously noted decomposition
${\bm{M}}_{\bm{\theta}}=\bm{C}\bm{C}^{T}$, we express $\bm{\Sigma}^{-1}$ as

$\displaystyle\bm{\Sigma}^{-1}=\frac{1}{\alpha}\bm{C}\bm{X}\bm{C}^{T}$

where

$\displaystyle\bm{X}=\alpha\bm{I}+\bm{C}^{T}\bm{\Psi}\bm{\Phi}^{2}\bm{\Psi}^{T}\bm{C}.$
(31)

This factorization implies that
$\bm{\Sigma}=\alpha\bm{C}^{-T}\bm{X}^{-1}\bm{C}^{-1}$ and facilitates easy
manipulation because ${\bm{X}}$ is a combination of diagonal and orthogonal
matrices. Applying the Sherman-Morrison-Woodbury formula to (31), we have

$\displaystyle\bm{X}^{-1}=\frac{1}{\alpha}\left(\bm{I}-\bm{C}^{T}\bm{\Psi}\bm{D}\bm{\Psi}^{T}\bm{C}\right)$

where $\bm{D}\in{\mathbb{R}}^{mN\times mN}$ is a diagonal matrix whose entries
are given by $\frac{\lambda_{i}}{\lambda_{i}+\alpha\rho_{j}}.$ Algebraic
simplifications yield

$\displaystyle\bm{\Sigma}={\bm{M}}_{\bm{\theta}}^{-1}-\bm{\Psi}\bm{D}\bm{\Psi}^{T}.$
(32)

The posterior covariance may be interpreted as the prior covariance
${\bm{M}}_{\bm{\theta}}^{-1}$ with its uncertainty reduced in the directions
of the columns of $\bm{\Psi}$ (29) with weight
$\frac{\lambda_{i}}{\lambda_{i}+\alpha\rho_{j}}$, a combination of the
eigenvalues from ${\bm{L}}$ and ${\bm{G}}$, along with the noise covariance
$\alpha$.

#### Compute the posterior mean

Recalling the form of the posterior mean in (26) and covariance in (32), and
computing the product

$\displaystyle{\bm{A}}^{T}{\bm{b}}=\sum\limits_{\ell=1}^{N}\left(\begin{array}[]{c}{\bm{y}}_{\ell}\\\
{\bm{y}}_{\ell}\otimes{\bm{M}}_{\bm{z}}{\bm{z}}_{\ell}\end{array}\right),$

we have

$\displaystyle\overline{{\bm{\theta}}}=\frac{1}{\alpha}\sum\limits_{\ell=1}^{N}\left({\bm{M}}_{\bm{\theta}}^{-1}-\bm{\Psi}\bm{D}\bm{\Psi}^{T}\right)\left(\begin{array}[]{c}{\bm{y}}_{\ell}\\\
{\bm{y}}_{\ell}\otimes{\bm{M}}_{\bm{z}}{\bm{z}}_{\ell}\end{array}\right).$

Using ${\bm{M}}_{\bm{\theta}}^{-1}={\bm{C}}^{-T}{\bm{C}}^{-1}$, recalling the
form of $\bm{\Psi}$’s columns (29), and writing matrix-vector products with
$\bm{\Psi}\bm{D}\bm{\Psi}^{T}$ in terms of linear solves involving ${\bm{L}}$,
we arrive at the posterior mean expression

$\displaystyle\overline{{\bm{\theta}}}=\frac{1}{\alpha}\sum\limits_{\ell=1}^{N}\left[\left(\begin{array}[]{c}a_{\ell}{\bm{u}}_{\ell}\\\
{\bm{u}}_{\ell}\otimes\zeta^{-2}({\bm{z}}_{\ell}-\tilde{{\bm{z}}})\end{array}\right)-\sum\limits_{i=1}^{N}b_{i,\ell}\left(\begin{array}[]{c}s_{i}\bm{u}_{i,\ell}\\\
\bm{u}_{i,\ell}\otimes\zeta^{-2}\bm{w}_{i}\end{array}\right)\right]$

involving the constants

$\displaystyle
a_{\ell}=1-\zeta^{-2}\tilde{{\bm{z}}}^{T}{\bm{M}}_{\bm{z}}({\bm{z}}_{\ell}-\tilde{{\bm{z}}})\qquad\text{and}\qquad
b_{i,\ell}=\zeta^{-2}({\bm{z}}_{\ell}-\tilde{{\bm{z}}})^{T}{\bm{M}}_{\bm{z}}\bm{Z}\bm{g}_{i}+({\bm{e}}^{T}\bm{g}_{i})a_{\ell}$

and the vectors

$\displaystyle\bm{u}_{\ell}={\bm{L}}^{-1}{\bm{y}}_{\ell}\qquad\text{and}\qquad\bm{u}_{i,\ell}=\left(\alpha{\bm{L}}+\lambda_{i}{\bm{I}}\right)^{-1}{\bm{u}}_{\ell}.$

We observe that $\overline{{\bm{\theta}}}$ depends on linear solves 444In
Section 6, we demonstrate how to efficiently invert ${\bm{L}}$ and
$\left(\alpha{\bm{L}}+\lambda_{i}{\bm{I}}\right)$ using a generalized singular
value decomposition when ${\bm{L}}$ is defined as the square of an elliptic
operator. involving ${\bm{L}}$ and $\alpha{\bm{L}}+\lambda_{i}{\bm{I}}$ where
the right hand sides arise from the data $\\{{\bm{y}}_{\ell}\\}_{\ell=1}^{N}$.
The posterior mean can be interpreted as a linear combination of the
discrepancy data preconditioned by the prior weighted by the noise covariance
$\alpha$ and the eigenpairs $(\lambda_{i},{\bm{g}}_{i})$ from the controller
data informed matrix ${\bm{G}}$.

###  4.4 Optimal solution update

To propagate $\overline{{\bm{\theta}}}$ through the post-optimality
sensitivity operator $\nabla_{{\bm{\theta}}}{\bm{F}}(\bm{0})$, we note the
Kronecker product structure of both ${\bm{B}}$ and $\overline{{\bm{\theta}}}$
which allows for computation of the post-optimality sensitivities using solves
in ${\mathbb{R}}^{m}$ and ${\mathbb{R}}^{n}$. The updated optimal solution is
given by $\tilde{{\bm{z}}}-{\bm{H}}^{-1}{\bm{B}}\overline{{\bm{\theta}}}$
where

$\displaystyle{\bm{B}}\overline{{\bm{\theta}}}$
$\displaystyle=\frac{1}{\alpha}\nabla_{{\bm{z}}}\tilde{{\bm{S}}}^{T}\nabla_{{\bm{u}},{\bm{u}}}{\bm{J}}\left[\sum\limits_{\ell=1}^{N}\left({\bm{u}}_{\ell}-\sum\limits_{i=1}^{N}b_{i,\ell}({\bm{e}}^{T}\bm{g}_{i})\bm{u}_{i,\ell}\right)\right]$
(34)
$\displaystyle+\frac{1}{\alpha}\sum\limits_{\ell=1}^{N}(\nabla_{\bm{u}}J{\bm{u}}_{\ell})\zeta^{-2}{\bm{M}}_{\bm{z}}({\bm{z}}_{\ell}-\tilde{{\bm{z}}})$
$\displaystyle-\frac{1}{\alpha}\sum\limits_{\ell=1}^{N}\sum\limits_{i=1}^{N}b_{i,\ell}(\nabla_{\bm{u}}{\bm{J}}\bm{u}_{i,\ell})\zeta^{-2}{\bm{M}}_{\bm{z}}\bm{w}_{i}.$

Equation (34) provides a systematic and interpretable combination of high-
fidelity data preconditioned by the prior in $\bm{u}_{\ell}$ and
$\bm{u}_{i,\ell}$, the low-fidelity model in $\tilde{\bm{S}}_{\bm{z}}^{T}$,
and the optimization objective in $\nabla_{\bm{u}}{\bm{J}}$ and
$\nabla_{{\bm{u}},{\bm{u}}}{\bm{J}}$.

##  5 Specifying the prior

The proposed approach requires specification of the matrix ${\bm{L}}$, the
scalar $\zeta$, and the scalar $\alpha$. This is necessary because the inverse
problem to estimate ${\bm{\delta}}$ is underdetermined and hence is posed in a
Bayesian formulation. These parameters play an important role imposing domain
expertise to constraint the discrepancy. This section explores a strategy to
choose them based on physical characteristics of the problem.

The noise covariance $\alpha>0$ normalizes the data misfit in the likelihood
thus weighting the data fit relative to the prior (25). Typically in Bayesian
inverse problems $\alpha$ is specified based on knowledge of the data
fidelity. The noise covariance in our framework corresponds to the
approximation error due to linearization of the discrepancy. As a rule of
thumb, $\alpha$ may be chosen based on the magnitude of the discrepancy data.
For instance, if the discrepancy data ${\bm{y}}_{\ell}$ has magnitude
$\mathcal{O}(1)$ then $\alpha=10^{-2}$ normalizes the difference
${\bm{y}}_{\ell}-{\bm{\delta}}({\bm{z}}_{\ell},{\bm{\theta}})$ by
$\sqrt{\alpha}=10^{-1}$ or equivalently $10\%$ of the discrepancy magnitude.
Comparing ${\bm{\delta}}({\bm{z}}_{\ell},\overline{{\bm{\theta}}})$ with
${\bm{y}}_{\ell}$ provides an easy check that the choice was reasonable.

Sampling the prior discrepancy aids in determining ${\bm{L}}$ and $\zeta$ by
analyzing the physical characteristics of the discrepancy samples. The matrix
${\bm{L}}$ imposes spatiotemporal characteristics (variance, correlation
length, boundary conditions, etc.) on the state space and the scalar $\zeta$
defines length scales on the controller space to affect how the discrepancy
varies with respect to ${\bm{z}}$. A sampling procedure is detailed below to
facilitate the specification of ${\bm{L}}$ and $\zeta$ independently of one
another.

Using the factorization ${\bm{M}}_{\bm{\theta}}^{-1}=\bm{C}^{-T}\bm{C}^{-1}$,
a prior sample of $\bm{\Theta}$ is given by $\bm{C}^{-T}\bm{\omega}_{p}$,
where $\bm{\omega}_{p}$ is a $p$ dimensional standard normal random vector.
Evaluating the model discrepancy at such samples yields

$\displaystyle{\bm{\delta}}({\bm{z}},\bm{C}^{-T}\bm{\omega}_{p})={\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{0}+\zeta^{-1}{\bm{L}}^{-\frac{1}{2}}\bm{\Omega}{\bm{M}}_{\bm{z}}^{\frac{1}{2}}({\bm{z}}-\tilde{{\bm{z}}})$

where $\bm{\omega}_{0}\in{\mathbb{R}}^{m}$ and
$\bm{\Omega}\in{\mathbb{R}}^{m\times n}$ have entries that are independent
identically distributed (i.i.d) samples from a standard normal distribution.

To avoid forming the dense matrix $\bm{\Omega}\in{\mathbb{R}}^{m\times n}$, we
evaluate ${\bm{\delta}}({\bm{z}},\bm{C}^{-T}\bm{\omega})$ at particular
${\bm{z}}$ values and visualize the resulting discrepancy, which is an element
of the state space. Let ${\bm{z}}_{r}\in{\mathbb{R}}^{n}$ denote a “reference
point” corresponding to a physically reasonable controller. We will evaluate
the discrepancy at a sequence of points
$\hat{{\bm{z}}}_{k}=\tilde{{\bm{z}}}+\frac{k}{K}({\bm{z}}_{r}-\tilde{{\bm{z}}})$,
$k=0,1,\dots,K$, where $K$ is an integer. Then for a fixed $\bm{\Theta}$
sample, or equivalently a fixed $\bm{\omega}_{p}$, we have

$\displaystyle{\bm{\delta}}(\hat{{\bm{z}}}_{k},\bm{C}^{-T}\bm{\omega}_{p})={\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{0}+\frac{k}{K}\zeta^{-1}{\bm{L}}^{-\frac{1}{2}}\bm{\Omega}{\bm{M}}_{\bm{z}}^{\frac{1}{2}}({\bm{z}}_{r}-\tilde{{\bm{z}}}).$

Rather than forming the large matrix $\bm{\Omega}\in{\mathbb{R}}^{m\times n}$
of i.i.d. random numbers, observe that
$\bm{\Omega}{\bm{M}}_{\bm{z}}^{\frac{1}{2}}({\bm{z}}_{r}-\tilde{{\bm{z}}})$ is
a Gaussian random vector in ${\mathbb{R}}^{m}$ whose entires are i.i.d.
Gaussians with mean $0$ and variance
$({\bm{z}}_{r}-\tilde{{\bm{z}}})^{T}{\bm{M}}_{\bm{z}}({\bm{z}}_{r}-\tilde{{\bm{z}}})$.
Thus to compute ${\bm{\delta}}(\hat{{\bm{z}}}_{k},\bm{C}^{-T}\bm{\omega}_{p})$
for each $k=0,1,\dots,K$, it is sufficient to apply ${\bm{L}}^{-\frac{1}{2}}$
to two random vectors $\bm{\omega}_{0},\bm{\omega}_{1}\in{\mathbb{R}}^{m}$,
sampled with i.i.d. standard normal entries, and evaluate

$\displaystyle{\bm{\delta}}(\hat{{\bm{z}}}_{k},\bm{C}^{-T}\bm{\omega}_{p})={\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{0}+\frac{k}{K}\zeta^{-1}\sqrt{({\bm{z}}_{r}-\tilde{{\bm{z}}})^{T}{\bm{M}}_{\bm{z}}({\bm{z}}_{r}-\tilde{{\bm{z}}})}{\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{1}.$

Our approach is summarized in Algorithm 1 where $S$ samples are generated by
looping over different standard normal random vectors
$\bm{\omega}_{0}^{s},\bm{\omega}_{1}^{s}$, $s=1,2,\dots,S$. Examination of
$\\{{\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{0}^{s}\\}_{s=1}^{S}$ and
$\\{{\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{1}^{s}\\}_{s=1}^{S}$ aids in tuning of
parameters in ${\bm{L}}$ so that the discrepancy realizations are physically
reasonable. Similarly, $\zeta$ may be tuned based on variability in
$\\{{\bm{\delta}}(\hat{{\bm{z}}}_{k},\bm{C}^{-T}\bm{\omega}_{p}^{s})\\}_{k=1}^{K}$
for each fixed $\bm{\omega}_{p}^{s}$, where
$\bm{\omega}_{p}^{s}\in{\mathbb{R}}^{p}$ is not computed but is used to denote
the sample corresponding to
$\bm{\omega}_{0}^{s},\bm{\omega}_{1}^{s}\in{\mathbb{R}}^{m}$.

Algorithm 1 Compute samples of the prior discrepancy.

1: Input: $S\in\mathbb{N},K\in\mathbb{N}$, and
${\bm{z}}_{r}\in{\mathbb{R}}^{n}$

2: Compute
$c=\sqrt{({\bm{z}}_{r}-\tilde{{\bm{z}}})^{T}{\bm{M}}_{\bm{z}}({\bm{z}}_{r}-\tilde{{\bm{z}}})}$

3: for $s=1,2,\dots,S$ do

4: Sample $\bm{\omega}_{0},\bm{\omega}_{1}\in{\mathbb{R}}^{m}$ with i.i.d.
standard normal entries

5: Compute ${\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{0}$ and
${\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{1}$

6: for $k=0,1,\dots,K$ do

7: Evaluate
${\bm{\delta}}(\hat{{\bm{z}}}_{k},\bm{C}^{-T}\bm{\omega}^{s})={\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{0}+\frac{k}{K}\zeta^{-1}c{\bm{L}}^{-\frac{1}{2}}\bm{\omega}_{1}$

8: end for

9: end for

##  6 Numerical results

We demonstrate our approach on an illustrative 1D control example, where model
discrepancy arises from a failure to include advection, and on an application
of fluid flow in 2D, where Stokes is used as the low-fidelity approximation of
Naiver-Stokes.

###  6.1 Illustrative example

A 1D steady state control problem is formulated that seeks an optimal source
to achieve a target state profile. The problem is formulated as

$\displaystyle\min_{z}\frac{1}{2}\int_{0}^{1}(\tilde{S}(z)(x)-T(x))^{2}dx+\frac{\beta}{2}\int_{0}^{1}z(x)^{2}dx$

where $T(x)=50-30(x-0.5)^{2}$ is the target, $z:[0,1]\to{\mathbb{R}}$ is a
source controller, $\beta=10$ is a regularization coefficient, and
$\tilde{S}(z)$ is the solution operator for the diffusion equation

$\displaystyle-\kappa u^{\prime\prime}=z\qquad$ $\displaystyle\text{on }(0,1)$
$\displaystyle\kappa u^{\prime}=hu\quad$ $\displaystyle\text{on }\\{0\\}$
$\displaystyle-\kappa u^{\prime}=hu\quad$ $\displaystyle\text{on }\\{1\\}$

with Robin boundary condition. The state and control spaces ${\mathcal{U}}$
and ${\mathcal{Z}}$ are $L^{2}(0,1)$. Our assumption is that the diffusion
model is missing physics, specifically it lacks advection. We consider the
high-fidelity model to be

$\displaystyle-\kappa u^{\prime\prime}+vu^{\prime}=z\qquad$
$\displaystyle\text{on }(0,1)$ $\displaystyle\kappa u^{\prime}=hu\quad$
$\displaystyle\text{on }\\{0\\}$ $\displaystyle-\kappa u^{\prime}=hu\quad$
$\displaystyle\text{on }\\{1\\}$

represented by solution operator $S(z)$. We take parameter values $\kappa=1$,
$v=0.5$, and $h=2$.

#### Optimization

After discretizing with linear finite elements on a uniform mesh of $200$
nodes, i.e. $m=n=200$, we solve both high and low-fidelity optimization
problems to illustrate the effect of model discrepancy, which is attributed to
the missing advection term. We observe a considerable difference in the
optimal sources, shown in Figure 4. Since advection moves the state from left
to right, including advection in the model results in a source skewed to the
left.

#### Prior model discrepancy

We define ${\bm{L}}=\bm{D}\bm{M}^{-1}\bm{D}$, where $\bm{D}$ is the
discretization of an elliptic operator $\gamma(-\epsilon\Delta+\mathcal{I})$,
equipped with zero Neumann boundary conditions, and $\bm{M}$ is the mass
matrix for the state discretization. The coefficients $\gamma$ and $\epsilon$
determine the magnitude and smoothness, respectively, of the prior discrepancy
realizations. The use of an elliptic operator ensures interpretability of the
hyper-parameters and computational convenience of efficient linear solves and
matrix factorizations for sampling [15, 39, 41].

Following Algorithm 1, we compute samples of the prior discrepancy to tune
$\gamma$ and $\epsilon$ so that the samples are physically consistent with the
magnitudes and smoothness inherent to this problem. The left panel of Figure 2
displays the samples corresponding to our choice of $\gamma=1$ and
$\epsilon=10^{-2}$. The observed data
${\bm{y}}_{\ell}={\bm{S}}({\bm{z}}_{\ell})-\tilde{{\bm{S}}}({\bm{z}}_{\ell})$
is shown in the right panel of Figure 2 to assess the discrepancy’s magnitude,
smoothness, and dependence on ${\bm{z}}$.

Figure 2: Left: 100 prior discrepancy samples (5 are colored for ease of
visualization) evaluated at $\tilde{{\bm{z}}}$; right: mean posterior
discrepancy evaluated at the two datapoints $\\{{\bm{z}}_{1},{\bm{z}}_{2}\\}$
compared with the discrepancy data $\\{{\bm{y}}_{1},{\bm{y}}_{2}\\}$.

The control length scale parameter $\zeta$ is determined by analyzing
variability of the prior samples with respect to ${\bm{z}}$, as described in
Section 5. Figure 3 shows $\\{\hat{{\bm{z}}}_{k}\\}_{k=0}^{11}$ in the
leftmost panel and three samples
$\\{{\bm{\delta}}(\hat{{\bm{z}}}_{k},\bm{C}^{-T}\bm{\omega}_{p}^{s})\\}_{k=0}^{11}$,
$s=1,2,3$, in the three rightmost panels. The scalar $\zeta$ may be tuned
according to the expected rate of variation of the discrepancy with respect to
${\bm{z}}$. In our case, $\zeta=2$.

Figure 3: Leftmost panel: sequence of controllers
$\hat{{\bm{z}}}_{k}=\tilde{{\bm{z}}}+\frac{k}{11}({\bm{z}}_{r}-\tilde{{\bm{z}}})$,
$k=0,1,\dots,11$; three rightmost panels: three prior discrepancy samples
evaluated at $\hat{{\bm{z}}}_{k}$, $k=0,1,\dots,11$.

#### Posterior model discrepancy

We evaluate the high-fidelity model twice to generate $N=2$ data pairs
$\\{{\bm{z}}_{\ell},{\bm{y}}_{\ell}\\}_{\ell=1}^{2}$ and take a noise variance
$\alpha=0.01$. The right panel of Figure 2 displays the mean posterior
discrepancy evaluated at ${\bm{z}}_{1}$ and ${\bm{z}}_{2}$. We observe that it
matches the data well. This is unsurprising given that
${\bm{S}}({\bm{z}})-\tilde{{\bm{S}}}({\bm{z}})$ is a linear function of
${\bm{z}}$ and $\alpha$ was taken small enough to encourage trust in the data.

#### Optimal solution update

Propagating the posterior model discrepancy mean through the post-optimality
sensitivity operator, we arrive at the updated optimal solution shown in
Figure 4. We observe that the updated optimal solution is much closer to the
high-fidelity solution ${\bm{z}}^{\star}$ than the low-fidelity solution
$\tilde{{\bm{z}}}$. Hence using the low-fidelity optimization problem along
with two high-fidelity forward solves, we are able to find a good
approximation of the high-fidelity optimal solution.

Figure 4: The low and high-fidelity solutions $\tilde{{\bm{z}}}$ and
${\bm{z}}^{\star}$ are shown by the blue and broken red lines, respectively.
The updated optimal solution is given by the black line.

###  6.2 Control of viscous fluid flow

In this subsection our approach is demonstrated on control of viscous fluid
flow. In particular, we consider the high-fidelity model to be governed by the
Naiver-Stokes equation and the low-fidelity model represented by the Stokes
equation. The goal is find the optimal distributed controller555We use the
notation $\bm{z}$ to denote the controller before discretization since it is
vector-valued. This should not be confused with the controller coordinates
$\bm{z}\in{\mathbb{R}}^{n}$ used throughout the article.
$\bm{z}=(z_{x},z_{y})$, defined on a subset of the domain, that minimizes the
vertical fluid flow. The low-fidelity optimization problem is

$\displaystyle\min_{z}\frac{1}{2}\int_{\chi}\tilde{v}_{y}(z)^{2}+\frac{\beta}{2}\int_{\Omega_{z}}(z_{x}^{2}+z_{y}^{2})$

where $\tilde{\bm{v}}(z)=(\tilde{v}_{x}(z),\tilde{v}_{y}(z))$ is the velocity
solution operator for the Stokes equation

$\displaystyle-\mu\nabla\bm{v}+\nabla p=\bm{g}+\bm{z}$ $\displaystyle\text{on
}\Omega$ $\displaystyle\nabla\cdot\bm{v}=0$ $\displaystyle\text{on }\Omega$

where $\Omega=(0,1)^{2}$, $\Omega_{z}=(0.2,.8)\times(0.1,0.4)$, and
$\chi=(0,1)\times(0,0.5)$. The controller is defined to be zero on
$\Omega\setminus\Omega_{z}$. The flow is driven by the effect of gravity
$\bm{g}=(0,9.81)$, viscosity $\mu=0.5$, and the inflow boundary conditions

$\displaystyle v_{x}(0,y)=6y(1-y)\qquad\text{and}\qquad v_{y}(x,1)=-2\sin(2\pi
x)^{2}.$

The state and control spaces ${\mathcal{U}}$ and ${\mathcal{Z}}$ are
$L^{2}(\Omega)$ and $L^{2}(\Omega_{z})$, respectively. We solve the
optimization problem with regularization coefficient $\beta=10^{-4}$. The
high-fidelity Naiver-Stokes model

$\displaystyle-\mu\nabla\bm{v}+(\bm{v}\cdot\nabla)\bm{v}+\nabla
p=\bm{g}+\bm{z}$ $\displaystyle\text{on }\Omega$
$\displaystyle\nabla\cdot\bm{v}=0$ $\displaystyle\text{on }\Omega$

includes the nonlinear convective term omitted in the Stokes equation.

#### Optimization

The three state variables $v_{x}$, $v_{y}$, and $p$ are discretized with
Taylor-Hood elements yielding $m=365622$ degrees of freedom. Figure 5 displays
the uncontrolled and optimally controlled solutions of the Stokes equation.
The inflow from the left and above, combined with the force of gravity,
produce a flow field with high downward velocity in the lower region of the
domain. The speed of this flow is minimized and the efficacy of the controller
(for the Stokes equation) is apparent from the center column of Figure 5.
However, the Naiver-Stokes solution evaluated at the optimal controller
exhibits faster flows.

Figure 5: Uncontrolled (top) and controlled (bottom) state solution for the
Stokes equation. The $x$-velocity $v_{x}$, $y$-velocity $v_{y}$, and pressure
$p$ are plotted from left to right. The uncontrolled $y$-velocity has a
magnitude around $6$ everywhere in $\chi$ whereas the controlled $y$-velocity
has a magnitude close to $0$.

#### Prior model discrepancy

Similar to the previous example, we define the prior covariance
${\bm{L}}^{-1}$ as the square of an inverse elliptic differential operator.
However, to mitigate the effect of high boundary variance in the prior we
follow [16] by imposing a Robin boundary condition which is optimized to
reduce the boundary effect. Because this condition is not sufficient to
completely eliminate the inflated variance near the boundary, a rescaling of
the covariance based on a point-wise variance estimate using $100$ samples is
also implemented.

The covariance is
${\bm{L}}^{-1}=\gamma^{2}(\bm{E}^{-\frac{1}{2}}\bm{D}\bm{M}^{-1}\bm{D}\bm{E}^{-\frac{1}{2}})^{-1}$,
where $\bm{D}$ is the discretization of the elliptic operator
$(-\epsilon\Delta+\mathcal{I})$, $\gamma=2$ and $\epsilon=0.2$, and $\bm{E}$
is a diagonal matrix that approximates the diagonal of
$(\bm{D}\bm{M}^{-1}\bm{D})^{-1}$, and $\bm{M}$ is the mass matrix for the
state discretization. We omit showing prior samples for conciseness, but
emphasize that the interpretably of the elliptic operator facilitates
specification of $\gamma$ and $\epsilon$ to ensure that the prior is
reflective of the physical characteristics of the discrepancy. We use $N=1$
high-fidelity model evaluation in this example. Close examination of (34)
shows that the controller length scale parameter $\zeta$ is immaterial if
$N=1$. The mean discrepancy does not vary with respect to ${\bm{z}}$ since
having only one high-fidelity evaluation is not sufficient data to fit the
variation.

To facilitate efficient computation, especially inversion of the shifted
linear systems $\alpha{\bm{L}}+\lambda_{i}{\bm{I}}$, we compute the truncated
generalized singular value decomposition of $\bm{D}^{-1}$ in the $\bm{M}$ and
$\bm{E}^{-1}$ inner product, which after algebraic manipulations gives the
eigenvalue decomposition of ${\bm{L}}^{-1}$. The subsequent linear solves
involving ${\bm{L}}$ may be efficiently approximated using the truncated
eigenvalue decomposition and the error in the approximation is controlled by
the rank, which we take to be $1000$, ensuring a relative truncation error of
$\mathcal{O}(10^{-3})$.

#### Posterior model discrepancy

We evaluate the Naiver-Stokes equation at the optimal Stokes controller
($N=1$) to generate the data pair $({\bm{z}}_{1},{\bm{y}}_{1})$ to fit the
model discrepancy. With a setting of $\alpha=0.01$, we solve the Bayesian
inverse problem to calibrate ${\bm{\delta}}({\bm{z}},{\bm{\theta}})$. Figure 6
displays the posterior mean of ${\bm{\delta}}$ evaluated at the controller
${\bm{z}}_{1}$ and shows a good data fit.

Figure 6: Observed data (top) and mean fit (bottom) for the state discrepancy.
The $x$-velocity $v_{x}$, $y$-velocity $v_{y}$, and pressure $p$ are plotted
from left to right.

#### Optimal solution update

Lastly, the mean from the posterior discrepancy is propagated through the
post-optimality sensitivity operator to produce the optimal solution update.
To illustrate the benefit of the controller update, Figure 7 displays the
optimal controller generated by solving the Stokes problem, the updated
optimal controller, and the optimal controller generated by solving the
control problem constrained by the Naiver-Stokes equation. Although our
assumption is that the latter controller is generally not available in
practice, we compute and display it here for comparison. Figure 7 demonstrates
that calibrating the discrepancy with $N=1$ high-fidelity forward solve
significantly improves the controller.

Figure 7: Controllers acting on the $x$-velocity component (top) and the
$y$-velocity component (bottom). The left column is the low-fidelity control
solution, the center column is the mean of the posterior optimal solution, and
the right column is the high-fidelity optimal solution computed by solving the
control problem constrained by the Naiver-Stokes model.

To further quantify the benefit of the controller update, Figure 8 displays
the vertical velocity $v_{y}$ of the Naiver-Stokes solution when evaluated at
the Stokes optimal controller (left), the updated optimal controller (center),
and the Naiver-Stokes optimal controller (right). The performance improvement
is significant given that it was achieved using only $N=1$ Naiver-Stokes solve
to update the Stokes solution. The extent of the improvement is quantified by
the value of the objective function, evaluated using the Naiver-Stokes model,
for each of the controllers (Table 1).

Figure 8: Vertical velocity state solutions corresponding to: Naiver-Stokes solve with low-fidelity control solution (left), Naiver-Stokes solve with updated control solution (center), and Naiver-Stokes with high-fidelity control solution (right). Controller | Objective Function Value
---|---
$\tilde{{\bm{z}}}$ | $.265$
$\tilde{{\bm{z}}}+\nabla_{{\bm{\theta}}}{\bm{F}}(\bm{0})\overline{{\bm{\theta}}}$ | $.113$
${\bm{z}}^{\star}$ | $.077$
Table 1: Value of the objective function $J(S(z),z)$, where $S(z)$ is the
Naiver-Stokes solution operator.

Having the viscosity $\mu=0.5$ in this example ensures that the Reynold number
is small. This justifies the Stokes approximation, though Figure 6
demonstrates that the difference in the Stokes and Naiver-Stokes state
solutions is non-trivial. We explored larger Reynolds numbers to stress test
the approach. The optimal solution update was less successful in these cases.
Nonetheless, we emphasize the significant optimal solution improvement in a
nontrivial $\mu=0.5$ case. Future work will explore approaches to quantify how
accurate the low-fidelity model must be, and improving on the optimal solution
updating strategy when it is not accurate enough.

##  7 Conclusion

It is common that high-fidelity models are too computationally intensive
and/or the code structures are not conducive to intrusive implementations of
efficient optimization algorithms, for instance, it may be difficult to
implement adjoints to enable efficient derivative computations. In such cases,
current practice is to construct low-fidelity models through either
simplifications of the equations from physics assumptions, or reduced order
modeling techniques to approximate the high-fidelity model at a lower
computational cost. The low-fidelity model is used to constrain an
optimization problem and its solution is assessed by evaluating the high-
fidelity model a small number of times. However, these high-fidelity solves
may be used for more than assessing the quality of the optimal solution. This
article provides a framework that systematically uses the high-fidelity data
to improve the optimal solution, without requiring any additional access to
the high-fidelity solver.

Our proposed approach is computationally scalable as a result of judiciously
manipulating linear algebra to enable closed form expressions for the
posterior. These efficiencies are possible by taking a wholistic perspective
that considers every aspect of the problem formulation (discrepancy
representation, Gaussian prior and likelihood, post-optimality sensitivity
operator, and the resulting matrix factorizations) to develop a mathematically
rigorous and computationally advantageous approach. Although the fundamental
assumption is a linear approximation of the model discrepancy to optimal
solution mapping, we demonstrate considerable improvements for nonlinear PDE-
constrained optimization problems with high-dimensional optimization
variables.

Since we constrain our analysis in a neighborhood of the low-fidelity
solution, there is no guarantee that the proposed approach will converge to
the solution of the high-fidelity optimization problem. If we have limited
high-fidelity model evaluations, as in our numerical results, then it is
unrealistic to expect convergence. This is a sacrifice made to afford an
approach that is pragmatic for large-scale applications where high-fidelity
simulation data may take hours or days to acquire and high-fidelity
derivatives are not accessible. Future work will explore an iterative approach
that uses high-fidelity data to update the solution and then takes a new
linear approximation about the updated solution. This may extend the
applicability of the approach for problems where high-fidelity data is limited
but may come in batches of evaluations.

There are many potential applications of our analysis including, but not
limited too, (1) using models in three spatial dimensions to improve optimal
solutions computed using models in two spatial dimensions, (2) using multi-
scale models to improve solutions computed using homogenized models, (3) using
coupled systems to improve optimal solutions computed using only a subset of
the system, and (4) using high-fidelity data from controlled experiments in
place of a high-fidelity model to improve design or control strategies. Other
mathematical questions regarding how and when to query the high-fidelity model
also remain and are the topic of ongoing research.

## Acknowledgements

This paper describes objective technical results and analysis. Any subjective
views or opinions that might be expressed in the paper do not necessarily
represent the views of the U.S. Department of Energy or the United States
Government. Sandia National Laboratories is 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-NA-0003525. SAND2022-14144 O.

## References

  * Agarwal and Biegler [2013] Agarwal, A., Biegler, L., 2013\.  A trust-region framework for constrained optimization using reduced order modeling.  Optimization and Engineering 14, 3–35.  doi:10.1007/s11081-011-9164-0.
  * Alexandrov et al. [1998] Alexandrov, N., Jr., J.D., Lewis, R., Torczon, V., 1998\.  A trust-region framework for managing the use of approximation models in optimization.  Structural Optimization 15, 16–23.  doi:10.1007/BF01197433.
  * Antil et al. [2018] Antil, H., Kouri, D.P., Lacasse, M., Ridzal, D. (Eds.), 2018\.  Frontiers in PDE-Constrained Optimization.  Springer.  doi:10.1007/978-1-4939-8636-1.
  * Arendt et al. [2012] Arendt, P.D., Apley, D.W., Chen, W., 2012.  Quantification of model uncertainty: Calibration, model discrepancy, and identifiability.  Journal of Mechanical Design 134, 1–12.  doi:10.1115/1.4007390.
  * Ascher and Haber [2001] Ascher, U.M., Haber, E., 2001\.  Grid refinement and scaling for distributed parameter estimation problems.  Inverse Problems 17, 571–590.  doi:10.1088/0266-5611/17/3/314.
  * Biegler et al. [2011] Biegler, L., Biros, G., Ghattas, O., Heinkenschloss, M., Keyes, D., Mallick, B., Marzouk, Y., Tenorio, L., van Bloemen Waanders, B., Willcox, K. (Eds.), 2011.  Large-Scale Inverse Problems and Quantification of Uncertainty.  John Wiley and Sons.  doi:10.1002/9780470685853.
  * Biegler et al. [2007] Biegler, L.T., Ghattas, O., Heinkenschloss, M., Keyes, D., van Bloemen Waanders, B. (Eds.), 2007.  Real-Time PDE-Constrained Optimization. volume 3.  SIAM Computational Science and Engineering.  doi:10.1137/1.9780898718935.
  * Biegler et al. [2003] Biegler, L.T., Ghattas, O., Heinkenschloss, M., van Bloemen Waanders, B. (Eds.), 2003.  Large-Scale PDE-Constrained Optimization. volume 30.  Springer-Verlag Lecture Notes in Computational Science and Engineering.  doi:10.1007/978-3-642-55508-4.
  * Biros and Ghattas [2005] Biros, G., Ghattas, O., 2005\.  Parallel Lagrange-Newton-Krylov-Schur methods for PDE-constrained optimization. Part I: the Krylov-Schur solver.  SIAM Journal on Scientific Computing 27, 687– 713.  doi:10.1137/S106482750241565X.
  * Bonnans and Shapiro [1998] Bonnans, J.F., Shapiro, A., 1998\.  Optimization problems with perturbations: A guided tour.  SIAM Review 40, 228–264.  doi:10.1137/S0036144596302644.
  * Borzi [2007] Borzi, A., 2007.  High-order discretization and multigrid solution of elliptic nonlinear constrained optimal control problems.  Journal of Computational and Applied Mathematics 200, 67–85.  doi:10.1016/j.cam.2005.12.023.
  * Brandes and Griesse [2007] Brandes, K., Griesse, R., 2007\.  Quantitative stability analysis of optimal solutions in PDE-constrained optimization.  Journal of Computational and Applied Mathematics 206, 908–926.  doi:10.1016/j.cam.2006.08.038.
  * Brynjarsdóttir and O’Hagan [2014] Brynjarsdóttir, J., O’Hagan, A., 2014\.  Learning about physical parameters: the importance of model discrepancy.  Inverse Problems 30, 1–24.  doi:10.1088/0266-5611/30/11/114007.
  * Bryson and Rumpfkeil [2018] Bryson, D.E., Rumpfkeil, M.P., 2018\.  Multifidelity Quasi-Newton method for design optimization.  AIAA Journal 56, 4074–4086.  doi:10.2514/1.J056840.
  * Bui-Thanh et al. [2013] Bui-Thanh, T., Ghattas, O., Martin, J., Stadler, G., 2013\.  A computational framework for infinite-dimensional Bayesian inverse problems. Part I: The linearized case, with applications to global seismic inversion.  SIAM Journal on Scientific Computing 35, A2494–A2523.  doi:10.1137/12089586X.
  * Daon and Stadler [2018] Daon, Y., Stadler, G., 2018\.  Mitigating the influence of the boundary on PDE-based covariance operators.  Inverse Problems and Imaging 12, 1083–1102.  doi:10.3934/ipi.2018045.
  * Fiacco and Ghaemi [1982] Fiacco, A.V., Ghaemi, A., 1982\.  Sensitivity analysis of a nonlinear structural design problem.  Computers and Operations Research 9, 29–55.  doi:10.1016/0305-0548(82)90005-3.
  * Gardner et al. [2021] Gardner, P., Rogers, T., Lord, C., Barthorpe, R., 2021\.  Learning model discrepancy: A Gaussian process and sampling-based approach.  Mechanical Systems and Signal Processing 152, 1–20.  doi:10.1016/j.ymssp.2020.107381.
  * Gorodetsky et al. [2020] Gorodetsky, A.A., Geraci, G., Eldred, M.S., Jakeman, J.D., 2020\.  A generalized approximate control variate framework for multifidelity uncertainty quantification.  Journal of Computational Physics 408, 1–29.  doi:10.1016/j.jcp.2020.109257.
  * Griesse [2004a] Griesse, R., 2004a.  Parametric sensitivity analysis in optimal control of a reaction-diffusion system – part II: practical methods and examples.  Optimization Methods and Software 19, 217–242.  doi:10.1080/10556780410001654250.
  * Griesse [2004b] Griesse, R., 2004b.  Parametric sensitivity analysis in optimal control of a reaction diffusion system. I. solution differentiability.  Numerical Functional Analysis and Optimization 25, 93–117.  doi:10.1081/NFA-120034120.
  * Haber and Ascher [2001] Haber, E., Ascher, U.M., 2001\.  Preconditioned all-at-once methods for large, sparse parameter estimation problems.  Inverse Problems 17, 1847–1864.  doi:10.1088/0266-5611/17/6/319.
  * Hart and van Bloemen Waanders [2023] Hart, J., van Bloemen Waanders, B., 2023\.  Enabling hyper-differential sensitivity analysis for ill-posed inverse problems.  Forthcoming in the SIAM Journal on Scientific Computing. arXiv:2106.11813 .
  * Hart et al. [2020] Hart, J., van Bloemen Waanders, B., Hertzog, R., 2020.  Hyper-differential sensitivity analysis of uncertain parameters in PDE-constrained optimization.  International Journal for Uncertainty Quantification 10, 225–248.  doi:10.1615/Int.J.UncertaintyQuantification.2020032480.
  * Hart et al. [2023] Hart, J., van Bloemen Waanders, B., Hood, L., Parish, J., 2023\.  Sensitivity driven experimental design to facilitate control of dynamical systems.  Journal of Optimization Theory and Applications doi:10.1007/s10957-023-02172-w.
  * Hazra and Schulz [2006] Hazra, S.B., Schulz, V., 2006\.  Simultaneous pseudo-timestepping for aerodynamic shape optimization problems with state constraints.  SIAM Journal on Scientific Computing 28, 1078–1099.  doi:10.1137/05062442X.
  * Higdon et al. [2008] Higdon, D., Gattiker, J., Williams, B., Rightley, M., 2008\.  Computer model calibration using high-dimensional output.  Journal of the American Statistical Association 103, 570–583.  doi:10.1198/016214507000000888.
  * Hintermuller and Vicente [2005] Hintermuller, M., Vicente, L.N., 2005\.  Space mapping for optimal control of partial differential equations.  SIAM J. Opt. 15, 1002–1025.  doi:10.1137/S105262340342907X.
  * Hinze et al. [2009] Hinze, M., Pinnau, R., Ulbrich, M., Ulbrich, S., 2009\.  Optimization with PDE Constraints.  Springer.  doi:10.1007/978-1-4020-8839-1.
  * Kennedy and O’Hagan [2001] Kennedy, M.C., O’Hagan, A., 2001\.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society 63, 425–464.  doi:10.1111/1467-9868.00294.
  * Laird et al. [2005] Laird, C.D., Biegler, L.T., van Bloemen Waanders, B., Bartlett, R.A., 2005\.  Time dependent contaminant source determination for municipal water networks using large scale optimization.  ASCE J. Water Res. Mgt. Plan. , 125–134.doi:10.1061/(ASCE)0733-9496(2005)131:2(125).
  * Ling et al. [2014] Ling, Y., Mullins, J., Mahadevan, S., 2014.  Selection of model discrepancy priors in Bayesian calibration.  Journal of Computational Physics 276, 665–680.  doi:10.1016/j.jcp.2014.08.005.
  * March and Willcox [2012] March, A., Willcox, K., 2012\.  Provably convergent multifidelity optimization algorithm not requiring high-fidelity derivatives.  AIAA Journal 50, 1079–1089.  doi:10.2514/1.J051125.
  * Maupin and Swiler [2020] Maupin, K.A., Swiler, L.P., 2020\.  Model discrepancy calibration across experimental settings.  Reliability Engineering & System Safety 200, 1–10.  doi:10.1016/j.ress.2020.106818.
  * Menhorn et al. [2020] Menhorn, F., Geraci, G., Seidl, D.T., Eldred, M.S., King, R., Bungartz, H.J., Marzouk, Y., 2020.  Higher moment multilevel estimators for optimization under uncertainty applied to wind plant design, in: AIAA SciTech Forum, pp. 1–16.  doi:10.2514/6.2020-1952.
  * Ng and Willcox [2014] Ng, L.W., Willcox, K.E., 2014\.  Multifidelity approaches for optimization under uncertainty.  International Journal for Numerical Methods in Engineering 100, 746–772.  doi:10.1002/nme.4761.
  * Oakley and O’Hagan [2002] Oakley, J., O’Hagan, A., 2002\.  Bayesian inference for the uncertainty distribution of computer model outputs.  Biometrika 89, 769–784.  URL: http://www.jstor.org/stable/4140536.
  * Peherstorfer et al. [2018] Peherstorfer, B., Willcox, K., Gunzburger, M., 2018.  Survey of multifidelity methods in uncertainty propagation, inference, and optimization.  SIAM Review 60, 550–591.  doi:10.1137/16M1082469.
  * Petra et al. [2014] Petra, N., Martin, J., Stadler, G., Ghattas, O., 2014\.  A computational framework for infinite-dimensional Bayesian inverse problems. Part II: Stochastic Newton MCMC with application to ice sheet flow inverse problems.  SIAM Journal on Scientific Computing 36, A1525–A1555.  doi:10.1137/130934805.
  * Saibaba et al. [2021] Saibaba, A.K., Hart, J., van Bloemen Waanders, B., 2021.  Randomized algorithms for generalized singular value decomposition with application to sensitivity analysis.  Numerical Linear Algebra with Applications 28, e2364.  doi:10.1002/nla.2364.
  * Stuart [2010] Stuart, A.M., 2010.  Inverse problems: A Bayesian perspective.  Acta Numerica 19, 451–559.  doi:10.1017/S0962492910000061.
  * Sunseri et al. [2022] Sunseri, I., Alexanderian, A., Hart, J., van Bloemen Waanders, B., 2022\.  Hyper-differential sensitivity analysis for nonlinear Bayesian inverse problems.  arXiv:2202.02219 .
  * Sunseri et al. [2020] Sunseri, I., Hart, J., van Bloemen Waanders, B., Alexanderian, A., 2020.  Hyper-differential sensitivity analysis for inverse problems constrained by partial differential equations.  Inverse Problems 36, 1–23.  doi:10.1088/1361-6420/abaf63.
  * Vogel [1999] Vogel, C.R., 1999.  Sparse matrix computations arising in distributed parameter identification.  SIAM Journal on Matrix Analysis and Applications , 1027–1037.doi:10.1137/S0895479897317703.
  * Vogel [2002] Vogel, C.R., 2002.  Computational Methods for Inverse Problems.  SIAM Frontiers in Applied Mathematics Series.  doi:10.1137/1.9780898717570.