Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling

1 Introduction

Model reduction constructs reduced models of large-scale systems of equations in a one-time high-cost offline phase and then uses the reduced models in an online phase to repeatedly compute accurate approximations of the full-model solutions with significantly reduced costs. In projection-based model reduction [46, 2, 6], a low-dimensional space is constructed that approximates well the high-dimensional solution space of the full model. Then, the full-model equations are projected onto the low-dimensional space and reduced solutions are computed by solving the projected equations. However, solutions of full models that describe transport-dominated behavior, e.g., convection-dominated flows, wave-type phenomena, shock propagation, typically exhibit high-dimensional features, which means that no low-dimensional space exists in which the full-model solutions can be approximated well; the Kolmogorov n-width is high. Thus, traditional model reduction methods typically fail for transport-dominated problems [34, 14, 51]. In this work, we exploit that transport-dominated problems typically have a rich structure that is local in nature, which we leverage to derive efficient reduced models.

Model reduction projects the full-model equations on a low-dimensional—reduced–space that is spanned by a set of basis vectors. There are many methods for constructing reduced spaces, including proper orthogonal decomposition (POD)

[7, 49], balanced truncation [29, 30], the reduced basis method [41, 53, 47, 24, 46], and interpolatory model reduction methods [23, 2]. In the following, we will be mostly concerned with full-model equations that are nonlinear in the states, where projecting the full-model equations onto a reduced space is typically insufficient to obtain a reduced model that is faster to solve than the full model, because the nonlinear terms entail computations that scale with the dimension of the full-model solution space. One remedy is empirical interpolation, which approximates nonlinear terms by evaluating them at a few, carefully selected interpolation points and approximating all other components via interpolation in a low-dimensional space [5, 22, 16]. We will build on the discrete empirical interpolation method (DEIM) [13, 15]

, which is the discrete counterpart of empirical interpolation. Note that there are other techniques for nonlinear model reduction, e.g., missing point estimation

[3]

and the Gauss-Newton with approximated tensors (GNAT) method

[12]. All these methods build on the assumption that there is a low-dimensional space in which the full-model solutions can be approximated well, which is violated in case of transport-dominated problems. However, the solutions of transport-dominated problems typically are low dimensional if considered locally in time, which we exploit by approximating the full-model solutions in local low-dimensional spaces that are adapted via low-rank basis updates over time [39, 62]. The basis updates are derived by querying the full model at selected points in the spatial domain. We derive an adaptive sampling scheme that selects where to query the full model to compute the basis updates. Our analysis shows that if the reduced-model residual is local in the spatial domain—which we observe for transport-dominated problems—then only few sampling points are necessary to adapt the local spaces with our adaptive sampling scheme, which makes our model reduction approach computationally efficient.

For reviewing the literature on model reduction for transport-dominated problems, we broadly distinguish between two lines of research. Literature that we categorize into the first line of research aims to transform the full model to recover a low-rank structure that can then be exploited with a reduced space. The work on transformations in the model reduction community seems to have originated from [34], where the transport dynamics are separated from the rest of the dynamics via freezing. In [44, 43], transport maps are constructed that reverse the transport and so recover a low-rank structure. Approaches have been developed that aim to find transformations numerically via, e.g., optimization [8]. Shifted POD [42] recovers the shift due to the transport and applies POD after having reversed the shift. The shift operator introduced in [42] is time dependent and separates different transports to be applicable to, e.g., multiple waves traveling at different wave speeds. In [56, 57], snapshots are transformed to better align them before the reduced bases are constructed. The second line of research on model reduction for transport-dominated problems constructs low-dimensional spaces that explicitly target transport-dominated problems. The work [14] formulates the greedy basis construction methods, developed in the reduced basis community [41, 53], in special norms that are better suited for transport-dominated problems. The work [1] constructs the basis via $L_{1}$ optimization and the authors of [52, 59] consider time-space discretizations. There is related literature on closure modeling [55, 18, 35, 36] where methods have been developed that construct bases with respect to different norms than the classical $L_{2}$ -optimal POD bases [25, 45]. Then, there are approaches that exploit local structure, such as the approach introduced in [51] that builds on the spatial locality of shock fronts. The work by Carlberg [11] builds on a special type of adaptation that enriches the reduced space during the online phase if an error indicator signals a high error. The dimension of the reduced space grows with each adaptation step, which means that the reduced model becomes computationally more expensive to solve over time. In [21], reduced bases are adapted over time via an auxiliary equation that describes the dynamics of the bases. Under certain conditions, the approach can be seen as an approximation of Lax pairs discussed in [26]. Dynamic low-rank approximation [27, 48, 33, 31, 32] adapts bases in time similar to our approach; however, the bases in dynamic low-rank approximation approaches are constructed so that they are valid for all parameters in a given parameter domain, which is problematic if the solution manifold contains high-dimensional features because the propagating coherent structure depends on parameters, as in, e.g., [51, Example 2.5]. In contrast, our approach is local in the parameter domain, additionally to being local in time and space. Furthermore, locality in the parameter domain allows our approach to derive basis updates from only few samples of the full model.

