The Discrete Empirical Interpolation Method: Canonical Structure and Formulation in Weighted Inner Product Spaces

TL;DR

This paper advances the Discrete Empirical Interpolation Method (DEIM) by detailing its canonical structure, tightening error bounds, and extending it to weighted inner product spaces for better physical property preservation.

Contribution

It provides a canonical characterization of DEIM, improves error bounds using QR factorization, and introduces weighted DEIM ($W$-DEIM) for applications in energy inner products.

Findings

Tighter error bounds for DEIM oblique projection.

Extension of DEIM to weighted inner product spaces.

Introduction of DGEIM as a generalized interpolation method.

Abstract

New contributions are offered to the theory and practice of the Discrete Empirical Interpolation Method (DEIM). These include a detailed characterization of the canonical structure; a substantial tightening of the error bound for the DEIM oblique projection, based on index selection via a strong rank revealing QR factorization; and an extension of the DEIM approximation to weighted inner products defined by a real symmetric positive-definite matrix $W$. The weighted DEIM ($W$-DEIM) can be deployed in the more general framework where the POD Galerkin projection is formulated in a discretization of a suitable energy inner product such that the Galerkin projection preserves important physical properties such as e.g. stability. Also, a special case of $W$-DEIM is introduced, which is DGEIM, a discrete version of the Generalized Empirical Interpolation Method that allows generalization of the interpolation via a dictionary of linear functionals.

Equations

$$\dot{x}(t)=Ax(t)+\mathbf{f}(x(t)),x(0)=x_0\in {\mathbb{R}}^m,$$

$$\dot{\hat{x}}(t)=V_k^{T}A{V}_k\hat{x}(t)+V_k^{T}A\overline{x}+V_k^T\mathbf{f}(\overline{x}+V_k\hat{x}(t)),\hat{x}(0)=V_k^T(x(0)-\overline{x}),$$

$$\dot{\hat{x}}(t)=V_k^{T}A{V}_k\hat{x}(t)+V_k^{T}A\overline{x}+V_k^T{U}_r{U}_r^T\mathbf{f}\left(\overline{x}+V_k\hat{x}(t)\right)+V_k^T({\mathbb{I}}_m-U_r{U}_r^T)\mathbf{f}\left(\overline{x}+V_k\hat{x}(t)\right),$$

$$\mathcal{S}\equiv \left(\begin{array}{ccc}{\mathbb{I}}_m(:,i_1)& \cdots & {\mathbb{I}}_m(:,i_r)\end{array}\right)\in {\mathbb{R}}^{m\times r},$$

$$\mathbb{D}\equiv U_r({\mathcal{S}}^T{U}_r{)}^{-1}{\mathcal{S}}^T.$$

$$\dot{\hat{x}}(t)\approx V_k^{T}A{V}_k\hat{x}(t)+V_k^{T}A\overline{x}+V_k^T\mathbb{D}\mathbf{f}\left(\overline{x}+V_k\hat{x}(t)\right)$$

$${\mathcal{S}}^T\mathbf{f}(\overline{x}+V_k\hat{x}(t))=\left(\begin{array}{c}{\varphi }_{i_1}({\overline{x}}_{i_1}+V_k(i_1,:)\hat{x}(t))\\ {\varphi }_{i_2}({\overline{x}}_{i_2}+V_k(i_2,:)\hat{x}(t))\\ ⋮\\ {\varphi }_{i_r}({\overline{x}}_{i_r}+V_k(i_r,:)\hat{x}(t))\end{array}\right)\equiv {\mathbf{f}}_{\mathcal{S}}({\mathcal{S}}^T\overline{x}+({\mathcal{S}}^T{V}_k)\hat{x}(t)),t=t_1,t_2,\dots$$

$$V_k^T\mathbb{D}\mathbf{f}(\overline{x}+V_k\hat{x}(t))=(V_k^T{U}_r)({\mathcal{S}}^T{U}_r{)}^{-1}{\mathbf{f}}_{\mathcal{S}}({\mathcal{S}}^T\overline{x}+({\mathcal{S}}^T{V}_k)\hat{x}(t))$$

$$\Vert f-\mathbb{D}f{\Vert }_2\le \kappa \Vert ({\mathbb{I}}_m-U_r{U}_r^T)f{\Vert }_2,\text{where}\kappa \equiv \Vert ({\mathcal{S}}^T{U}_r{)}^{-1}{\Vert }_2.$$

$$F\equiv \left(\begin{array}{ccc}\mathbf{f}(x(t_1))& \cdots & \mathbf{f}(x(t_n))\end{array}\right)\in {\mathbb{R}}^{m\times n}$$

$$A\Pi =Q\left(\begin{array}{cc}{R}_{11}& R_{12}\\ 0& R_{22}\end{array}\right),$$

$${\sigma }_j(R_{11})$$

$${\sigma }_{r+j}(A)$$

$$\frac{{\sigma }_j(A)}{\sqrt{1+{\eta }^{2}r(n-r)}}$$

$${\sigma }_j(R_{22})$$

$$\left|{\left(R_{11}^{-1}{R}_{12}\right)}_{ij}\right|\le \eta ,1\le i\le r,1\le j\le n-r.$$

$$A\Pi =Q\left(\begin{array}{cc}{R}_{11}& R_{12}\end{array}\right),$$

$$\frac{{\sigma }_j(A)}{\sqrt{1+{\eta }^{2}m(n-m)}}$$

$$\frac{1}{\sqrt{1+{\eta }^{2}r(m-r)}}\le {\sigma }_j({\mathcal{S}}^{T}V)\le 1,1\le j\le r,$$

$$1\le \Vert ({\mathcal{S}}^{T}V{)}^{-1}{\Vert }_2\le \sqrt{1+{\eta }^{2}r(m-r)}.$$

$$V^T\left(\begin{array}{cc}{\Pi }_1& {\Pi }_2\end{array}\right)=Q\left(\begin{array}{cc}{R}_{11}& R_{12}\end{array}\right),$$

$$\frac{1}{\sqrt{1+{\eta }^{2}r(m-r)}}\le {\sigma }_j(R_{11})\le 1,1\le j\le r.$$

$$\kappa \le \sqrt{1+{\eta }^{2}r(m-r)}.$$

$$\mathbb{D}=U_r({\mathcal{S}}^T{U}_r{)}^{-1}{\mathcal{S}}^T=({\mathcal{P}}_{\mathcal{S}}{\mathcal{P}}_{U_r}{)}^{†},$$

$$({\mathcal{P}}_{\mathcal{S}}{\mathcal{P}}_{U_r}{)}^{†}=U_r({\mathcal{S}}^T{U}_r{)}^{†}{\mathcal{S}}^T.$$

$${\mathcal{S}}^T(\mathbb{D}f)={\mathcal{S}}^T{U}_r({\mathcal{S}}^T{U}_r{)}^{†}{\mathcal{S}}^{T}f={\mathcal{S}}^{T}f.$$

$${\mathcal{S}}^T(\mathbb{D}f)={\mathcal{S}}^T{U}_r({\mathcal{S}}^T{U}_r{)}^{†}{\mathcal{S}}^{T}f={\mathcal{P}}_{\mathcal{X}}({\mathcal{S}}^{T}f),\mathcal{X}=\mathcal{R}({\mathcal{S}}^T{U}_r).$$

$$1={\sigma }_1=\dots ={\sigma }_{\ell }>{\sigma }_{\ell +1}\ge \dots \ge {\sigma }_{\ell +p}>{\sigma }_{\ell +p+1}=\dots ={\sigma }_{\min (r,s)}=0.$$

$$\mathbb{D}=({\mathcal{P}}_{\mathcal{S}}{\mathcal{P}}_{U_r}{)}^{†}=Z\left(\begin{array}{ccc}{\mathbb{I}}_{\ell }& & \\ & \underset{i=1}{\overset{p}{⨁}}{T}_i& \\ & & 0\end{array}\right)Z^T,T_i=\left(\begin{array}{cc}1& 0\\ \tan {\psi }_{\ell +i}& 0\end{array}\right).$$

$${\mathcal{P}}_{\mathcal{S}}$$

$${\mathcal{P}}_{U_r}$$

$$({\mathcal{P}}_{\mathcal{S}}{\mathcal{P}}_{U_r}{)}^{†}=Z\left(\begin{array}{ccc}{\mathbb{I}}_{\ell }& & \\ & \underset{i=1}{\overset{p}{⨁}}(J_i{\Psi }_i{)}^{†}& \\ & & 0\end{array}\right)Z^T.$$

Full text

1 Introduction

Suppose we want to run numerical simulations of a physical reality described by a set of ordinary differential equations (ODEs)

[TABLE]

where $A\in\mathbb{R}^{m\times m}$, $\mathbf{f}:\mathbb{R}^{m}\longrightarrow\mathbb{R}^{m}$. Often, of interest is $y(t)=Cx(t)$, with some given $p\times m$ matrix $C$. Such a system of ODEs can arise from discretization of a spatial differential operator in time dependent PDEs (e.g. method of lines), e.g. for the purposes of prediction and/or control, or optimization with respect to a set of parameters. In a parameter dependent case we have $A=A(\mu)$, $x=x(t;\mu)$, $x_{0}=x_{0}(\mu)$, and $\mathbf{f}(\cdot;\mu)$ is also parameter dependent, where the parameter $\mu$, that may carry e.g. information on material properties, is from a parameter domain $\mathcal{P}\subset\mathbb{R}^{d}$, $d\geq 1$.111To keep the notation simple, we suppress the explicit parameter dependence until numerical experiments in §5. The function $\mathbf{f}(\cdot;\cdot)$ is in general assumed to be nonlinear. Large dimension $m$ (say, $m>10^{5}$) makes the task computationally intractable for multiple query problems, and one is forced to devise and use a reduced order system that emulates (1).

