Felix Black, Philipp Schulze, Benjamin Unger 11institutetext: Institute of
Mathematics, Technische Universität Berlin, Berlin, Germany
11email<EMAIL_ADDRESS>22institutetext: Stuttgart Center
for Simulation Science, Universität Stuttgart, Stuttgart, Germany
22email<EMAIL_ADDRESS>
# Decomposition of flow data via gradient-based transport optimization

Felix Black 11 Philipp Schulze 11 Benjamin Unger 22

###### Abstract

We study an optimization problem related to the approximation of given data by
a linear combination of transformed modes. In the simplest case, the
optimization problem reduces to a minimization problem well-studied in the
context of proper orthogonal decomposition. Allowing transformed modes in the
approximation renders this approach particularly useful to compress data with
transported quantities, which are prevalent in many flow applications. We
prove the existence of a solution to the infinite-dimensional optimization
problem. Towards a numerical implementation, we compute the gradient of the
cost functional and derive a suitable discretization in time and space. We
demonstrate the theoretical findings with three challenging numerical
examples.

###### keywords:

nonlinear model order reduction, transport-dominated phenomena, transformed
modes, gradient-based optimization

##  1 Introduction

Projection-based model order reduction (MOR) typically relies on the fact that
the solution manifold of a (parametrized) differential equation can be
approximately embedded in a low-dimensional linear subspace. The best subspace
of a given dimension, where best is understood as the minimal worst-case
approximation error, is characterized by the Kolmogorov $n$-widths [12]. In
practice, the minimizing subspace for the $n$-widths is difficult to compute.
Instead, one relies on the proper orthogonal decomposition (POD) [10], which
is typically combined with a greedy-search within the parameter domain, to get
an approximate solution. In more detail, for given parameters
$\mu_{\sigma}\in\mathbb{M}$ ($\sigma=1,\ldots,\ell$), associated data samples
$z(t,x;\mu_{\sigma})$ with time variable $t\in\mathbb{T}\vcentcolon=[0,T]$,
space variable $x\in\Omega\subseteq\mathbb{R}^{d}$, and desired dimension
$r\in\mathbb{N}$ of the low-dimensional subspace, POD determines orthonormal
basis functions of a low-dimensional subspace solving the minimization problem

$\left\\{\begin{aligned}
&\min\frac{1}{2}\sum_{\sigma=1}^{\ell}\int_{0}^{T}\bigg{\|}z(t,x;\mu_{\sigma})-\sum_{i=1}^{r}\alpha_{i}(t;\mu_{\sigma})\varphi_{i}(x)\bigg{\|}^{2}\,\mathrm{d}t\\\
&\mathrm{with}~{}\alpha_{i}(t;\mu_{\sigma})\vcentcolon=\langle
z(t,x;\mu_{\sigma}),\varphi_{i}(x)\rangle\quad\text{for
}i=1,\ldots,r,\,\sigma=1,\ldots,\ell,\\\
&\mathrm{s.t.}~{}\langle\varphi_{i},\varphi_{j}\rangle=\delta_{ij}\quad\text{for
}i,j=1,\ldots,r.\end{aligned}\right.$ (1)

If the $n$-widths, respectively the Hankel singular values for linear
dynamical systems [31], decay fast, then one can expect to construct an
effective reduced-order model (ROM) able to approximate the full dynamics with
a small approximation error. Although one can show exponential decay for a
large class of problems [16], it is well-known, see for instance [6, 9], that
the decay of the $n$-widths for flow problems is typically slow, thus
conspiring against MOR.

To remedy this issue prevalent in transport-dominated phenomena, several
strategies have been proposed in the literature. We refer to [21, 30, 25, 6,
20, 22, 29, 32, 13] to name just a few.

One promising approach, introduced in [24, 28, 14] and formalized in [3, 5],
is to replace the POD minimization problem with

$\min\frac{1}{2}\sum_{\sigma=1}^{\ell}\int_{0}^{T}\bigg{\|}z(t,x;\mu_{\sigma})-\sum_{i=1}^{r}\alpha_{i}(t;\mu_{\sigma})\varphi_{i}(x-p_{i}(t;\mu_{\sigma}))\bigg{\|}^{2}\,\mathrm{d}t,$
(2)

thus accounting explicitly for the transportation of quantities throughout the
spatial domain. Consequently, the linear subspace in the Kolmogorov $n$-widths
is replaced with a subspace able to adapt itself to the solution over time,
hence rendering this a nonlinear approach. Note that in contrast to the POD
minimization problem (1), we do not require the modes to be orthogonal to each
other. This is due to the fact that in the setting of (2), we would need to
require orthogonality of $\varphi_{i}(x-p_{i}(t;\mu_{\sigma}))$ and
$\varphi_{j}(x-p_{j}(t;\mu_{\sigma}))$ for all $i\neq j$, $t\in[0,T]$, and
$\sigma=1,\ldots,\ell$, which is in general not a reasonable assumption, cf.
[3, Ex. 4.4] for an illustrative example.

In the past years, there have been some attempts of solving discretized
versions of (2) or related minimization problems. In [24], the authors propose
a heuristic iterative method for computing a decomposition of a given snapshot
matrix by an approximation ansatz as in (2). The numerical experiments
indicate promising results, but it is not clear in which situations the
proposed method actually determines an optimal solution. Another heuristic has
been recently proposed in [4, sec. 5.2.1] and it is applied to snapshot data
of a wildland fire simulation. This method is based on a decomposition of the
snapshot matrix and involves a small number of singular value decompositions
without requiring an iterative procedure. The numerical results presented in
[4] demonstrate the effectiveness of this approach, but it is in general not
optimal in the sense of the minimization problem (2). In contrast, the method
introduced in [28] directly solves a fully discretized version of (2) by
determining optimal modes $\varphi$ and coefficients $\alpha$, but assumes the
paths $p$ to be given or determined in a pre-processing step. A similar
optimization problem has been proposed in [23] and aims at approximating the
snapshot matrix by a sum of matrices representing different reference frames
while achieving a fast singular value decay in each of the reference frames.
The corresponding cost function is shown to be an upper bound for a fully
discretized version of (2) and the cost functions coincide for the special
case that only one reference frame is considered, i.e., if all ansatz
functions in (2) are shifted by the same amount. Again, the paths are not
considered as part of the optimization problem, but instead determined in a
pre-processing step via peak or front tracking. On the contrary, the authors
in [18] focus on determining optimal paths, whereas the determination of
optimal ansatz functions and coefficients is not addressed. Moreover, to
simplify the optimization problem, the paths are sought within a low-
dimensional subspace consisting of predefined time-dependent library
functions. As in the case of the other mentioned works, also the authors in
[18] consider a fully discrete optimization problem.

We conclude that a gradient-based algorithm for the full optimization problem
(2) is currently not available. Besides, a rigorous proof showing that (2) has
a solution is missing in the literature. A notable exception is provided in
[3, Thm. 4.6], albeit under the assumption that the path variables $p_{i}(t)$
are known a-priori. In this paper we aim to close this gap. Our main
contributions are the following:

  1. 1.

We show in Theorem 3.1 the existence of a minimizing solution for the
optimization problem (6), which generalizes the minimization problem (2).
Afterward, we reformulate the constrained minimization problem (6) as
unconstrained problem (9) by adding appropriate penalty terms and conclude
from Theorem 3.1 that also the unconstrained problem has a solution, cf.
Corollary 3.5. In addition, Theorem 3.7 details that the solution of the
unconstrained problem converges to the solution of the constraint problem.

  2. 2.

We compute the gradient of the unconstrained problem in Theorem 3.10, which
enables the use of gradient-based methods to solve (2). In this context, a
remarkable finding is that the paths have to be sufficiently smooth (e.g. in
$H^{1}(0,T)$), since otherwise some directional derivatives of the cost
functional with respect to the paths may not exist, cf. Example 3.13.

  3. 3.

We discuss the discretization of the gradient in space and time in section 4
and explicitly compute the path-dependent inner products for the shift
operator with periodic boundary conditions in Example 4.1. Finally, the
effectiveness of gradient-based optimization is demonstrated for several
examples in section 5.

#### Notation

We denote the space of real $m\times n$ matrices by $\mathbb{R}^{m\times n}$
and the transpose of a matrix $A$ is written as $A^{\top}$. Furthermore, for a
vector with $n$ entries all equal to one we use the symbol $\mathbf{1}_{n}$.
Besides, for abbreviating diagonal and blockdiagonal matrices we use

$\operatorname{diag}(a_{1},\ldots,a_{n})\vcentcolon=\begin{bmatrix}a_{1}&&\\\
&\ddots&\\\
&&a_{n}\end{bmatrix},\quad\mathrm{blkdiag}(A_{1},\ldots,A_{n})\vcentcolon=\begin{bmatrix}A_{1}&&\\\
&\ddots&\\\ &&A_{n}\end{bmatrix},$

respectively, where $a_{1},\ldots,a_{n}$ are scalars and $A_{1},\ldots,A_{n}$
matrices of arbitrary size. For the Kronecker product of two matrices $A$ and
$B$ we write $A\otimes B$. The space of square-integrable functions mapping
from an interval $(a,b)$ to a Banach space $\mathscr{X}$ is denoted by
$L^{2}(a,b;\mathscr{X})$ and, similarly, the space of essentially bounded
measurable functions by $L^{\infty}(a,b;\mathscr{X})$. Furthermore, we use
$H^{1}(a,b;\mathscr{X})$ for the Sobolev subspace of functions in
$L^{2}(a,b;\mathscr{X})$ possessing also a weak derivative in
$L^{2}(a,b;\mathscr{X})$. The corresponding subspace consisting of
$H^{1}(a,b;\mathscr{X})$ functions whose values at the boundaries $a$ and $b$
coincide is denoted by $H^{1}_{\mathrm{per}}(a,b;\mathscr{X})$. Besides, for
the space of continuous functions from $[a,b]$ to $\mathscr{X}$ we use the
symbol $C([a,b];\mathscr{X})$. For the special case $\mathscr{X}=\mathbb{R}$,
we omit the last argument, i.e., we write, for instance, $L^{2}(a,b)$ instead
of $L^{2}(a,b;\mathbb{R})$.

##  2 Preliminaries and problem formulation

To formalize the optimization problem (2), we introduce the following spaces
and notation. Consider a real Hilbert space
$\left(\mathscr{X},\left<\cdot,\cdot\right>_{\mathscr{X}}\right)$ with induced
norm $\|\cdot\|_{\mathscr{X}}$, and let $\mathscr{Y}$ denote a dense subspace
of $\mathscr{X}$ that itself is a reflexive Banach space with norm
$\|\cdot\|_{\mathscr{Y}}$. Our standing assumption is that we are minimizing
the mean-squared distance to the data $z\in L^{2}(0,T;\mathscr{Y})$ in the
Bochner space $L^{2}(0,T;\mathscr{X})$ with the additional requirement that
the modes are elements of $\mathscr{Y}$.

To formalize the meaning of $\varphi_{i}(x-p_{i}(t))$ in (2), we follow the
notation in [3] and introduce a family of linear and bounded operators
$\mathcal{T}_{i}\colon\mathscr{P}_{i}\to\mathscr{B}(\mathscr{X})$ with real,
finite-dimensional vector space $\mathscr{P}_{i}$, for which we postulate the
following properties, taken from [3, Ass. 4.1].

###### Assumption

For every $i=1,\ldots,r$, every $\varphi_{i}\in\mathscr{Y}$, and every
$p_{i}\in\mathscr{P}_{i}$, the operator $\mathcal{T}_{i}(p_{i})$ is
$\mathscr{Y}$-invariant, i.e.,
$\mathcal{T}_{i}(p_{i})\mathscr{Y}\subseteq\mathscr{Y}$, and the mapping

$\mathscr{P}_{i}\to\mathscr{X},\qquad
p_{i}\mapsto\mathcal{T}_{i}(p_{i})\varphi_{i}$

is continuous.

A particular example for such a family of operators is given by the shift
operator with periodic boundary conditions, see for instance [3, Ex. 5.2]. For
further examples we refer to [2, 13].

For the ease of presentation, we restrict ourselves to the case
$\mathscr{P}_{i}=\mathbb{R}$, and emphasize that all results can be
generalized to $\mathscr{P}_{i}=\mathbb{R}^{m_{i}}$ for some
$m_{i}\in\mathbb{N}$. For

$\mathscr{Z}\vcentcolon=L^{2}(0,T;\mathbb{R}^{r})\times
H^{1}(0,T;\mathbb{R}^{r})\times\mathscr{Y}^{r}$ (3)

let us define the cost functional

$J\colon\mathscr{Z}\to\mathbb{R},\qquad\left(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}\right)\mapsto\frac{1}{2}\left\|z-\sum_{i=1}^{r}\alpha_{i}\mathcal{T}_{i}\left(p_{i}\right)\varphi_{i}\right\|_{L^{2}(0,T;\mathscr{X})}^{2}$
(4)

