Oblique projection for scalable rank-adaptive reduced-order modelling of nonlinear stochastic partial differential equations with time-dependent bases

(a) Setup

Consider the nonlinear stochastic PDE given by

augmented with appropriate initial and boundary conditions. In the above equation, $v=v(x,t;\mathit{\xi })$, $x$ is the spatial coordinate, $\mathit{\xi }\in {\mathbb{R}}^d$ are the set of random parameters, $t$ is the time and $f(v;x,t,\mathit{\xi })$ includes the nonlinear spatial differential operators. We assume generic nonlinear PDEs, where the nonlinearity of $f$ versus $v$ may be non-polynomial, e.g. exponential, fractional, etc. For the sake of simplicity in the exposition, we consider a collocation/strong-form discretization of equation (2.1) in $x$ and $\mathit{\xi }$. Because of the simplicity of the resulting discrete system, this choice facilitates an uncluttered illustration of the main contribution of this paper, which is focused on the efficient low-rank approximation of nonlinear MDEs. However, the presented methodology can also be applied to other types of discretizations, for example, weak form discretizations (finite element, etc.). Examples of collocation/strong-form discretizations in the spatial domain are Fourier/polynomial spectral collocation schemes or finite-difference discretizations. Example collocation schemes in the random domain include the probabilistic collocation method (PCM) [32] or any Monte Carlo-type sampling methods [33–35]. Applying any of the above schemes to equation (2.1) leads to the following nonlinear MDE: where $I=[0,T_f]$ denotes the time interval, $\mathbf{\text{V}}(t):I\to {\mathbb{R}}^{n\times s}$ is a matrix with $n$ rows corresponding to collocation points in the spatial domain and $s$ columns corresponding to collocation/sampling points of the parameters $\mathit{\xi }$ and $\mathcal{F}(t,\mathbf{\text{V}}):I\times {\mathbb{R}}^{n\times s}\to {\mathbb{R}}^{n\times s}$ is obtained by discretizing $f(v;x,t,\mathit{\xi })$ in $x$ and $\mathit{\xi }$. Equation (2.2) is augmented with appropriate initial conditions, i.e. $\mathbf{\text{V}}(t_0)={\mathbf{\text{V}}}_0$. We also assume that boundary conditions are already incorporated into equation (2.2), which can be accomplished in a number of ways, for example by using weak treatment of the boundary conditions [36]. For the remainder of this paper, we will refer to equation (2.2) as the FOM, which will be used as the ground truth for evaluating the performance of the proposed methodology. For the problems targeted in this work, we assume $n>s$ without loss ofgenerality.

The presented methodology is limited to explicit time integration schemes. For the computational complexity analysis, we consider sparse discretization schemes for spatial discretization, which means that each row is dependent on $p_a$ rows, where $p_a\ll n$. The majority of discretization schemes, e.g. finite difference, finite volume, finite element, spectral element, result in sparse row dependence. As a result, the computational cost of computing each column of FOM (equation (2.2)) is $\mathcal{O}(n)$ and the cost of solving MDE (2.2) for all $s$ columns scales with $\mathcal{O}(ns)$. We also note that the presented methodology is not limited to sparse spatial discretizations and can be applied to dense discretizations as well. See remark 2.10 for more details.

(b) Preliminaries

In this section, we present some of the definitions of matrix manifolds, tangent spaces, orthogonal and oblique projections and CUR decomposition.

Definition 2.1. (Low-rank matrix manifolds)

The low-rank matrix manifold ${\mathcal{M}}_r$ is defined as the set

$${\mathcal{M}}_r=\{\widehat{\mathbf{\text{V}}}\in {\mathbb{R}}^{n\times s}:\hspace{0.17em}\text{rank}(\widehat{\mathbf{\text{V}}})=r\},$$

of matrices of fixed rank $r$. Any member of the set ${\mathcal{M}}_r$ is denoted by a hat symbol $(\widehat{\hspace{0.17em}\hspace{0.17em}})$, e.g. $\widehat{\mathbf{\text{V}}}$.

Any member of ${\mathcal{M}}_r$ may be represented by $\widehat{\mathbf{\text{V}}}=\mathbf{\text{U}}\mathit{\Sigma }{\mathit{Y}}^T$, where $\mathbf{\text{U}}\in {\mathbb{R}}^{n\times r}$ and $\mathbf{\text{Y}}\in {\mathbb{R}}^{s\times r}$ are a set of orthonormal columns and $\mathit{\Sigma }\in {\mathbb{R}}^{r\times r}$ is a rank-$r$ matrix. The rank-$r$ matrix $\widehat{\mathbf{\text{V}}}$ may also be represented via the multiplication of two matrices, i.e. $\widehat{\mathbf{\text{V}}}=\mathbf{\text{U}}{\mathit{Y}}^T$, where $\mathbf{\text{U}}\in {\mathbb{R}}^{n\times r}$ and $\mathbf{\text{Y}}\in {\mathbb{R}}^{s\times r}$ have full column rank.

Definition 2.2. (Tangent space)

The tangent space of manifold ${\mathcal{M}}_r$ at $\widehat{\mathbf{\text{V}}}$, represented with the decomposition of $\widehat{\mathbf{\text{V}}}=\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T$, is the set of matrices in the form of Koch & Lubich [2]:

$${\mathcal{T}}_{\widehat{\mathbf{\text{V}}}}{\mathcal{M}}_r=\{\delta \mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T+\mathbf{\text{U}}\delta \mathit{\Sigma }{\mathbf{\text{Y}}}^T+\mathbf{\text{U}}\mathit{\Sigma }\delta {\mathbf{\text{Y}}}^T:\hspace{0.17em}\delta {\mathbf{\text{U}}}^T\mathbf{\text{U}}=\mathbf{\text{0}}\text{and}\delta {\mathbf{\text{Y}}}^T\mathbf{\text{Y}}=\mathbf{\text{0}}\},$$

where $\delta \mathbf{\text{U}}\in {\mathbb{R}}^{n\times r}$ and $\delta \mathbf{\text{Y}}\in {\mathbb{R}}^{s\times r}$.

Definition 2.3. (Orthogonal projection onto the tangent space)

The orthogonal projection of matrix $\mathbf{\text{W}}\in {\mathbb{R}}^{n\times s}$ onto the tangent space of manifold ${\mathcal{M}}_r$ at $\widehat{\mathbf{\text{V}}}$, represented with the decomposition of $\widehat{\mathbf{\text{V}}}=\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T$, is given by Koch & Lubich [2, lemma 4.1]:

$${\mathcal{P}}_{{\mathcal{T}}_{\widehat{\mathbf{\text{V}}}}}(\mathbf{\text{W}})=\mathbf{\text{U}}{\mathbf{\text{U}}}^T\mathbf{\text{W}}+\mathbf{\text{W}}\mathbf{\text{Y}}{\mathbf{\text{Y}}}^T-\mathbf{\text{U}}{\mathbf{\text{U}}}^T\mathbf{\text{W}}\mathbf{\text{Y}}{\mathbf{\text{Y}}}^T.$$2.3

In the above projection, $\mathbf{\text{U}}{\mathbf{\text{U}}}^T$ and $\mathbf{\text{Y}}{\mathbf{\text{Y}}}^T$ are orthogonal projections onto spaces spanned by the columns of $\mathbf{\text{U}}$ and $\mathbf{\text{Y}}$. We denote these orthogonal projections with:

where the symbol $∟$ indicates orthogonal projection. In the following, we define oblique projectors. We first explain the notation that is used in this section. Let $\mathbf{\text{U}}\in {\mathbb{R}}^{n\times r}$ and $\mathbf{\text{Y}}\in {\mathbb{R}}^{s\times r}$ be matrices whose columns are orthonormal and let $\mathbf{\text{p}}=[p_1,p_2,\dots ,p_{r^{\prime }}]\in {\mathbb{N}}^{r^{\prime }}$ and $\mathbf{\text{s}}=[s_1,s_2,\dots ,s_{r^{\prime }}]\in {\mathbb{N}}^{r^{\prime }}$ be vectors containing row and column indices, where the number of indices can be greater than or equal to the dimension of the subspaces spanned by $\mathbf{\text{U}}$ and $\mathbf{\text{Y}}$, i.e. $r^{\prime }\ge r$. Also, $r^{\prime }\le n$ for row indices and $r^{\prime }\le s$ for column indices. We use MATLAB indexing where $\mathbf{\text{V}}(\mathbf{\text{p}},:)\in {\mathbb{R}}^{r^{\prime }\times s}$ selects all columns at the $\mathbf{\text{p}}$ rows, and $\mathbf{\text{V}}(:,\mathbf{\text{s}})\in {\mathbb{R}}^{n\times r^{\prime }}$ selects all rows at the $\mathbf{\text{s}}$ columns of the matrix $\mathbf{\text{V}}$. We also use the indexing matrices, $\mathbf{\text{P}}={\mathbf{\text{I}}}_n(:,\mathbf{\text{p}})\in {\mathbb{R}}^{n\times r^{\prime }}$ and $\mathbf{\text{S}}={\mathbf{\text{I}}}_s(:,\mathbf{\text{s}})\in {\mathbb{R}}^{s\times r^{\prime }}$, where ${\mathbf{\text{I}}}_n$ and ${\mathbf{\text{I}}}_s$ are identity matrices of size $n\times n$ and $s\times s$, respectively. It is easy to verify that ${\mathbf{\text{P}}}^T\mathbf{\text{U}}\equiv \mathbf{\text{U}}(\mathbf{\text{p}},:)$ and ${\mathbf{\text{S}}}^T\mathbf{\text{Y}}\equiv \mathbf{\text{Y}}(\mathbf{\text{s}},:)$. Let ${(\hspace{0.17em})}^{†}$ denote the Moore–Penrose pseudoinverse of a matrix, i.e. ${\mathbf{\text{A}}}^{†}={({\mathbf{\text{A}}}^T\mathbf{\text{A}})}^{-1}{\mathbf{\text{A}}}^T$.

Definition 2.4. (Oblique projection)

Let $\mathbf{\text{U}}\in {\mathbb{R}}^{n\times r}$ and $\mathbf{\text{Y}}\in {\mathbb{R}}^{s\times r}$ be orthonormal matrices and let $\mathbf{\text{p}}\in {\mathbb{N}}^{r^{\prime }}$ and $\mathbf{\text{s}}\in {\mathbb{N}}^{r^{\prime }}$ be sets of distinct row and column indices, respectively. Oblique projectors onto Ran($\mathbf{\text{U}}$) and Ran($\mathbf{\text{Y}}$) are defined as [37]

$${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}=\mathbf{\text{U}}{({\mathbf{\text{P}}}^T\mathbf{\text{U}})}^{†}{\mathbf{\text{P}}}^T\text{and}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}=\mathbf{\text{S}}{({\mathbf{\text{Y}}}^T\mathbf{\text{S}})}^{†}{\mathbf{\text{Y}}}^T,$$2.5

provided ${({\mathbf{\text{P}}}^T\mathbf{\text{U}})}^T({\mathbf{\text{P}}}^T\mathbf{\text{U}})\in {\mathbb{R}}^{r\times r}$ and ${({\mathbf{\text{Y}}}^T\mathbf{\text{S}})}^T({\mathbf{\text{Y}}}^T\mathbf{\text{S}})\in {\mathbb{R}}^{r\times r}$ are invertible.

In the above definition, the symbol $\mathrm{\angle }$ distinguishes these projectors from orthogonal projectors. For a given matrix $\mathbf{\text{A}}\in {\mathbb{R}}^{n\times s}$, ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ operates on the left side of the matrix and ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$ operates on the right side of the matrix. It is easy to verify that ${({\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})}^2={\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ and ${({\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})}^2={\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$. It is also easy to verify that the oblique projection of $\mathbf{\text{A}}$ belongs to the manifold of rank-$r$ matrices, i.e. ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{A}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}\in {\mathcal{M}}_r$.

For the special case of $r^{\prime }=r$, the oblique projectors ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ and ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$ are also interpolatory projectors. In this case, the oblique projectors become ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}=\mathbf{\text{U}}{({\mathbf{\text{P}}}^T\mathbf{\text{U}})}^{-1}{\mathbf{\text{P}}}^T$ and ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}=\mathbf{\text{S}}{({\mathbf{\text{Y}}}^T\mathbf{\text{S}})}^{-1}{\mathbf{\text{Y}}}^T$. Unlike orthogonal projection or a general oblique projection, the interpolatory projection is guaranteed to match the original matrix at the selected rows and columns, i.e.