This work is organized as follows. Section 2 sets up the problem and gives the problem formulation. Section 3 demonstrates the local structure in transport-dominated problems and proposes an approach to adaptively sample the full model to construct basis updates. Numerical results in Section 4 on benchmark examples and a model combustor of a single-element rocket engine demonstrate that our approach achieves significant speedups and is applicable to a wide range of problems. In particular, the numerical results indicate that our approach faithfully approximates interactions between propagating coherent structures traveling at different speeds. Concluding remarks are provided in Section 5.

2 Preliminaries

We briefly discusses model reduction with empirical interpolation in Section 2.1 and then demonstrate on a toy example in Section 2.2 why these traditional model reduction methods fail for problems exhibiting transport-dominated phenomena.

2.1 Model reduction with empirical interpolation

Consider the system of discretized equations that is obtained after discretizing a partial differential equation (PDE) in space and time

q_{k} (μ) = f (q_{k + 1} (μ), μ), k = 1, \dots, K,

(1)

where $q_{k} (μ) \in R^{N}$ is the $N$ -dimensional state at time $k$ and parameter $μ \in D$ , with parameter domain $D$ . The number of time steps is $K \in N$ . The function $f : R^{N} \to R^{N}$ describes the operators of the discretized PDE and typically is nonlinear in the state $q_{k + 1} (μ)$ . The time discretization is implicit in time, which means that at each time step $k = 1, \dots, K$ , a potentially nonlinear, large-scale system of equations has to be solved, e.g., with Newton’s method.

To derive a reduced model of the full model (1) with empirical interpolation [5, 13, 15], consider the trajectory at parameter $μ \in D$

Q (μ) = [q_{1} (μ), \dots, q_{K} (μ)] \in R^{N \times K},

which is the matrix with the states $q_{1} (μ), \dots, q_{K} (μ)$ as columns. Let the columns of $U = [u_{1}, \dots, u_{n}] \in R^{N \times n}$ be the POD basis of dimension $n ≪ N$ obtained from the snapshot matrix

Q = [Q (μ_{1}), \dots, Q (μ_{M})] \in R^{N \times M K},

with parameters $μ_{1}, \dots, μ_{M} \in D$ . The space spanned by the columns of $U$ is denoted as $U \subset R^{N}$ and is a subspace of $R^{N}$

. The critical assumption of traditional model reduction is that the singular values of

Q

decay fast so that only few basis vectors are necessary to approximate well the columns of

Q

in the corresponding space

U

. Following QDEIM, introduced in [15], select the interpolation points

p_{1}, \dots, p_{n} \in {1, \dots, N}

and define the corresponding interpolation points matrix

P = [e_{p_{1}}, \dots, e_{p_{n}}] \in R^{N \times n}

, where

e_{p_{i}} \in R^{N}

is the

p_{i}

-th canonical unit vector with entry 1 at the

p_{i}

-th component and entry 0 at all other components. Note that all of the following discussion can be extended in a straightforward way to oversampled empirical interpolation (ODEIM) [38]. Define

~ f (q (μ); μ) = (P^{T} U)^{- 1} P^{T} f (q (μ); μ),

so that $U ~ f (q (μ); μ)$ is the DEIM approximation of $f (q (μ); μ)$ at state $q (μ) \in R^{N}$ and parameter $μ \in D$ . Note that computing $P^{T} f (q (μ); μ)$ typically requires evaluating $f (q (μ); μ)$ at the $n$ interpolation points $p_{1}, \dots, p_{n}$ only, see [5, 13, 15]. We overload the notation $~ f$ in the following so that if we have a reduced state $~ q (μ) \in R^{n}$ , then $~ f (~ q (μ); μ) = (P^{T} U)^{- 1} P^{T} f (U ~ q (μ); μ)$ . The reduced model corresponding to $~ f$ is

{~ q}_{k} (μ) = ~ f ({~ q}_{k + 1} (μ); μ), k = 1, \dots, K,

(2)

with the reduced trajectory $~ Q (μ) = [{~ q}_{1} (μ), \dots, {~ q}_{K} (μ)] \in R^{n \times K}$ . Note that we approximate the state and the nonlinear function in the same space $U$ , which is in contrast to the original use of DEIM in [13] and similar to model reduction via missing point estimation [3]. Once a reduced model (2) is constructed in the offline phase, it is solved in the online phase. The one-time high costs of constructing the reduced model are compensated by approximating the full-model solutions with reduced-model solutions for a large number of parameters online, see, e.g., [46, 2, 6, 40] for details on the wide range of outer-loop and many-query applications where such an offline/online splitting is beneficial.

Figure 1: Advection equation: Plot (a) shows the initial condition and the direction of the transport. Plot (b) shows the slow decay of the singular values of the snapshots.


(a) initial condition	(b) singular values

2.2 Problem formulation

It has been observed that states of problems with transport-dominated behavior can require DEIM (reduced) spaces $U$ with high dimensions, see, e.g., [34, 14, 51]. The following toy example illustrates this by demonstrating the slow decay of the singular values of the snapshot matrix of solutions of the advection equation. Let $Ω = [- 1, 1] \subset R$ be the spatial domain and consider the advection equation

\partial_{t} q (x, t) + μ \partial_{x} q (x, t) = 0, x \in Ω,

(3)

with periodic boundary conditions $q (1, t) = q (- 1, t)$ and time $t \in [0, \infty)$ . Set the initial condition to

q (x, 0) = \frac{1}{\sqrt{π 0.02}} exp (- \frac{x^{2}}{0.0002}),