and for $C>0$ the space

$\mathcal{A}_{C}\vcentcolon=\left\\{(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\in\mathscr{Z}\,\left|\,\max\left\\{\|\varphi_{i}\|_{\mathscr{Y}},\|\alpha_{i}\|_{L^{2}\left(0,T\right)},\|p_{i}\|_{H^{1}\left(0,T\right)}\right\\}\leq
C\right.\right\\},$ (5)

where we use the notation $\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{r})$
to denote the coefficients of $\boldsymbol{\alpha}$ and analogously for
$\boldsymbol{p}$ and $\boldsymbol{\varphi}$. To ensure that the norm in (4) is
defined, we invoke the following assumption, which is for instance satisfied
if the family of operators $\mathcal{T}_{i}(\cdot)$ is uniformly bounded, cf.
[3, Lem 4.2].

###### Assumption

For every $i=1,\ldots,r$, $\alpha_{i}\in L^{2}(0,T)$, $p_{i}\in H^{1}(0,T)$,
and every $\varphi_{i}\in\mathscr{Y}$, we assume

$\alpha_{i}\mathcal{T}_{i}(p_{i})\varphi_{i}\in L^{2}(0,T;\mathscr{X}).$

With these preparations, the constrained minimization problem that we are
interested in takes the form

$\displaystyle\min_{\left(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}\right)}J\left(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}\right),\quad\text{s.\,t.
}\left(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}\right)\in\mathcal{A}_{C}.$
(6)

Before we proceed with our main results, let us make the following remarks:

  * •

To simplify the notation, we have implicitly set $\ell=1$ in (2), thus
assuming a single data sample. We emphasize that it is straightforward to
generalize all results to $\ell>1$.

  * •

The restriction of the optimization parameters to the admissible set stated in
(5) is used for proving the existence of a minimizing solution, cf. Theorem
3.1, and helps to circumvent the problem that the cost functional in (4) is
not coercive. For instance, since the transformed modes
$\mathcal{T}_{i}\left(p_{i}\right)\varphi_{i}$ may be linearly dependent, the
corresponding coefficients may become unbounded, even though the value of the
cost functional remains constant. The restriction to an admissible set of
bounded functions as in (5) allows to avoid such difficulties. Moreover, also
from an application point of view, we note that it is usually reasonable to
restrict to bounded optimization parameters, since when decomposing flow data,
for example, we are usually not interested in unbounded coefficients,
discontinuous paths, or modes which are less regular than the given flow data.

##  3 Main results

As first main result, we discuss the existence of a solution for the
optimization problem (6), thus generalizing [3, Thm. 4.6] to include the path
variables.

######  Theorem 3.1.

Assume that the reflexive Banach space $\mathscr{Y}$ is compactly embedded
into $\mathscr{X}$, and let $z\in L^{2}(0,T;\mathscr{Y})$. Furthermore, let
the family of transformation operators satisfiy Assumptions 2 and 2. Then the
constraint minimization problem (6) has a solution for every $C>0$.

######  Proof 3.2.

The proof follows along the lines of the proof of [3, Thm. 4.6], with slight
modifications to account for the optimization with respect to the path
variables. Let $C>0$. We first observe that the optimization problem possesses
a finite infimum $J^{\star}\geq 0$. This follows directly from $J\geq 0$ and
$(0,0,0)\in\mathcal{A}_{C}$. We may thus choose a sequence
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})_{k\in\mathbb{N}}\in\mathcal{A}_{C}$
satisfying

$\lim_{k\to\infty}J(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})=J^{\star}.$

Additionally, we have

$\left\|(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})\right\|_{L^{2}(0,T;\mathbb{R}^{r})\times
H^{1}(0,T;\mathbb{R}^{r})\times\mathscr{Y}^{r}}^{2}\\\
=\|\boldsymbol{\alpha}^{k}\|_{L^{2}(0,T;\mathbb{R}^{r})}^{2}+\|\boldsymbol{p}^{k}\|_{H^{1}(0,T;\mathbb{R}^{r})}^{2}+\|\boldsymbol{\varphi}^{k}\|_{\mathscr{Y}^{r}}^{2}\leq
3rC^{2}$

for all $k\in\mathbb{N}$, such that the Eberlein-S̆muljan theorem [33, Thm.
21.D] ensures the existence of a weakly convergent subsequence
$(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})_{n\in\mathbb{N}}\subseteq\mathcal{A}_{C}$
with weak limit
$(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})\in\mathcal{A}_{C}$,
cf. [33, Prop. 21.23 (c)]. Due to the compact embeddings
$\mathscr{Y}\hookrightarrow\mathscr{X}$ and
$H^{1}(0,T;\mathbb{R}^{r})\hookrightarrow L^{2}(0,T;\mathbb{R}^{r})$, we
conclude that $(\boldsymbol{\varphi}^{k_{n}})_{n\in\mathbb{N}}$ and
$(\boldsymbol{p}^{k_{n}})_{n\in\mathbb{N}}$ converge strongly in $\mathscr{X}$
and $L^{2}(0,T;\mathbb{R}^{r})$ to $\boldsymbol{\varphi}^{\star}$ and
$\boldsymbol{p}^{\star}$, respectively, cf. [33, Prop. 21.35]. Using [27, Thm.
3.12], we conclude the existence of a subsequence, for which we use the same
indexing, such that $(\boldsymbol{p}^{k_{n}})_{n\in\mathbb{N}}$ converges
pointwise to $\boldsymbol{p}^{\star}$ for almost all $t\in(0,T)$.

For the next part of the proof, we introduce the mapping

$\displaystyle\beta\colon\mathscr{Z}\to
L^{2}(0,T;\mathscr{X}),\qquad(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\mapsto\sum_{i=1}^{r}\alpha_{i}\mathcal{T}_{i}(p_{i})\varphi_{i}$
(7)

with $\mathscr{Z}$ as defined in (3) and notice

$J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})=\frac{1}{2}\|z-\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\|_{L^{2}(0,T;\mathscr{X})}^{2}.$

If
$\beta(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})$
converges weakly to
$\beta(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})$,
then the weak sequential lower semicontinuity of the norm, see for instance
[33, Prop. 21.23 (c)], implies that
$(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})$
is a minimizer of $J$. It thus remains to show that
$\beta(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})$
converges weakly to
$\beta(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})$.

To this end, we observe that

$\displaystyle\|\mathcal{T}_{i}(p_{i}^{k_{n}}(t))\varphi_{i}^{k_{n}}-\mathcal{T}_{i}(p_{i}^{\star}(t))\varphi^{\star}\|_{\mathscr{X}}$
$\displaystyle\leq\|\mathcal{T}_{i}(p_{i}^{k_{n}}(t))\varphi_{i}^{k_{n}}-\mathcal{T}_{i}(p_{i}^{\star}(t))\varphi_{i}^{k_{n}}\|_{\mathscr{X}}+\|\mathcal{T}_{i}(p_{i}^{\star}(t))\varphi_{i}^{k_{n}}-\mathcal{T}_{i}(p_{i}^{\star}(t))\varphi_{i}^{\star}\|_{\mathscr{X}}$

together with Assumption 2 and the strong convergence of
$(\varphi_{i}^{k_{n}})_{n\in\mathbb{N}}$ in $\mathscr{X}$ implies

$\displaystyle\|\mathcal{T}_{i}(p_{i}^{k_{n}}(t))\varphi_{i}^{k_{n}}-\mathcal{T}_{i}(p_{i}^{\star}(t))\varphi^{\star}\|_{\mathscr{X}}\to
0\qquad\text{for }n\to\infty$

for $i=1,\ldots,r$ and almost all $t\in(0,T)$. Let $f\in
L^{2}(0,T;\mathscr{X})$. Then clearly

$\langle
f(t),\mathcal{T}_{i}(p_{i}^{k_{n}}(t))\varphi_{i}^{k_{n}}\rangle_{\mathscr{X}}\to\langle
f(t),\mathcal{T}_{i}(p_{i}^{\star}(t))\varphi_{i}^{\star}\rangle_{\mathscr{X}}\qquad\text{for
}n\to\infty$

for $i=1,\ldots,r$ and almost all $t\in(0,T)$ such that [33, Prop 21.23 (j)]
implies

$\displaystyle\sum_{i=1}^{r}\left\langle\alpha_{i}^{k_{n}},\langle
f,\mathcal{T}_{i}(p_{i}^{k_{n}})\varphi_{i}^{k_{n}}\rangle_{\mathscr{X}}\right\rangle_{L^{2}(0,T)}\to\sum_{i=1}^{r}\left\langle\alpha_{i}^{\star},\langle
f,\mathcal{T}_{i}(p_{i}^{\star})\varphi_{i}^{\star}\rangle_{\mathscr{X}}\right\rangle_{L^{2}(0,T)}$

for $n\to\infty$ and thus

$\displaystyle\beta(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})\rightharpoonup\beta(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})\qquad\text{for
}n\to\infty,$

which completes the proof. ∎

For numerical methods, it may be easier to work with unconstrained
optimization problems. To this end, we use a penalty method, see for instance
[15, Cha. 13.1], i.e., we add the constraint equation with a penalty parameter
to the cost functional. In more detail, we assume for $C>0$ a penalty
functional