In a projection based model order reduction, one constructs a suitable low dimensional subspace $\mathcal{V}_{k}$ as the range of an $m\times k$ orthonormal matrix $V_{k}$ ($V_{k}^{T}V_{k}=\mathbb{I}_{k}$) and seeks an approximation of the form $x(t)\approx\overline{x}+V_{k}\widehat{x}(t)$, $\widehat{x}\in\mathbb{R}^{k}$. The solution $x(t)$ is stored at a set of discrete times (also known as snapshots) and $\overline{x}$ is the average over the snapshots. The matrix $V_{k}$ can be, e.g., the POD basis of the $k$ leading left singular vectors of the centered snapshots $x(t_{i})-\overline{x}$, computed at the discrete times $t_{i}$ from high resolution numerical simulations in the off-line phase; possibly over a parameter grid. It is assumed that $k\ll m$. By enforcing the orthogonality of the residual and the space $\mathcal{V}_{k}$, one obtains Galerkin projection of the original problem

[TABLE]

where $A_{k}=V_{k}^{T}AV_{k}$ is $k\times k$, $V_{k}^{T}A\overline{x}\in\mathbb{R}^{k}$, but the projected nonlinear forcing term $V_{k}^{T}\mathbf{f}(\overline{x}+V_{k}\widehat{x}(t))$ still involves the dimension $m$, in computing $\widetilde{x}(t)=\overline{x}+V_{k}\widehat{x}(t)$ and $f=\mathbf{f}(\widetilde{x}(t))$ (at a sequence of discrete values $t=t_{i}$), as well as in computing $V_{k}^{T}f$. For large $m$, this carries substantial computational effort and heavy memory traffic.

The Discrete Empirical Interpolation (DEIM) [15] method provides a way to alleviate these burdens, and to efficiently approximate $\mathbf{f}(\cdot)$ from a learned subspace. DEIM originates in the Empirical Interpolation Method (EIM) [24], [5], [42] and it uses the reduced basis provided by the POD. For a related discrete version of EIM see [26]. Here, for the reader’s convenience, we first briefly review the main steps of DEIM approximation and its error estimate, and then we place it in the more general concepts of GEIM and PBDW.

1.1 DEIM

Suppose we have empirically determined an $r$-dimensional subspace $\mathcal{U}_{r}$ as the range of an orthonormal $U_{r}$ such that $U_{r}U_{r}^{T}\mathbf{f}(\overline{x}+V_{k}\widehat{x}(t))\approx\mathbf{f}(\overline{x}+V_{k}\widehat{x}(t))$. This can be done e.g. by the POD, which will determine a suitable dimension $r$ from the decay of the singular values of the matrix of snapshots. The tacit assumption is that $\mathbf{f}(\cdot)$ is from a set of functions $\mathcal{F}$ with small Kolmogorov $r$-width [35], [36, Chapter 6], [42]. Inserting the orthogonal projection $U_{r}U_{r}^{T}$ into (2) gives

[TABLE]

where $\mathbb{I}_{m}\in\mathbb{R}^{m\times m}$ denotes the identity, and the last term (the POD error, as seen from $\mathcal{V}_{k}$) can be neglected. However, this still does not solve the problem of computational complexity because it requires all $m$ components of $\mathbf{f}\left(\overline{x}+V_{k}\widehat{x}(t)\right)$, and the matrix vector product $U_{r}^{T}\mathbf{f}\left(\overline{x}+V_{k}\widehat{x}(t)\right)$ takes $\mathcal{O}(mr)$ operations for every time point $t=t_{i}$.

The DEIM [15] trick is to select a submatrix of the $m\times m$ identity $\mathbb{I}_{m}$,

[TABLE]

and to replace the orthogonal projector $U_{r}U_{r}^{T}$ by the oblique projector

[TABLE]

Note that $\mathbb{D}$ has an interpolating property at the $r$ selected coordinates, $\mathscr{S}^{T}\mathbb{D}f=\mathscr{S}^{T}f$. The alternative for (3) is thus

[TABLE]

where in the matrix product $V_{k}^{T}\mathbb{D}$, the factor $V_{k}^{T}U_{r}(\mathscr{S}^{T}U_{r})^{-1}$ can be pre-computed in the off-line phase. Obviously, important is only the component of the error $(\mathbb{I}_{m}-\mathbb{D})\mathbf{f}(\overline{x}+V_{k}\widehat{x}(t))$ that lies in $\mathcal{V}_{k}$.

The on-line computation $\mathscr{S}^{T}\mathbf{f}(\overline{x}+V_{k}\widehat{x}(t))$ at any particular $t$ involves only $r$ values $\mathbf{f}_{i_{j}}(\overline{x}+V_{k}\widehat{x}(t))$, $j=1,\ldots,r$. If $\mathbf{f}$ is defined at a vector $x=(x_{i})_{i=1}^{m}$ component-wise as222For a general nonlinear $\mathbf{f}(x)=(\varphi_{1}(x_{\mathcal{I}_{1}}),\varphi_{2}(x_{\mathcal{I}_{2}}),\ldots,\varphi_{m}(x_{\mathcal{I}_{m}}))^{T}$, where $x_{\mathcal{I}_{j}}$ ($\mathcal{I}_{j}\subseteq\{1,\ldots,m\}$) denotes a sub-array of $x$ needed to evaluate $\varphi_{j}(x)$, the situation is more complicated, see [15, §3.5]. $\mathbf{f}(x)=(\phi_{1}(x_{1}),\phi_{2}(x_{2}),\ldots,\phi_{m}(x_{m}))^{T}$ then

[TABLE]

and the computational complexity of

[TABLE]

becomes independent of the dimension $m$, once the time independent matrices are precomputed in the off-line phase.333In the sequel, for the sake of simplicity, we do not include centering of the snapshots. This tremendously reduces both the flop count and the memory traffic in the (on-line) simulation.

The error of the DEIM oblique projection can be bounded in the Euclidean norm by that of the orthogonal projector,

[TABLE]

The condition number $\kappa$ determines the quality of the approximation, and satisfies $\kappa\leq\mathcal{O}\left(m^{(r-1)/2}\right)/\|U_{r}(:,1)\|_{\infty}$ [15]. In practical situations, however, this bound is pessimistic and $\kappa$ is much lower. (Using the concept of maximal volume [34], [23], it can be shown that there exists a strategy such that $\kappa\leq\sqrt{1+r(m-r)}$.)

1.1.1 Variations and generalizations

DEIM has been successfully deployed in many applications, and tuned for better performance, giving rise to the localized DEIM [48], unassembled DEIM (UDEIM) [60], [59], matrix DEIM [66], [46], nonnegative DEIM (NNDEIM) [1], and Q-DEIM [19]. The latter is an orthogonal variant of DEIM, which can be efficiently implemented with high-performance libraries such as LAPACK [2] and ScaLAPACK [9]. Furthermore, Q-DEIM admits a better condition number bound, $\kappa\leq\sqrt{m-r+1}\mathcal{O}(2^{r})$; it allows randomized sampling; and it can work with only a subset of the rows of $U_{r}$ for computing selection matrices $\mathscr{S}$ while keeping $\kappa$ moderate.

1.2 GEIM and PBDW

In many applications, the functions’ values may not be available through point evaluation because, e.g., they are from a class that does not contain continuous functions, there is no analytical expression, or they may be noisy sensor data (measurements) obtained by weighted averaging. In those cases, point-wise interpolation may not be possible, nor even desirable – for a most illuminating discussion see [40]. This motivated the development of a generalization of EIM, GEIM (Generalized Empirical Interpolation Method), which replaces point interpolation by more general evaluation functionals selected from a dictionary; see [45, Chapters 4, 5] and [39], [40], [41]

These ideas have been further extended in the Parametrized-Background Data-Weak approach to data assimilation (PBDW) [43]. PBDW is an elaborate data assimilation scheme whose weak formulation naturally fits variational framework for (parametrized) PDEs, and facilitates error estimates (both a priori and a posteriori) with the capability to identify optimal observation functionals. Additional insights and analysis of PBDW with respect to noise in the data, and an updating strategy for many-query scenarios, are provided in [44]; the optimality of the approximation is established in [6]. Furthermore, [6] contains a multi-space extension.

In the context of empirical interpolation, PBDW allows more (generalized) approximation positions than the cardinality of the POD basis, thus calling for least squares approximation. In particular, it contains GEIM as a special (one-space) case.

1.3 Proper inner-product space structure

In applications in engineering and applied sciences the solution of (1) represents an approximation of a function from an appropriate function space, that is subject to governing equations that describe a physical reality. The quality of an approximation is then naturally measured in an appropriate (interpretable) metric in that space.

For instance, in many applications the natural ambient space is (weighted) $L^{2}(\Omega)$, $L^{2}(\Omega)=\{f:\Omega\longrightarrow\mathbb{R}\;:\;\int_{\Omega}|f(x)|^{2}\rho(x)dx<\infty\},$ with the Hilbert space structure generated by the inner product $(f,g)_{L^{2}(\Omega)}=\int_{\Omega}f(x)g(x)\rho(x)dx,$ and with the corresponding induced norm $\|f\|_{L^{2}(\Omega)}=\sqrt{(f,f)_{L^{2}(\Omega)}}$. Both the weight function $\rho(\cdot)$ and a quadrature formula in the course of constructing a discrete (finite $m$-dimensional) framework yield a weighted inner product in $\mathbb{R}^{m}$, $(u,v)_{W}=v^{T}Wu$, where $W$ is the corresponding symmetric positive definite matrix. Then the natural framework for devising e.g. a POD approximation [63, §1.2] is given by the Hilbert space structure of $(\cdot,\cdot)_{W}$. Further, for the equations of e.g. compressible fluid flow, Galerkin projection in an $(\cdot,\cdot)_{L^{2}(\Omega)}$ inner product may not preserve the underlying physics, such as energy conservation or stability, see e.g. [53], [27, §3.4.3]. Different inner products (with corresponding norms) may yield substantially different results, see e.g. [21], [58], [31]. In model order reduction, for instance, a Galerkin projection may be naturally defined in a Lyapunov inner product, generated by the positive definite solution $W$ of a Lyapunov matrix equation, see e.g. [32, §6.1], [54], [52], [27, §5.4.3]. For further examples and in-depth discussion see [4] [33], [12], [69, 68], [51], [47], [20]. It should be clear that the use of a weighted inner product in the POD-DEIM framework does not guarantee the stability of the reduced system, unless the DEIM is additionally adapted to a particular structure. An excellent example of energy stable DEIM approximation is the NNDEIM [1].