which is the probability density function of a normal random variable with standard deviation 0.01 and mean 0, see Figure

1a. Discretize (3) with a second-order upwind scheme and

N = 8192

inner grid points in the spatial domain and time step size

δ t = 10^{- 6}

and end time

T = 0.08

. The singular values of the trajectory

Q (μ)

are plotted in Figure 1b for

μ = 10

. According to the decay of the singular values, a DEIM space of more than

n = 200

dimensions is necessary to approximate the trajectory

Q (μ)

with a projection error below

10^{- 15}

in the Euclidean norm, which is a rather slow decay compared to the fast decay observed in many other problems [6]. This example demonstrates that efficient reduced models of transport-dominated problems will have to exploit structure beyond the classical low-rank structure that traditional model reduction methods rely on.

Figure 2: Advection equation: The plot in (a) shows that the singular values of a local trajectory decay orders of magnitude faster than the singular values of a trajectory that is global in time in this example, which we exploit via online adaptive basis updates. Plot (b) indicates that the squared residual of DEIM approximations of states of transport-dominated problems decays rapidly. We will show that this fast decay of the residual means that basis updates can be derived from only few components—samples—of the residual.


(a) local low-rank structure	(b) decay of residual (local coherence)

3 Exploiting local structure via online adaptive basis updates

We propose AADEIM (adaptive bases and adaptive sampling) to exploit local structure to construct online adaptive reduced models of full models that exhibit transport-dominated behavior. We focus on two types of local structure: Local low-rankness that we exploit via online adaptive bases and local coherence that enables updating the bases computationally efficient via few samples from the full model. Section 3.1 describes local low-rankness and local coherence in more detail. Section 3.2 discusses the adaptation of the DEIM basis to exploit local low-rank structure and Section 3.3 shows that only few samples are necessary to derive basis updates if the adapted DEIM spaces are locally coherent. Algorithm 1 in Section 3.4 summarizes our approach and provides practical considerations.

In most of this section, we drop the dependence on the parameter $μ$ of the function $f$ , the state $q_{k}$ at time step $k$ , and the trajectory $Q$ , as well as their reduced counterparts. Parametrization of our approach AADEIM is discussed in Section 3.4.

3.1 Local structure in transport-dominated problems

Consider the toy example introduced in Section 2.2. Let $w \in N$ be a window size and consider the local trajectory $Q_{k} = [q_{k - w + 1}, \dots, q_{k}] \in R^{N \times w}$ at time step $k$ , which consists of the $w$ states from time steps $k - w + 1$ to time step $k$ . Figure 2a compares the singular values of the local trajectory $Q_{k}$ to the singular values of the global trajectory $Q$ for $w = 500$ and $k = 1$ . The results indicate that the singular values of the local trajectory $Q_{k}$ decay orders of magnitude faster than the singular values of the global trajectory $Q$ . We call this behavior—that the singular values of local trajectories decay fast while the singular values of the global trajectory decay slowly—local low-rank structure. By adapting the DEIM spaces, we will exploit this local low-rank structure in Section 3.2. Similar observations about local low-rank structure are exploited in [27, 48, 33, 31, 32, 21, 37], cf. Section 1. In Appendix A, we analyze analytically an example to demonstrate its low-rank structure.

Let us now consider the local DEIM space $U_{k}$ of dimension $n = 3$ and the corresponding local interpolation points matrix $P_{k}$ obtained from the local trajectory $Q_{k}$ . The residual of approximating the state $q_{k}$ at time step $k = 500$ with DEIM in $U_{k}$ is

r_{k} = q_{k} - U_{k} (P_{k}^{T} U_{k})^{- 1} P_{k}^{T} q_{k} .

Figure 2b shows the decay of the sorted squared components of $r_{k}$ . The decay is fast, which means that there is a high residual locally only, i.e., the residual at a few components dominates whereas at most of the $N = 8192$ components of $r_{k}$ the residual is low. Intuitively, such a fast decay of the residual means that the basis matrix $U_{k}$ of the local DEIM space needs to be corrected at a few components only. We will show that such a fast decay of the residual can be the result of a local coherence structure of the local DEIM spaces, which we will make precise and will exploit with adaptive sampling in Section 3.3.

3.2 Exploiting local low-rank structure: Basis updates

To exploit local low-rank structure as described in Section 3.1, we adapt the DEIM spaces and the DEIM interpolation points during the time steps $k = 1, \dots, K$ in the online phase. The adaptation is initialized with the DEIM basis matrix $U_{1}$ and interpolation points matrix $P_{1}$ at time step $k = 1$ . Then, at each time step $k = 1, \dots, K$ , the DEIM basis matrix $U_{k}$ is adapted via an additive low-rank update to $U_{k + 1}$ . The update is based on ADEIM [39]. The interpolation points matrix $P_{k + 1}$ is derived with QDEIM from the adapted basis matrix $U_{k + 1}$ . In the following description, the DEIM interpolant is adapted at each time step $k = 1, \dots, K$ for ease of exposition, even though all of the following directly applies to situations where the adaptation is performed at selected time steps only, e.g., every other time step or via a criterion that decides adaptively when to update the DEIM interpolant.

3.2.1 Adaptation with ADEIM