$\Lambda_{C}\colon\mathscr{Z}\to\mathbb{R}$ (8)

with the following properties available.

###### Assumption 3.3

The penalty functional (8) is continuous, weakly sequentially lower
semicontinuous, non-negative and has the following properties:

  * •

We have
$\Lambda_{C}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})=0$ if,
and only if,
$(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\in\mathcal{A}_{C}$.

  * •

For any sequence
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})$ with

$\max\\{\|\boldsymbol{\alpha}^{k}\|_{L^{2}(0,T;\mathbb{R}^{r})},\|\boldsymbol{p}^{k}\|_{H^{1}(0,T;\mathbb{R}^{r})},\|\boldsymbol{\varphi}^{k}\|_{\mathscr{Y}^{r}}\\}\to\infty\qquad\text{for
}k\to\infty,$

we have
$\Lambda_{C}(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})\to\infty$
for $k\to\infty$.

######  Example 3.4.

The penalty functional

$\displaystyle\Lambda_{C}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\vcentcolon=\sum_{i=1}^{r}\max\\{0,\max\\{\|\alpha_{i}\|_{L^{2}(0,T)},\left\|p_{i}\right\|_{H^{1}(0,T)},\left\|\varphi_{i}\right\|_{\mathscr{Y}}\\}-C\\}$

satisfies Assumption 3.3.

The penalized cost functional is then given as

$\widetilde{J}_{C}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi},\lambda)\vcentcolon=J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})+\lambda\Lambda_{C}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$

with penalty coefficient $\lambda>0$. The associated (unconstrained)
minimization problem is thus given by

$\displaystyle\min_{(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\in\mathscr{Z}}\widetilde{J}_{C}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi},\lambda)$
(9)

with $\mathscr{Z}$ as defined in (3) and given $\lambda>0$.

######  Corollary 3.5.

Let the assumptions of Theorem 3.1 and Assumption 3.3 be satisfied. Then for
any $C>0$ and any $\lambda>0$ the optimization problem (9) has a solution.

######  Proof 3.6.

Similarly as in the proof of Theorem 3.1, we conclude the existence of a
finite infimum, such that we can choose a minimizing sequence
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})\in\mathscr{Z}$.
Due to Assumption 3.3, we deduce that
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})_{k\in\mathbb{N}}$
is bounded in $\mathscr{Z}$, i.e., there exists some $\widetilde{C}>0$ such
that
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})\in\mathcal{A}_{\widetilde{C}}$
for all $k\in\mathbb{N}$. The remaining proof thus follows along the lines of
the proof of Theorem 3.1. ∎

######  Theorem 3.7.

Let $(\lambda^{k})_{k\in\mathbb{N}}\subseteq\mathbb{R}$ denote a non-
decreasing sequence of positive numbers with
$\lim_{k\to\infty}\lambda^{k}=\infty$, and let $C>0$. For $k\in\mathbb{N}$,
let
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})\in\mathscr{Z}$
denote a solution of (9) with penalty parameter $\lambda^{k}$. If the
assumptions of Corollary 3.5 are satisfied, then any limit point of
$(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})_{k\in\mathbb{N}}$
is a solution of (6).

######  Proof 3.8.

The proof follows along the lines of the proof of the main theorem in [15,
Cha. 13.1]. Let
$(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})\in\mathcal{A}_{C}$
denote a minimizer of (6) with minimum $J^{\star}$. Then for every
$k\in\mathbb{N}$ we have

$\displaystyle\widetilde{J}_{C}(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k},\lambda^{k})$
$\displaystyle\leq\widetilde{J}_{C}(\boldsymbol{\alpha}^{k+1},\boldsymbol{p}^{k+1},\boldsymbol{\varphi}^{k+1},\lambda^{k})$
$\displaystyle\leq\widetilde{J}_{C}(\boldsymbol{\alpha}^{k+1},\boldsymbol{p}^{k+1},\boldsymbol{\varphi}^{k+1},\lambda^{k+1})$

and

$\displaystyle
J(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k})\leq\widetilde{J}_{C}(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k},\lambda^{k})\leq\widetilde{J}_{C}(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star},\lambda^{k})={J}(\boldsymbol{\alpha}^{\star},\boldsymbol{p}^{\star},\boldsymbol{\varphi}^{\star})=J^{\star}.$

Thus
$(\widetilde{J}_{C}(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k},\lambda^{k}))_{k\in\mathbb{N}}$
is a monotone sequence bounded above by $J^{\star}$. We thus set

$\widetilde{J}_{C}^{\star}\vcentcolon=\lim_{k\to\infty}\widetilde{J}_{C}(\boldsymbol{\alpha}^{k},\boldsymbol{p}^{k},\boldsymbol{\varphi}^{k},\lambda^{k})\leq
J^{\star}.$ (10)

Let
$(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})$
denote a convergent subsequence with limit
$(\boldsymbol{\alpha}^{\dagger},\boldsymbol{p}^{\dagger},\boldsymbol{\varphi}^{\dagger})$
and set

$J^{\dagger}\vcentcolon=\lim_{n\to\infty}J(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})=J(\boldsymbol{\alpha}^{\dagger},\boldsymbol{p}^{\dagger},\boldsymbol{\varphi}^{\dagger}),$
(11)

using the continuity of $J$. Subtracting (10) from (11) yields

$\lim_{n\to\infty}\lambda^{k_{n}}\Lambda_{C}(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})=\widetilde{J}_{C}^{\star}-J^{\dagger}.$

Assumption 3.3 and $\lambda^{k_{n}}\to\infty$ for $n\to\infty$ together with
the continuity of $\Lambda_{C}$ thus implies

$\Lambda_{C}(\boldsymbol{\alpha}^{\dagger},\boldsymbol{p}^{\dagger},\boldsymbol{\varphi}^{\dagger})=\lim_{n\to\infty}\Lambda_{C}(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})=0,$

showing
$(\boldsymbol{\alpha}^{\dagger},\boldsymbol{p}^{\dagger},\boldsymbol{\varphi}^{\dagger})\in\mathcal{A}_{C}$.
We conclude

$\displaystyle
J^{\dagger}=\lim_{n\to\infty}J(\boldsymbol{\alpha}^{k_{n}},\boldsymbol{p}^{k_{n}},\boldsymbol{\varphi}^{k_{n}})\leq
J^{\star},$

which completes the proof. ∎

Although (9) is an unconstrained optimization problem, we still have to choose
a suitable constant $C>0$ for the admissible set. Let us emphasize that the
proofs of Theorem 3.1 and Corollary 3.5 heavily depend on the fact that we
have bounded sequences, which is the main reason for the constant $C>0$ in the
admissible set (5). However, we observed faster convergence in our numerical
experiments when considering the unconstrained minimization problem without
penalization. For this reason and the sake of a concise presentation, we
consider in the following only the unconstrained optimization problem (9) with
penalty parameter $\lambda=0$. Nevertheless, we emphasize that adding the
derivatives of the penalty terms to the gradient formulas is straightforward
as long as the partial Fréchet derivatives of $\Lambda_{C}$ are available.

To solve the optimization problem (9) with penalty parameter $\lambda=0$
numerically, we employ a gradient-based algorithm and thus have to compute the
gradient of the objective function (4). It is easy to see that the directional
derivatives of $J$ with respect to the coefficient function
$\boldsymbol{\alpha}\in L^{2}(0,T;\mathbb{R}^{r})$ and the modes
$\boldsymbol{\varphi}\in\mathscr{Y}^{r}$ in directions $\boldsymbol{d}\in
L^{2}(0,T;\mathbb{R}^{r})$ and $\boldsymbol{h}\in\mathscr{Y}^{r}$,
respectively, are given by

$\displaystyle\partial_{\boldsymbol{\alpha},\boldsymbol{d}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle=\sum_{i=1}^{r}\left\langle\sum_{j=1}^{r}\alpha_{j}\mathcal{T}_{j}(p_{j})\varphi_{j}-z,d_{i}\mathcal{T}_{i}(p_{i})\varphi_{i}\right\rangle_{\\!\\!\\!L^{2}(0,T;\mathscr{X})},$
(12a)
$\displaystyle\partial_{\boldsymbol{\varphi},\boldsymbol{h}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle=\sum_{i=1}^{r}\left\langle\sum_{j=1}^{r}\alpha_{j}\mathcal{T}_{j}(p_{j})\varphi_{j}-z,\alpha_{i}\mathcal{T}_{i}(p_{i})h_{i}\right\rangle_{\\!\\!\\!L^{2}(0,T;\mathscr{X})}.$
(12b) The situation is slightly different for the partial derivative with
respect to the path variable. First of all, we have to ensure that the
transformed modes are differentiable (with respect to the path variable),
i.e., we have to evoke the following assumption.

###### Assumption 3.9

For every $\varphi_{i}\in\mathscr{Y}$ and every $i=1,\ldots,r$, the mapping

$\displaystyle\mathbb{R}\to\mathscr{X},\qquad
p_{i}\mapsto\mathcal{T}_{i}(p_{i})\varphi_{i},$

is continuously differentiable with derivatives in $\mathscr{X}$. For
$p_{i}\in\mathbb{R}$ we denote the derivative by $\frac{\partial}{\partial
p_{i}}\mathcal{T}_{i}\left(p_{i}\right)\varphi_{i}\in\mathscr{X}$ and assume
$\alpha_{i}\frac{\partial}{\partial p_{i}}\mathcal{T}_{i}(p_{i})\varphi_{i}\in
L^{2}(0,T;\mathscr{X})$ for all $\alpha_{i}\in L^{2}(0,T)$ and all $p_{i}\in
H^{1}(0,T)$.

In this case, the directional derivative in direction $\boldsymbol{g}\in
H^{1}(0,T;\mathbb{R}^{r})$ is given as

$\partial_{\boldsymbol{p},\boldsymbol{g}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})=\sum_{i=1}^{r}\left\langle\sum_{j=1}^{r}\alpha_{j}\mathcal{T}_{j}(p_{j})\varphi_{j}-z,\alpha_{i}\left[\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\varphi_{i}\right]g_{i}\right\rangle_{\\!\\!\\!L^{2}(0,T;\mathscr{X})}.$
(12c)

Note that the Sobolev embedding theorems, see for instance [33, Thm.
21.A.(d)], imply $g_{i}\in C([0,T])\subseteq L^{\infty}(0,T)$, such that (12c)
is defined.

######  Theorem 3.10.

Let the transformation operators satisfy Assumptions 2, 2, and 3.9. Let
$(\boldsymbol{\alpha},\boldsymbol{\varphi},\boldsymbol{p})\in\mathscr{Z}$ and
assume

$\displaystyle\mathcal{T}_{i}(p_{i})\varphi_{i}$ $\displaystyle\in
L^{\infty}(0,T;\mathscr{X}),$ (13a)
$\displaystyle\alpha_{i}\|\mathcal{T}_{i}(p_{i})\|$ $\displaystyle\in
L^{2}(0,T)$ (13b)

for $i=1,\ldots,r$, then the partial Fréchet derivatives of the cost
functional $J$ (defined in (4)) with respect to the coefficients, paths, and
modes at
$(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\in\mathscr{Z}$ are
given by

$\displaystyle\partial_{\boldsymbol{\alpha}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})(\boldsymbol{d})$
$\displaystyle\vcentcolon=\partial_{\boldsymbol{\alpha},\boldsymbol{d}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),$
$\displaystyle\forall\boldsymbol{d}\in L^{2}(0,T;\mathbb{R}^{r}),$ (14a)
$\displaystyle\partial_{\boldsymbol{p}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})(\boldsymbol{g})$
$\displaystyle\vcentcolon=\partial_{\boldsymbol{p},\boldsymbol{g}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),$
$\displaystyle\forall\boldsymbol{g}\in H^{1}(0,T;\mathbb{R}^{r}),$ (14b)
$\displaystyle\partial_{\boldsymbol{\varphi}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})(\boldsymbol{h})$
$\displaystyle\vcentcolon=\partial_{\boldsymbol{\varphi},\boldsymbol{h}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),$
$\displaystyle\forall\boldsymbol{h}\in\mathscr{Y}^{r},$ (14c)