The use of a proper inner product is implicitly assumed in the abstract framework of PBDW, including the special case of GEIM. The resulting numerical realization of the proper inner product results in the discrete $(\cdot,\cdot)_{W}$ inner product. From the numerical point of view, this is not a mere change to another inner product, as the condition number of $W$ becomes an important factor both in the theoretical projection error bound and in the computation in finite precision arithmetic. Hence, it seems natural and important to revise the numerical implementation of DEIM oblique projection, to place it in the wider context of PBDW, and to ensure its robustness independent of the possibly high condition number of the weight matrix $W$.

1.4 Scaling of variables

We discuss difficulties due to scaling issues in the practical computation of a POD basis and construction of a DEIM projection, and argue that, when appropriate, the DEIM projection must be weighted in a consistent manner with the POD basis.

Scaling issues discussed here arise from two sources. First, when unknowns $x_{i}(t)$ represent different physical quantities, such as velocity and pressure, and the numerical values of one of them, say pressure, can dominate all others by several orders of magnitude. Second, a single variable can vary over a wide range. In both scenarios, the components of $\mathbf{f}(x(t))$ may vary over several orders of magnitude, so that the matrix of nonlinear snapshots

[TABLE]

has graded rows, with widely varying norms.

Let us try to understand the computational ramifications. Suppose the rows of $F=\left(\begin{smallmatrix}B\cr s\end{smallmatrix}\right)$ are permuted so that $B$ contains the rows with large norm, $s$ the rows with small norm, so that in the Frobenius norm $\|B\|_{F}\gg\|s\|_{F}$. Typically $m\gg n$, and let the thin SVD be $F=U\Sigma K^{T}$, where $U$ is $m\times n$ orthonormal, $\Sigma$ is diagonal and $K$ is an orthogonal matrix. An economical way to compute $U_{r}$, often used in practice, is to first compute the eigenvalue decomposition $G\equiv F^{T}F=K\Sigma^{2}K^{T}$. Then choose a suitable $r$, compute $U_{r}=FK(:,1:r)\Sigma(1:r,1:r)^{-1}$, and apply a Gram-Schmidt correction to improve the numerical orthonormality of $U_{r}$. Since $F^{T}F=B^{T}B+s^{T}s\approx B^{T}B$, in this procedure the contribution of $s$ to the computation of $K$ and $U_{r}$ is marginal and the subdominant variables are almost invisible.

Further, the POD basis may inherit the graded structure of $F$. Assume, for the purpose of demonstration, that the dominant $r$ singular values of $F$ are nearly equal and much larger than the remaining, subdominant, singular values. Rearranging the thin SVD $F=U\Sigma K^{T}$, where $K$ is an orthogonal matrix, shows that the row norms of $U\Sigma=FK$ are the same as those of $F$. Furthermore, the row norms of the matrix of the leading $r$ singular vectors $U_{r}$ are distributed like the corresponding row norms of $F$. The indices corresponding to dominant variables have dominant rows in $U_{r}$, which creates difficulties for the representation of subdominant variables.

Moreover, the DEIM [15] and Q-DEIM [19] are based on greedy algorithms that try to identify an $r\times r$ submatrix of $U_{r}$ of maximal volume, thus preferring row indices corresponding to dominant variables and ignoring the others. The resulting small approximation error in the Euclidean norm is misleading, though. Without prior scaling, relevant and informative subdominant variables are unnecessarily suppressed.444Recall the discussion in §1.3.

Finally, it should be pointed out that strongly graded $F$ poses intrinsic computational difficulties for any algorithm for computing $U_{r}$. Even the backward error $\delta F$ that corresponds to the numerical computation, and which is small in the sense that $\|\delta F\|_{F}/\|F\|_{F}$ is small, may wipe out the information on the subdominant variables. The corresponding entries of the left singular vectors $u_{k}$, $k=1,\ldots,r$, are computed possibly with large relative error, as numerical methods in general compute the singular vectors with error such that $\|\delta u_{k}\|_{2}$ is appropriately bounded by the machine roundoff times a condition number [64], [56, V.4]. Tiny components of $u_{k}$ are usually computed with large relative error.

1.5 Contributions and overview of the paper

Our contributions to the theory and practice of empirical interpolation methods in the framework of PBDW approximations (in particular, EIM and GEIM) are towards numerical linear algebra and matrix theory; the goal is to setup a more general algorithmic schemes and principles for development of numerical methods with sharp error bounds, and for their successful software implementations and applications in scientific computing.

In Section 2, we present a substantial improvement of the bound on the condition number $\kappa$ in (5). The selection operator $\mathscr{S}$ is based on local maximal volume approach [34, 23], implemented via a strong rank-revealing QR decomposition [25]; the resulting DEIM condition number is $\kappa\leq\sqrt{1+\eta^{2}\,r(m-r)}$, with tunable parameter $\eta\geq 1$. In §3, we present a canonical form for the DEIM projector, which is based on the well known structure of oblique projections. It provides better understanding of the structure of DEIM and its approximation error. In §4 we introduce and give detailed analysis of the weighted DEIM ($W$-DEIM) which naturally applies in the situations discussed in §1.3, §1.4. The goal is to establish a universal framework for DEIM projections in weighted inner product spaces, where inner products are induced by positive definite matrices $W$ of various origins and with various interpretations.

We present several algorithms for computing the $W$-DEIM approximation. The algorithms come in two flavors depending on whether generalized or pointwise interpolation is desired. When generalized interpolation is to be used, we present different algorithms depending on whether $W$ is dense or sparse. When pointwise interpolation is used, our analysis shows that the condition number of $W$ plays a role in the error analysis. To mitigate this issue, we present an algorithm for which the spectral condition number $\kappa_{2}(W)=\|W\|_{2}\|W^{-1}\|_{2}$ enters the error bound, up to the factor of $\sqrt{m}$, as $\sqrt{\min_{D=\mathrm{diag}}\kappa_{2}(DWD)}$.

$W$-DEIM can be considered as a numerical realization of disretized (one-space) PBDW that includes, as a special case, a discrete version of the generalized empirical interpolation method (GEIM). In §5, we corroborate the results with numerical examples.

2 Nearly optimal subset selection

We review strong rank revealing QR methods for matrices with at least as many rows as columns (§2.1); and present an extension to matrices with fewer rows than columns and apply it to matrices with orthonormal rows (§2.2). This yields a new DEIM selection with superior error bound.

2.1 Tall and skinny matrices

For ${A}\in\mathbb{R}^{m\times n}$ with $m\geq n$, and a target rank $r<n$, a QR factorization with column pivoting computes

[TABLE]

where ${\Pi}\in\mathbb{R}^{n\times n}$ is a permutation; ${Q}\in\mathbb{R}^{m\times m}$ is an orthogonal matrix; ${R}_{11}\in\mathbb{R}^{r\times r}$ and ${R}_{22}\in\mathbb{R}^{(m-r)\times(n-r)}$ are upper triangular; and ${R}_{12}\in\mathbb{R}^{r\times(n-r)}$.

Let $\sigma_{1}(A)\geq\cdots\geq\sigma_{n}(A)\geq 0$ be the singular values of $A$. Singular value interlacing [22, Corollary 8.6.2] implies for the non-increasingly ordered singular values $\sigma_{j}({R}_{11})$ and $\sigma_{j}({R}_{22})$ of the diagonal blocks ${R}_{11}$ and ${R}_{22}$, respectively,

[TABLE]

So-called rank-revealing QR (RRQR) factorizations [14, 25] try to make the singular values of ${R}_{11}$ as large as possible, and those of ${R}_{22}$ as small as possible. In particular, the strong RRQR (sRRQR) factorization [25, Algorithm 4] with tuning parameter $\eta\geq 1$ computes a triangular matrix $R$ whose diagonal blocks have singular values within, essentially, a polynomial factor (in $n$ and $r$) of the singular values of ${A}$,

[TABLE]

and whose off-diagonal block is bounded by

[TABLE]

For $\eta>1$ the sRRQR factorization can be computed in $\mathcal{O}\left((m+n\log_{\eta}{n})n^{2}\right)$ arithmetic operations [25, Section 4.4]. Recommended values for $\eta$ are small fractional powers of $n$ [25, Section 4.4], such as $\eta=10\sqrt{n}$ [25, Section 6], which result in a $\mathcal{O}(mn^{2})$ time complexity.

The traditional Businger-Golub QR with column pivoting [11], [22, Algorithm 5.4.1] often achieves the above bounds in practice, but fails spectacularly on contrived examples such as the Kahan matrix [30], [25, Section 6]. Sometimes, the failure is caused by the software implementation of a RRQR factorization; for details see [18].

2.2 Short and fat matrices

The sRRQR factorization can be adapted to matrices with fewer rows than columns, $m<n$, and of full row rank, to select a well-conditioned $m\times m$ submatrix.

For ${A}\in\mathbb{R}^{m\times n}$ with $m\leq n$, and target rank $r=m$, a QR factorization with column pivoting computes

[TABLE]

where ${\Pi}\in\mathbb{R}^{n\times n}$ is a permutation matrix, ${Q}\in\mathbb{R}^{m\times m}$ is an orthogonal matrix, ${R}_{11}\in\mathbb{R}^{m\times m}$ is upper triangular, and ${R}_{12}\in\mathbb{R}^{m\times(n-m)}$.

A sRRQR factorization, in particular, is computed with a simplified version of [25, Algorithm 4]. A column of $R_{11}$ is swapped with one in $R_{12}$ until $\left|\left({R}_{11}^{-1}{R}_{12}\right)_{i,j}\right|\leq\eta$, $1\leq i\leq m,1\leq j\leq n-m$. From [10, Lemma 3.1] follows

[TABLE]

Given a matrix $V$ with $r$ orthonormal columns, this algorithm can be used to select a well conditioned $r\times r$ submatrix.

Lemma 1**.**

Let ${V}\in\mathbb{R}^{m\times r}$ with ${V}^{T}{V}={I}_{r}$. Applying [25, Algorithm 4] with target rank $r$ and tuning parameter $\eta\geq 1$ to $V^{T}$ gives a submatrix $\mathscr{S}\in\mathbb{R}^{m\times r}$ of $\mathbb{I}_{m}$ with

[TABLE]

and

[TABLE]

Proof.

Applying [25, Algorithm 4] to $V^{T}$ gives

[TABLE]