The other extreme is when all the rows or columns are selected, i.e. $r^{\prime }=n$ or $r^{\prime }=s$. Take, for example, the projector ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ when $r^{\prime }=n$. In this case, ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ becomes the same as the orthogonal projector, i.e. ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\equiv {\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{∟}$. To show this, first note that ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ is invariant with respect to the ordering of the row indices ($\mathbf{\text{p}}$) and when $r^{\prime }=n$, $\mathbf{\text{p}}$ can be taken to be: $\mathbf{\text{p}}=[1,2,\dots ,n]$. In this case, $\mathbf{\text{P}}\equiv {\mathbf{\text{I}}}_n$. Therefore: where we have used the orthonormality condition of $\mathbf{\text{U}}$, i.e. ${\mathbf{\text{U}}}^T\mathbf{\text{U}}={\mathbf{\text{I}}}_r$, where ${\mathbf{\text{I}}}_r$ is the $r\times r$ identity matrix. The analogous relationship exists for ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$, when $r^{\prime }=s$. In the following, we define CUR decompositions, which are closely related to the oblique projections.

Definition 2.5. (CUR decomposition)

A CUR decomposition of matrix $\mathbf{\text{V}}$ is a rank-$r$ approximation of $\mathbf{\text{V}}$ in the form of $\mathbf{\text{V}}\approx \mathbf{\text{C}}\mathbf{\text{U}}\mathbf{\text{R}}$, where $\mathbf{\text{C}}\in {\mathbb{R}}^{n\times r}$ and $\mathbf{\text{R}}\in {\mathbb{R}}^{r\times s}$ are actual columns and rows of matrix $\mathbf{\text{V}}$, i.e. $\mathbf{\text{C}}=\mathbf{\text{V}}(:,\mathbf{\text{s}})$ and $\mathbf{\text{R}}=\mathbf{\text{V}}(\mathbf{\text{p}},:)$. The matrix $\mathbf{\text{U}}\in {\mathbb{R}}^{r\times r}$ is computed such that $\mathbf{\text{C}}\mathbf{\text{U}}\mathbf{\text{R}}$ is a good approximation to $\mathbf{\text{V}}$. The CUR of matrix $\mathbf{\text{V}}$ is denoted with $\mathtt{\text{CUR}}(\mathbf{\text{V}})$.

Here, the matrices $\mathbf{\text{C}}$, $\mathbf{\text{U}}$ and $\mathbf{\text{R}}$ are different from the matrices defined in previous sections. Different CUR decompositions can be obtained for the same matrix depending on two factors: (i) the selection of columns and rows and (ii) the method used to compute the matrix $\mathbf{\text{U}}$. For more details on CUR decompositions, we refer the reader to Mahoney & Drineas [38]. Finally, it is easy to verify that $\mathtt{\text{CUR}}(\mathbf{\text{V}})\in {\mathcal{M}}_r$. The connection between CUR decomposition and oblique projections is shown in §2g.

(c) Time-continuous variational principle

The central idea behind TDB-based low-rank approximation is that the bases evolve optimally to minimize the residual due to low-rank approximation error. The residual is obtained by substituting an SVD-like low-rank approximation into the FOM so that $\mathbf{\text{V}}(t)$ is closely approximated by the rank-$r$ matrix

where $\mathbf{\text{U}}(t)\in {\mathbb{R}}^{n\times r}$ is a time-dependent orthonormal spatial basis for the column space, $\mathbf{\text{Y}}(t)\in {\mathbb{R}}^{s\times r}$ is a time-dependent orthonormal parametric basis for the row space, $\mathit{\Sigma }(t)\in {\mathbb{R}}^{r\times r}$ is, in general, a full matrix and $r\ll \min (n,s)$ is the rank of the approximation.

Because this is a low-rank approximation, it cannot satisfy the FOM exactly and there will be a residual equal to:

This residual is minimized via the first-order optimality conditions of the variational principle given by subject to orthonormality constraints on $\mathbf{\text{U}}$ and $\mathbf{\text{Y}}$. Since the above variational principle involves the time-continuous equation (i.e. no temporal discretization is applied), the idea is to minimize the instantaneous residual by optimally updating $\mathbf{\text{U}}$, $\mathit{\Sigma }$ and $\mathbf{\text{Y}}$ in time. Therefore, we refer to this as the time-continuous variational principle. As indicated in [2,7], the optimality conditions of equation (2.8) lead to closed-form evolution equations for $\mathbf{\text{U}}$, $\mathit{\Sigma }$ and $\mathbf{\text{Y}}$: where $\mathbf{\text{F}}\in {\mathbb{R}}^{n\times s}$ is a matrix defined as $\mathbf{\text{F}}=\mathcal{F}(t,\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T)$, and $\mathbf{\text{I}}$ is the identity matrix of appropriate dimensions. The above variational principle is the same as the Dirac–Frenkel time-dependent variational principle in the quantum chemistry literature [8] or the dynamical low-rank approximation (DLRA) [2]. In [2], the geometry of the tangent space, ${\mathcal{T}}_{\widehat{\mathbf{\text{V}}}}{\mathcal{M}}_r$, was exploited to solve the constrained residual minimization problem given by equation (2.8). In this setting, the residual, $\mathcal{I}(\dot{\widehat{\mathbf{\text{V}}}})=||\dot{\widehat{\mathbf{\text{V}}}}-\mathcal{F}(t,\widehat{\mathbf{\text{V}}})|{|}_F$, is minimized with the constraint that $\dot{\widehat{\mathbf{\text{V}}}}(t)\in {\mathcal{T}}_{\widehat{\mathbf{\text{V}}}}{\mathcal{M}}_r$. The solution to the above minimization problem is obtained by where ${\mathcal{P}}_{{\mathcal{T}}_{\widehat{\mathbf{\text{V}}}}}$ is the orthogonal projection onto the tangent space ${\mathcal{T}}_{\widehat{\mathbf{\text{V}}}}$ at $\widehat{\mathbf{\text{V}}}=\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T$ [2, Lemma 4.1]. It is easy to show that equations (2.9a)–(2.9c) can be recovered from equation (2.10).

As was shown in [39], it is possible to derive a similar variational principle for the DO decomposition, $\widehat{\mathbf{\text{V}}}(t)={\mathbf{\text{U}}}_{\mathrm{DO}}(t){\mathbf{\text{Y}}}_{\mathrm{DO}}^T(t)$, whose optimality conditions are constrained to the orthonormality of the spatial modes, ${\mathbf{\text{U}}}_{\mathrm{DO}}^T{\mathbf{\text{U}}}_{\mathrm{DO}}=\mathbf{\text{I}}$, via the DO condition, ${\dot{\mathbf{\text{U}}}}_{\mathrm{DO}}^T{\mathbf{\text{U}}}_{\mathrm{DO}}=\mathbf{\text{0}}$. However, for the sake of simplicity and unlike the original DO formulation presented in [1], an evolution equation for the mean field is not derived. Without loss of generality, the low-rank DO evolution equations become

and where $\mathbf{\text{C}}={\mathbf{\text{Y}}}_{\mathrm{DO}}^T{\mathbf{\text{Y}}}_{\mathrm{DO}}$ is the low-rank correlation matrix. Note that the low-rank approximation based on DO is equivalent to equation (2.6), i.e. ${\mathbf{\text{U}}}_{\mathrm{DO}}{\mathbf{\text{Y}}}_{\mathrm{DO}}^T=\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T$. Similarly, the BO decomposition, $\widehat{\mathbf{\text{V}}}(t)={\mathbf{\text{U}}}_{\mathrm{BO}}(t){\mathbf{\text{Y}}}_{\mathrm{BO}}^T(t)$, which is subject to BO conditions, ${\mathbf{\text{U}}}_{\mathrm{BO}}^T{\mathbf{\text{U}}}_{\mathrm{BO}}=\text{diag}({\lambda }_1,\dots ,{\lambda }_r)$ and ${\mathbf{\text{Y}}}_{\mathrm{BO}}^T{\mathbf{\text{Y}}}_{\mathrm{BO}}=\mathbf{\text{I}}$, is also identical to DO and equation (2.6). As it was shown in [11], one can derive matrix differential equations that transform the factorization $\{\mathbf{\text{U}},\mathit{\Sigma },\mathbf{\text{Y}}\}$ to $\{{\mathbf{\text{U}}}_{\mathrm{DO}},{\mathbf{\text{Y}}}_{\mathrm{DO}}\}$ or $\{{\mathbf{\text{U}}}_{\mathrm{BO}},{\mathbf{\text{Y}}}_{\mathrm{BO}}\}$. The equivalence of DO and BO formulations was shown in [12]. Using the DO/BO terminology, equations (2.9a)–(2.9c) have both DO and BO conditions, i.e. the DO conditions for $\mathbf{\text{U}}$ and $\mathbf{\text{Y}}$: ${\dot{\mathbf{\text{U}}}}^T\mathbf{\text{U}}=\mathbf{\text{0}}$ and ${\dot{\mathbf{\text{Y}}}}^T\mathbf{\text{Y}}=\mathbf{\text{0}}$ as well as bi-orthonormality conditions: ${\mathbf{\text{U}}}^T\mathbf{\text{U}}=\mathbf{\text{I}}$ and ${\mathbf{\text{Y}}}^T\mathbf{\text{Y}}=\mathbf{\text{I}}$. Despite their equivalence, these three factorizations have different numerical performances in the presence of small singular values. As was shown in [11], equations (2.9a)–(2.9c) outperform both DO and BO.

Despite the potential of equations (2.9a)–(2.9c) to significantly reduce the computational cost of solving massive MDEs like equation (2.2), there are still a number of outstanding challenges for most practical problems of interest. As highlighted in the Introduction, computing $\mathbf{\text{F}}=\mathcal{F}(t,\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T)$ requires $\mathcal{O}(ns)$ operations that scale with the size of the FOM. This involves applying the nonlinear map ($\mathcal{F}$) on every column of the matrix $\widehat{\mathbf{\text{V}}}=\mathbf{\text{U}}\mathit{\Sigma }{\mathbf{\text{Y}}}^T$. While it is possible to achieve $\mathcal{O}(n+s)$ for the special cases of homogeneous linear and quadratic nonlinear $\mathcal{F}$, this comes at the expense of a highly intrusive process, that requires a careful term-by-term treatment of the right side of equations (2.9a)–(2.9c) [21, Appendix B]. Furthermore, solving equations (2.9b) and (2.9c) becomes unstable when $\mathit{\Sigma }$ is singular or near singular. This is particularly problematic because it is often necessary to retain very small singular values in order to have an accurate approximation.

While the low-rank approximation based on TDBs can be cast in different, yet equivalent formulations, we have chosen equations (2.9a)–(2.9c) over DO/BO/DDO decompositions to highlight the underlying challenges. Since DO/BO/DDO decompositions exhibit all of the above challenges, addressing these challenges in the context of equations (2.9a)–(2.9c) automatically addresses the DO/BO/DDO challenges as well.

(d) Time-discrete variational principle

To address the challenges of low-rank approximations based on TDB using the time-continuous variational principle, we consider a time-discrete variational principle for rank-adaptive matrix approximations, which has recently been applied in [40] and also [25,41] in the context of tensors. To this end, consider an explicit Runge–Kutta temporal discretization of equation (2.2):

where $\mathrm{\Delta }t$ is the step size and $\overline{\mathbf{\text{F}}}$ is obtained via an explicit Euler or Runge–Kutta scheme. For example, the first-order explicit Euler method is given by: $\overline{\mathbf{\text{F}}}=\mathcal{F}(t^{k-1},{\widehat{\mathbf{\text{V}}}}^{k-1})$. In the above equation, it is important to note that ${\mathbf{\text{V}}}^k$ is not the FOM solution since the right-hand side is computed using the low-rank state from the previous time step. Despite using the rank-$r$ ${\widehat{\mathbf{\text{V}}}}^{k-1}$ in equation (2.12), ${\mathbf{\text{V}}}^k$ will not be a rank-$r$ matrix, i.e. ${\mathbf{\text{V}}}^k\notin {\mathcal{M}}_r$. Excluding some rare exceptions, taking one step according to equation (2.12) will put ${\mathbf{\text{V}}}^k$ off the rank-$r$ manifold. Therefore, to solve the MDE while remaining on ${\mathcal{M}}_r$, a rank truncation is needed to map the solution back onto the rank-$r$ manifold at each time step. In other words, we need to approximate ${\mathbf{\text{V}}}^k$ with a rank-$r$ matrix, ${\widehat{\mathbf{\text{V}}}}^k$, such that where ${\mathbf{\text{R}}}^k$ is the low-rank approximation error.

The time-discrete variational principle can be stated as finding the best ${\widehat{\mathbf{\text{V}}}}^k\in {\mathcal{M}}_r$ such that the Frobenius norm of the residual is minimized [25]:

The solution of the above residual minimization scheme is obtained via where $\mathtt{\text{SVD}}({\mathbf{\text{V}}}^k)={\mathbf{\text{U}}}_{\mathrm{best}}^k{\mathit{\Sigma }}_{\mathrm{best}}^k{\mathbf{\text{Y}}}_{\mathrm{best}}^{k^T}$ is the rank-$r$ truncated SVD of matrix ${\mathbf{\text{V}}}^k$, where ${\mathbf{\text{U}}}_{\mathrm{best}}^k\in {\mathbb{R}}^{n\times r}$ and ${\mathbf{\text{Y}}}_{\mathrm{best}}^k\in {\mathbb{R}}^{s\times r}$ are the matrices of the first $r$ left and right singular vectors of ${\mathbf{\text{V}}}^k$, respectively and ${\mathit{\Sigma }}_{\mathrm{best}}^k\in {\mathbb{R}}^{r\times r}$ is the matrix of singular values.

An important advantage of equation (2.15) over equations (2.9a)–(2.9c) is that the time advancement according to equation (2.15) does not become singular in the presence of small singular values. While this solves the issue of ill-conditioning, computing equation (2.15) at each iteration of the time stepping scheme is cost prohibitive. This computational cost is due to two sources: (i) computing the nonlinear map $\mathcal{F}(t^{k-1},{\widehat{\mathbf{\text{V}}}}^{k-1})$ to obtain ${\mathbf{\text{V}}}^k$ and (ii) computing $\mathtt{\text{SVD}}({\mathbf{\text{V}}}^k)$. The cost of (i) alone makes the solution of the time-discrete variational principle as expensive as the FOM, i.e. $\mathcal{O}(ns)$. Besides the flops cost associated with computing ${\mathbf{\text{V}}}^k$, the memory cost of storing ${\mathbf{\text{V}}}^k$ is prohibitive for most realistic applications. On the other hand, computing the exact SVD of ${\mathbf{\text{V}}}^k$ scales with $\min \{\mathcal{O}(n^3),\mathcal{O}(s^3)\}$. While this cost is potentially alleviated by fast algorithms for approximating the SVD, e.g. randomized SVD [42] or incremental QR [37], for general nonlinearities in $\mathcal{F}$, (i) is unavoidable. This ultimately leads to a computational cost that exceeds that of the FOM.

(e) Low-rank approximation via an oblique projection

To overcome these challenges, we present an oblique projection scheme that enables a cost-effective approximation to the rank-$r$ $\mathtt{\text{SVD}}({\mathbf{\text{V}}}^k)$. Before presenting our methodology, we provide a geometric interpretation of $\mathtt{\text{SVD}}({\mathbf{\text{V}}}^k)$. In particular, $\mathtt{\text{SVD}}({\mathbf{\text{V}}}^k)$ can be interpreted as an orthogonal projection onto the column and row space of ${\widehat{\mathbf{\text{V}}}}_{\mathrm{best}}^k$ as shown below:

where ${\mathbf{\text{P}}}_{{\mathbf{\text{U}}}_{\mathrm{best}}^k}^{∟}={\mathbf{\text{U}}}_{\mathrm{best}}^k{\mathbf{\text{U}}}_{\mathrm{best}}^{k^T}$ and ${\mathbf{\text{P}}}_{{\mathbf{\text{Y}}}_{\mathrm{best}}^k}^{∟}={\mathbf{\text{Y}}}_{\mathrm{best}}^k{\mathbf{\text{Y}}}_{\mathrm{best}}^{k^T}$ are orthogonal projections onto the column and row space of ${\widehat{\mathbf{\text{V}}}}_{\mathrm{best}}^k$, respectively. The approximation ${\widehat{\mathbf{\text{V}}}}_{\mathrm{best}}^k$ is the optimal rank-$r$ approximation of ${\mathbf{\text{V}}}^k$, however as mentioned in the previous section, computing ${\widehat{\mathbf{\text{V}}}}_{\mathrm{best}}^k$ is more expensive than solving the FOM.

In the following, we present a methodology that computes an accurate approximation to ${\widehat{\mathbf{\text{V}}}}_{\mathrm{best}}^k$ in a cost-effective manner. From the geometric perspective, our approach is to use an oblique projection onto a set of rank-$r$ orthonormal column (${\mathbf{\text{U}}}^k$) and row (${\mathbf{\text{Y}}}^k$) subspaces:

The above equation is analogous to equation (2.16) where the orthogonal projections are replaced with oblique projections. While orthogonal projection requires access to the entire ${\widehat{\mathbf{\text{V}}}}^{k-1}+\mathrm{\Delta }t\overline{\mathbf{\text{F}}}$ matrix, the oblique projectors can be designed to require the computation of ${\widehat{\mathbf{\text{V}}}}^{k-1}+\mathrm{\Delta }t\overline{\mathbf{\text{F}}}$ at $\mathcal{O}(r)$ columns and rows. This geometric perspective is depicted in figure 1a.

• Download figure
• Open in new tab
• Download PowerPoint

From the matrix decomposition point of view, the above approximation may be represented via a CUR decomposition:

where CUR represents the algorithmic implementation of a CUR decomposition. See figure 1c.

From the residual minimization perspective, the SVD can be viewed as a Galerkin projection where $||{\mathbf{\text{R}}}^k|{|}_F$ is minimized. On the other hand, the presented approach based on interpolatory projection (a special case of oblique projection) can be viewed as a collocated scheme where the residual is set to zero at $r$ strategically selected rows and columns of the residual matrix, ${\mathbf{\text{R}}}^k$, in equation (2.13). To this end, we present an algorithm to set ${\mathbf{\text{R}}}^k(\mathbf{\text{p}},:)=\mathbf{\text{0}}$ and ${\mathbf{\text{R}}}^k(:,\mathbf{\text{s}})=\mathbf{\text{0}}$, where $\mathbf{\text{p}}\in {\mathbb{N}}^r$ and $\mathbf{\text{s}}\in {\mathbb{N}}^r$ are vectors containing the row and column indices at which the residual is set to zero. This simply requires ${\widehat{\mathbf{\text{V}}}}^k(\mathbf{\text{p}},:)={\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ and ${\widehat{\mathbf{\text{V}}}}^k(:,\mathbf{\text{s}})={\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$. See figure 1b. In §2h, we consider oblique projection for the general case when $r^{\prime }>r$ (definition 2.4), where the residual at the selected rows and columns is not guaranteed to be zero.

Although the approach we will present is equivalent to the oblique projection of equation (2.17), ${\mathbf{\text{U}}}^k$ and ${\mathbf{\text{Y}}}^k$ are the unknown column and row bases of ${\widehat{\mathbf{\text{V}}}}^k$ at the current time step. Therefore, equation (2.17) cannot be readily used for the computation of ${\widehat{\mathbf{\text{V}}}}^k$. In the following, we present a methodology to compute these bases (and subsequently ${\widehat{\mathbf{\text{V}}}}^k$) by strategically sampling $r$ columns and rows of ${\mathbf{\text{V}}}^k$. While there are many possible choices for the indices $\mathbf{\text{p}}$ and $\mathbf{\text{s}}$, selecting these points should be done in a principled manner, to ensure the residual at all points remains small. To compute these points, we use the discrete empirical interpolation method (DEIM) [43] which has been shown to provide near-optimal sampling points for computing CUR matrix decompositions [37]. A similar approach was recently applied in [21] to accelerate the computation of equations (2.9a)–(2.9c), by only sampling $\mathbf{\text{F}}$ at a small number of rows and columns. However, the approach presented in [21] still suffers from the issue of ill-conditioning.

To compute the DEIM points, the rank-$r$ SVD (or an approximation) is required [37]. Since we do not have access to the rank-$r$ SVD at the current time step, $k$, we use the approximation of the SVD from the previous time step, ${\widehat{\mathbf{\text{V}}}}^{k-1}={\mathbf{\text{U}}}^{k-1}{\mathit{\Sigma }}^{k-1}\mathbf{\text{Y}}{{\hspace{0.17em}}^{k-1}}^T$, to compute the DEIM points. The initial approximation is ideally obtained from the FOM initial condition as the rank-$r$ $\mathtt{\text{SVD}}({\mathbf{\text{V}}}_0)$. The algorithm for computing ${\widehat{\mathbf{\text{V}}}}^k$ using interpolation is as follows:

(i)	Compute the sampling indices, $\mathbf{\text{p}}\leftarrow \mathtt{DEIM}({\mathbf{\text{U}}}^{k-1})$, and $\mathbf{\text{s}}\leftarrow \mathtt{DEIM}({\mathbf{\text{Y}}}^{k-1})$, in parallel.
(ii)	Compute ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ and ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ by taking one step according to equation (2.12) at the selected rows and columns, in parallel.
(iii)	Compute $\mathbf{\text{Q}}\in {\mathbb{R}}^{n\times r}$ as the orthonormal basis for the range of ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ by QR decomposition such that ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})=\mathbf{\text{Q}}\mathbf{\text{R}}$, where $\mathbf{\text{R}}\in {\mathbb{R}}^{r\times r}$.
(iv)	Interpolate every column of ${\mathbf{\text{V}}}^k$ onto the orthonormal basis $\mathbf{\text{Q}}$ at sparse indices $\mathbf{\text{p}}$: $$\mathbf{\text{Z}}=\mathbf{\text{Q}}{(\mathbf{\text{p}},:)}^{-1}{\mathbf{\text{V}}}^k(\mathbf{\text{p}},:),$$2.19 where $\mathbf{\text{Z}}\in {\mathbb{R}}^{r\times s}$ is the matrix of interpolation coefficients such that $\mathbf{\text{Q}}\mathbf{\text{Z}}$ interpolates ${\mathbf{\text{V}}}^k$ onto the basis $\mathbf{\text{Q}}$ at the interpolation points indexed by $\mathbf{\text{p}}$.
(v)	Compute the SVD of $\mathbf{\text{Z}}$ so that $$\mathbf{\text{Z}}={\mathbf{\text{U}}}_{\mathbf{\text{Z}}}{\mathit{\Sigma }}^k{{\mathbf{\text{Y}}}^k}^T,$$2.20 where ${\mathbf{\text{U}}}_{\mathbf{\text{Z}}}\in {\mathbb{R}}^{r\times r}$, ${\mathit{\Sigma }}^k\in {\mathbb{R}}^{r\times r}$ and ${\mathbf{\text{Y}}}^k\in {\mathbb{R}}^{s\times r}$.
(vi)	Compute ${\mathbf{\text{U}}}^k\in {\mathbb{R}}^{n\times r}$ as the in-subspace rotation: $${\mathbf{\text{U}}}^k=\mathbf{\text{Q}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}.$$2.21

In Step (i), the details of the DEIM algorithm can be found in [43, algorithm 1]. A DEIM algorithm based on the QR factorization, a.k.a QDEIM, may also be used [44,45]. Both DEIM and QDEIM are sparse selection algorithms and they perform comparably in the cases considered in this paper. We explain here how the above algorithm addresses the three challenges mentioned in the Introduction (1).

(i)	Computational efficiency: The above procedure returns the updated low-rank approximation ${\widehat{\mathbf{\text{V}}}}^k=\mathbf{\text{QZ}}={\mathbf{\text{U}}}^k{\mathit{\Sigma }}^k{{\mathbf{\text{Y}}}^k}^T$, and only requires sampling ${\mathbf{\text{V}}}^k$ at $r$ rows and columns. This alone significantly reduces both the required number of flops and memory, compared to computing the entire ${\mathbf{\text{V}}}^k$. Furthermore, instead of directly computing the SVD of the $n\times s$ matrix ${\mathbf{\text{V}}}^k$, we only require computing the QR of the $n\times r$ matrix ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$, and the SVD of the $r\times s$ matrix $\mathbf{\text{Z}}$. This reduces the computational cost to $\mathcal{O}(s+n)$ for $r\ll s$ and $r\ll n$. Moreover, in most practical applications, computing ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ is the costliest part of the algorithm, which requires solving $s$ samples of the FOM. However, since these samples are independent of each other, the columns of ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ can be computed in parallel. Similarly, each row of ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ can be computed in parallel.
(ii)	Intrusiveness: While this significantly reduces the computational burden, perhaps an equally important outcome is the minimally intrusive nature of the above approach. For example, when the columns of ${\mathbf{\text{V}}}^k$ are independent, e.g. random samples, ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ can be computed by directly applying equation (2.12) to the selected columns of the low-rank approximation from the previous time step. This effectively allows for existing numerical implementations of equation (2.12) to be used as a black box for computing ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$. The non-intrusive column sampling in the presented algorithm is the counterpart of solving equation (2.9b) in DLRA and equation (2.11a) in DO. However, equations (2.9b) and (2.11a) require deriving and implementing new PDEs, whereas the presented algorithm allows an existing deterministic solver to be used in a black box fashion, in which a suitable column space basis is extracted. On the other hand, the rows of ${\mathbf{\text{V}}}^k$ are in general dependent, based on a known map for the chosen spatial discretization scheme, e.g. sparse discretizations like finite difference, finite element, or dense discretization schemes, e.g. global spectral methods. Therefore, computing ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ does require specific knowledge of the governing equations, namely the discretized differential operators. Based on the discretization scheme, one can determine a set of adjacent points, ${\mathbf{\text{p}}}_a$, that are required for computing the derivatives at the points specified by $\mathbf{\text{p}}$. While this introduces an added layer of complexity, this is much less intrusive than deriving and implementing reduced-order operators for each term in the governing equations; which we emphasize again, is only feasible for homogeneous linear or quadratic nonlinear equations. In the present work, that bottleneck is removed, regardless of the type of nonlinearity.
(iii)	Ill-conditioning: The presented algorithm is robust in the presence of small or zero singular values. First note that the inversion of the matrix of singular values is not required in the presented algorithm. In fact, the conditioning of the algorithm depends on $\mathbf{\text{Q}}(\mathbf{\text{p}},:)$ and $\mathbf{\text{Y}}(\mathbf{\text{s}},:)$, and the DEIM algorithm ensures that these two matrices are well-conditioned. To illustrate this point, let us consider the case of overapproximation where the rank of ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ is $r_1<r$. In this case, equations (2.9a)–(2.9c) and equations (2.11a)–(2.11b) cannot be advanced because $\mathit{\Sigma }\in {\mathbb{R}}^{r\times r}$ and $\mathbf{\text{C}}\in {\mathbb{R}}^{r\times r}$ will be singular, i.e. rank$(\mathit{\Sigma })=$rank$(\mathbf{\text{C}})=r_1<r$. On the other hand, despite ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ being rank-deficient, $\mathbf{\text{Q}}$ will still be a full rank matrix in the presented algorithm. While there is no guarantee that a subset of rows of $\mathbf{\text{Q}}$, i.e. $\mathbf{\text{Q}}(\mathbf{\text{p}},:)$ is well-conditioned, the DEIM is a greedy algorithm that is designed to keep $||\mathbf{\text{Q}}{(\mathbf{\text{p}},:)}^{-1}||$ as small as possible in a near-optimal fashion. In §2h, we show that oversampling further improves the condition number of the presented algorithm, and in theorem 2.8, we show that $||\mathbf{\text{Y}}{(\mathbf{\text{s}},:)}^{-1}||$ plays an equally important role in maintaining a well-conditioned algorithm.

As we will show in §2g, the low-rank approximation computed above is equivalent to a CUR matrix decomposition that interpolates ${\mathbf{\text{V}}}^k$ at the selected rows and columns. Therefore, we refer to the above procedure as the TDB-CUR algorithm.

(f) Computing ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$

Up until this point, we have considered ${\mathbf{\text{V}}}^k={\widehat{\mathbf{\text{V}}}}^{k-1}+\mathrm{\Delta }t\overline{\mathbf{\text{F}}}$ to be an $n\times s$ matrix resulting from an explicit Runge–Kutta temporal discretization of equation (2.12). We showed that sparse row and column measurements, ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ and ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$, could be used to efficiently compute an approximation to the rank-$r$ SVD of ${\mathbf{\text{V}}}^k$. While ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$ is straightforward to compute for independent random samples, as discussed in §2e, computing ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ depends on a set of adjacent points, ${\mathbf{\text{p}}}_a$, according to the spatial discretization scheme. As a result, for higher-order integration schemes, special care must be taken in the computation of ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$. To demonstrate this, we consider the second-order explicit Runge–Kutta scheme where

After determining the row indices, $\mathbf{\text{p}}$ and ${\mathbf{\text{p}}}_a$, ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ can be computed as follows:

(i)	Compute ${\widehat{\mathbf{\text{V}}}}^{k-1}([\mathbf{\text{p}},{\mathbf{\text{p}}}_a],:)={\mathbf{\text{U}}}^{k-1}([\mathbf{\text{p}},{\mathbf{\text{p}}}_a],:){\mathit{\Sigma }}^{k-1}{{\mathbf{\text{Y}}}^{k-1}}^T$.
(ii)	Compute the first stage ${\mathbf{\text{F}}}_1=\mathcal{F}(t^{k-1},{\widehat{\mathbf{\text{V}}}}^{k-1})$ at the $\mathbf{\text{p}}$ rows as $${\mathbf{\text{F}}}_1(\mathbf{\text{p}},:)=\mathcal{F}(t^{k-1},{\widehat{\mathbf{\text{V}}}}^{k-1}([\mathbf{\text{p}},{\mathbf{\text{p}}}_a],:)).$$ Note, if the explicit Euler method is used, $\overline{\mathbf{\text{F}}}={\mathbf{\text{F}}}_1$, and no additional steps are required. Simply compute ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)={\widehat{\mathbf{\text{V}}}}^{k-1}+\mathrm{\Delta }t{\mathbf{\text{F}}}_1(\mathbf{\text{p}},:)$. If a higher-order scheme is used, proceed with the following steps.
(iii)	The final stage of the second-order integration scheme requires taking a half step to evaluate $\mathcal{F}$ at the midpoint: $${\mathbf{\text{F}}}_2=\mathcal{F}(t^{k-1}+\frac{1}{2}\mathrm{\Delta }t,{\widehat{\mathbf{\text{V}}}}^{k-1}+\frac{1}{2}\mathrm{\Delta }t{\mathbf{\text{F}}}_1).$$ Note that for the second-order scheme, $\overline{\mathbf{\text{F}}}={\mathbf{\text{F}}}_2$. Here, we require ${\mathbf{\text{F}}}_2(\mathbf{\text{p}},:)$, given by $${\mathbf{\text{F}}}_2(\mathbf{\text{p}},:)=\mathcal{F}(t^{k-1}+\frac{1}{2}\mathrm{\Delta }t,{\widehat{\mathbf{\text{V}}}}^{k-1}([\mathbf{\text{p}},{\mathbf{\text{p}}}_a],:)+\frac{1}{2}\mathrm{\Delta }t{\mathbf{\text{F}}}_1([\mathbf{\text{p}},{\mathbf{\text{p}}}_a],:)).$$ Note that we now require ${\mathbf{\text{F}}}_1({\mathbf{\text{p}}}_a,:)$ to evaluate the above expression. While this can be computed according to Step (ii), where ${\mathbf{\text{p}}}_a$ will have its own set of adjacent points ${\mathbf{\text{p}}}_{aa}$, this process quickly gets out of hand, especially as more stages are added to the integration scheme. As a result, for higher-order schemes, the efficiency afforded by the presented algorithm will deteriorate, and the resulting implementation will become increasingly complex. To overcome these challenges, we instead compute the low-rank approximation ${\widehat{\mathbf{\text{F}}}}_i\approx {\mathbf{\text{F}}}_i$, using the sparse row and column measurements, ${\mathbf{\text{F}}}_i(\mathbf{\text{p}},:)$ and ${\mathbf{\text{F}}}_i(:,\mathbf{\text{s}})$, which are already required for computing ${\mathbf{\text{V}}}^k(\mathbf{\text{p}},:)$ and ${\mathbf{\text{V}}}^k(:,\mathbf{\text{s}})$. Here, the subscript denotes the $i\mathrm{th}$ stage of the integration scheme. The first step is to compute ${\mathbf{\text{U}}}_{{\mathbf{\text{F}}}_i}$ as an orthonormal basis for the $\mathrm{Ran}({\mathbf{\text{F}}}_i(:,\mathbf{\text{s}}))$, using QR. Next, compute the oblique projection of ${\mathbf{\text{F}}}_i$ onto ${\mathbf{\text{U}}}_{{\mathbf{\text{F}}}_i}$, such that ${\widehat{\mathbf{\text{F}}}}_i={\mathbf{\text{U}}}_{{\mathbf{\text{F}}}_i}{\mathbf{\text{U}}}_{{\mathbf{\text{F}}}_i}{(\mathbf{\text{p}},:)}^{-1}{\mathbf{\text{F}}}_i(\mathbf{\text{p}},:)$. Using this low-rank approximation, ${\mathbf{\text{F}}}_i({\mathbf{\text{p}}}_a,:)$ is readily approximated by ${\widehat{\mathbf{\text{F}}}}_i({\mathbf{\text{p}}}_a,:)={\mathbf{\text{U}}}_{{\mathbf{\text{F}}}_i}({\mathbf{\text{p}}}_a,:){\mathbf{\text{U}}}_{{\mathbf{\text{F}}}_i}{(\mathbf{\text{p}},:)}^{-1}{\mathbf{\text{F}}}_i(\mathbf{\text{p}},:)$. Although we have considered the second-order Runge–Kutta method in the example above, this approach is easily extended to higher-order Runge–Kutta methods. It is straightforward to show that the above procedure is equivalent to a CUR decomposition of matrix ${\mathbf{\text{F}}}_i$, similar to our previous work [21].