with directional derivatives as defined in (12).

######  Proof 3.11.

It suffices to show that $J$ is partially Fréchet differentiable with respect
to the coefficients, paths, and modes. Let
$(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),(\boldsymbol{d},\boldsymbol{g},\boldsymbol{h})\in\mathscr{Z}$.
Using (13a) we obtain

$\displaystyle J(\boldsymbol{\alpha}$
$\displaystyle+\boldsymbol{d},\boldsymbol{p},\boldsymbol{\varphi})-J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})-\partial_{\boldsymbol{\alpha},\boldsymbol{d}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})=\frac{1}{2}\left\|\sum_{i=1}^{r}d_{i}\mathcal{T}_{i}(p_{i})\varphi_{i}\right\|_{L^{2}(0,T;\mathscr{X})}^{2}$
$\displaystyle\leq\frac{1}{2}\max_{i=1,\ldots,r}\|\mathcal{T}_{i}(p_{i})\varphi_{i}\|_{L^{\infty}(0,T;\mathscr{X})}^{2}\left(\sum_{i=1}^{r}\|d_{i}\|_{L^{2}(0,T)}\right)^{2}$
$\displaystyle\leq\frac{r^{2}}{2}\max_{i=1,\ldots,r}\|\mathcal{T}_{i}(p_{i})\varphi_{i}\|_{L^{\infty}(0,T;\mathscr{X})}^{2}\|\boldsymbol{d}\|_{L^{2}(0,T;\mathbb{R}^{r})}^{2}$

and thus

$\displaystyle\lim_{\|\boldsymbol{d}\|_{L^{2}(0,T;\mathbb{R}^{r})}\to
0}\frac{\left|J(\boldsymbol{\alpha}+\boldsymbol{d},\boldsymbol{p},\boldsymbol{\varphi})-J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})-\partial_{\boldsymbol{\alpha},\boldsymbol{d}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\right|}{\|\boldsymbol{d}\|_{L^{2}(0,T;\mathbb{R}^{r})}}=0.$

We conclude that $J$ is Fréchet differentiable with respect to the
coefficients with Fréchet derivative as in (14a). For the partial derivative
with respect to the modes we obtain

$\displaystyle J(\boldsymbol{\alpha}$
$\displaystyle,\boldsymbol{p},\boldsymbol{\varphi}+\boldsymbol{h})-J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})-\partial_{\boldsymbol{\varphi},\boldsymbol{h}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})=\frac{1}{2}\left\|\sum_{i=1}^{r}\alpha_{i}\mathcal{T}_{i}(p_{i})h_{i}\right\|_{L^{2}(0,T;\mathscr{X})}^{2}$
$\displaystyle\leq\frac{1}{2}\int_{0}^{T}\left(\sum_{i=1}^{r}|\alpha_{i}(t)|\|\mathcal{T}_{i}(p_{i}(t))\|\|h_{i}\|_{\mathscr{X}}\right)^{2}\,\mathrm{d}t$
$\displaystyle\leq\frac{\|\boldsymbol{h}\|_{\mathscr{X}^{r}}^{2}}{2}\int_{0}^{T}\left(\sum_{i=1}^{r}|\alpha_{i}(t)|\|\mathcal{T}_{i}(p_{i}(t))\|\right)^{2}\,\mathrm{d}t.$

Using (13b), we observe that the integral is finite. Similarly as before, we
thus conclude that $J$ is Fréchet differentiable with respect to the modes
with Fréchet derivative as in (14c). We conclude our proof for the partial
derivative with respect to the path variable. Note that the Sobolev embedding
theorem [1, Thm 4.12, Part I, Case A] implies that the Sobolev space
$H^{1}(0,T)$ is continuously embedded into the space $L^{\infty}(0,T)$, i.e.,
there exists a constant $\gamma>0$ independent of $g_{i}$, such that
$\|g_{i}\|_{L^{\infty}(0,T)}\leq\gamma\|g_{i}\|_{H^{1}(0,T)}$. We define

$f_{i}(p_{i},\varphi_{i},g_{i})\vcentcolon=\mathcal{T}_{i}(p_{i}+g_{i})\varphi_{i}-\mathcal{T}_{i}(p_{i})\varphi_{i}-\left[\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\right]g_{i}.$

For $g_{i}\equiv 0$ we have $f_{i}(p_{i},\varphi_{i},g_{i})=0$ for almost all
$t\in(0,T)$. For $\|g_{i}\|_{H^{1}(0,T)}\neq 0$, let us define
$\widehat{\mathbb{T}}_{i}\vcentcolon=\\{t\in(0,T)\mid g_{i}(t)\neq 0\\}$. Then

$\frac{\left\langle
z,\alpha_{i}f_{i}(p_{i},\varphi_{i},g_{i})\right\rangle_{L^{2}(0,T;\mathscr{X})}}{\|g_{i}\|_{H^{1}(0,T)}}\leq\gamma\int_{\widehat{\mathbb{T}}_{i}}\alpha_{i}(t)\left\langle
z(t),\frac{f_{i}(p_{i}(t),\varphi_{i},g_{i}(t))}{|g_{i}(t)|}\right\rangle_{\mathscr{X}}\,\mathrm{d}t.$

From Assumption 3.9 we conclude

$\lim_{\|g_{i}\|_{H^{1}(0,T)}\to 0}\frac{\left\langle
z,\alpha_{i}f_{i}(p_{i},\varphi_{i},g_{i})\right\rangle_{L^{2}(0,T;\mathscr{X})}}{\|g_{i}\|_{H^{1}(0,T)}}=0,$

and thus

$\lim_{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}\to
0}\frac{\sum_{i=1}^{r}\left\langle
z,\alpha_{i}f_{i}(p_{i},\varphi_{i},g_{i})\right\rangle_{L^{2}(0,T;\mathscr{X})}}{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}}=0.$

Furthermore, using $\beta$ as defined in (7), we obtain

$\displaystyle\tfrac{1}{2}\|\beta(\boldsymbol{\alpha},\boldsymbol{p}$
$\displaystyle+\boldsymbol{g},\boldsymbol{\varphi})\|_{L^{2}(0,T;\mathscr{X})}^{2}-\tfrac{1}{2}\left\|\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\right\|_{L^{2}(0,T;\mathscr{X})}^{2}$
$\displaystyle\phantom{=}\qquad-\sum_{j=1}^{r}\left\langle\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),\alpha_{j}\left[\tfrac{\partial}{\partial
p_{j}}\mathcal{T}_{j}(p_{j})\varphi_{j}\right]g_{j}\right\rangle_{L^{2}(0,T;\mathscr{X})}$
$\displaystyle=\tfrac{1}{2}\left\|\beta(\boldsymbol{\alpha},\boldsymbol{p}+\boldsymbol{g},\boldsymbol{\varphi})-\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})+\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\right\|_{L^{2}(0,T;\mathscr{X})}^{2}$
$\displaystyle\phantom{=}\qquad-\tfrac{1}{2}\left\|\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\right\|_{L^{2}(0,T;\mathscr{X})}^{2}$
$\displaystyle\phantom{=}\qquad-\sum_{j=1}^{r}\left\langle\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),\alpha_{j}\left[\tfrac{\partial}{\partial
p_{j}}\mathcal{T}_{j}(p_{j})\varphi_{j}\right]g_{j}\right\rangle_{L^{2}(0,T;\mathscr{X})}$
$\displaystyle=\tfrac{1}{2}\|\beta(\boldsymbol{\alpha},\boldsymbol{p}+\boldsymbol{g},\boldsymbol{\varphi})-\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\|_{L^{2}(0,T;\mathscr{X})}^{2}$
$\displaystyle\phantom{=}\qquad+\sum_{j=1}^{r}\big{\langle}\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),\alpha_{j}f_{j}(p_{j},\varphi_{j},g_{j})\big{\rangle}_{L^{2}(0,T;\mathscr{X})}.$

Similarly as before, we obtain

$\displaystyle\lim_{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}\to
0}\frac{1}{2}\frac{\|\beta(\boldsymbol{\alpha},\boldsymbol{p}+\boldsymbol{g},\boldsymbol{\varphi})-\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\|_{L^{2}(0,T;\mathscr{X})}^{2}}{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}}=0,$
$\displaystyle\lim_{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}\to
0}\frac{\big{\langle}\beta(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),\sum_{j=1}^{r}\alpha_{j}f_{j}(p_{j},\varphi_{j},g_{j})\big{\rangle}_{L^{2}(0,T;\mathscr{X})}}{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}}=0.$

Combining the previous results, we thus infer

$\displaystyle\lim_{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}\to
0}\frac{|J(\boldsymbol{\alpha},\boldsymbol{p}+\boldsymbol{g},\boldsymbol{\varphi})-J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})-\partial_{\boldsymbol{p},\boldsymbol{g}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})|}{\|\boldsymbol{g}\|_{H^{1}(0,T;\mathbb{R}^{r})}}=0,$

which concludes the proof. ∎

######  Remark 3.12.

If the family of transformation operators is uniformly bounded, i.e., there
exists some $\overline{C}>0$ such that

$\|\mathcal{T}_{i}(p_{i})\|\leq\overline{C}\qquad\text{for all
}p_{i}\in\mathbb{R},$

then it is easy to see that condition (13) is satisfied. Note that in this
case Assumption 2 is also satisfied, cf. [3, Lem 4.2]. An example for such a
family of operators is (again) the periodic shift operator.

Let us emphasize that it is essential for the directional derivative
$\partial_{\boldsymbol{p},\boldsymbol{g}}J$ to have the path variable and
associated directions in $H^{1}(0,T;\mathbb{R}^{r})$. The following example
details that if we take a direction in $L^{2}(0,T;\mathbb{R}^{r})$, then the
directional derivative may not be finite.

######  Example 3.13.