Let $k$ be the current time step and $U_{k}$ and $P_{k}$ the DEIM basis matrix and the DEIM interpolation points matrix, respectively. Consider the matrix $F_{k} \in R^{N \times w}$ and the coefficient matrix $C_{k} = (P_{k}^{T} U_{k})^{- 1} P_{k}^{T} F_{k}$ . The residual of the DEIM approximation of the columns of $F_{k}$ is $R_{k} = U_{k} C_{k} - F_{k}$ . Let now $S_{k} = [e_{s_{1}^{(k)}}, \dots, e_{s_{m}^{(k)}}] \in R^{N \times m}$ be the sampling points matrix corresponding to the $s_{1}^{(k)}, \dots, s_{m}^{(k)} \in {1, \dots, N}$ samplings points with $m > n$ . ADEIM [39] adapts the DEIM basis matrix $U_{k}$ to $U_{k + 1}$ with a low-rank update $α_{k} β_{k}^{T} \in R^{N \times n}$

U_{k + 1} = U_{k} + α_{k} β_{k}^{T},

with $α_{k} \in R^{N \times r}, β_{k} \in R^{w \times r}$ with rank $r \in N$ . The ADEIM update $α_{k} β_{k}^{T}$ is the rank- $r$ matrix that minimizes the residual $U_{k + 1} C_{k} - F_{k}$ at the sampling points $S_{k}$ in the Frobenius norm, which means that the ADEIM update $α_{k} β_{k}^{T}$ minimizes

∥ S_{k}^{T} ((U_{k} + α_{k} β_{k}^{T}) C_{k} - F_{k}) ∥_{F}^{2},

see [39] for details on how to compute $α_{k} β_{k}^{T}$ and the computational costs of computing the update.

Critical for the adaptation is the matrix $F_{k}$ , because ADEIM adapts the space $U_{k}$ such that the residual of approximating the columns of $F_{k}$ is minimized at the sampling points. A potentially good choice for the columns of $F_{k}$ would be the states $q_{k - w + 1}, \dots, q_{k}$ of the full model at time steps $k - w + 1, \dots, k$ ; however, the states of the full model are unavailable because their availability would mean the full model has been solved. Instead, we set the columns of $F_{k}$ as follows. Let ${˘ S}_{k} \in R^{N \times (N - m)}$ be the complementary sampling points matrix derived from the points ${1, \dots, N} ∖ {s_{1}^{(k)}, \dots, s_{m}^{(k)}}$ that have not been selected as sampling points. Let further ${^q}_{k} \in R^{N}$ be the vector with

S_{k}^{T} {^q}_{k} = S_{k}^{T} f (U_{k} {~ q}_{k}), {˘ S}_{k}^{T} {^q}_{k} = {˘ S}_{k}^{T} U_{k} (S_{k}^{T} U_{k})^{+} S_{k}^{T} f (U_{k} {~ q}_{k}),

(4)

which means that the components of ${^q}_{k}$ corresponding to the sampling points in $S_{k}$ are equal to the components of $f (U_{k} {~ q}_{k})$ and all other components are approximated via DEIM given by $U_{k}$ and $S_{k}$ . Note that $m > n$ and therefore the Moore-Penrose inverse $(S_{k}^{T} U_{k})^{+}$ is used in (4), instead of the inverse when $U_{k}$ and $P_{k}$ are used. Note further that ${~ q}_{k}$ is the reduced state at time $k$ and that $f$ is the function that defines the full model (1). The matrix $F_{k}$ that we use in the following for adaptation is

F_{k} = [{^q}_{k - w + 1}, \dots, {^q}_{k}] .

The vectors ${^q}_{k - w + 1}, \dots, {^q}_{k}$ serve as surrogates of the full-model states $q_{k - w}, \dots, q_{k - 1}$ to which the DEIM space is adapted, because the full-model states $q_{k + 1}$ satisfy $q_{k} = f (q_{k + 1})$ for $k = 1, \dots, K$ , see equation (1). Computing a vector ${^q}_{k}$ requires evaluating the full-model function $f$ at the sampling points $S_{k}$ .

3.2.2 Analysis of the error of the adapted space

We now provide an analysis of the ADEIM adaptation. Let $U$ and $¯ U$ be two $n$ -dimensional subspaces of $R^{N}$ . We measure the distance between $U$ and $¯ U$ as

d (¯ U, U) = ∥ ¯ U - U U^{T} ¯ U ∥_{F}^{2},

(5)

where $U$ and $¯ U$ are orthonormal basis matrices of $U$ and $¯ U$ , respectively. The distance $d (¯ U, U)$ is symmetric and invariant under orthogonal basis transformations, which is true because $d (¯ U, U) = n - ∥ U^{T} ¯ U ∥_{F}^{2}$ holds. The $n$ -dimensional subspace of $R^{N}$ to which we want to adapt at iteration $k$ is denoted as ${¯ U}_{k + 1}$ and the adapted space is $U_{k + 1}$ . The following lemma quantifies the reduction of the residual after an ADEIM update and establishes Proposition 7 that bounds the error $d ({¯ U}_{k + 1}, U_{k + 1})$ of the adapted space $U_{k + 1}$ with respect to the space ${¯ U}_{k + 1}$ .