(g) Equivalence to a CUR decomposition and oblique projection

Before presenting the details of our methodology in §2e, we discuss how the presented approach can be understood as an oblique projection (equation (2.17)) or alternatively as a CUR decomposition (equation (2.18)). In this section, we show that (i) the presented algorithm is equivalent to a CUR decomposition (theorem 2.6) and (ii) the matrix ${\widehat{\mathbf{\text{V}}}}^k$ is obtained via an oblique projection, which requires access to only the selected rows and columns of ${\mathbf{\text{V}}}^k$ (theorem 2.7). In theorems 2.6 and 2.7, $\mathbf{\text{P}}$ and $\mathbf{\text{S}}$ are matrices of size $n\times r$ and $s\times r$, respectively. We drop the time step $k$ for simplicity.

Theorem 2.6.

Let $\hat{\mathbf{V}}=\mathbf{Q}\mathbf{Z}$ be the low-rank approximation of $\mathbf{V}$ computed according to the TDB-CUR algorithm. Then: (i) $\hat{\mathbf{V}}=\mathbf{Q}\mathbf{Z}$ is equivalent to the CUR factorization given by $(\mathbf{V}\mathbf{S}){({\mathbf{P}}^T\mathbf{V}\mathbf{S})}^{-1}({\mathbf{P}}^{\mathrm{T}}\mathbf{V})$. (ii) The low-rank approximation is exact at the selected rows and columns, i.e. ${\mathbf{P}}^T\hat{\mathbf{V}}={\mathbf{P}}^T\mathbf{V}$ and $\hat{\mathbf{V}}\mathbf{S}=\mathbf{V}\mathbf{S}$.