Consider the shift operator $\mathcal{T}(p)\varphi=\varphi(\cdot-p)$ with
periodic embedding into the spaces $\mathscr{X}\vcentcolon=L^{2}(0,2\pi)$ and
$\mathscr{Y}\vcentcolon=H^{1}_{\mathrm{per}}(0,2\pi)$, cf. [3, Ex. 4.3 and
5.12]. It is well-known, that the shift operator is a semi-group with
generator $-\tfrac{\partial}{\partial x}$, see for instance [8, Sec. II.2.10].
Let $z(t,x)=t^{-1/3}\cos(x)$, $r=1$, $p_{1}\equiv 0$, and
$\varphi_{1}(x)=\sin(x)$. Then for any $\alpha_{1},g_{1}\in L^{2}(0,T)$, we
obtain

$\displaystyle\partial_{\boldsymbol{p},\boldsymbol{g}}J(\alpha_{1},p_{1},\varphi_{1})$
$\displaystyle=-\left\langle\alpha_{1}\varphi_{1}-z,\alpha_{1}\tfrac{\partial}{\partial
x}\varphi_{1}g_{1}\right\rangle_{L^{2}(0,T;\mathscr{X})}$
$\displaystyle=\langle z,\alpha_{1}\tfrac{\partial}{\partial
x}\varphi_{1}g_{1}\rangle_{L^{2}(0,T;\mathscr{X})}=\|\cos\|^{2}_{\mathscr{X}}\int_{0}^{T}t^{-1/3}\alpha_{1}(t)g_{1}(t)\,\mathrm{d}t.$

We notice that for $\alpha_{1}(t)=g_{1}(t)=t^{-1/3}$ we have
$\alpha_{1},g_{1}\in L^{2}(0,T)$ but the product $t^{-1/3}\alpha_{1}g_{1}$ is
not in $L^{1}(0,T)$. We conclude
$\partial_{\boldsymbol{p},\boldsymbol{g}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\not\in\mathbb{R}$.

######  Remark 3.14.

To ensure $\boldsymbol{p}\in H^{1}(0,T;\mathbb{R}^{r})$ during a (numerical)
optimization, we may choose a suitable low-dimensional subspace with
continuously differentiable basis functions, such as the space of polynomials
with given maximal degree. The associated gradient is easily computed from
Theorem 3.10 via the chain rule. Besides the reduced computational cost, such
an approach yields an interpretable representation for the wave speeds. We
refer to [18] for a similar idea in a fully discretized setting.

##  4 Discretization

Towards a numerical implementation, we derive discretized versions of the
partial derivatives from Theorem 3.10. To shorten notation, we introduce for
$(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\in\mathscr{Z}$ and
$\boldsymbol{h}\in\mathscr{Y}^{r}$ the quantities

$\displaystyle v_{i}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\bigg{\langle}\sum_{j=1}^{r}\alpha_{j}\mathcal{T}_{j}(p_{j})\varphi_{j}-z,\mathcal{T}_{i}(p_{i})\varphi_{i}\bigg{\rangle}_{\\!\\!\\!\mathscr{X}},$
(15a)
$\displaystyle\xi_{i}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\bigg{\langle}\sum_{j=1}^{r}\alpha_{j}\mathcal{T}_{j}(p_{j})\varphi_{j}-z,\alpha_{i}\left[\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\varphi_{i}\right]\bigg{\rangle}_{\\!\\!\\!\mathscr{X}},$
(15b)
$\displaystyle\mu_{i}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi},\boldsymbol{h})$
$\displaystyle\vcentcolon=\bigg{\langle}\sum_{j=1}^{r}\alpha_{j}\mathcal{T}_{j}(p_{j})\varphi_{j}-z,\alpha_{i}\mathcal{T}_{i}(p_{i})h_{i}\bigg{\rangle}_{\\!\\!\\!\mathscr{X}},$
(15c)

for $i=1,\ldots,r$.

We start our exposition with the discretization with respect to time. To this
end, consider a time grid $0=t_{0}<t_{1}<\ldots<t_{m}=T$ and associated
quadrature rule defined by weights $\omega_{\ell}\geq 0$ for
$\ell=0,\ldots,m$. The approximation of the directional derivative of $J$ with
respect to $\boldsymbol{\alpha}$ is thus given by

$\displaystyle\partial_{\boldsymbol{\alpha},\boldsymbol{d}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle=\sum_{i=1}^{r}\int_{0}^{T}v_{i}(\boldsymbol{\alpha}(t),\boldsymbol{p}(t),\boldsymbol{\varphi})d_{i}(t)\,\mathrm{d}t$
$\displaystyle\approx\sum_{i=1}^{r}\sum_{k=0}^{m}w_{k}v_{i}(\boldsymbol{\alpha}(t_{k}),\boldsymbol{p}(t_{k}),\boldsymbol{\varphi})d_{i}(t_{k})$
$\displaystyle=(\boldsymbol{v}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}))^{\top}\left(I_{r}\otimes
W\right)\boldsymbol{d}^{m},$

where $\otimes$ is the Kronecker product and

$\displaystyle W$
$\displaystyle\vcentcolon=\mathrm{diag}\left(w_{0},\ldots,w_{m}\right)\in\mathbb{R}^{(m+1)\times(m+1)},$
$\displaystyle\boldsymbol{v}_{i}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\begin{bmatrix}v_{i}(\boldsymbol{\alpha}(t_{0}),\boldsymbol{p}(t_{0}),\boldsymbol{\varphi})&\cdots
v_{i}(\boldsymbol{\alpha}(t_{m}),\boldsymbol{p}(t_{m}),\boldsymbol{\varphi})\end{bmatrix}^{\top}\in\mathbb{R}^{m+1},$
$\displaystyle\boldsymbol{v}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\begin{bmatrix}\boldsymbol{v}_{1}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})^{\top}&\cdots&\boldsymbol{v}_{r}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{r(m+1)},$
$\displaystyle\boldsymbol{d}_{i}^{m}$
$\displaystyle\vcentcolon=\begin{bmatrix}d_{i}(t_{0})&\cdots
d_{i}(t_{m})\end{bmatrix}^{\top}\in\mathbb{R}^{m+1},$
$\displaystyle\boldsymbol{d}^{m}$
$\displaystyle\vcentcolon=\begin{bmatrix}(\boldsymbol{d}_{1}^{m})^{\top}&\cdots&(\boldsymbol{d}_{r}^{m})^{\top}\end{bmatrix}^{\top}\in\mathbb{R}^{r(m+1)}.$

The time-discrete approximation of the partial derivative is thus given by

$\partial_{\boldsymbol{\alpha}}J^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\vcentcolon=(\boldsymbol{v}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}))^{\top}\left(I_{r}\otimes
W\right)\in\mathbb{R}^{1\times r(m+1)}.$

Analogously, the time-discrete approximation of the partial derivative of $J$
with respect to the path variables, i.e.,
$\partial_{\boldsymbol{p}}J(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$,
is given by

$\partial_{\boldsymbol{p}}J^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\vcentcolon=(\boldsymbol{\xi}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}))^{\top}\left(I_{r}\otimes
W\right)\in\mathbb{R}^{1\times r(m+1)},$

with
$\boldsymbol{\xi}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
defined analogously as
$\boldsymbol{v}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$.
In the same fashion, we obtain the time-discrete approximation for the
directional derivative with respect to the mode variables as

$\displaystyle\partial_{\boldsymbol{\varphi},\boldsymbol{h}}J^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})\vcentcolon=(\boldsymbol{\mu}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi},\boldsymbol{h}))^{\top}\left(I_{r}\otimes
W\right)\mathbf{1}_{r(m+1)}\in\mathbb{R},$

where we denote by $\mathbf{1}_{r(m+1)}\in\mathbb{R}^{r(m+1)}$ the vector with
all entries equal to $1$, and $\boldsymbol{\mu}^{m}$ defined analogously as
$\boldsymbol{v}^{m}$.

For the spatial discretization, let $\mathscr{Y}_{n}$ denote an
$n$-dimensional subspace of $\mathscr{Y}$ with basis functions
$\psi_{1},\ldots,\psi_{n}\in\mathscr{Y}$. Let us define for $i,j=1,\ldots,r$
and $p_{i},p_{j}\in\mathbb{R}$ the matrices
$M_{i,j},N_{i,j},F_{i},G_{i}\in\mathbb{R}^{n\times n}$ via

$\displaystyle\left[M_{i,j}(p_{i},p_{j})\right]_{k,\ell}$
$\displaystyle\vcentcolon=\langle\mathcal{T}_{j}(p_{j})\psi_{k},\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}},$
(16a) $\displaystyle\left[N_{i,j}(p_{i},p_{j})\right]_{k,\ell}$
$\displaystyle\vcentcolon=\langle\mathcal{T}_{j}(p_{j})\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}},$ (16b)
$\displaystyle\left[F_{i}(p_{i})\right]_{k,\ell}$
$\displaystyle\vcentcolon=\langle\psi_{k},\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}},$
(16c) $\displaystyle\left[G_{i}(p_{i})\right]_{k,\ell}$
$\displaystyle\vcentcolon=\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}},$ (16d)

for $k,\ell=1,\ldots,r$. For the data $z\in L^{2}(0,T;\mathscr{Y})$, the modes
$\varphi_{i}\in\mathscr{Y}$, and directions $h_{i}\in\mathscr{Y}$, we consider
the approximations

$z(t)\approx\sum_{\ell=1}^{n}\widehat{z}_{\ell}(t)\psi_{\ell},\qquad\varphi_{i}\approx\sum_{\ell=1}^{n}\widehat{\varphi}_{i,\ell}\psi_{\ell},\qquad
h_{i}\approx\sum_{\ell=1}^{n}\widehat{h}_{i,\ell}\psi_{\ell},$

with

$\displaystyle\boldsymbol{\widehat{z}}(t)\vcentcolon=[\widehat{z}_{1}(t)\;\cdots\;\widehat{z}_{n}(t)]^{\top}\in\mathbb{R}^{n},\qquad\qquad\boldsymbol{\widehat{\varphi}}_{i}\vcentcolon=[\widehat{\varphi}_{i,1}\;\cdots\;\widehat{\varphi}_{i,n}]^{\top}\in\mathbb{R}^{n},$
$\displaystyle\boldsymbol{\widehat{h}}_{i}\vcentcolon=[\widehat{h}_{i,1}\;\cdots\;\widehat{h}_{i,n}]^{\top}\in\mathbb{R}^{n}.$

With these preparations, we obtain the spatial approximation of the inner
products in (15) as

$\displaystyle\widehat{v}_{i}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\bigg{(}\sum_{j=1}^{r}\alpha_{j}\boldsymbol{\widehat{\varphi}}_{j}^{\top}M_{i,j}(p_{i},p_{j})-\widehat{\boldsymbol{z}}^{\top}F_{i}(p_{i})\bigg{)}\boldsymbol{\widehat{\varphi}}_{i}\in\mathbb{R},$
$\displaystyle\widehat{\xi}_{i}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\alpha_{i}\bigg{(}\sum_{j=1}^{r}\boldsymbol{\widehat{\varphi}}_{j}^{\top}N_{i,j}(p_{i},p_{j})-\widehat{\boldsymbol{z}}^{\top}G_{i}(p_{i})\bigg{)}\boldsymbol{\widehat{\varphi}}_{i}\in\mathbb{R},$
$\displaystyle\widehat{\mu}_{i}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\alpha_{i}\bigg{(}\sum_{j=1}^{r}\alpha_{j}\boldsymbol{\widehat{\varphi}}_{j}^{\top}M_{i,j}(p_{i},p_{j})-\widehat{\boldsymbol{z}}^{\top}F_{i}(p_{i})\bigg{)}\in\mathbb{R}^{1\times
n}.$