Let $C_{k} = (P_{k}^{T} U_{k})^{- 1} P_{k}^{T} F_{k}$ be the coefficient matrix and $R_{k} = U_{k} C_{k} - F_{k}$ the residual matrix. Let further $S_{k}$ be the sampling points matrix and $U_{k + 1} = U_{k} + α_{k} β_{k}^{T}$ the adapted basis matrix of the adapted space $U_{k + 1}$ . Let $¯ r$ be the rank of $S_{k}^{T} R_{k}$ and let $r \in N$ with $r \leq ¯ r$ be the rank of the update $α_{k} β_{k}^{T}$ . Then,

{∥ ∥ S_{k}^{T} (U_{k + 1} C_{k} - F_{k}) ∥ ∥}_{F}^{2} = ∥ S_{k}^{T} R_{k} ∥_{F}^{2} - r \sum i = 1 σ_{i}^{2},

where $σ_{1} \geq σ_{2} \geq \dots \geq σ_{¯ r} > 0$ are the singular values of $S_{k}^{T} R_{k}$ .

Proof.

Follows from [39, Lemma 3.5] by transforming the generalized symmetric positive definite eigenproblem into the symmetric eigenproblem with matrix $(S_{k}^{T} R_{k})^{T} (S_{k}^{T} R_{k})$

, which then is equivalent to computing the singular value decomposition of

S_{k}^{T} R_{k}

so that the squared singular values of

S_{k}^{T} R_{k}

are equal to the eigenvalues of the generalized eigenproblem. This shows this lemma by using the last identity of

[39, proof of Lemma 3.5]. ∎

Consider the same setup as in Lemma 3.2.2 and additionally assume that $F_{k}$ has rank $n$ and its columns are in the $n$ -dimensional space ${¯ U}_{k + 1}$ . Let further $U_{k + 1}$ be the adapted space derived with the rank- $r$ ADEIM update and sampling points matrix $S_{k}$ . Then, it holds

d (U_{k + 1}, {¯ U}_{k + 1}) \leq \frac{ρ_{k}^{2}}{σ_{min}^{2} (F_{k})}

(6)

with

ρ_{k}^{2} = ∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2} + ¯ r \sum i = r + 1 σ_{i}^{2},

(7)

where ${˘ S}_{k}$ denotes the complementary sampling matrix of $S_{k}$ and $σ_{min} (F_{k})$ is the minimal singular value of $F_{k}$ .

Proof.

First, note that the matrix ${˘ S}_{k}^{T} α_{k} β_{k}^{T}$ has zero entries only, because the ADEIM update $α_{k} β_{k}^{T}$ updates only rows of $U_{k}$ corresponding to the sampling points $S_{k}$ , see [39, Lemma 3.5]. Thus, it holds

{˘ S}_{k}^{T} U_{k + 1} = {˘ S}_{k}^{T} U_{k} .

(8)

Consider now the norm of the residual with respect to the adapted space $U_{k + 1}$

$∥ U_{k + 1} C_{k} - F_{k} ∥_{F}^{2} =$	$∥ S_{k}^{T} (U_{k + 1} C_{k} - F_{k}) ∥_{F}^{2} + ∥ {˘ S}_{k}^{T} (U_{k + 1} C_{k} - F_{k}) ∥_{F}^{2}$	(9)
$=$	$∥ S_{k}^{T} (U_{k} C_{k} - F_{k}) ∥_{F}^{2} - r \sum i = 1 σ_{i}^{2} + ∥ {˘ S}_{k}^{T} (U_{k + 1} C_{k} - F_{k}) ∥_{F}^{2}$	(10)
$=$	$∥ S_{k}^{T} (U_{k} C_{k} - F_{k}) ∥_{F}^{2} - r \sum i = 1 σ_{i}^{2} + ∥ {˘ S}_{k}^{T} (U_{k} C_{k} - F_{k}) ∥_{F}^{2}$	(11)
$=$	$¯ r \sum i = r + 1 σ_{i}^{2} + ∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2} .$	(12)

Equation (9) follows because ${˘ S}_{k}$ is the complementary sampling points matrix of $S_{k}$ . Equation (10) follows from Lemma 3.2.2. Equation (11) holds because of (8). Equation (12) uses $R_{k} = U_{k} C_{k} - F_{k}$ and that $\sum_{i = 1}^{r} σ_{i}^{2} = ∥ S_{k}^{T} R_{k} ∥_{F}^{2}$ . Consider now the projection of the columns of $F_{k}$ onto $U_{k + 1}$ and observe that

∥ F_{k} - U_{k + 1} U_{k + 1}^{T} F_{k} ∥_{F}^{2} \leq ∥ U_{k + 1} C_{k} - F_{k} ∥_{F}^{2} .

(13)

Note that the basis matrix $U_{k + 1}$ obtained with the ADEIM update is not necessarily orthonormal and therefore it might be necessary to consider the oblique projection $U_{k + 1} (U_{k + 1}^{T} U_{k + 1})^{- 1} U_{k + 1}^{T} F_{k}$ onto the column span $U_{k + 1}$ of $U_{k + 1}$ in (13); however, the projection error in the Frobenius norm is the same for both projections. With (12) and the definition of $ρ_{k}^{2}$ follows

∥ F_{k} - U_{k + 1} U_{k + 1}^{T} F_{k} ∥_{F}^{2} \leq ρ_{k}^{2} .

(14)