Proof.

(i)

According to the TDB-CUR algorithm, $\mathbf{\text{Q}}$ is a basis for the $\mathrm{Ran}(\mathbf{\text{V}}(:,\mathbf{\text{s}}))$. Therefore, $\mathbf{\text{V}}(:,\mathbf{\text{s}})=\mathbf{\text{VS}}={\mathbf{\text{QQ}}}^T\mathbf{\text{V}}\mathbf{\text{S}}$, and it follows that ${\mathbf{\text{P}}}^T\mathbf{\text{VS}}={\mathbf{\text{P}}}^T{\mathbf{\text{QQ}}}^T\mathbf{\text{VS}}$. Substituting this result into the CUR factorization gives $$(\mathbf{\text{VS}}){({\mathbf{\text{P}}}^T\mathbf{\text{VS}})}^{-1}({\mathbf{\text{P}}}^T\mathbf{\text{V}})={\mathbf{\text{QQ}}}^T\mathbf{\text{VS}}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}}{\mathbf{\text{Q}}}^T\mathbf{\text{VS}})}^{-1}{\mathbf{\text{P}}}^T\mathbf{\text{V}}.$$ Rearranging the above expression gives the desired result $$(\mathbf{\text{VS}}){({\mathbf{\text{P}}}^T\mathbf{\text{VS}})}^{-1}({\mathbf{\text{P}}}^T\mathbf{\text{V}})={\mathbf{\text{QQ}}}^T\mathbf{\text{VS}}{({\mathbf{\text{Q}}}^T\mathbf{\text{VS}})}^{-1}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}})}^{-1}{\mathbf{\text{P}}}^T\mathbf{\text{V}}=\mathbf{\text{QZ}}=\widehat{\mathbf{\text{V}}},$$ where we have used $\mathbf{\text{Z}}={({\mathbf{\text{P}}}^T\mathbf{\text{Q}})}^{-1}{\mathbf{\text{P}}}^T\mathbf{\text{V}}=\mathbf{\text{Q}}{(\mathbf{\text{p}},:)}^{-1}\mathbf{\text{V}}(\mathbf{\text{p}},:)$, from equation (2.19).