We thus obtain the space- and time-discretized partial derivatives as

$\displaystyle\partial_{\boldsymbol{\alpha}}\widehat{J}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=(\widehat{\boldsymbol{v}}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}))^{\top}\left(I_{r}\otimes
W\right)\in\mathbb{R}^{1\times(m+1)r},$
$\displaystyle\partial_{\boldsymbol{p}}\widehat{J}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=(\widehat{\boldsymbol{\xi}}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}))^{\top}\left(I_{r}\otimes
W\right)\in\mathbb{R}^{1\times(m+1)r},$
$\displaystyle\partial_{\boldsymbol{\varphi}}\widehat{J}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi})$
$\displaystyle\vcentcolon=\mathbf{1}_{r(m+1)}^{\top}\left(I_{r}\otimes
W\right)\mathcal{M}\in\mathbb{R}^{1\times nr},$

with
$\mathcal{M}\vcentcolon=\mathrm{blkdiag}(\widehat{\mu}_{1}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}),\ldots,\widehat{\mu}_{r}^{m}(\boldsymbol{\alpha},\boldsymbol{p},\boldsymbol{\varphi}))$.

We conclude this section with a specific computation of the quantities
depending on the inner products for the periodic shift operator and $P_{1}$
finite elements.

######  Example 4.1.

Let us assume we have a one-dimensional domain $\Omega=(0,1)$ and a
corresponding equidistant grid of step size $h\vcentcolon=\tfrac{1}{n}$. We
discretize $\mathscr{Y}=H^{1}_{\mathrm{per}}(\Omega)$ via periodic $P_{1}$
finite element functions. For $\mathscr{X}=L^{2}(0,1)$ and shift operator with
periodic embedding, we observe

$M_{i,j}(p_{i},p_{j})=F_{i}(p_{i}-p_{j})\qquad\text{and}\qquad
N_{i,j}(p_{i},p_{j})=G_{i}(p_{i}-p_{j}).$

For $p_{i}=qh+\tilde{p}_{i}$ with $q\in\mathbb{Z}$ and $\tilde{p}_{i}\in[0,h)$
we obtain

For further details, including the computation of
$\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}}$, we refer to
Appendix A.

##  5 Numerical examples

For our numerical examples, we use an equidistant time grid
$t_{i}\vcentcolon=i\tau$ with step size $\tau>0$. The weights for the time
integration are chosen based on the trapezoidal rule. For the discretization
in space, we also use an equidistant grid and follow Example 4.1 for
approximating the inner products occurring in the cost functional and the
gradient. The optimization itself is carried out with the MATLAB®package
GRANSO with default settings, see [7]. The GRANSO algorithm is based on a
quasi-Newton solver and can handle non-convex, non-smooth, constrained
optimization problems. For computing the relative $L^{2}$ errors of the
approximations, we discretized the corresponding integrals by the trapezoidal
rule.

For notational convenience, we assumed so far that there is exactly one mode
per transformation operator. In practice, it is often more reasonable to
cluster the modes into different reference frames, see, for instance, [3, sec.
7.1]. Thus, we use the clustered approximation ansatz

$z\approx\sum_{\rho=1}^{q}\sum_{i=1}^{r_{\rho}}\alpha_{\rho,i}\mathcal{T}_{\rho}\left(p_{\rho}\right)\varphi_{\rho,i}$
(18)

for the following numerical experiments and emphasize that this only requires
a minor and straightforward modification of the gradient. We denote the
approximation based on our optimization results with sPOD, not to be confused
with the original sPOD algorithm from [24]. Furthermore, we use dashed lines
in the plots to display the (optimized) path variables.

###  5.1 Viscous Burgers’ equation

We consider the one-dimensional viscous Burgers’ equation

$\displaystyle\tfrac{\partial}{\partial t}z(t,x)$
$\displaystyle=\tfrac{1}{\mathrm{Re}}\tfrac{\partial^{2}}{\partial
x^{2}}z(t,x)-z(t,x)\tfrac{\partial}{\partial x}z(t,x),$
$\displaystyle(t,x)\in(0,2)\times(0,1),$ (19)

and, following [17], use the analytical solution

$\displaystyle
z(t,x)=\frac{x}{t+1}\bigg{(}1+\sqrt{\tfrac{t+1}{\exp\left(\tfrac{\mathrm{Re}}{8}\right)}}\exp\left(\mathrm{Re}\tfrac{x^{2}}{4t+4}\right)\\!\\!\\!\bigg{)}^{\\!\\!-1}$

with Reynolds number $\mathrm{Re}=1000$ for our experiment. We test our
algorithm with data obtained from the analytical solution, evaluated on a grid
with $100$ equidistant intervals in space and time, respectively. We optimize
for an approximation with a single frame, i.e., $q=1$ in (18), and $r_{1}=r=2$
modes, and supply the first snapshot and the zero vector as starting values
for the modes. We initialize the corresponding coefficients as a constant
function with value $1$. For the path, we start with a straight line given by
$p(t)=\frac{37}{200}t$, evaluated at the time grid points. The results are
depicted in Figure 1, detailing that already with $r=2$, an accurate
approximation with a relative $L^{2}$ error of less than
$\mathrm{3}\mathrm{\char 37\relax}$ can be achieved, while the POD
approximation is not able to reproduce the shock front.

$x$$t$ (a) snapshot data

$x$$t$ (b) sPOD $r=2$

$x$$t$ (c) POD $r=2$

Figure 1: Burgers’ equation – comparison of original data and approximations

Besides, we compare the relative $L^{2}$ errors of the approximations for
different mode numbers. For the initialization of the optimization algorithm,
we use the first $r$ snapshots. The coefficients and the path for the
optimization are initialized as before. The resulting errors are presented in
Table 1, detailing the superior approximation capabilities of our method for
this test case.

Table 1: Burgers’ equation – comparison of relative $L^{2}$ errors $r$ | sPOD | POD
---|---|---
$1$ |  $1.217\text{\times}{10}^{-1}$ | $4.499\text{\times}{10}^{-1}$
$2$ |  $2.887\text{\times}{10}^{-2}$ | $2.853\text{\times}{10}^{-1}$
$3$ |  $1.312\text{\times}{10}^{-2}$ | $2.102\text{\times}{10}^{-1}$
$4$ |  $8.406\text{\times}{10}^{-3}$ | $1.654\text{\times}{10}^{-1}$
$5$ |  $6.565\text{\times}{10}^{-3}$ | $1.347\text{\times}{10}^{-1}$

###  5.2 Nonlinear Schrödinger equation

In this section, we consider the nonlinear Schrödinger equation

$\displaystyle\mathrm{i}\tfrac{\partial}{\partial t}z(x,t)$
$\displaystyle=-\tfrac{1}{2}\tfrac{\partial^{2}}{\partial
x^{2}}z(t,x)+\kappa|z(t,x)|^{2}z(t,x),$ $\displaystyle z(0,t)$
$\displaystyle=2\operatorname{sech}(x+7)\exp(2\mathrm{i}x)+2\operatorname{sech}(x-7)\exp(-2\mathrm{i}x),$

as presented in [18, Example 4]. We compute a solution using the code from
[18] on a uniform grid of $501\times 1024$ points on the domain
$\mathbb{T}\times\Omega=[0,2\pi]\times[-15,15]$. The absolute value of the
numerical solution is presented in Figure 2(a).

We initialize our algorithm by assuming an approximation with two frames, each
with two modes. As starting value for the modes, we use

$\displaystyle\varphi_{\rho,1}(x)=2\operatorname{sech}(x-(-1)^{\rho}7)\exp(-(-1)^{\rho}2\mathrm{i}x)\qquad\text{and}\qquad\varphi_{\rho,2}\equiv
0$

for $\rho=1,2$. The coefficients are initialized as constants, with value $1$
at each time point. For the initial paths, we use $p_{1}(t)=2t$ and
$p_{2}(t)=-2t$. The corresponding approximation and the absolute error are
presented in Figures 2(b) and 2(c). We observe that the error results mainly
from the complicated wave dynamics in the middle of the spatial and time
domain, whereas the transported wave profiles are captured accurately. Let us
emphasize that the error is very localized such that it can be captured with
only a few additional POD modes.

$x$$t$ (a) snapshot data

$x$$t$ (b) sPOD $r=4$

$x$$t$ (c) absolute error

Figure 2: Nonlinear Schrödinger equation – original data, approximation, and
absolute error for initial paths $p_{1}(t)=2t$ and $p_{2}(t)=-2t$

We notice that the optimizer does not keep the linear path over the whole time
domain. Instead, as depicted in Figure 2(b), in the middle of the
computational domain, the paths jump between the wavefronts. Inspecting the
snapshot matrix in the co-moving frame along the path $p_{1}(t)=2t$ in Figure
3 provides a possible explanation: the vertical wavefront features an offset
after the two waves have crossed. The optimizer needs to account for this
offset, which explains the jump.

$x$$t$ Figure 3: Nonlinear Schrödinger equation – transformation to the co-
moving frame along $p_{1}(t)=2t$

Let us emphasize that with a different initialization for the path variables,
the optimizer finds another local minimum with a similar approximation
quality. Using piecewise linear paths as depicted in Figure 4(a), we observe
that the resulting optimized path smoothes out the edges of the initial path
in the middle of the domain (cf. Figure 4(c)) and does not feature any jumps.
It is smooth and tracks the wavefronts as if the waves reflect off each other.

$x$$t$ (a) initial paths

$x$$t$ (b) sPOD $r=4$

$x$$t$ (c) zoom: initial vs. optimized path

Figure 4: Nonlinear Schrödinger equation – original data and approximation for
a piecewise linear initialization of the path variables as depicted in Figure
4(a)

###  5.3 FitzHugh–Nagumo wave train

We follow [11] and consider the FitzHugh–Nagumo model given by

$\displaystyle\tfrac{\partial}{\partial t}u(t,x)$
$\displaystyle=\nu\tfrac{\partial^{2}}{\partial
x^{2}}u(t,x)-v(t,x)+u(t,x)(1-u(t,x))(u(t,x)-a),$ (20)
$\displaystyle\tfrac{\partial}{\partial t}v(t,x)$
$\displaystyle=\epsilon(bu(t,x)-v(t,x)),$

with spatial domain $(0,500)$ and time interval $(0,1000)$. The partial
differential equation (20) is closed by periodic boundary conditions and the
initial condition

$u(0,x)=\tfrac{1}{2}\left(1+\sin\left(\tfrac{\pi}{50}x\right)\right),\quad
v(0,x)=\tfrac{1}{2}\left(1+\cos\left(\tfrac{\pi}{50}x\right)\right).$