Let now ${¯ U}_{k + 1}$ be an orthonormal basis matrix of ${¯ U}_{k + 1}$ . Since ${¯ U}_{k + 1}$ is spanned by the columns of $F_{k}$ , there exists a full-rank matrix ${~ F}_{k} \in R^{n \times w}$ such that $F_{k} = {¯ U}_{k + 1} {~ F}_{k}$ . We obtain

$∥ F_{k} - U_{k + 1} U_{k + 1}^{T} F_{k} ∥_{F}^{2} =$	$∥ ({¯ U}_{k + 1} - U_{k + 1} U_{k + 1}^{T} {¯ U}_{k + 1}) {~ F}_{k} ∥_{F}^{2}$	(15)
$\geq$	$∥ {¯ U}_{k + 1} - U_{k + 1} U_{k + 1}^{T} {¯ U}_{k + 1} ∥_{F}^{2} σ_{min}^{2} ({~ F}_{k})$
$=$	$∥ {¯ U}_{k + 1} - U_{k + 1} U_{k + 1}^{T} {¯ U}_{k + 1} ∥_{F}^{2} σ_{min}^{2} (F_{k}),$

where we used $σ_{min} ({~ F}_{k}) = σ_{min} (F_{k})$ , which holds because ${¯ U}_{k + 1}$ is orthonormal. The matrix $F_{k}$ has rank $n$ and therefore the smallest singular value $σ_{min} (F_{k}) > 0$ is positive. Combining (15) with (14) and the definition of $d ({¯ U}_{k + 1}, U_{k + 1})$ shows (6).

∎

3.3 Exploiting local coherence: Adaptive sampling

Proposition 7 shows that the choice of the sampling points $S_{k}$ influences the bound of the error $d ({¯ U}_{k + 1}, U_{k + 1})$ of the adapted space $U_{k + 1}$ . In this section, we derive an adaptive sampling strategy that minimizes the upper bound derived in Proposition 7 in case of full-rank ADEIM updates. Then, we show that if our adaptive sampling strategy is used, the error $d ({¯ U}_{k + 1}, U_{k + 1})$ decays at least as fast with the number of sampling points $m$ as the norm of the DEIM residual, which in turn means that only few sampling points are necessary to adapt the basis if the residual decays fast.

3.3.1 Adaptive sampling strategy based on residual

Following Proposition 7 and the decay factor $ρ_{k}$ defined in (7), we select the sampling points so that $∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2}$ is minimized. Note that $σ_{min} (F_{k})$ in (6) is independent of the sampling points $S_{k}$ and therefore it is sufficient to derive sampling points that lead to a small decay factor $ρ_{k}$ .

Consider the component-wise residual

r_{i} = ∥ R_{k}^{T} e_{i} ∥_{2}^{2}

(16)

and let $j_{1}, \dots, j_{N}$ be an ordering of $1, \dots, N$ such that $r_{j_{1}} \geq r_{j_{2}} \geq \dots \geq r_{j_{N}}$ . Then, select the first $j_{1}, \dots, j_{m}$ components as sampling points and form the corresponding sampling points matrix

S_{k} = [e_{j_{1}}, \dots, e_{j_{m}}] .

(17)

If a full-rank ADEIM update is applied, i.e., $¯ r = r$ in Proposition 7, then this choice of sampling points is optimal in the sense that the bound $ρ_{k}$ is minimized.

We now show that a fast decay in the residual implies a fast decay in the error bound of $d ({¯ U}_{k + 1}, U_{k + 1})$ .

Consider the same setup as in Proposition 7. Let $j_{1}, \dots, j_{N}$ be an ordering of ${1, \dots, N}$ such that the component-wise residual defined in (16) decays as

r_{j_{i}} \leq c_{1} e^{- c_{2} i^{a}},

(18)

with rate $a > 1$ and constants $c_{1}, c_{2} > 0$ . The error of the adapted space $U_{k + 1}$ is bounded as

d ({¯ U}_{k + 1}, U_{k + 1}) \leq c_{3} e^{- c_{2} m^{a}},

(19)

if a full-rank ADEIM update is applied and $m$ sampling points are selected via the adaptive sampling (17). The constant $c_{3} = c_{1} / ((1 - e^{- c_{2}}) σ_{min}^{2} (F_{k}))$ is independent of $m$ . In particular, setting the number of sampling points to

m \geq min {N, {(- \frac{1}{c_{2}} log (\frac{ϵ}{c_{3}}))}^{1 / a}}

(20)

guarantees $d ({¯ U}_{k + 1}, U_{k + 1}) \leq ϵ$ for a threshold $ϵ > 0$ .

Proof.

In case of a full-rank update, the factor $ρ_{k}^{2}$ in (6) is

ρ_{k}^{2} = ∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2},

as shown in Proposition 7. Then, we obtain with (18) and the adaptive sampling (17)

	$∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2} =$	$N \sum i = m + 1 ∥ R_{k}^{T} e_{j_{i}} ∥_{2}^{2} \leq c_{1} N \sum i = m + 1 e^{- c_{2} i^{a}} \leq c_{1} \infty \sum i = m + 1 e^{- c_{2} i^{a}}$
	$\leq$	$c_{1} \infty \sum i = 0 e^{- c_{2} (i + (m + 1))^{a}} \leq c_{1} e^{- c_{2} (m + 1)^{a}} \infty \sum i = 0 e^{- c_{2} i^{a}},$