where $Q\in\mathbb{R}^{r\times r}$ is an orthogonal matrix; $R_{11}\in\mathbb{R}^{r\times r}$ is upper triangular; and $\begin{pmatrix}\Pi_{1}&\Pi_{2}\end{pmatrix}\in\mathbb{R}^{m\times m}$ is a permutation matrix with $\Pi_{1}\in\mathbb{R}^{m\times r}$.

Since ${V}$ has $r$ orthonormal columns, $\sigma_{j}(V)=1$, $1\leq j\leq r$. From (6) follows

[TABLE]

Set $\mathscr{S}=\Pi_{1}$, so the first block column equals ${V}^{T}\mathscr{S}={V}^{T}{\Pi}_{1}={Q}{R}_{11}$. Since ${Q}$ is an orthogonal matrix, ${V}^{T}\mathscr{S}$ has the same singular values as ${R}_{11}$. ∎

Lemma 1, applied with $V=U_{r}$, implies a tremendous improvement for the error of the oblique projector $\mathbb{D}$ in (5). If $\mathscr{S}$ is computed from a sRRQR factorization of the transposed POD basis $U_{r}^{T}$, the condition number is bounded by

[TABLE]

3 Canonical representation of $\mathbb{D}$

The DEIM operator is an oblique projection and, as such, it possesses certain canonical structure that is revealed in an appropriately chosen basis. In this section we derive representation of the DEIM projection operator in a particular basis, in order to gain better understanding of the effectiveness of DEIM. As already mentioned in §1.2, the PBDW framework [43] allows selecting $s\geq r$ approximation points, and we will proceed with the general case of rectangular $\mathscr{S}^{T}U_{r}$. The oversampling has been successfully used in the related context of missing point estimation, see [3] [49], [70], [71].

We adopt the following notation. Let $\mathscr{S}\in\mathbb{R}^{m\times s}$ be a selection of $s$ columns of the identity $\mathbb{I}_{m}$ and let $U_{r}\in\mathbb{R}^{m\times r}$. Define the orthogonal projectors $\mathcal{P}_{\mathscr{S}}=\mathscr{S}\mathscr{S}^{T}$ and $\mathcal{P}_{U_{r}}=U_{r}U_{r}^{T}$ onto $\mathcal{R}(\mathscr{S})$ and $\mathcal{R}(U_{r})$, respectively.

3.1 Generalization of oblique DEIM

We first derive a representation of the DEIM projection in terms of $\mathcal{P}_{\mathscr{S}}$ and $\mathcal{P}_{U_{r}}$. Suppose $\mathscr{S}$ and $U_{r}$ have full column rank, then the DEIM projector $\mathbb{D}$ can be written as [7, Theorem 2.2.3]

[TABLE]

where the superscript $\dagger$ denotes the Moore-Penrose inverse. Note that the expression $(\mathcal{P}_{\mathscr{S}}\mathcal{P}_{U_{r}})^{\dagger}$ does not require existence of the inverse $(\mathscr{S}^{T}U_{r})^{-1}$; in fact it does not even require $\mathscr{S}$ and $U_{r}$ to have the same number of columns, or the same rank.

We now consider the case that $\mathscr{S}^{T}U_{r}\in\mathbb{R}^{s\times r}$ is a rectangular matrix where $s\neq r$. In this case, one can check (e.g., using the SVD of $\mathscr{S}^{T}U_{r}$) that it holds

[TABLE]

This observation leads to a general definition of the DEIM projection as $\mathbb{D}=(\mathcal{P}_{\mathscr{S}}\mathcal{P}_{U_{r}})^{\dagger}$ which is valid when $\mathscr{S}$ has different number of columns as $U_{r}$, and different rank. We now investigate whether this generalization retains the properties of interpolation ($\mathscr{S}^{T}\mathbb{D}f=\mathscr{S}^{T}f$) and projection ( $\mathbb{D}\mathcal{P}_{U_{r}}=\mathcal{P}_{U_{r}}$).

With the observation $\mathop{\operator@font rank}\nolimits(\mathbb{D})=\mathop{\operator@font rank}\nolimits(\mathscr{S}^{T}U_{r})=\min\{s,r\}$, suppose that $s\neq r$ and split the analysis into two cases.

Case $\mathop{\operator@font rank}\nolimits(\mathbb{D})=s<r$

1. (a)

The interpolation property still holds, i.e.,

[TABLE]

The reason for this is because $\mathscr{S}^{T}U_{r}$ has full row rank, and $(\mathscr{S}^{T}U_{r})^{\dagger}$ is a right multiplicative inverse. 2. (b)

On the other hand, the projection property is lost, i.e., $\mathbb{D}\mathcal{P}_{U_{r}}\neq\mathcal{P}_{U_{r}}$. However, $\mathbb{D}$ is still a projector, $\mathbb{D}^{2}=\mathbb{D}$. To find its range, let $W_{s}$ be the leading right $s$ singular vectors of $\mathscr{S}^{T}U_{r}$. Then the DEIM projection operator. $\mathbb{D}\mathcal{P}_{U_{r}}=\mathcal{P}_{V_{s}}$, where $V_{s}=U_{r}W_{s}$ spans an $s$–dimensional subspace of $\mathcal{R}(U_{r})$. Therefore, $\mathbb{D}$ is a projector onto $\mathcal{R}(V_{s})\subset\mathcal{R}(U_{r})$. 2. 2.

Case $\mathop{\operator@font rank}\nolimits(\mathbb{D})=r<s$

1. (a)

The interpolation property does not hold, i.e., $\mathscr{S}^{T}(\mathbb{D}f)\neq\mathscr{S}^{T}f$. This is because $(\mathscr{S}^{T}U_{r})^{\dagger}$ is no longer a right multiplicative inverse. However, $\mathscr{S}^{T}(\mathbb{D}f)$ is the least square projection of $\mathscr{S}^{T}f$ onto the range of $\mathscr{S}^{T}U_{r}$. To see this

[TABLE] 2. (b)

In this case $\mathbb{D}\mathcal{P}_{U_{r}}=\mathcal{P}_{U_{r}}$ since $(\mathscr{S}^{T}U_{r})^{\dagger}$ is a left multiplicative inverse of $\mathscr{S}^{T}U_{r}$.

As can be seen above, when the DEIM operator is generalized to the setting $s\neq r$ only the projection property or the interpolation property is retained but not both simultaneously. For related developments, see  [43], [71], [13].

3.2 Canonical structure of $\mathbb{D}$

We present the following theorem that sheds light onto the canonical structure of the DEIM operator $\mathbb{D}$.

Theorem 2**.**

Let $U_{r}\in\mathbb{R}^{m\times r}$ and $\mathscr{S}\in\mathbb{R}^{m\times s}$ have orthonormal columns, and $\mathbb{D}=U_{r}(\mathscr{S}^{T}U_{r})^{\dagger}\mathscr{S}^{T}$ and assume that $1\leq r,s\leq m$. Let $\ell\equiv\mathrm{dim}(\mathcal{R}(\mathscr{S})\bigcap\mathcal{R}(U_{r}))$, set $p\equiv\mathrm{rank}(\mathbb{D})-\ell$, and let the singular values $\sigma_{i}=\cos\psi_{i}$ of $\mathscr{S}^{T}U_{r}$ be ordered as

[TABLE]

(Here $0<\psi_{\ell+1}\leq\ldots\leq\psi_{\ell+p}<\pi/2$ are the acute principal angles between the ranges of $\mathscr{S}$ and $U_{r}$.)

(i)* There exists an orthogonal $m\times m$ matrix $Z$ such that the matrix $\mathbb{D}$ can be represented as*

[TABLE]

Here the $\mathbf{0}$ block is of size $m-\ell-2p$.

(ii)* The DEIM projector $\mathbb{D}$ satisfies $\|\mathbb{D}\|_{2}=1/\cos\psi_{\ell+p}$. If, in addition, $\mathbb{D}\neq\mathbf{0}$ and $\mathbb{D}\neq\mathbb{I}_{m}$, then $\|\mathbb{D}\|_{2}=\|\mathbb{I}_{m}-\mathbb{D}\|_{2}=1/\cos\psi_{\ell+p}$.*

Proof.

The above representation follows immediately from the canonical representation of a pair of orthogonal projectors [65]. In a particularly constructed orthonormal basis given by the columns of $Z$, the two projectors have the following matrix representations:

[TABLE]

with $\psi_{\ell+i}$’s as stated in the theorem, and $D_{s}$, $D_{u}$ are diagonal matrices with diagonal entries [math] or $1$ and such that $D_{s}D_{u}=\mathbf{0}$. Note that each $(D_{s})_{ii}=1$ ($(D_{u})_{ii}=1$) corresponds to a direction in the range of $\mathscr{S}$ ($U_{r}$) orthogonal to the entire range of $U_{r}$ ($\mathscr{S}$). In the special case when $\mathscr{S}^{T}U_{r}$ is invertible, $D_{s}=D_{u}=\mathbf{0}$.

The expression for $\mathbb{D}$ is obtained by multiplying the representations in (12) and (13), and taking the pseudoinverse. It follows that

[TABLE]

A direct evaluation shows that

[TABLE]

From the canonical representation (11) each block $T_{i}$ has the norm

[TABLE]

Therefore, it also follows that $\|\mathbb{D}\|_{2}=1/\cos\psi_{\ell+p}$. From (11) we can also derive the canonical form of $\mathbb{I}_{m}-\mathbb{D}$:

[TABLE]

The $\mathbf{0}$ block has dimensions $\ell$, whereas the identity block has dimensions ${m-\ell-2p}.$ When $\mathbb{D}\neq\mathbf{0},\mathbb{I}_{m}$, from [29, Corollary 5.2] and [57] it follows that $\|\mathbb{D}\|_{2}=\|\mathbb{I}_{m}-\mathbb{D}\|_{2}$. ∎

The novelty and the importance of Theorem 2 are in the interpretation in the DEIM setting, allowing for a deeper understanding of the structure of the DEIM projection and its error. For related usage of canonical angles between subspaces [8], see the construction of the favorable bases in [6].

Remark 3**.**