(ii)

Using the above result, $\widehat{\mathbf{\text{V}}}=(\mathbf{\text{VS}}){({\mathbf{\text{P}}}^T\mathbf{\text{VS}})}^{-1}({\mathbf{\text{P}}}^T\mathbf{\text{V}})$, we show the selected rows of $\widehat{\mathbf{\text{V}}}$ are exact, i.e. $\widehat{\mathbf{\text{V}}}(\mathbf{\text{p}},:)=\mathbf{\text{V}}(\mathbf{\text{p}},:)$: $${\mathbf{\text{P}}}^T\widehat{\mathbf{\text{V}}}=({\mathbf{\text{P}}}^T\mathbf{\text{VS}}){({\mathbf{\text{P}}}^T\mathbf{\text{VS}})}^{-1}({\mathbf{\text{P}}}^T\mathbf{\text{V}})={\mathbf{\text{P}}}^T\mathbf{\text{V}}.$$ Similarly for the columns, $$\widehat{\mathbf{\text{V}}}\mathbf{\text{S}}=(\mathbf{\text{VS}}){({\mathbf{\text{P}}}^T\mathbf{\text{VS}})}^{-1}({\mathbf{\text{P}}}^T\mathbf{\text{VS}})=\mathbf{\text{VS}}.$$ This completes the proof.