For the parameter values, we choose $\nu=1$, $a=-0.1$, $\epsilon=0.05$, and
$b=0.3$. The spatial discretization of (20) is performed via a central sixth-
order finite-difference scheme with mesh width $h=0.5$ and for the time
integration we use MATLAB®’s ode45 function based on a time grid with step
size $1$. The corresponding numerical solution for the variable $u$ is
depicted in Figure 5(a).

For the optimization we consider only the data of the variable $u$ and use an
approximation with one reference frame to account for the traveling wave
train. Furthermore, we reduce the computational complexity by considering the
optimization problem only in terms of the path, whereas the coefficients and
modes are computed in each iteration via a truncated singular value
decomposition of the snapshot matrix shifted into the co-moving reference
frame. Here we exploit that the periodic shift operator is isometric, such
that we can solve the optimization problem via classical POD with transformed
data, cf. [3, Thm. 4.8]. As starting value for the path, we choose a linear
function in $t$ with slope $1.04$, which we determined by inspecting the first
and the last snapshot of the original data. The corresponding approximation
obtained from the optimization procedure is depicted in Figure 5(b). As
reference approximation, we consider a POD approximation with the same number
of modes in Figure 5(c). The corresponding total relative errors are $15\%$
for the approximation based on shifted modes and $31\%$ for the POD
approximation. We note that in contrast to the Burgers test case considered in
section 5.1, the traveling wave train can be better approximated by POD due to
the lack of a traveling shock wave. Correspondingly, the difference between
the POD approximation and the approximation based on shifted modes for the
considered FitzHugh–Nagumo test case is less striking than the one observed
for the example in section 5.1.

$x$$t$ (a) snapshot data

$x$$t$ (b) sPOD $r=4$

$x$$t$ (c) POD $r=4$

Figure 5: FitzHugh–Nagumo model – comparison of original data and
approximations

##  6 Summary

In this paper, we analyze the problem of determining an optimal approximation
of given snapshot data by a linear combination of dynamically transformed
modes. This data compression can, for instance, be used for model order
reduction of transport-dominated systems, see, for example, [3] for an
approach that makes use of such decompositions to construct dynamical ROMs via
projection. As optimization parameters, we consider the modes, the
corresponding coefficients or amplitudes, and the so-called path variables,
which parameterize the coordinate transforms applied to the modes. We first
show that the considered infinite-dimensional optimization problem possesses a
minimizing solution if the admissible set is constrained such that the
optimization parameters are norm bounded. Afterward, we derive a corresponding
unconstrained optimization problem by adding an appropriate penalization term
and show that the unconstrained problem also has a solution. Furthermore, we
demonstrate that if the penalization coefficient tends to infinity, each limit
point of the corresponding sequence of minimizers is a solution to the
original constrained optimization problem. To derive a gradient-based
optimization procedure, we compute the partial Fréchet derivatives of the
unconstrained cost functional and discuss their space and time discretization.
Finally, we apply the optimization procedure to some numerical test cases and
observe that the optimized decompositions are significantly more accurate than
corresponding approximations obtained by the classical proper orthogonal
decomposition with the same number of modes.

After full discretization, the optimization problem still features a large
number of optimization parameters scaling with the number of grid points in
space and time. Thus, an interesting future research direction is to
investigate approaches for reducing the computational complexity of the
optimization procedure, for instance, by using multigrid optimization
techniques [19], or by making use of low-dimensional parametrizations of the
optimization parameters, cf. Remark 3.14. Such a parametrization seems to be
especially promising for reducing the complexity in the path variables since
our numerical experiments revealed that the optimization procedure is
sensitive with respect to the paths. Furthermore, let us emphasize that,
although we have only discussed applications in a one-dimensional spatial
domain with a periodic shift operator, our framework is not restricted to this
case. Thus, another promising direction for the future is to explore the
applicability to problems with higher-dimensional spatial domains using
different transformation operators, see e.g. [13, 26, 29] for some
contributions in this direction.

#### Acknowledgments

The work of F. Black is supported by the Deutsche Forschungsgemeinschaft (DFG,
German Research Foundation) Collaborative Research Center (CRC) 1029
_Substantial efficiency increase in gas turbines through direct use of coupled
unsteady combustion and flow dynamics_ , project number 200291049. P. Schulze
acknowledges funding by the DFG CRC Transregio 154 _Mathematical Modelling,
Simulation and Optimization Using the Example of Gas Networks_ , project
number 239904186. B. Unger acknowledges funding from the DFG under Germany’s
Excellence Strategy – EXC 2075 – 390740016 and is thankful for support by the
Stuttgart Center for Simulation Science (SimTech).

## References

  * [1] R. A. Adams and J. J. F. Fournier.  Sobolev Spaces.  Elsevier, Oxford, UK, second edition, 2003.
  * [2] W. J. Beyn and V. Thümmler.  Freezing solutions of equivariant evolution equations.  SIAM J. Appl. Dyn. Syst., 3(2):85–116, 2004.
  * [3] F. Black, P. Schulze, and B. Unger.  Projection-based model reduction with dynamically transformed modes.  ESAIM: Math. Model. Numer. Anal., 54(6):2011–2043, 2020.
  * [4] F. Black, P. Schulze, and B. Unger.  Efficient wildland fire simulation via nonlinear model order reduction.  ArXiv preprint 2106.11381, 2021.
  * [5] F. Black, P. Schulze, and B. Unger.  Model order reduction with dynamically transformed modes for the wave equation.  PAMM, 20(1):e202000321, 2021.
  * [6] N. Cagniart, Y. Maday, and B. Stamm.  Model order reduction for problems with large convection effects.  In B. N. Chetverushkin, W. Fitzgibbon, Y. A. Kuznetsov, P. Neittaanmäki, J. Periaux, and O. Pironneau, editors, Contributions to Partial Differential Equations and Applications, pages 131–150. Springer, Cham, Switzerland, 2019.
  * [7] F. E. Curtis, T. Mitchell, and M. L. Overton.  A BFGS-SQP method for nonsmooth, nonconvex, constrained optimization and its evaluation using relative minimization profiles.  Optim. Methods Softw., 32(1):148–181, 2017.
  * [8] K. J. Engel and R. Nagel.  One-Parameter Semigroups for Linear Evolution Equations.  Springer, New York, NY, USA, 2000.
  * [9] C. Greif and K. Urban.  Decay of the Kolmogorov $n$-width for wave problems.  Appl. Math. Lett., 96:216–222, 2019.
  * [10] M. Gubisch and S. Volkwein.  Chapter 1: Proper orthogonal decomposition for linear-quadratic optimal control.  In P. Benner, A. Cohen, M. Ohlberger, and K. Willcox, editors, Model Reduction and Approximation, pages 3–63. SIAM, Philadelphia, PA, USA, 2017\.
  * [11] J. Koch.  Data-driven modeling of nonlinear traveling waves.  Chaos, 31(4):043128, 2021.
  * [12] A. Kolmogoroff.  Über die beste Annäherung von Funktionen einer gegebenen Funktionenklasse.  Ann. Math., 37(1):107–110, 1936.
  * [13] P. Krah, M. Sroka, and J. Reiss.  Model order reduction of combustion processes with complex front dynamics.  In F. J. Vermolen and C. Vuik, editors, Numerical Mathematics and Advanced Applications ENUMATH 2019, pages 803–811, Cham, Switzerland, 2021\. Springer.
  * [14] K. Lee and K. T. Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  J. Comput. Phys., page 108973, 2020.
  * [15] D. G. Luenberger and Y. Ye.  Linear and Nonlinear Programming.  Springer New York, NY, USA, fourth edition, 2016.
  * [16] Y. Maday, A. T. Patera, and G. Turinici.  A priori convergence theory for reduced-basis approximations of single-parameter elliptic partial differential equations.  J. Sci. Comput., 17(1):437–446, 2002.
  * [17] R. Mauli, B. Lusch, and P. Balaprakash.  Reduced-order modeling of advection-dominated systems with recurrent neural networks and convolutional autoencoders.  Phys. Fluids, 33:037106, 2021.
  * [18] A. Mendible, S. L. Brunton, A. Y. Aravkin, W. Lowrie, and J. N. Kutz.  Dimensionality reduction and reduced-order modeling for traveling wave physics.  Theor. Comput. Fluid Dyn., 34:385–400, 2020.
  * [19] S. G. Nash.  A multigrid approach to discretized optimization problems.  Optim. Methods Softw., 14(1–2):99–116, 2000.
  * [20] M. Nonino, F. Ballarin, G. Rozza, and Y. Maday.  Overcoming slowly decaying Kolmogorov $n$-width by transport maps: application to model order reduction of fluid dynamics and fluid structure interaction problems.  ArXiv preprint 1911.06598, 2019.
  * [21] M. Ohlberger and S. Rave.  Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing.  C. R. Acad. Sci. Paris, 351(23–24):901–906, 2013.
  * [22] B. Peherstorfer.  Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling.  SIAM J. Sci. Comput., 42(5):A2803–A2836, 2020.
  * [23] J. Reiss.  Optimization-based modal decomposition for systems with multiple transports.  SIAM J. Sci. Comput., 43(3):A2079–A2101, 2021.
  * [24] J. Reiss, P. Schulze, J. Sesterhenn, and V. Mehrmann.  The shifted proper orthogonal decomposition: a mode decomposition for multiple transport phenomena.  SIAM J. Sci. Comput., 40(3):A1322–A1344, 2018.
  * [25] D. Rim, S. Moe, and R. J. LeVeque.  Transport reversal for model reduction of hyperbolic partial differential equations.  SIAM/ASA J. Uncertain. Quantif., 6(1):118–150, 2018.
  * [26] D. Rim, B. Peherstorfer, and K. T. Mandli.  Manifold approximations via transported subspaces: model reduction for transport-dominated problems.  ArXiv preprint 1912.13024v2, 2020.
  * [27] W. Rudin.  Real and Complex Analysis.  McGraw Hill, New York, NY, USA, third edition, 1987.
  * [28] P. Schulze, J. Reiss, and V. Mehrmann.  Model reduction for a pulsed detonation combuster via shifted proper orthogonal decomposition.  In R. King, editor, Active Flow and Combustion Control 2018, pages 271–286. Springer, Cham, Switzerland, 2019.
  * [29] T. Taddei.  A registration method for model order reduction: data compression and geometry reduction.  SIAM J. Sci. Comput., 42(2):A997–A1027, 2020.
  * [30] T. Taddei, S. Perotto, and A. Quarteroni.  Reduced basis techniques for nonlinear conservation laws.  ESAIM Math. Model. Numer. Anal., 49(3):787–814, 2015.
  * [31] B. Unger and S. Gugercin.  Kolmogorov $n$-widths for linear dynamical systems.  Adv. Comput. Math., 45(5-6):2273–2286, 2019.
  * [32] G. Welper.  Transformed snapshot interpolation with high resolution transforms.  SIAM J. Sci. Comput., 42(4):A2037–A2061, 2020.
  * [33] E. Zeidler.  Nonlinear Functional Analysis and its Applications IIa: Linear Monotone Operators.  Springer, New York, NY, USA, 1990.

##  Appendix A Shifted Inner Products of Hat Functions