Let $\mathscr{S}^{T}U_{r}$ be invertible and $\mathbb{D}=U_{r}(\mathscr{S}^{T}U_{r})^{-1}\mathscr{S}^{T}$. If $f\in\mathcal{R}(U_{r})$ then both the DEIM error and the orthogonal projection error are zero, as $\mathbb{D}f=\mathcal{P}_{U_{r}}f=f$. In the case $f\neq\mathcal{P}_{U_{r}}f$, write $f-\mathbb{D}f=(\mathbb{I}_{m}-\mathcal{P}_{U_{r}})f+(\mathcal{P}_{U_{r}}-\mathbb{D})f$; verify that the summands are orthogonal, apply Pythagoras’ theorem to get

[TABLE]

Since $f\notin\mathcal{R}(U_{r})$, we can factor out $\|(\mathbb{I}_{m}-\mathcal{P}_{U_{r}})f\|^{2}$ to get

[TABLE]

(This is illustrated graphically in Figure 1.) Next, introduce the partition of $f$, represented in the basis $Z$, as follows:

[TABLE]

Now, straightforward computation for each $i=1,\ldots,p$ reveals that

[TABLE]

Together this gives

[TABLE]

Since $\mathscr{S}^{T}U_{r}$ is invertible, from the proof of Theorem 2, we have $D_{u}=\mathbf{0}$, and therefore

[TABLE]

Since $f\notin\mathcal{R}(U_{r})$, we can divide throughout by $\|f-\mathcal{P}_{U_{r}}f\|_{2}^{2}$ to obtain the inequality

[TABLE]

Combining this inequality with the relation for $\|\mathbb{D}f-\mathcal{P}_{U_{r}}f\|_{2}^{2}$ into (14) gives

[TABLE]

Therefore, $\kappa^{\prime}\leq\|\mathbb{D}\|_{2}=1/\cos\psi_{\ell+p}$. This result, of course, reproduces the bound (5). However, the analysis shows that a tighter condition number $\kappa^{\prime}$ can be obtained by considering how the contributions of the error are weighted in the principal directions identified in Theorem 2. **

3.3 Connection to CS decomposition

The structure of $\mathbb{D}$ can also be analyzed using the Cosine–Sine (CS) decomposition [55]. Assume for simplicity that the rows of $U_{r}$ are ordered so that $\mathscr{S}=\mathbb{I}_{m}(:,1:r)$. If this is not the case, we work with $\Pi^{T}\mathbb{D}\Pi$, where $\Pi$ is a permutation matrix. Assume that $\mathscr{S}^{T}U_{r}$ is invertible and therefore, the DEIM operator is $\mathbb{D}=U_{r}(\mathscr{S}^{T}U_{r})^{-1}\mathscr{S}^{T}$. Further, let $\mathscr{S}_{\perp}=\mathbb{I}_{m}(:,r+1:m)$.

With these assumptions, $U_{r}$ has the CS decomposition

[TABLE]

Here $\Omega_{1},\Gamma\in\mathbb{R}^{r\times r}$ and $\Omega_{2}\in\mathbb{R}^{(m-r)\times(m-r)}$ are orthogonal matrices and

[TABLE]

We can therefore represent $\mathbb{D}$ as

[TABLE]

where $\mathrm{Tan}\Psi=\Omega_{2}\mathrm{Sin}\Psi(\mathrm{Cos}\Psi)^{-1}\Omega_{1}^{T}=\Omega_{2}\mathrm{diag}(\tan\psi_{i})_{i=1}^{r}\Omega_{1}^{T}$. Similarly, we have

[TABLE]

and we (again) see that $\|\mathbb{D}\|_{2}=\|\mathbb{I}_{m}-\mathbb{D}\|_{2}=\sqrt{1+\|\mathrm{Tan}\Psi\|_{2}^{2}}.$ For further insights on the tangents between subspaces, see e.g., [67].

4 Weighted DEIM

As discussed in §1.3, the discrete analogue of a (generalized) interpolatory projection based approximation must be constructed within an appropriate weighted inner product, and the selection of the interpolation indices must ensure sharp error bounds. In particular, care must be taken to control how the condition number of the positive definite weight matrix $W$ influences the projection error, expressed in the $W$-weighted norm $\|u\|_{W}=\sqrt{u^{T}Wu}$. In this section, we address this issue and propose two new algorithms for $W$-weighted variants of DEIM.

To set the scene and to introduce notation, in §4.1 we recall the weighted POD. In §4.2, we propose $W$-DEIM oblique projection that relies on a more general form of the selection operator, and in the numerical realization uses $W$ implicitly through its Cholesky factor. In this case, although the pointwise interpolation is lost, the more general interpolation condition in the sense of GEIM holds true. In §4.6 and §4.6.1 we propose alternative methods for point selection in the weighted setting that allow for pointwise interpolation; however, the resulting approximation error bounds depend on the condition number of $W$ or on the condition number of optimally scaled $W$.

4.1 Setting the scene

Let $W\in\mathbb{R}^{m\times m}$ be symmetric positive definite, and define the weighted inner product for $u,v\in\mathbb{R}^{m}$ by $(u,v)_{W}\equiv v^{T}Wu.$ Let $W=LL^{T}$ be a factorization where the nonsingular matrix $L$ is a Cholesky factor or the positive definite square root $L=W^{1/2}$. The original problem might give rise to a nonsingular matrix $L$, so that the weight matrix $W=LL^{T}$ is then given implicitly by its factor $L$. Recall that any two square “Cholesky” factors of $W$ are related by an orthogonal matrix $Q$, so that $W^{1/2}=LQ$ [28, page 67, Exercise (x)].

Remark 4**.**

In the weighted norm $\|u\|_{W}\equiv\sqrt{(u,u)_{W}}=\sqrt{u^{T}Wu}=\|L^{T}u\|_{2},$ the induced operator norm of an $M\in\mathbb{R}^{m\times m}$ equals

[TABLE]

Further, in the $W$-inner product space, the adjoint of $M$ is $M^{[T]}\equiv W^{-1}M^{T}W,$ where $M^{T}$ is the transpose of $M$. **

The POD basis with respect to $(\cdot,\cdot)_{W}$ is determined by the 3-step procedure in Algorithm 1. For more details see [63]. For the sake of simplicity, we do not include centering of the snapshots matrix $Y$.

Algorithm 1 computes a matrix $\widehat{U}$ whose columns are $W$-orthonormal, i.e., $\widehat{U}^{T}W\widehat{U}=\mathbb{I}_{r},$ and the POD projection in the weighted inner product space is represented by

[TABLE]

Note that $\widehat{{\mathcal{P}}}_{\widehat{U}}^{2}=\widehat{{\mathcal{P}}}_{\widehat{U}}$ and that $\widehat{{\mathcal{P}}}_{\widehat{U}}^{[T]}=\widehat{{\mathcal{P}}}_{\widehat{U}}.$ In fact, $Y=\widehat{U}\Sigma V^{T}$ is a GSVD [62] of $Y$.

Remark 5**.**

For $\mathbb{R}^{m\times r}\ni\widehat{U}:(\mathbb{R}^{r},(\cdot,\cdot)_{2})\longrightarrow(\mathbb{R}^{m},(\cdot,\cdot)_{W})$, the adjoint matrix in the two inner products is, by definition, given as $\widehat{U}^{<T>}=\widehat{U}^{T}W$. Hence, $\widehat{U}^{<T>}\widehat{U}=\mathbb{I}_{r}$ and we can write the $W$-orthogonal projector (15) conveniently in the usual form as $\widehat{{\mathcal{P}}}_{\widehat{U}}=\widehat{U}\widehat{U}^{<T>}$. Recalling the discussion from §1.3, the projected problem (2) is then computed in the sense of $(\cdot,\cdot)_{W}$. **

4.2 $W$-DEIM

Once a discrete inner product $(\cdot,\cdot)_{W}$ has been chosen to capture the geometric framework (e.g., for Petrov-Galerkin projection, POD), one needs to define an appropriate DEIM projection operator in this weighted setting. Furthermore, the resulting quantities are now measured in the weighted norm $\|x\|_{W}$. To that end, using the notation introduced in Remark 5, we define a $W$-DEIM projector as follows.

Definition 6**.**

Let $\widehat{U}\in\mathbb{R}^{m\times r}$ be $W$-orthogonal. With a full column rank generalized selection operator $\mathbb{S}\in\mathbb{R}^{m\times s}$ (where $s\geq r$), define a weighted $W$-DEIM projector

[TABLE]

In the above definition, in addition to the use of a more general inner product, we also allow for tall rectangular $\mathbb{S}^{T}W\widehat{U}{\in\mathbb{C}^{s\times r}}$. The only constraint is that $\mathbb{S}^{T}W\widehat{U}$ has full column rank. However, in practice, we will use the square nonsingular case.

For the moment, we leave the (generalized) selection operator $\mathbb{S}$ unspecified, and we remark that the columns of $\mathbb{S}$ need not be the columns of the identity matrix. 555In fact, one can also allow full row rank to obtain a further variation of the DEIM projection as discussed in §3.1, but we omit this for the sake of brevity. As in the case of DEIM, the matrix $\mathbb{D}$ is an oblique projector, i.e., it satisfies $\mathbb{D}^{2}=\mathbb{D}.$

The following proposition is a recast of [71, Proposition 2.1] to the $\|\cdot\|_{W}$ norm.

Proposition 7**.**

Let $\mathbb{D}$ be as in Definition 6 and let $\mathbb{S}^{T}W\widehat{U}$ have full column rank. Then

[TABLE]

Proof.

Since $\mathbb{S}^{T}W\widehat{U}$ has full column rank, $(\mathbb{S}^{T}W\widehat{U})^{\dagger}$ is a left inverse, so that $\mathbb{D}\widehat{{\mathcal{P}}}_{\widehat{U}}=\widehat{{\mathcal{P}}}_{\widehat{U}},$ hence $(\mathbb{I}_{m}-\mathbb{D})\widehat{{\mathcal{P}}}_{U}=0$. Consequently for any vector $f\in\mathbb{R}^{m}$

[TABLE]

Since $\mathbb{D}$ is non-trivial projector ($\mathbb{D}\neq\mathbf{0}$, $\mathbb{D}\neq\mathbb{I}_{m}$) it holds that $\|\mathbb{D}\|_{W}=\|\mathbb{I}_{m}-\mathbb{D}\|_{W}$, and (17) follows. ∎