where the last inequality holds because $i \geq 0, a > 1$ and thus $(i + (m + 1))^{a} \geq i^{a} + (m + 1)^{a}$ . Using that $a > 1$ , we obtain $i^{a} \geq i$ and $(m + 1)^{a} \geq m^{a}$ and thus

∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2} \leq c_{1} e^{- c_{2} (m + 1)^{a}} \infty \sum i = 0 e^{- c_{2} i^{a}} \leq c_{1} e^{- c_{2} m^{a}} \infty \sum i = 0 e^{- c_{2} i} = c_{1} e^{- c_{2} m^{a}} \frac{1}{1 - e^{- c_{2}}} .

(21)

Set $c_{3} = c_{1} / ((1 - e^{- c_{2}}) σ_{min}^{2} (F_{k}))$ to obtain (19) with Proposition 7. Setting $m$ as in (20) shows $d ({¯ U}_{k + 1}, U_{k + 1}) \leq ϵ$ . ∎

3.3.2 Local coherence

Following, e.g., [9], define the coherence of a space $U$ as

γ (U) = \frac{N}{n} max i = 1, \dots, N ∥ U^{T} e_{i} ∥_{2}^{2},

where $U$ is an orthonormal basis of $U$ and $e_{i}$ is the $i$ -th canonical unit vector. Define further the local coherence as

γ_{i} (U) = \frac{N}{n} ∥ U^{T} e_{i} ∥_{2}^{2},

for $i = 1, \dots, N$ , see, e.g., [28, 60]. We refer to [10, 9, 28, 60] for details.

We now show that the rate of the decay of the local coherence of the current DEIM space $U_{k}$ and of the space ${¯ U}_{k + 1}$ to which we want to adapt is reflected in the decay of the bound of the error $d ({¯ U}_{k + 1}, U_{k + 1})$ of the adapted space $U_{k + 1}$ with respect to the number of sampling points $m$ . Thus, we now relate the decay of the error $d ({¯ U}_{k + 1}, U_{k + 1})$ of the adapted space to a property of the spaces $U_{k}$ and ${¯ U}_{k + 1}$ , namely their decay of the local coherence.

Let $j_{1}, \dots, j_{N}$ be an ordering of ${1, \dots, N}$ such that

γ_{j_{i}} (U_{k}) \leq c_{4} exp (- c_{5} i^{a}), γ_{j_{i}} ({¯ U}_{k + 1}) \leq {¯ c}_{4} exp (- {¯ c}_{5} i^{¯ a})

(22)

holds for $i = 1, \dots, N$ , with $c_{4}, {¯ c}_{4}, c_{5}, {¯ c}_{5} > 0$ and $a, ¯ a > 1$ . Let the columns of $F_{k}$ be in ${¯ U}_{k + 1}$ . Define $C_{k} = (P_{k}^{T} U_{k})^{- 1} P_{k}^{T} F_{k}$ , then for the residual $R_{k} = U_{k} C_{k} - F_{k}$ holds

(23)

with a constant $c_{7} > 0$ that is independent of $a, ¯ a, c_{5}, {¯ c}_{5}$ .

Proof.

Denote with $U_{k}^{(i)}$ , ${¯ U}_{k + 1}^{(i)}$ , and $F_{k}^{(i)}$ the $i$ -th row of $U$ , ${¯ U}_{k + 1}$ , and $F_{k}$ , respectively. The matrix ${¯ U}_{k + 1}$ is an orthonormal basis matrix of ${¯ U}_{k + 1}$ . Let further ${~ F}_{k} \in R^{n \times w}$ be a matrix such that $F_{k} = {¯ U}_{k + 1} {~ F}_{k}$ . Note that $∥ F_{k} ∥_{2} = ∥ {~ F}_{k} ∥_{2}$ , because ${¯ U}_{k + 1}$ is orthonormal. Then, with $C_{k} = (P_{k}^{T} U_{k})^{- 1} P_{k}^{T} F_{k}$ , follows

	$∥ R_{k}^{T} e_{j_{i}} ∥_{2}^{2} =$	$∥ F_{k}^{(j_{i})} - U_{k}^{(j_{i})} C_{k} ∥_{2}^{2}$
	$=$	$∥ F_{k}^{(j_{i})} ∥_{2}^{2} - 2 F_{k}^{(j_{i})} (U_{k}^{(j_{i})} C_{k})^{T} + ∥ U_{k}^{(j_{i})} C_{k} ∥_{2}^{2}$
	$\leq$	$∥ F_{k}^{(j_{i})} ∥_{2}^{2} + 2 ∥ F_{k}^{(j_{i})} ∥_{2} ∥ U_{k}^{(j_{i})} ∥_{2} ∥ C_{k} ∥_{2} + ∥ U_{k}^{(j_{i})} ∥_{2}^{2} ∥ C_{k} ∥_{2}^{2} .$

Now make the following estimate

∥ F_{k}^{(j_{i})} ∥_{2}^{2} \leq ∥ {¯ U}_{k + 1}^{(j_{i})} ∥_{2}^{2} ∥ {~ F}_{k} ∥_{2}^{2} = ∥ {¯ U}_{k + 1}^{(j_{i})} ∥_{2}^{2} ∥ F_{k} ∥_{2}^{2},

because $∥ F_{k} ∥_{2} = ∥ {~ F}_{k} ∥_{2}$ . Further, we have