Now we show that $\widehat{\mathbf{\text{V}}}$ is an oblique projection of $\mathbf{\text{V}}$ onto the selected columns and rows of $\mathbf{\text{V}}$. In particular, the oblique projector involved is an interpolatory projector. For the sake of brevity, we drop the superscript $k$ in the following.

Theorem 2.7.

Let $\hat{\mathbf{V}}=\mathbf{Q}\mathbf{Z}$ be the low-rank approximation of $\mathbf{V}$ computed according to the TDB-CUR algorithm. Then $\hat{\mathbf{V}}={\mathbf{P}}_{\mathbf{U}}^{\mathrm{\angle }}\mathbf{V}{\mathbf{P}}_{\mathbf{Y}}^{\mathrm{\angle }}$ where ${\mathbf{P}}_{\mathbf{U}}^{\mathrm{\angle }}$ and ${\mathbf{P}}_{\mathbf{Y}}^{\mathrm{\angle }}$ are oblique projectors onto Ran($\mathbf{U}$) and Ran($\mathbf{Y}$), respectively, according to equation (2.5).

Proof.

We first show that ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ can be represented versus $\mathbf{\text{Q}}$ as the interpolation basis. To this end, replacing $\mathbf{\text{U}}=\mathbf{\text{Q}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}$ in the definition of ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ results in:

$${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}=\mathbf{\text{U}}{({\mathbf{\text{P}}}^T\mathbf{\text{U}})}^{-1}{\mathbf{\text{P}}}^T=\mathbf{\text{Q}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}})}^{-1}{\mathbf{\text{P}}}^T=\mathbf{\text{Q}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}^{-1}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}})}^{-1}{\mathbf{\text{P}}}^T=\mathbf{\text{Q}}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}})}^{-1}{\mathbf{\text{P}}}^T.$$

where we have used the fact that ${\mathbf{\text{U}}}_{\mathbf{\text{Z}}}$ is a square orthonormal matrix and therefore, ${\mathbf{\text{U}}}_{\mathbf{\text{Z}}}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}^{-1}=\mathbf{\text{I}}$. Similarly, ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$ can be represented versus ${\mathbf{\text{Z}}}^T$ as the interpolation basis by replacing ${\mathbf{\text{Y}}}^T={\mathit{\Sigma }}^{-1}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}^{-1}\mathbf{\text{Z}}$ in ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$:

$${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}=\mathbf{\text{S}}{({\mathbf{\text{Y}}}^T\mathbf{\text{S}})}^{-1}{\mathbf{\text{Y}}}^T=\mathbf{\text{S}}{({\mathit{\Sigma }}^{-1}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}^{-1}\mathbf{\text{Z}}\mathbf{\text{S}})}^{-1}{\mathit{\Sigma }}^{-1}{\mathbf{\text{U}}}_{\mathbf{\text{Z}}}^{-1}\mathbf{\text{Z}}=\mathbf{\text{S}}{(\mathbf{\text{Z}}\mathbf{\text{S}})}^{-1}\mathbf{\text{Z}}.$$

Using these projection operators we have