The condition number that amplifies the POD projection error $\|f-\widehat{{\mathcal{P}}}_{\widehat{U}}f\|_{W}$ is the weighted norm $\|\mathbb{D}\|_{W}$. A naive application of the result in Remark 4 suggests the bound $\|\mathbb{D}\|_{W}\leq\sqrt{\kappa(W)}\|\mathbb{D}\|_{2}$. That is, the condition number of the inner product matrix $W$ could potentially amplify the $W$-DEIM projection error. However, by a clever choice of $\mathbb{S}$ we can eliminate the factor $\sqrt{\kappa(W)}$.

Definition 8**.**

Let the weighted selection operator $\mathbb{S}$ and the corresponding $W$-DEIM projector $\mathbb{D}$, respectively, be defined as

[TABLE]

where $\mathscr{S}$ is an $m\times s$ index selection operator ($s$ selected columns of the identity $\mathbb{I}_{m}$, $s\geq r$).

Note that while $\mathbb{S}$ is possibly dense, $\mathscr{S}$ is a sparse matrix. We now present a result that quantifies the condition number $\|\mathbb{D}\|_{W}$ for the specific choice of selection operator $\mathbb{S}$.

Proposition 9**.**

Let $\mathbb{S}$ and $\mathbb{D}$ be defined as in (19). Then $\mathbb{S}^{T}W\mathbb{S}=\mathbb{I}_{k}$ and $\|D\|_{W}=\|(\mathscr{S}^{T}U_{r})^{\dagger}\|_{2}$.

Proof.

Recall that $L^{T}\widehat{U}=U_{r}$ and by (19), $\mathbb{S}^{T}L=\mathscr{S}^{T}$. Following Remark 4, straightforward computation yields

[TABLE]

where, by (19), $\mathbb{S}^{T}L=\mathscr{S}^{T}$, and thus $\|\mathbb{D}\|_{W}=\|U_{r}(\mathscr{S}^{T}U_{r})^{\dagger}\mathscr{S}^{T}\|_{2}=\|(\mathscr{S}^{T}U_{r})^{\dagger}\|_{2}$. ∎

Therefore, with this choice of $\mathbb{S}$, the condition number of $\|W\|_{2}$ does not explicitly appear in the bounds. However, the dependence on $W$ is implicitly contained in the matrix $U_{r}$ of the left singular vectors, and in the definition of $\widehat{U}$.

In §4.6 we present alternative choices for the Selection Operator $\mathbb{S}$ which can ensure pointwise interpolation.

Remark 10**.**

To obtain the canonical structure of $W$-DEIM, one follows the derivation from §3, properly adapted to the structure induced by $(\cdot,\cdot)_{W}$. **

4.3 How to choose $\mathscr{S}$

Recall that $\mathscr{S}$ contains carefully chosen columns of $\mathbb{I}_{m}$. The index selection to determine the columns of $\mathscr{S}$ can be computed using the original DEIM approach [15]. Another approach, Q-DEIM proposed in [19], uses a rank revealing QR factorization [11], implemented in high performance software libraries such as LAPACK [2] and ScaLAPACK [9].

However, in this paper, we adopt the strong Rank Revealing QR (sRRQR) factorization [25, Algorithm 4]. We present a result that characterizes the error of $W$-DEIM

Theorem 11**.**

Applying sRRQR [25, Algorithm 4] to $U_{r}$ produces an index selection operator $\mathscr{S}$ whose $W$-DEIM* projection error satisfies*

[TABLE]

Proof.

Combining (17) and Proposition 9 gives

[TABLE]

Since $U_{r}$ has orthonormal columns sRRQR [25, Algorithm 4] gives a selection operator $\mathscr{S}\in\mathbb{R}^{m\times r}$ such that $\mathscr{S}^{T}U_{r}$ is invertible. Applying Lemma 1 to bound $\|(\mathscr{S}^{T}U_{r})^{-1}\|_{2}$ gives the desired result. ∎

The importance of this result is that the point selection can also be applied in the weighted inner product case, and the resulting error bound similar as the DEIM bound in §2.

4.4 On the interpolating property and its generalization

Recall that the original DEIM formulation allows pointwise interpolation $\mathscr{S}^{T}\mathbb{D}f=\mathscr{S}^{T}f$, i.e., the projection $\mathbb{D}f$ and $f$ match exactly for a set of indices $i_{1},\dots,i_{r}$ determined by the columns of $\mathscr{S}$. In the case of $W$-DEIM, the following interpolation properties hold.

Proposition 12**.**

Let $\mathbb{S}^{T}W\widehat{U}$ be invertible and let $\mathbb{D}$ be as in Definition 6. Then $\mathbb{S}^{T}W\mathbb{D}f=\mathbb{S}^{T}Wf$.

This can be readily verified; since $\mathbb{S}^{T}W\widehat{U}$ is invertible then

[TABLE]

With the choice $\mathbb{S}=L^{-T}\mathscr{S}$, Proposition 12 simplifies to

[TABLE]

Hence, $W$-DEIM cannot in general interpolate $f\in\mathbb{R}^{m}$ at the selected indices $f_{i_{j}}=\phi_{i_{j}}(x_{i_{j}})$, $j=1,\ldots,r$. An exception to this is the case that $W$ has diagonal entries, see §4.5.1 for details. But, in many applications the discretized functions values may not be available through point evaluation either because there is no analytical expression or they may be sensor data corrupted by noise. In those cases, pointwise interpolation may not be possible, nor desirable – for a most illuminating discussion see [40].

4.4.1 DGEIM

The DEIM is a realization of the discrete version of the Empirical Interpolation Method (EIM) [5] in which interpolation was handled by only using pointwise function evaluation. In the same way we can interpret the interpolation condition (22) as a discrete version of GEIM, DGEIM, as a particular case of $W$-DEIM.

To this end, consider a more general concept of interpolation using a family of linear functionals, see [17, Chapter 11]. Introduce in (22) a column partition of $L=\begin{pmatrix}\ell_{1}&\dots&\ell_{m}\end{pmatrix}$ and rewrite it as

[TABLE]

If we interpret $\ell_{i}\in\mathbb{R}^{m}$ as the discretized Riesz representation of a given linear functional, then (22) interpolates the desired function $f$ at selected functionals. (The point interpolation corresponds to using the point evaluation functional, $(\ell_{i})_{j}=W_{ji}=W_{ij}=(\ell_{j})_{i}=\delta_{ij}$, where $\delta_{ij}$ is the Kronecker delta.)

4.5 How to ensure sparse selection

The original DEIM approximation was computationally efficient because it only required evaluating a small number of components of the vector $f$. However, in the computation of $\mathbb{D}f$, the factor $\mathscr{S}^{T}L^{T}f$ may, in the worst case, require many, or possibly all, components of $f$. This might make $W$-DEIM computationally inefficient. It is clear that the selection is sparse when the matrix $L$ is sparse.

The analysis is subdivided into three different cases. When the weighting matrix is sparse, or diagonal, the Cholesky factor $L$ is also sparse. When $W$ is sparse, reordering the matrix may lead to sparse factors $L$. On the other hand, if $W$ is dense, we must resort to an inexact sparse factorization. These cases are discussed below.

4.5.1 Diagonal weighting matrix $W$

If $W=\mathrm{diag}(w_{i})_{i=1}^{m}$, then

[TABLE]

and the computation of $\mathbb{D}f=W^{-1/2}U_{r}(\mathscr{S}^{T}U_{r})^{-1}\mathscr{S}^{T}{W}^{1/2}f$ requires only the indices $i_{1},\ldots,i_{r}$ of $f$ selected by $\mathscr{S}$. Furthermore, in this case the interpolation condition (22) simplifies to

[TABLE]

i.e., $\mathbb{D}$ is an interpolating projection.

4.5.2 Sparse weighting matrix $W$

In some cases, the matrix $W$ that defines a discrete inner product is large and sparse, and possibly contains additional block structure, see e.g., [32, §5.4]. Examples of sparse weighting matrices are discussed in the section on numerical experiments (Section 5).

When $W$ is sparse, one can take advantage of sparse factorization techniques to compute a pivoted factorization $\Pi^{T}W\Pi={L}_{s}{L}_{s}^{T}$, where the permutation matrix $\Pi$ is determined to produce a sparse Cholesky factor ${L}_{s}$. (In fact, the permutation matrix $\Pi$ has the additional benefit of making ${L}_{s}$ well conditioned for inversion by trying to improve diagonal dominance.) Then we factor $W=LL^{T}$ with $L=\Pi{L}_{s}$, and we have

[TABLE]

Since $\mathscr{S}^{T}(L_{s}^{T}\Pi^{T})$ will select only a small portion of the rows of a sparse matrix $L_{s}^{T}\Pi^{T}$, the product $\mathscr{S}^{T}L_{s}^{T}\Pi^{T}f$ is expect to require only relatively small number of the entries of $f$. An efficient implementation of this procedure would deploy the data structure and algorithms from the sparse matrices technology.

We now see an advantage of pure algebraic selection of the interpolation indices, as featured in the Q-DEIM version of the method [19]. In Q-DEIM, the index selection is computed by a rank revealing (column) pivoting in the QR factorization of the $r\times m$ matrix $U_{r}^{T}$, where $r\ll m$. The role of pivoting is to select an $r\times r$ submatrix of $U_{r}$ with small inverse. Hence, as argued in [19], it might be possible to find such a submatrix without having to touch all rows of $U_{r}$.

One possible way to improve sparsity is to lock certain columns of $U_{r}^{T}$ (whose indices correspond to non-sparse rows of $L_{s}^{T}$) and exclude them from the pivot selection. Since $m\gg r$, it is very likely that even with some columns of $U_{r}^{T}$ excluded, the selection will perform well. In fact, pivoting in the QR factorization can be modified to prefer indices that correspond to most sparse rows of $L_{s}^{T}$.

4.5.3 General dense positive definite $W$

In the most difficult case, the natural inner product is defined with large dimensional dense positive definite $W$ that is also difficult to compute. For instance, as mentioned in §1.3, $W$ can be the Gramian obtained by solving a large scale Lyapunov equation, or replaced by an empirical approximation based on the method of snapshots.