∥ C_{k} ∥_{2} \leq ∥ (P_{k}^{T} U_{k})^{- 1} ∥_{2} ∥ F_{k} ∥_{2} .

Set now ${¯ c}_{6} = n / N {¯ c}_{4}$ and $c_{6} = n / N c_{4} ∥ (P_{k}^{T} U_{k})^{- 1} ∥_{2}^{2}$ to bound

∥ R_{k}^{T} e_{j_{i}} ∥_{2}^{2} \leq ∥ F_{k} ∥_{2}^{2} ({¯ c}_{6} e^{- {¯ c}_{5} i^{a}} + 2 \sqrt{{¯ c}_{6} c_{6}} e^{- \frac{}{{¯ c}_{5}} 2 i^{¯ a} - \frac{c_{5}}{2} i^{a}} + c_{6} e^{- c_{5} i^{a}}) .

Now set $c_{7} = c_{6} + {¯ c}_{6} + 2 \sqrt{c_{6} {¯ c}_{6}}$ to obtain

∥ R_{k}^{T} e_{j_{i}} ∥_{2}^{2} \leq c_{7} ∥ F_{k} ∥_{2}^{2} e^{- min {c_{5}, {¯ c}_{5}} i^{min {a, ¯ a}}},

which shows the proposition. ∎

Since we consider subspaces of finite-dimensional spaces only, i.e., $N$ is finite, the bounds in (22) hold for any subspace by increasing the constants and choosing the rate close to 1. However, Lemma 23 still is meaningful because it shows that the constants and rates that appear in (22) are obtained in the bound of the residual in (23) as well. Thus, the local coherence structure is directly reflected in the decay of the DEIM residual.

The following proposition combines the local coherence structure exploited in Lemma 23 and the adaptive sampling of Proposition 20 to derive the number of sampling points $m$ that are required to achieve $d ({¯ U}_{k + 1}, U_{k + 1}) \leq ϵ$ in probability.

Assume that Lemma 23 applies and consider the same setting as in Proposition 7 except that $F_{k} = {¯ U}_{k + 1} {~ F}_{k}$ with ${~ F}_{k}$ being an $n \times n$ matrix with independent and identically distributed standard Gaussian entries. Then, $d ({¯ U}_{k + 1}, U_{k + 1}) \leq ϵ$ with probability at least $1 - δ$ if a full-rank update is applied and

m \geq min ⎧ ⎨ ⎩ N, {(- \frac{1}{min {c_{5}, {¯ c}_{5}}} log (\frac{ϵ δ^{2}}{c_{8}}))}^{1 / min {a, ¯ a}} ⎫ ⎬ ⎭,

(24)

where $c_{8}$ is a constant independent of $m, a$ , and $¯ a$ .

Proof.

Since Lemma 23 applies, we have with (21) that

∥ {˘ S}_{k}^{T} R_{k} ∥_{F}^{2} \leq \frac{c_{7} ∥ F_{k} ∥_{2}^{2}}{1 - e^{- min {c_{5}, {¯ c}_{5}}}} e^{- min {c_{5},}

Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling

Rank-adaptive structure-preserving reduced basis methods for Hamiltonian systems

An adaptive, training-free reduced-order model for convection-dominated problems based on hybrid snapshots

Manifold Approximations via Transported Subspaces: Model reduction for transport-dominated problems

Numerical Analysis of Automodel Solutions for Superdiffusive Transport

Mapping of coherent structures in parameterized flows by learning optimal transportation with Gaussian models

A Cost-Efficient Space-Time Adaptive Algorithm for Coupled Flow and Transport

Reduced models with nonlinear approximations of latent dynamics for model premixed flame problems

1 Introduction

2 Preliminaries

2.1 Model reduction with empirical interpolation

2.2 Problem formulation

3 Exploiting local structure via online adaptive basis updates

3.1 Local structure in transport-dominated problems

3.2 Exploiting local low-rank structure: Basis updates

3.2.1 Adaptation with ADEIM

3.2.2 Analysis of the error of the adapted space

Proof.

Proof.

3.3 Exploiting local coherence: Adaptive sampling

3.3.1 Adaptive sampling strategy based on residual

Proof.

3.3.2 Local coherence

Proof.

Proof.

Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling

Related Research

Rank-adaptive structure-preserving reduced basis methods for Hamiltonian systems

An adaptive, training-free reduced-order model for convection-dominated problems based on hybrid snapshots

Manifold Approximations via Transported Subspaces: Model reduction for transport-dominated problems

Numerical Analysis of Automodel Solutions for Superdiffusive Transport

Mapping of coherent structures in parameterized flows by learning optimal transportation with Gaussian models

A Cost-Efficient Space-Time Adaptive Algorithm for Coupled Flow and Transport

Reduced models with nonlinear approximations of latent dynamics for model premixed flame problems

1 Introduction

2 Preliminaries

2.1 Model reduction with empirical interpolation

2.2 Problem formulation

3 Exploiting local structure via online adaptive basis updates

3.1 Local structure in transport-dominated problems

3.2 Exploiting local low-rank structure: Basis updates

3.2.1 Adaptation with ADEIM

3.2.2 Analysis of the error of the adapted space

Proof.

Proof.

3.3 Exploiting local coherence: Adaptive sampling

3.3.1 Adaptive sampling strategy based on residual

Proof.

3.3.2 Local coherence

Proof.

Proof.