$${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}=\mathbf{\text{Q}}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}})}^{-1}{\mathbf{\text{P}}}^T\mathbf{\text{V}}\mathbf{\text{S}}{(\mathbf{\text{Z}}\mathbf{\text{S}})}^{-1}\mathbf{\text{Z}}.$$2.22

Using the results of theorem 2.6, Part (ii), we have: $\mathbf{\text{V}}(\mathbf{\text{p}},\mathbf{\text{s}})=\widehat{\mathbf{\text{V}}}(\mathbf{\text{p}},\mathbf{\text{s}})$. Therefore:

$${\mathbf{\text{P}}}^T\mathbf{\text{V}}\mathbf{\text{S}}=\mathbf{\text{V}}(\mathbf{\text{p}},\mathbf{\text{s}})=\widehat{\mathbf{\text{V}}}(\mathbf{\text{p}},\mathbf{\text{s}})=\mathbf{\text{Q}}(\mathbf{\text{p}},:)\mathbf{\text{Z}}(:,\mathbf{\text{s}})={\mathbf{\text{P}}}^T\mathbf{\text{Q}}\mathbf{\text{Z}}\mathbf{\text{S}}.$$

Using this result in equation (2.22) yields:

$${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}=\mathbf{\text{Q}}{({\mathbf{\text{P}}}^T\mathbf{\text{Q}})}^{-1}{\mathbf{\text{P}}}^T\mathbf{\text{Q}}\mathbf{\text{Z}}\mathbf{\text{S}}{(\mathbf{\text{Z}}\mathbf{\text{S}})}^{-1}\mathbf{\text{Z}}=\mathbf{\text{Q}}\mathbf{\text{Z}}=\widehat{\mathbf{\text{V}}}.$$

This result completes the proof.

In the following theorem, we show that the oblique projection error is bounded by an error factor multiplied by the maximum of orthogonal projection errors onto $\mathbf{\text{U}}$ or $\mathbf{\text{Y}}$. We follow a similar procedure to that used in [37], however, in [37] the CUR is computed based on orthogonal projections onto the selected columns and rows, whereas in the presented TDB-CUR algorithm, oblique projectors are used. Without loss of generality, we consider a generic oblique projection, where the indexing matrices $\mathbf{\text{P}}$ and $\mathbf{\text{S}}$ are of size $n\times r^{\prime }$ and $s\times r^{\prime }$, respectively, and in general, $r^{\prime }\ge r$ (see definition 2.4 for details). In the following, we use the second norm ($||\sim ||\equiv ||\sim |{|}_2$).

Theorem 2.8.

Let ${\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}$ and ${\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}$ be oblique projectors according to definition 2.4 and let $\mathbf{U}\in {\mathbb{R}}^{n\times r}$ and $\mathbf{Y}\in {\mathbb{R}}^{s\times r}$ be a set of orthonormal matrices, i.e. ${\mathbf{U}}^T\mathbf{U}=\mathbf{I}$ and ${\mathbf{Y}}^T\mathbf{Y}=\mathbf{I}$. Let ${ϵ}_f\ge 0$ be given by: ${ϵ}_f=\text{min}\{{\eta }_p(1+{\eta }_s),{\eta }_s(1+{\eta }_p)\}-1$, where ${\eta }_p=||{({\mathbf{P}}^{\mathbf{T}}\mathbf{U})}^{†}||$ and ${\eta }_s=||{({\mathbf{S}}^T\mathbf{Y})}^{†}||$ and ${\hat{\sigma }}_{r+1}=\text{max}\{||\mathbf{V}-\mathbf{U}{\mathbf{U}}^T\mathbf{V}||,||\mathbf{V}-\mathbf{V}\mathbf{Y}{\mathbf{Y}}^T||\}$. Then the error of the oblique projection is bounded by

$$||\mathbf{\text{V}}-{\mathbf{\text{P}}}_{\mathbf{U}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{Y}}^{\mathrm{\angle }}||\le (1+{ϵ}_f){\hat{\sigma }}_{r+1}.$$2.23

Proof.

First note that ${ϵ}_f\ge 0$ because ${\eta }_p\ge 1$ and ${\eta }_q\ge 1$. The error matrix can be written as

$$\begin{array}{rl}\mathbf{\text{V}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}& \hspace{0.17em}=(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{V}}+{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}\\ & \hspace{0.17em}=(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{V}}+{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})\end{array}$$

where $\mathbf{\text{I}}$ is the identity matrix of appropriate size. Also, $\mathcal{P}\mathbf{\text{U}}=\mathbf{\text{U}}{({\mathbf{\text{P}}}^T\mathbf{\text{U}})}^{-1}{\mathbf{\text{P}}}^T\mathbf{\text{U}}=\mathbf{\text{U}}$. Therefore,$(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{U}}=\mathbf{\text{0}}$. Similarly, ${\mathbf{\text{Y}}}^T(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})=\mathbf{\text{0}}$. Therefore,

$$\begin{array}{rl}||\mathbf{\text{V}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||& \hspace{0.17em}\le ||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{V}}||+||{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})||\\ & \hspace{0.17em}=||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})(\mathbf{\text{V}}-\mathbf{\text{U}}{\mathbf{\text{U}}}^T\mathbf{\text{V}})||+||{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}(\mathbf{\text{V}}-\mathbf{\text{V}}\mathbf{\text{Y}}{\mathbf{\text{Y}}}^T)(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})||\\ & \hspace{0.17em}\le (||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})||+||{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}||||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})||){\hat{\sigma }}_{r+1}\\ & \hspace{0.17em}={\eta }_p(1+{\eta }_s){\hat{\sigma }}_{r+1}.\end{array}$$

In the above inequality, we have made use of the fact that $||\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}||=||{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}||={\eta }_p$ and $||\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||=||{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||={\eta }_s$ as long as the projectors are neither null nor the identity [46]. In the second line of the above inequality, we have made use of $(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{U}}=\mathbf{\text{0}}$ and ${\mathbf{\text{Y}}}^T(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})=\mathbf{\text{0}}$. Similarly, it is possible to express the error matrix as

$$\begin{array}{rl}\mathbf{\text{V}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}& \hspace{0.17em}=\mathbf{\text{V}}(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})+\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}\\ & \hspace{0.17em}=\mathbf{\text{V}}(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})+(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}.\end{array}$$

Therefore, another error bound can be obtained as

$$\begin{array}{rl}||\mathbf{\text{V}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||& \hspace{0.17em}\le ||\mathbf{\text{V}}(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})||+||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})\mathbf{\text{V}}{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||\\ & \hspace{0.17em}=||(\mathbf{\text{V}}-\mathbf{\text{V}}\mathbf{\text{Y}}{\mathbf{\text{Y}}}^T)(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})||+||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }})(\mathbf{\text{V}}-\mathbf{\text{U}}{\mathbf{\text{U}}}^T\mathbf{\text{V}}){\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||\\ & \hspace{0.17em}\le (||(\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }})||+||\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{U}}}^{\mathrm{\angle }}||||{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||){\hat{\sigma }}_{r+1}\\ & \hspace{0.17em}={\eta }_s(1+{\eta }_p){\hat{\sigma }}_{r+1},\end{array}$$

where $||\mathbf{\text{I}}-{\mathbf{\text{P}}}_{\mathbf{\text{Y}}}^{\mathrm{\angle }}||={\eta }_s$ is used. Combining the above two inequalities yields inequality (2.23).

In the above error bound, when $\mathbf{\text{U}}$ and $\mathbf{\text{Y}}$ are the $r$ most dominant exact left and right singular vectors of $\mathbf{\text{V}}$, then ${\hat{\sigma }}_{r+1}={\sigma }_{r+1}$, where ${\sigma }_{r+1}$ is the $r+1$th singular value of $\mathbf{\text{V}}$, since

In that case, $1+{ϵ}_f$ is the error factor of the CUR decomposition when compared against the optimal rank-$r$ reduction error obtained by SVD. As demonstrated in our numerical experiments, the TDB-CUR algorithm closely approximates the rank-$r$ SVD approximation of $\mathbf{\text{V}}$.