If computational complexity requires enforcing sparsity of the selection operator, then we can resort to inexact sparse factorization of the form $\Pi^{T}W\Pi+\delta W=\widetilde{L_{s}}\widetilde{L_{s}}^{T}$, i.e., we compute $W\approx(\Pi\widetilde{L_{s}})(\Pi\widetilde{L_{s}})^{T}$. The resulting approximation has the backward error $\Delta W=\Pi\delta W\Pi^{T}$ as a result of a thresholding strategy to produce the sparse factor $\widetilde{L_{s}}$. We mention two possibilities here. The incomplete Cholesky factorization is one candidate, see e.g. [37]. The matrix $W$ can also be sparsified by zeroing entries $W_{ij}$ if e.g., $|W_{ij}|/\sqrt{W_{ii}W_{jj}}$ is below some threshold.

Let us identify, for simplicity, $W\equiv\Pi^{T}W\Pi=LL^{T}$, so that $W+\delta W=\widetilde{L_{s}}\widetilde{L_{s}}^{T}$. Set $\widetilde{\mathbb{D}}=\widehat{U}(\mathscr{S}^{T}U_{r})^{-1}\mathscr{S}^{T}\widetilde{L_{s}}$. Then

[TABLE]

Now, from $\widetilde{\mathbb{D}}f=\mathbb{D}f+(\widetilde{\mathbb{D}}-\mathbb{D})f$ we have

[TABLE]

One can also justify using the sparsified weighting matrix in a backward sense, i.e. using $W+\delta W$ as the generator of the inner product. This line of reasoning via the incomplete factorization requires further analysis which we defer to our future work. Of course, in the case of dense $W$, saving the work in evaluating $f$ by the generalized interpolation (23) is nearly impossible as it may require too many entries to be practical. In that case, one can resort to point-wise interpolation that we discuss next.

4.6 Pointwise-interpolating $W$-DEIM

Note that in the formula for the $W$-DEIM projection in Definition 6 there is a certain freedom in choosing $\mathbb{S}$. The key in our formulation is indeed that we have left it as an adaptable device. In the case of the original DEIM with $W=\mathbb{I}_{m}$, $\mathbb{S}\equiv\mathscr{S}$ is a submatrix of $\mathbb{I}_{m}$, resulting in more efficient computation of the projection [15]. If a generalized interpolation of the type (22) and (23) is desired, then $\mathbb{S}^{T}=\mathscr{S}^{T}L^{-1}$ as in (19) in Proposition 9 will accomplish the task.

On the other hand, if we want point-wise interpolation

[TABLE]

also in the weighted case with a general positive definite $W$, then this can be obtained using the following definition.

Definition 13**.**

Let the weighted selection operator $\mathbb{S}$ and the corresponding $W$-DEIM projector $\mathbb{D}$, respectively, be defined as

[TABLE]

Here $\widehat{U}$ is $W$-orthogonal and $\mathscr{S}$ has columns from the identity matrix $\mathbb{I}_{m}$

Note that the relations (24), $\mathbb{D}\widehat{{\mathcal{P}}}_{\widehat{U}}=\widehat{{\mathcal{P}}}_{\widehat{U}},$ and the error estimate (17) still apply. However, now, the condition number $\|\mathbb{D}\|_{W}$ will depend on the specific choice of $\mathscr{S}$. We now show how to pick the indices that determine the columns of $\mathscr{S}$.

The algorithm proceeds as follows. First, as in Algorithm 1, a thin generalized SVD [62] of the $m\times n_{s}$ snapshot matrix $Y$ is computed and truncated to obtain low rank approximation $Y\approx\widehat{U}\widehat{\Sigma}\widehat{V}^{T}$, where $\widehat{V}^{T}\widehat{V}=\mathbb{I}_{r}$ and $\widehat{U}\in\mathbb{R}^{m\times r}$ is $W$-orthonormal, i.e., $\widehat{U}^{T}W\widehat{U}=\mathbb{I}_{r}$.

Then, the thin QR of $\widehat{U}=Q_{\widehat{U}}R_{\widehat{U}}$ is computed, and strong RRQR is applied to $Q_{\widehat{U}}^{T}$, to obtain the selection operator $\mathscr{S}$ (whose columns come from the $m\times m$ identity matrix). Finally, we set $\mathbb{S}\equiv W^{-1}\mathscr{S}$. This procedure is summarized in Algorithm 2, where the first two steps are implemented as in Algorithm 1. The corresponding error bound is given in Theorem 14.

Theorem 14**.**

Assume that the DEIM projection operator $\mathbb{D}$ is defined as in Algorithm 2. Then

[TABLE]

Proof.

Note that $\|\mathbb{D}\|_{W}=\|L^{T}\mathbb{D}L^{-T}\|_{2}\leq\sqrt{\kappa_{2}(W)}\|\mathbb{D}\|_{2}.$ We now bound $\|\mathbb{D}\|_{2}$. Consider the thin QR of $\widehat{U}=Q_{\widehat{U}}R_{\widehat{U}}$, where $R_{\widehat{U}}$ must be nonsingular. Then

[TABLE]

Since $Q_{\widehat{U}}$ and $\mathscr{S}$ have orthonormal columns, $\|\mathbb{D}\|_{2}=\|(\mathscr{S}^{T}Q_{\widehat{U}})^{-1}\|_{2}$. The rest of the proof is similar to Theorem 11. ∎

4.6.1 Scaling invariant error bound

Note that, compared to Theorem 11, the error bound (26) has an additional factor of $\sqrt{\kappa_{2}(W)}$. For highly ill-conditioned matrices $W$, this considerably inflates the error bound and possibly the actual error as well. It is instructive to see how a simple trick can improve this undesirable situation.

Let $\Delta=\mathrm{diag}(\sqrt{W_{ii}})_{i=1}^{m}$ and $W_{s}=\Delta^{-1}W\Delta^{-1}$; note that this scaling ensures $(W_{s})_{ii}=1$ for all $i=1,\dots,m$. It is well known (see [61]) that this diagonal equilibration nearly minimizes the spectral condition number over all diagonal scalings,

[TABLE]

where $\mathcal{D}^{m}$ is the space of diagonal $m\times m$ matrices. The task is to eliminate the scaling factor $\Delta$ from the bound on $\|\mathbb{D}\|_{W}$ (by the use of a different subset selection) and to replace $\sqrt{\kappa_{2}(W)}$ with $\sqrt{\kappa_{2}(W_{s})}$ – which can be a substantial improvement for certain applications of interest. To that end, we must examine how $W$ influences the structure of $\widehat{U}$, and interweave assembling of $\widehat{U}$ with the construction of the DEIM selection operator. The selection operator is $\mathbb{S}^{T}=\mathscr{S}^{T}W^{-1}$, as in Algorithm 2.

We use the expression for the weighted POD basis $\widehat{U}$ as in Algorithm 1, i.e. $\widehat{U}=L^{-T}U_{r}$, where $W=LL^{T}$ and $U_{r}^{T}U_{r}=\mathbb{I}_{r}$. If we define $L_{s}=\Delta^{-1}L$, then $W_{s}=L_{s}L_{s}^{T}$; $L_{s}$ has rows of unit Euclidean length, and, since $\mathbb{D}=\widehat{U}(\mathscr{S}^{T}\widehat{U})^{-1}\mathscr{S}$,

[TABLE]

The key observation is that $\mathscr{S}^{T}\Delta^{-1}=\widehat{\Delta}^{-1}\mathscr{S}^{T}$, where $\widehat{\Delta}$ is a diagonal matrix with the vector $\mathscr{S}^{T}\Delta$ on its diagonal. This cancels out $\Delta$,

[TABLE]

Let now $L_{s}^{-T}U_{r}=Q_{\widehat{U}}R_{s}$ be the QR factorization. (Note that $L_{s}^{-T}U_{r}=\Delta\widehat{U}$.) Then

[TABLE]

and we conclude that DEIM selection using $Q_{\widehat{U}}$ yields the desired bound

[TABLE]

These considerations are summarized in Algorithm 3 and Theorem 15.

Theorem 15**.**

Assume that the DEIM projection operator $\mathbb{D}$ is defined as in Algorithm 3. Then

[TABLE]

Remark 16**.**

It follows from (27) that the DEIM projection error bound (28) that applies to Algorithm 3 is never much worse ($\sqrt{\kappa_{2}(W_{s})}\leq\sqrt{m}\sqrt{\kappa_{2}(W)}$) and it is potentially substantially better666Take e.g. diagonal and highly ill-conditioned $W$. ($\sqrt{\kappa_{2}(W_{s})}\ll\sqrt{\kappa_{2}(W)}$) than the estimate (26) that holds for Algorithm 2. Although the two algorithms determine $\mathscr{S}$ from different orthonormal matrices, the factor $\sqrt{1+\eta^{2}r(m-r)}$ is the same, because of the property of the sRRQR. **

Remark 17**.**

In both Algorithm 3 and Algorithm 2, the sRRQR and computation of $\mathscr{S}$ can be replaced with the Q-DEIM selection [19], which is more efficient, essentially nearly as robust, but with weaker theoretical bound. However, the weaker upper bound on $\kappa$ is unlikely to make a substantial difference in practical computations, and both algorithms can be implemented using Q-DEIM. **

Remark 18**.**

For better numerical properties, the Cholesky factorization can be computed with pivoting, $\Pi^{T}W\Pi=LL^{T}$, i.e. $W=(\Pi L)(\Pi L)^{T}$, and we can easily modify Algorithm 3 to work implicitly with $\Pi L$ instead of $L$. **

Remark 19**.**

Note that the computation in Line 1. of Algorithm 3 can be rephrased as the GSVD of $Y$, $Y=U_{Y}\Sigma V^{T}$, where $U_{Y}=L^{-T}U$ is $W$-orthogonal, $U_{Y}^{T}WU_{Y}=\mathbb{I}_{m}$; see Algorithm 2. Then the matrix $\widehat{U}$ optionally returned by Algorithm 3 is $\widehat{U}=U_{Y}(:,1:r)=L^{-T}U_{r}$. Since $L_{s}^{-T}=\Delta L^{-T}$, the matrix $L_{s}^{-T}U_{r}$ in Line 4. can be expressed as $L_{s}^{-T}U_{r}=\Delta L^{-T}U_{r}=\Delta\widehat{U}$. **

5 Numerical Examples