In this section, we explicitly compute the path-dependent inner products for a
particular example. Let us assume a one-dimensional domain $\Omega=(0,1)$,
which we discretize with an equidistant grid with step size
$h\vcentcolon=\tfrac{1}{n}$ for some $n\in\mathbb{N}$. We consider the spaces
$\mathscr{X}=L^{2}(\Omega)$ and $\mathscr{Y}=H^{1}_{\mathrm{per}}(\Omega)$.
For the discretization of $\mathscr{Y}$ we choose periodic $P_{1}$ finite
elements, which are given as

$\displaystyle\psi_{\ell}(x)=\begin{cases}\frac{x-(\ell-1)h}{h}&\text{for
}x\in((\ell-1)h,\ell h],\\\ \frac{(\ell+1)h-x}{h}&\text{for }x\in(\ell
h,(\ell+1)h),\\\ 0&\text{otherwise},\end{cases}\quad\text{for
}\ell=1,\ldots,n-1,$ and
$\displaystyle\psi_{n}(x)=\begin{cases}\frac{x-(n-1)h}{h}&\text{for
}x\in((n-1)h,nh],\\\ \frac{h-x}{h}&\text{for }x\in(0,h),\\\
0&\text{otherwise},\end{cases}$

with derivatives

$\displaystyle\psi^{\prime}_{k}=\begin{cases}\frac{1}{h}&\text{for
}x\in((k-1)h,kh),\\\ -\frac{1}{h}&\text{for }x\in(kh,(k+1)h),\\\
0&\text{otherwise},\end{cases}\quad\text{for }k=1,\ldots,n-1,$ and
$\displaystyle\psi^{\prime}_{n}=\begin{cases}\frac{1}{h}&\text{for
}x\in((n-1)h,nh),\\\ -\frac{1}{h}&\text{for }x\in(0,h),\\\
0&\text{otherwise}.\end{cases}$

For the family of transformation operators we choose the shift operator with
periodic embedding, cf. Example 3.13, i.e.,
$\mathcal{T}_{i}(p_{i})\varphi_{i}=\varphi_{i}(\cdot-p_{i})$.

######  Lemma A.1.

Consider the settings as described above. Then for $i,j=1,\ldots,r$,
$\varphi_{i}\in\mathscr{Y}$, and $p_{i}\in\mathbb{R}$ we have

$\displaystyle\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\varphi_{i}=-\mathcal{T}_{i}(p_{i})\tfrac{\partial}{\partial
x}\varphi_{i},\qquad M_{i,j}(p_{i},p_{j})=F_{i}(p_{i}-p_{j}),$ $\displaystyle
N_{i,j}(p_{i},p_{j})=G_{i}(p_{i}-p_{j}),$

with $M_{i,j}$, $N_{i,j}$, $F_{i}$, and $G_{i}$ defined in (16).

######  Proof A.2.

The statement is a mere consequence of the fact that the shift operator is
unitary and a semi-group and $-\tfrac{\partial}{\partial x}$ is its generator,
see for instance [8, Sec. II.2.10]. ∎

It is thus sufficient to compute the inner products for $F_{i}$ and $G_{i}$.
We first observe that for $p_{i}=qh$ for some $q\in\mathbb{Z}$ we obtain

$\displaystyle\mathcal{T}_{i}(p_{i})\psi_{\ell}=\psi_{\ell}(\cdot-
qh)=\psi_{\ell+q}=\mathcal{T}_{i}(0)\psi_{\ell+q}$

with the understanding that for $\ell+q\not\in\\{1,\ldots,n\\}$ we use
$((\ell+q-1)\mod n)+1$ instead. It is thus sufficient to compute the inner
products for $p_{i}\in[0,h)$.

######  Lemma A.3.

Consider the setting as described above. Then for $p_{i}\in[0,h)$ we obtain

and

$\displaystyle\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle=\begin{cases}-\tfrac{1}{h^{2}}\left(2p_{i}h-\tfrac{3}{2}p_{i}^{2}\right),&\text{if
}\ell=k,\\\ -\tfrac{1}{2h^{2}}(h-p_{i})^{2},&\text{if }\ell=k+1,\\\
-\tfrac{1}{h^{2}}\left(2p_{i}^{2}-2hp_{i}-\tfrac{1}{2}(h-p_{i})^{2}\right),&\text{if
}\ell=k-1,\\\ \tfrac{1}{2h^{2}}p_{i}^{2},&\text{if }\ell=k-2,\\\
0,&\text{otherwise.}\end{cases}$

######  Proof A.4.

We immediately obtain
$\langle\psi_{k},\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}}=0$ and
$\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{\ell}\rangle_{\mathscr{X}}=0$ for
$k=1,\ldots,n$ and $\ell\not\in\\{k-2,k-1,k,k+1\\}$. For $k=1,\ldots,n$ we
compute

$\displaystyle\langle\psi_{k}$
$\displaystyle,\mathcal{T}_{i}(p_{i})\psi_{k}\rangle_{\mathscr{X}}=\int_{0}^{1}\psi_{k}(x)\psi_{k}(x-p_{i})\,\mathrm{d}x$
$\displaystyle=\int_{(k-1)h+p_{i}}^{kh}\psi_{k}(x)\psi_{k}(x-p_{i})\,\mathrm{d}x+\int_{kh}^{kh+p_{i}}\psi_{k}(x)\psi_{k}(x-p_{i})\,\mathrm{d}x$
$\displaystyle\phantom{=}\quad+\int_{kh+p_{i}}^{(k+1)h}\psi_{k}(x)\psi_{k}(x-p_{i})\,\mathrm{d}x$
$\displaystyle=\int_{0}^{h-p_{i}}\tfrac{x(x+p_{i})}{h^{2}}\,\mathrm{d}x+\int_{0}^{p_{i}}\tfrac{(x+h-p_{i})(h-x)}{h^{2}}\,\mathrm{d}x-\int_{p_{i}-h}^{0}\tfrac{(p_{i}-x)x}{h^{2}}\,\mathrm{d}x$
$\displaystyle=\tfrac{1}{h^{2}}\left(\tfrac{2}{3}(h-p_{i})^{3}+p_{i}(h-p_{i})^{2}+p_{i}h(h-p_{i})+\tfrac{1}{6}p_{i}^{3}\right).$

Similarly, we obtain

$\displaystyle\langle\psi_{k},\mathcal{T}_{i}(p_{i})\psi_{k+1}\rangle_{\mathscr{X}}$
$\displaystyle=\int_{kh+p_{i}}^{(k+1)h}\psi_{k}(x)\psi_{k+1}(x-p_{i})\,\mathrm{d}x$
$\displaystyle=\tfrac{1}{h^{2}}\int_{0}^{h-p_{i}}(h-p_{i}-x)x\,\mathrm{d}x=\tfrac{1}{6h^{2}}(h-p_{i})^{3}.$

Furthermore, we have

$\displaystyle\langle\psi_{k}$
$\displaystyle,\mathcal{T}_{i}(p_{i})\psi_{k-1}\rangle_{\mathscr{X}}$
$\displaystyle=\int_{(k-1)h}^{(k-1)h+p_{i}}\psi_{k}(x)\psi_{k-1}(x-p_{i})\,\mathrm{d}x+\int_{(k-1)h+p_{i}}^{kh}\psi_{k}(x)\psi_{k-1}(x-p_{i})\,\mathrm{d}x$
$\displaystyle\phantom{=}\qquad+\int_{kh}^{kh+p_{i}}\psi_{k}(x)\psi_{k-1}(x-p_{i})\,\mathrm{d}x$
$\displaystyle=\int_{0}^{p_{i}}\tfrac{(x+h-p_{i})x}{h^{2}}\,\mathrm{d}x+\int_{0}^{h-p_{i}}\tfrac{(h-x)(x+p_{i})}{h^{2}}\,\mathrm{d}x-\int_{-p_{i}}^{0}\tfrac{x(h-p_{i}-x)}{h^{2}}\,\mathrm{d}x$
$\displaystyle=\tfrac{1}{h^{2}}\left(\tfrac{1}{6}(h-p_{i})^{3}-\tfrac{1}{3}p_{i}^{3}+h^{2}p_{i}\right).$

In addition, we obtain

$\displaystyle\langle\psi_{k},\mathcal{T}_{i}(p_{i})\psi_{k-2}\rangle_{\mathscr{X}}$
$\displaystyle=\int_{(k-1)h}^{(k-1)h+p_{i}}\psi_{k}(x)\psi_{k-2}(x-p_{i})\,\mathrm{d}x$
$\displaystyle=\tfrac{1}{h^{2}}\int_{0}^{p_{i}}(p_{i}-x)x\,\mathrm{d}x=\tfrac{1}{6h^{2}}p_{i}^{3}.$

For the derivatives, we obtain for suitable $k\in\\{1,\ldots,n\\}$

$\displaystyle\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{k}\rangle_{\mathscr{X}}$
$\displaystyle=-\int_{0}^{1}\psi_{k}(x)\psi_{k}^{\prime}(x-p_{i})\,\mathrm{d}x$
$\displaystyle=-\int_{p_{i}}^{h}\tfrac{x}{h^{2}}\,\mathrm{d}x+\int_{-h}^{p_{i}-h}\tfrac{x}{h^{2}}\,\mathrm{d}x-\int_{p_{i}-h}^{0}\tfrac{x}{h^{2}}\,\mathrm{d}x$
$\displaystyle=-\tfrac{1}{h^{2}}\left(2p_{i}h-\tfrac{3}{2}p_{i}^{2}\right),$
$\displaystyle\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{k+1}\rangle_{\mathscr{X}}$
$\displaystyle=-\int_{kh+p_{i}}^{(k+1)h}\psi_{k}(x)\psi_{k}^{\prime}(x-p_{i})\,\mathrm{d}x=\int_{p_{i}-h}^{0}\tfrac{x}{h^{2}}\,\mathrm{d}x$
$\displaystyle=-\tfrac{1}{2h^{2}}(h-p_{i})^{2},$
$\displaystyle\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{k-1}\rangle_{\mathscr{X}}$
$\displaystyle=\int_{0}^{p_{i}}-\tfrac{x}{h^{2}}\,\mathrm{d}x+\int_{p_{i}}^{h}\tfrac{x}{h^{2}}\,\mathrm{d}x-\int_{-h}^{p_{i}-h}\tfrac{x}{h^{2}}\,\mathrm{d}x$
$\displaystyle=-\tfrac{1}{h^{2}}\left(2p_{i}^{2}-2hp_{i}-\tfrac{1}{2}(h-p_{i})^{2}\right),$
$\displaystyle\langle\psi_{k},\tfrac{\partial}{\partial
p_{i}}\mathcal{T}_{i}(p_{i})\psi_{k-2}\rangle_{\mathscr{X}}$
$\displaystyle=\int_{(k-1)h}^{(k-1)h+p_{i}}\tfrac{x-(k-1)h}{h^{2}}\,\mathrm{d}x=\tfrac{1}{h^{2}}\int_{0}^{p_{i}}x\,\mathrm{d}x=\tfrac{1}{2h^{2}}p_{i}^{2},$

which concludes the proof.∎