In this section, we show numerical examples that highlight the benefits of our proposed algorithms.

5.1 Example 1

This example is based on [19, Example 3.1]. In this example we study the performance of sRRQR [25, Algorithm 4] for subset selection compared to the DEIM approach [15] and Q–DEIM [19]. Therefore, we let the weighting matrix $W=\mathbb{I}_{m}$. Let

[TABLE]

The snapshot set is generated by taking $40$ evenly spaced values of $\mu$ and $n=10,000$ evenly spaced points in time. The snapshots are collected in a matrix of size $10000\times 40$, the thin SVD of this matrix is computed and the left singular vectors corresponding to the first $34$ modes are used to define $U_{r}$.

To test the interpolation accuracy, we compute its value using the DEIM approximation at $200$ evenly spaced points in the $\mu$-domain. Three different subset selection procedures were used: DEIM, Pivoted QR labeled Q–DEIM, and sRRQR. In each case, we report the relative error defined as

[TABLE]

The results of the comparison are provided in Figure 2.

We observe that while all three methods are very accurate, Q–DEIM and sRRQR are much more accurate compared to DEIM for this example. Furthermore, from the right plot in Figure 2, we see that sRRQR is more accurate compared to both Q–DEIM and sRRQR. In practice, the performance of sRRQR is very similar to Q–DEIM, except for some adversarial cases in which Q–DEIM can fail spectacularly. In the subsequent examples, we use sRRQR for subset selection.

5.2 Example 2

Our next example is inspired by the Nonlinear RC-Ladder circuit, which is a standard benchmark problem for model reduction (see, for example [16, Section 6]). The underlying model is given by a dynamical system of the form

[TABLE]

where $g(x)=\exp(40x)+x-1$ and $u(t)=\exp(-t)$ and $N=1000$. The diagonal matrix $D$ is chosen to have entries

[TABLE]

The diagonal matrix $D$ induces the norm $\|\cdot\|_{D}$ and the relative error between the full and the reduced order models are measured in this norm.

The dynamical system is simulated over $t=[0,7]$ seconds and $2000$ snapshots of the dynamical system and the nonlinear function are collected with equidistant time steps. Based on the decay rate of the snapshots, we vary the number of basis vectors from $5$ to $40$. The relative error is defined to be

[TABLE]

where $x(t)$ is the solution of the dynamical system at time $t$, whereas $\hat{x}(t)$ is the reduced order approximation at the same time. The relative error as a function of number of retained basis vectors is plotted in left panel of Figure 3. On the right, the reconstruction of the first component of the dynamical system $x_{1}(t)$ is shown; here $k=40$ basis vectors were retained. As can be seen, the reconstruction error is low and the $W$-DEIM, indeed, approximates the large-scale dynamical system accurately.

5.3 Example 3

This example is inspired by [48, Section 2.3]. The spatial domain is taken to be $\Omega=[0,1]^{2}$ and the parameter domain is $\mathcal{D}=[0,1]^{2}$. We define a function $g:\Omega\times\mathcal{D}\rightarrow\mathbb{R}$ which satisfies

[TABLE]

where $h(z;\mu)=((1-z)-(0.99\cdot\mu-1))^{2}$. The function that is to be interpolated is

[TABLE]

Depending on the parameter $\mu$, it has a sharp peak in one of the four corners of $\Omega$. The function is discretized on a $100\times 100$ grid in $\Omega$, and parameter samples are drawn from a $25\times 25$ equispaced grid in $\mathcal{D}$. These $625$ snapshots are used to construct the DEIM approximation. We choose three different weighting functions: $W_{1}$ is the identity matrix, $W_{2}$ is the weighting matrix corresponding to the $L^{2}(\Omega)$ inner product, and $W_{3}$ is the weighting matrix corresponding to the $H^{1}(\Omega)$ inner product.

We then compute the average relative error over a test sample corresponding to a $11\times 11$ equispaced grid in $\mathcal{D}$. The relative error is defined to be

[TABLE]

The POD basis is computed using Algorithm 1, whereas the subset selection is done using sRRQR [25, Algorithm 4]. The results of the interpolation errors as a function of number of DEIM interpolation points retained, is displayed in the left panel of Figure 4. On the right hand panel of the same figure, we display the error constants $\|\mathbb{D}\|_{W}=\|(\mathscr{S}^{T}U_{r})^{-1}\|$. As can be seen, although the error constants increase with increasing number of basis vectors, the overall interpolation error decreases resulting an effective approximation.

5.4 Example 4

This is a continuation of Example 3. We use the same setup as before; however, we compare the different algorithms for $W$-DEIM. In ‘Method 1’ we use Algorithm 1 to generate the POD basis, while the subset selection is done using sRRQR [25, Algorithm 4]. The error constant for this method is $\eta_{1}\equiv\|(\mathscr{S}^{T}U_{r})^{-1}\|_{2}$. In ‘Method 2’ we use Algorithm 2 with error constant $\eta_{2}\equiv\sqrt{\kappa_{2}(W)}\|(\mathscr{S}^{T}Q_{\widehat{U}})^{-1}\|_{2}$ and in ‘Method 3’ we use Algorithm 3 with error constant $\eta_{3}\equiv\sqrt{\kappa_{2}(W_{s})}\|(\mathscr{S}^{T}Q_{\widehat{U}})^{-1}\|_{2}$.

In Algorithm 2, the GSVD of the snapshot matrix w.r.t. the weighting matrix $W$ was computed as follows. First, the weighted QR was computed using [38, Algorithm 2] to obtain $Y=Q_{Y}R_{Y}$. Note that $Q_{Y}^{T}WQ_{Y}=\mathbb{I}_{n_{s}}$. Then the SVD of $R_{Y}$ is computed as $R_{Y}=U_{R}\Sigma V^{T}$. We obtain the GSVD of $Y={U}_{Y}\Sigma V^{T}$, where now ${U}_{Y}=Q_{Y}U_{R}$.

For a given weighting matrix, we define the relative error as

[TABLE]

The DEIM operators $\mathbb{D}_{j}$ correspond to the different Methods described above. In Figure 5 we plot the relative error using the DEIM approximation and error constants; here the weighting matrix $W=W_{2}$ corresponds to the $L^{2}(\Omega)$ inner product. As can be seen, the overall interpolation error from all three methods are comparable. However, the error constants for Method 2 are highest as expected, since it involves $\sqrt{\kappa_{2}(W)}$. In Figure 6 we repeat the same experiment; however, the weighting matrix $W=W_{3}$ corresponds to the $H^{1}(\Omega)$ inner product. Our conclusions are similar to the previous weighting matrix. Note here that $W_{3}$ is more ill-conditioned than $W_{2}$ and furthermore, for $W=W_{3}$ we have that $\kappa_{2}(W)\approx\kappa_{2}(W_{s})$. Therefore, the difference between the error constants for Methods 2 and 3 is very small.

In conclusion, for the application at hand, all three $W$-DEIM methods produce comparable results. Methods 2 and 3 maybe desirable if factorization of $W$ is computationally expensive, or even infeasible.

5.5 Example 5

In this example, we consider a parameterized PDE based on [50, Section 8.5]. Consider the following parameterized PDE form defined on domain $\Omega=[0,1]^{2}$ with boundary $\partial\Omega$

[TABLE]

Here $\boldsymbol{\mu}=[\mu_{1},\mu_{2},\mu_{3}]$, and $\mathbf{n}$ is the normal vector. The wind velocity $\mathbf{b}(\mu_{1})$ is taken to be as $\boldsymbol{b}=[\cos\mu_{1},\sin\mu_{1}]$, which is a constant in space but depends nonlinearly on the parameter $\mu_{1}$. The source term $s(\boldsymbol{\mu})$ has the form of a Gaussian function centered at $(\mu_{2},\mu_{3})$ and spread $0.25$

[TABLE]

The goal of this problem is to construct a reduced order model for the solution $u(\mathbf{x};\boldsymbol{\mu})$ in the domain $\Omega$ over the range of parameters $\mu_{1}\in[0,2\pi]$, $\mu_{2}\in[0.2,0.8]$ and $\mu_{3}\in[0.15,0.35]$. A POD based approach is used to reduce the model of the parameterized PDE with DEIM/WDEIM approximation for the source term.

As the weighting matrix $W$, we choose the arising from the discrete representation of the $H_{1}(\Omega)$ inner product. For constructing the WPOD and WDEIM bases, we first generated a training set of parameters $\boldsymbol{\mu}$ of $1000$ points generated by Latin Hypercube sampling; then the source term and the solution of the PDE is computed at each training point $\boldsymbol{\mu}$. The maximum dimension for the WPOD and WDEIM bases were chosen to be $20$ and $24$ respectively based on the decay of the singular values. From the same snapshot set we also compute bases for the POD and DEIM with dimensions $20$ and $24$ respectively. For both approaches, we use the PQR for computing for point selections. We report the errors used by both approaches in Figure 7. The errors were averaged over $10$ different randomly generate samples in the parameter range. In the left panel, we compare the error in the DEIM approximation and WDEIM approximations of the source term $s(\mathbf{x};\boldsymbol{\mu})$, whereas in the right panel, we consider the error in the solution of $u(\mathbf{x};\boldsymbol{\mu})$ over the same test samples. For the right panel, the dimension of the DEIM/W-DEIM was chosen to be $24$ and the dimension of the POD/WPOD basis dimension was chosen to be $28$. All the errors were computed with the weighted norm $\|\cdot\|_{W}$. We see that the error in our approach (WPOD-WDEIM) is comparable with that of the POD-DEIM approach, whereas the error in the WDEIM approach is slightly better than the error in the error in the DEIM approach.

6 Conclusions

The main contributions of this work are: (i) it defines a new index selection operator, based on strong rank revealing QR factorization, that nearly attains the optimal DEIM projection error bound; (ii) it facilitates the understanding of the canonical structure of the DEIM projection; (iii) it establishes a core numerical linear algebra framework for the DEIM projection in weighted inner product spaces; (iv) it defines a discrete version of the Generalized Empirical Interpolation Method (GEIM). We believe that these will be useful for further development of the DEIM idea and its applications in scientific computing.

Acknowledgements

We are indebted to Ilse Ipsen and in particular to the two anonymous referees for constructive criticism and many suggestions that have improved the presentation of the paper.