# A DEIM Tucker Tensor Cross Algorithm and its Application to Dynamical Low-
Rank Approximation

Behzad Ghahremani, Hessam Babaee∗ Department of Mechanical Engineering and
Materials Science, University of Pittsburgh, 3700 O’Hara Street, Pittsburgh,
PA, 15213, USA
* Corresponding Author<EMAIL_ADDRESS>
###### Abstract

We introduce a Tucker tensor cross approximation method that constructs a low-
rank representation of a $d$-dimensional tensor by sparsely sampling its
fibers. These fibers are selected using the discrete empirical interpolation
method (DEIM). Our proposed algorithm is referred to as DEIM fiber sampling
(DEIM-FS). For a rank-$r$ approximation of an $\mathcal{O}(N^{d})$ tensor,
DEIM-FS requires access to only $dNr^{d-1}$ tensor entries, a requirement that
scales linearly with the tensor size along each mode. We demonstrate that
DEIM-FS achieves an approximation accuracy close to the Tucker-tensor
approximation obtained via higher-order singular value decomposition at a
significantly reduced cost. We also present DEIM-FS (iterative) that does not
require access to singular vectors of the target tensor unfolding and can be
viewed as a black-box Tucker tensor algorithm. We employ DEIM-FS to reduce the
computational cost associated with solving nonlinear tensor differential
equations (TDEs) using dynamical low-rank approximation (DLRA). The
computational cost of solving DLRA equations can become prohibitive when the
exact rank of the right-hand side tensor is large. This issue arises in many
TDEs, especially in cases involving non-polynomial nonlinearities, where the
right-hand side tensor has full rank. This necessitates the storage and
computation of tensors of size $\mathcal{O}(N^{d})$. We show that DEIM-FS
results in significant computational savings for DLRA by constructing a low-
rank Tucker approximation of the right-hand side tensor on the fly. Another
advantage of using DEIM-FS is to significantly simplify the implementation of
DLRA equations, irrespective of the type of TDEs. We demonstrate the
efficiency of the algorithm through several examples including solving high-
dimensional partial differential equations.

###### keywords:

Cross approximation, dynamical low-rank approximation, time-dependent bases,
Tucker tensor

††journal: Computer Methods in Applied Mechanics and Engineering

A novel cross Tucker tensor algorithm is introduced aimed at constructing
near-optimal low-rank Tucker tensor models by sampling a relatively small
number of elements from the target tensor.

The algorithm strategically samples fibers using the discrete empirical
interpolation method.

An iterative cross algorithm is introduced, which operates without the need
for accessing singular vectors of the tensor unfolding. It can be regarded as
a black-box algorithm for constructing Tucker tensor models.

The cross algorithm is utilized to decrease both the computational cost and
memory requirements when solving high-dimensional nonlinear tensor
differential equations using dynamical low-rank approximation.

##  1 Introduction

Multi-dimensional tensors play a critical role in many applications in science
and engineering [1]. However, performing computational tasks involving high-
dimensional tensors or even storing them suffers from the curse of
dimensionality: for a tensor of size $N_{1}\times N_{2}\times\dots\times
N_{d}$, the number of its elements increases exponentially as $d$ grows,
resulting in $\mathcal{O}(N^{d})$ elements, where $N=N_{i}$, $i=1,\dots,d$ is
assumed. Various tensor low-rank approximations have been developed to
mitigate this issue by leveraging multi-dimensional correlations [2]. These
dimension reduction techniques aim to decrease the total number of tensor
elements while allowing a controllable loss of accuracy. Some of the most
common tensor low-rank approximation schemes include Tucker tensor
decomposition [3], CANDECOMP/PARAFAC (CP) [4], hierarchical Tucker tensor
decomposition [5], and tensor train decomposition [6].

The Tucker tensor low-rank approximation reduces the total number of elements
of a $d$-dimensional tensor from $N^{d}$ to $r^{d}+rdN$, where $r\ll N$
represents the rank of the unfolded tensors along each mode. The Tucker tensor
low-rank approximation is the building block for many tensor low-rank
approximations and it finds applications in various fields, including computer
vision, deep neural networks, data mining, numerical analysis, neuroscience,
and more. For an excellent review of these applications, refer to [7].

One of the applications of the Tucker tensor decomposition is to reduce the
computational cost of solving multi-dimensional partial differential equations
(PDEs) using dynamical low-rank approximation (DLRA) [8]. Discretizing these
PDEs in all dimensions except time results in tensor differential equations
(TDEs) in the form of $d\mathcal{V}/dt=\mathscr{F}(\mathcal{V})$, where
$\mathcal{V}\in\mathbb{R}^{N_{1}\times N_{2}\times\dots\times N_{d}}$ and
$\mathscr{F}(\mathcal{V})$ is the right had side tensor of the same size as
$\mathcal{V}$. Example applications include the Schrodinger equation [9], the
Fokker-Planck equation [10], the Boltzmann transport equation [11], and the
Hamilton-Jacobi-Bellman equations [12], among others. Solving these TDEs, even
in moderate dimensions ($3<d<7$), using traditional numerical methods such as
finite difference and finite element methods, encounters the issue of the
curse of dimensionality [13]. DLRA mitigates this issue by solving TDEs on the
manifold of low-rank Tucker tensors, in which explicit evolution equations for
a low-rank Tucker model are obtained. DLRA for Tucker tensor form is the
extension of DLRA for matrix differential equations [14], which has found many
diverse applications such as kinetics [15, 16], linear sensitivity analyses
[17] and species transport equations in turbulent combustion [18].

Reducing the computational cost of solving DLRA equations is the primary
motivating application for the developments presented in this work. The DLRA
evolution equations involve the $\mathscr{F}(\mathcal{V})$ tensor. When the
exact rank of $\mathscr{F}(\mathcal{V})$ is high, the computational cost of
DLRA increases. This occurs for linear TDEs with a large number of right-hand
side terms and for TDEs with high-order polynomial nonlinearities. In the case
of TDEs with nonpolynomial nonlinearities, such as exponential or fractional
nonlinearity, the computational cost of solving the DLRA evolution equations
exceeds that of the full-order model (FOM). This increased cost arises in TDEs
with nonpolynomial nonlinearity because $\mathscr{F}(\mathcal{V})$ is a full-
rank tensor even when $\mathcal{V}$ is low rank. Explicit formation of this
tensor necessitates memory and floating-point operation costs similar to those
of the FOM, i.e., $\mathcal{O}(N^{d})$.

One solution to mitigate the issue of the cost of DLRA for nonlinear TDEs is
to approximate $\mathscr{F}(\mathcal{V})$ with another Tucker tensor
approximation. This idea was recently successfully employed for solving
nonlinear matrix differential equations (MDEs) on low-rank matrix manifolds
[19, 20]. As demonstrated in Section 2.8, once $\mathscr{F}(\mathcal{V})$ is
approximated in the Tucker form, it can be efficiently incorporated into the
DLRA evolution equations.

To maintain the computational advantages of DLRA, any competitive Tucker
tensor decomposition algorithm should satisfy stringent accuracy-versus-cost
criteria: (i) the algorithm should be fast—preferably depending linearly on
$N$ in terms of floating-point operations (flops) and memory requirements;
(ii) the algorithm should be accurate, as the error introduced by the Tucker
low-rank approximation of the right-hand side tensor should not exceed the
DLRA errors, which are typically very small, often within the range of
$\mathcal{O}(10^{-6})$ to $\mathcal{O}(10^{-10})$ in relative errors,
depending on the DLRA rank and the local temporal integration error. To this
end, we review various techniques for computing Tucker tensor approximation
and evaluate them based on the aforementioned criteria.

Determining the _best_ Tucker tensor low-rank approximation of a tensor lacks
a known closed solution. The higher-order singular value decomposition (HOSVD)
is a reliable approach for computing a near-optimal Tucker tensor
decomposition [21]. The computational cost of computing the Tucker tensor
decomposition using HOSVD scales at least linearly with the total number of
elements of a tensor, i.e., $\mathcal{O}(N^{d})$. For example, a six-
dimensional probability density function with $N=100$ grid points in each
dimension results in a tensor with one trillion ($10^{12}$) elements. This
cost is prohibitive for many applications — certainly for DLRA — when $N^{d}$
is very large. The computational expense stems from two primary sources:

  1. 1.

Computing the SVD: The HOSVD algorithm requires performing the SVD of large
matrices to compute the factor matrices. Specifically, it involves computing
the SVD of $d$ matrices of size $N\times N^{d-1}$ obtained by unfolding the
tensor along its $d$ modes. The core is then computed via an orthogonal
projection of the tensor onto factor matrices. The computational cost of
performing $d$ SVD scales with $\mathcal{O}(dN^{(d+1)})$.

  2. 2.

Data access: The HOSVD requires access to all elements of the tensor for
computing the factor matrices and the core tensor. When dealing with very
large values of $N^{d}$, it might be impossible to hold the entire tensor in
memory due to resource limitations. Consequently, multiple loading of the
entire data in chunks from the disk becomes necessary, which is significantly
slower compared to loading from memory. For DLRA, computing any tensor element
requires applying $\mathscr{F}(\sim)$ to $\mathcal{V}$, implying that data
access incurs flops costs in addition to the memory requirements.

Numerous algorithms have been proposed to tackle these issues. For reducing
the SVD cost, various approaches exist; for instance, randomized HOSVD [22],
higher-order interpolatory decomposition [23], and sequentially truncated
HOSVD [24]. However, all of these algorithms require access to all tensor
entries. In contrast, cross algorithms tackle both of these issues
simultaneously.

Cross algorithms are powerful techniques for constructing low-rank matrix and
tensor approximations. The first cross approximation was proposed for matrix
low-rank approximations [25]. Cross algorithms are also known as
_pseudoskeleton_ or _CUR decomposition_. The simplest cross algorithm is an
interpolatory CUR matrix decomposition, wherein a rank-$r$ matrix is
constructed by sampling only $r$ columns and rows of the matrix. The clear
advantage of the interpolatory CUR algorithm is that the low-rank
approximation does not require all the entries of the matrix. However, the
accuracy of the low-rank approximation obtained via CUR critically depends on
the selection of rows and columns. As shown in [25], the accuracy of the CUR
approximation depends on the determinant of the intersection submatrix, which
is referred to as _matrix volume_. In particular, the columns and rows should
be selected such that the matrix volume is maximized. Since this selection
problem is NP-hard, several heuristic algorithms have been proposed including
Maxvol [26, 27], Cross2D [28, 29], leverage score [30], and discrete empirical
interpolation method (DEIM) [31] algorithms.

The cross algorithms have been extended to various tensor low-rank
approximations, including Tucker tensor decomposition [32], tensor train
decomposition [33], and hierarchical Tucker tensor decomposition [34]. Similar
to the interpolatory CUR algorithm, these algorithms do not require access to
all entries of the tensor. We refer the readers to [35] for an excellent
review of various tensor cross approximations.

The cross approximation, introduced in [32], is known as fiber sampling Tucker
decomposition (FSTD). FSTD is an elegant extension of matrix CUR decomposition
to Tucker tensor approximations. However, as we demonstrate, FSTD becomes ill-
conditioned as rank increases and therefore, the approximation error cannot be
reduced to machine precision. The issue of ill-conditioning is similarly
encountered in matrix CUR decompositions, where the intersection submatrix
becomes singular as rank increases. This issue is mitigated by applying the QR
factorization to either the selected columns or rows [33]. As a result, FSTD
is not well-suited for DLRA due to the stringent accuracy-versus-cost
requirements that DLRA demands.

In this paper, we present a novel Tucker cross algorithm that addresses the
aforementioned challenges. In particular, the contributions of this paper are
the following:

  1. 1.

We introduce DEIM-FS — a Tucker tensor cross algorithm by sampling $r^{d-1}$
fibers along each mode of a tensor. Therefore, it requires access to
$dr^{d-1}N$ entries of the tensor. The fiber selection (FS) is guided by the
DEIM algorithm. The DEIM algorithm requires access to the exact or approximate
singular vectors of the tensor unfolding along each mode. As we will discuss
in this paper, DLRA is one such application, where the approximate singular
vectors are available at no additional cost. We demonstrate that the low-rank
approximation error of the DEIM-FS algorithm is comparable to that achieved by
HOSVD.

  2. 2.

We present DEIM-FS (iterative) for the problems where the exact or approximate
singular vectors are not available. DEIM-FS (iterative) starts with a random
guess for fibers and iteratively applies DEIM-FS until convergence. In
practice, a small number of iterations are needed. As a result, DEIM-FS
(iterative) can be regarded as a _black-box_ Tucker cross algorithm.

  3. 3.

We present DLRA-DEIM-FS to reduce the computational cost of solving DLRA
evolution equations for nonlinear TDEs. DLRA-DEIM-FS constructs a low-rank
Tucker tensor approximation of the right-hand side of the TDE by applying
DEIM-FS to $\mathscr{F}(\mathcal{V})$. We augment the DEIM-FS with _rank-
adaptivity_ , where the Tucker tensor rank is adjusted on the fly to meet an
error threshold criterion.

The paper is organized as follows. The methodologies are discussed in Section
2, demonstrations and results are presented in Section 3, and the conclusions
follow in Section 4.

##  2 Methodology

###  2.1 Definitions and notations

We first introduce the notation used for vectors, matrices, and tensors.
Vectors are denoted in bold lowercase letters (e.g. $\mathbf{a}$), matrices
are denoted by bold uppercase letters (e.g. $\mathbf{A}$), and tensors by
uppercase calligraphic letters (e.g. $\mathcal{F}$). The symbol
$\bigtimes_{n}$ is used to denote the $n$-mode product. The $n$-mode product
of a tensor $\mathcal{F}\in\mathbb{R}^{N_{1}\times N_{2}\times\dots\times
N_{d}}$ with a matrix $\mathbf{M}\in\mathbb{R}^{J\times N_{n}}$ is obtained by
$\mathcal{F}\bigtimes_{n}\mathbf{M}$ and is of size $N_{1}\times\dots\times
N_{n-1}\times J\times N_{n+1}\times\dots\times N_{d}$. We denote the unfolding
of tensor $\mathcal{F}$ along its $n$-th mode with $\mathcal{F}_{(n)}$.
Unfolding a tensor involves reshaping its elements in such a way that results
in a matrix instead of a tensor. For instance, a tensor of size $3\times
5\times 6$ can be unfolded along the second axis as a matrix of size $5\times
18$ [7]. The Frobenius norm of a tensor is shown by $\|\mathcal{F}\|_{F}$ and
is defined as:

$\|\mathcal{F}\|_{F}=\sqrt{\sum_{n_{1}=1}^{N_{1}}\sum_{n_{2}=1}^{N_{2}}\dots\sum_{n_{d}=1}^{N_{d}}a^{2}_{n_{1}n_{2}\dots
n_{d}}},$ (1)

where $a_{n_{1}n_{2}\dots n_{d}}$ are the entries of tensor $\mathcal{F}$. We
use typewriter font to denote algorithms, e.g., SVD or DEIM-FS. We use the
MATLAB indexing notation where $\mathbf{A}(\mathbf{p},:)$ selects all columns
at the $\mathbf{p}$ rows and $\mathbf{A}(:,\mathbf{p})$ selects all rows at
the $\mathbf{p}$ columns of matrix $\mathbf{A}$ and
$\mathbf{p}=[p_{1},p_{2},\dots,p_{q}]$ is the integer vector containing the
selected indices. We also use MATLAB notation for computing the SVD of a
matrix. For example, consider for $\mathbf{A}\in\mathbb{R}^{m\times n}$. Then
$[\mathbf{U},\mathbf{\Sigma},\mathbf{V}]=\texttt{SVD}(\mathbf{A},r)$ means
computing the SVD of $\mathbf{A}$ and truncating at rank $r$, where $r<m$ and
$r<n$, $\mathbf{U}\in\mathbb{R}^{n\times r}$ is the matrix of left singular
vectors, $\mathbf{\Sigma}\in\mathbb{R}^{r\times r}$ is the matrix of singular
values and $\mathbf{V}\in\mathbb{R}^{m\times r}$ is the matrix of right
singular vectors. If any of the singular matrices are not needed, the symbol
($\sim$) is used. For example,
$[\mathbf{U},\sim,\sim]=\texttt{SVD}(\mathbf{A},r)$ returns only the first $r$
left singular vectors.

###  2.2 Dynamical low-rank approximation for Tucker tensors

As mentioned in the Introduction, DLRA is the primary motivation for the
developments in this paper. The DLRA formulation for Tucker tensors is used to
solve high-dimensional PDEs on the manifold of low-rank Tucker tensors [8]. In
the following, we briefly introduce DLRA and explain the computational cost
issues related to DLRA.

We consider a general PDE given by:

$\frac{\partial v(\mathbf{x},t)}{\partial
t}=f\bigl{(}v(\mathbf{x},t)\bigr{)},$ (2)

augmented with appropriate initial and boundary conditions. Here
$\mathbf{x}=(x_{1},x_{2},\dots,x_{d})\in\mathbb{R}^{d}$, $t$ is time, $f$ is a
general nonlinear differential operator, and $d$ is the dimension of the
problem. We consider the differential operators in $\mathbf{x}$ being
discretized using a method of line. We discretize the differential operators
of Eq. (2) in $\mathbf{x}$ using a method of lines, which results in the
following TDE:

$\frac{d\mathcal{V}}{dt}=\mathscr{F}\bigl{(}\mathcal{V}\bigr{)},$ (3)

where $\mathcal{V}(t)\in\mathbb{R}^{N_{1}\times N_{2}\times\dots\times N_{d}}$
is the solution tensor and $\mathscr{F}$ is the discrete representation of
$f$. Here, $N_{1},N_{2},\dots,N_{d}$ are the number of the discretized points
along each mode of the tensor. We refer to Eq. 3 as the FOM.

Discretizing Eq. (2) via classical methods such as finite-difference and
spectral methods results in a system where the degrees of freedom grow
exponentially fast as the number of dimensions grows. One approach to mitigate
this issue is to solve Eq. (3) on a manifold of low-rank Tucker tensors. To
this end, consider a low-rank Tucker tensor approximation of $\mathcal{V}$ [8]
as shown below:

$\mathcal{V}(t)\approx\hat{\mathcal{V}}(t)=\mathcal{S}(t)\times_{1}\mathbf{U}^{(1)}(t)\times_{2}\mathbf{U}^{(2)}(t)\dots\times_{d}\mathbf{U}^{(d)}(t),$
(4)

where $\times_{n}$ is tensor mode product,
$\mathcal{S}\in\mathbb{R}^{r_{1}\times r_{2}\times\dots\times r_{d}}$ is the
core tensor, $\mathbf{U}^{(i)}\in\mathbb{R}^{N_{i}\times r_{i}}$ are the
orthonormal time-dependent bases or factor matrices along the corresponding
mode of the tensor, i.e., ${\mathbf{U}^{(i)}}^{T}\mathbf{U}^{(i)}=\mathbf{I}$,
where $\mathbf{I}$ is the identity matrix and $r_{i}\ll N_{i},i=1,2,\dots,d$
are the rank along each tensor mode. Substituting the low-rank approximation
given by Eq. (4) into Eq. (3) results in a residual equal to:

$\mathcal{R}(\dot{\mathcal{S}},\dot{\mathbf{U}}^{(1)},\dot{\mathbf{U}}^{(2)},\dots,\dot{\mathbf{U}}^{(d)})=\\\
\Biggl{\|}\frac{d~{}\bigl{(}\mathcal{S}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}~{}\dots\times_{d}\mathbf{U}^{(d)}\bigr{)}}{dt}-\mathscr{F}\bigl{(}\mathcal{S}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}~{}\dots\times_{d}\mathbf{U}^{(d)}\bigr{)}\Bigg{\|}_{F}^{2}.$
(5)

The evolution equations for orthonormal bases and the core tensor are obtained
by minimizing the above residual subject to the orthonormality constraints of
the bases, which results in the following evolution equations for the core and
the factor matrices [8]:

$\displaystyle\dot{\mathcal{S}}$
$\displaystyle=\mathcal{F}~{}\bigtimes_{i=1}^{d}~{}{\mathbf{U}^{(i)}}^{T},$
(6a) $\displaystyle\dot{\mathbf{U}}^{(i)}$
$\displaystyle=\Bigl{(}\mathbf{I}-\mathbf{U}^{(i)}\mathbf{U}^{(i)^{T}}\Bigr{)}~{}\Bigl{[}\mathcal{F}\bigtimes_{k\neq
i}{\mathbf{U}^{(k)}}^{T}\Bigr{]}_{(i)}~{}\mathcal{S}_{(i)}^{\dagger},$ (6b)

where $\mathbf{I}$ is the identity matrix,
$\mathcal{S}_{(i)}^{\dagger}=\mathcal{S}_{(i)}^{T}\bigl{(}\mathcal{S}_{(i)}\mathcal{S}_{(i)}^{T})^{-1}$
is the pseudo-inverse, $\mathcal{F}\in\mathbb{R}^{N_{1}\times
N_{2}\times\dots\times N_{d}}$ is a tensor defined as
$\mathcal{F}=\mathscr{F}\bigl{(}\mathcal{S}\times_{1}\mathbf{U}^{(1)}\times_{2}\mathbf{U}^{(2)}~{}\dots\times_{d}\mathbf{U}^{(d)}\bigr{)}$.
For further details on the derivation of the evolution equations, refer to
[8]. Eqs. (6a) and (6b) represent the DLRA evolution equations in the Tucker
tensor form. For recent developments related to a stable time integration
scheme and rank adaptivity of DLRA equations, see [36, 37].

The computational advantage of the DLRA evolution equations over the FOM (Eq.
3) is that in Eqs. (6a) and (6b), the solution is sought in terms of the
factor matrices ($\mathbf{U}^{(i)}$) and the core tensor ($\mathcal{S}$)
instead of the full-dimensional tensor ($\mathcal{V}(t)$). The memory
requirement for storing $\\{\mathcal{S},\mathbf{U}^{(i)}\\}$ is
$\mathcal{O}(r^{d})+\mathcal{O}(rdN)$. However, solving $\mathcal{V}(t)$ using
the FOM requires $\mathcal{O}(N^{d})$ memory. For simplicity in computational
complexity analysis, we assume $r=r_{1}=r_{2}=\ldots=r_{d}$ and
$N=N_{1}=N_{2}=\ldots=N_{d}$.

The computational savings of the DLRA in memory and floating-point operations
(flops) are lost when dealing with TDEs featuring general nonlinearity. In
these cases, the right-hand side (RHS) tensor $\mathcal{F}$ is full rank,
necessitating its computation and storage in memory. The computational cost of
computing $\mathcal{F}$ scales at $\mathcal{O}(N^{d})$, mirroring the
computational complexity of solving the FOM. Even in linear TDEs, a
substantial number of terms on the right-hand side may result in a large exact
rank for $\mathcal{F}$. Therefore, the issue of high computational cost is not
only limited to general nonlinearity, but the issue arises for problems in
which the exact rank of $\mathcal{F}$ is large. Consequently, evaluating
$\mathcal{F}$ remains the primary computational bottleneck for DLRA. To
address this challenge, we introduce an efficient and innovative cross
algorithm that constructs a low-rank Tucker tensor approximation of
$\mathcal{F}$ by selectively sampling a few fibers of $\mathcal{F}$ using the
DEIM algorithm.

In the next section, we first provide the utility of the DEIM algorithm for
nonlinear reduced-order modeling. The cross algorithm presented in this paper
is inspired by our previous work [19], where we developed a CUR algorithm for
DLRA of nonlinear MDEs. In Section 2.4, we provide a brief overview of the
algorithm presented in [19].

###  2.3 DEIM for low-rank approximation of vector differential equations

In reduced-order modeling based on proper orthogonal decomposition (POD),
analogous challenges arise due to general nonlinearity. Consider the FOM given
by: $d\mathbf{v}/dt=f(\mathbf{v})$, where $\mathbf{v}\in\mathbb{R}^{N}$
represents the state vector, and
$f(\mathbf{v}):\mathbb{R}^{N}\rightarrow\mathbb{R}^{N}$ is the RHS vector. Let
$\mathbf{U}\in\mathbb{R}^{N\times r}$ denote the matrix of orthonormal POD
modes, and $\mathbf{y}\in\mathbb{R}^{r}$ be the POD coefficients, such that
$\hat{\mathbf{v}}=\mathbf{U}\mathbf{y}\in\mathbb{R}^{N}$ is the POD
approximation of $\mathbf{v}$. Here, $r\ll N$ represents the number of POD
modes. The ROM is obtained via Galerkin projection:
$d\mathbf{y}/dt=\mathbf{U}^{T}f(\mathbf{U}\mathbf{y})$. If $f(\sim)$ involves
a polynomial nonlinearity of degree $p$, it becomes feasible to compute
$\mathbf{U}^{T}f(\mathbf{U}\mathbf{y})$ with a computational cost of at least
$\mathcal{O}(r^{p})$. This implies the potential avoidance of forming the
vector $f(\mathbf{U}\mathbf{y})\in\mathbb{R}^{N}$, thus preventing the
computational cost of solving the FOM from scaling with the FOM size ($N$)
[38]. However, in scenarios where $f(\sim)$ has high-order polynomial
nonlinearity (i.e., a large $p$), the computational cost of solving the ROM
can become considerable. Furthermore, in cases where $f(\sim)$ has non-
polynomial nonlinearity, the explicit formation of the vector
$f(\mathbf{U}\mathbf{y})$ becomes necessary. This results in the loss of
computational savings offered by the POD-ROM, as computing
$f(\mathbf{U}\mathbf{y})$ requires $\mathcal{O}(N)$ operations—equivalent to
solving the FOM.

One computationally efficient remedy is to interpolate the
$f(\mathbf{U}\mathbf{y})$ onto a low-rank basis using the DEIM algorithm [39].
This involves sampling $f(\mathbf{U}\mathbf{y})$ at only a few strategically
selected points. To explain this algorithm, let
$\mathbf{f}=f(\mathbf{U}\mathbf{y})$ and $\mathbf{U}_{f}\in\mathbb{R}^{N\times
r_{f}}$ be a low-rank basis for the vector $\mathbf{f}$, where $r_{f}$ is the
number of POD modes for the vector $\mathbf{f}$. The basis $\mathbf{U}_{f}$ is
computed in the offline stage as the left singular vectors of the RHS snapshot
matrix. The DEIM algorithm [39, Algorithm 1] yields a set of near-optimal
sampling points for interpolation of vector $\mathbf{f}$ onto
$\mathbf{U}_{f}$:

$\mathbf{p}=\texttt{DEIM}(\mathbf{U}_{f}),$ (7)

where $\mathbf{p}=[p_{1},p_{2},\dots,p_{r}]$ is the vector of sampling point
indices. For convenience, the DEIM algorithm is provided in Appendix 1. The
vector $\mathbf{f}$ can be interpolated onto $\mathbf{U}_{f}$ using:

$\hat{\mathbf{f}}=\mathbf{U}_{f}\mathbf{U}_{f}(\mathbf{p},:)^{-1}\mathbf{f}(\mathbf{p}).$
(8)

It is easy to verify that $\hat{\mathbf{f}}$ and $\mathbf{f}$ are equal to
each other at interpolation points, i.e.,
$\hat{\mathbf{f}}(\mathbf{p})=\mathbf{f}(\mathbf{p})$. Incorporating Eq. 8
into the POD-ROM results in:
$d\mathbf{y}/dt=\mathbf{U}^{T}\mathbf{U}_{f}\mathbf{U}_{f}(\mathbf{p},:)^{-1}\mathbf{f}(\mathbf{p})$.
The small matrix
$\mathbf{U}^{T}\mathbf{U}_{f}\mathbf{U}_{f}(\mathbf{p},:)^{-1}\in\mathbb{R}^{r\times
r_{f}}$ can be computed and stored in the offline stage. The key advantage of
using the DEIM algorithm in POD-ROM is that it requires evaluating
$f(\mathbf{U}\mathbf{y})$ at a only small number of points ($r_{f}$)-
regardless of the type of nonlinearity of $f(\sim)$. As shown in [39], the
approximation error of $\hat{\mathbf{f}}$ is bounded by the best approximation
error by a magnification factor:

$\|\hat{\mathbf{f}}-\mathbf{f}\|\leq\eta\|(\mathbf{I}-\mathbf{U}_{f}^{T}\mathbf{U}_{f})\mathbf{f}\|$
(9)

where $\eta=\|\mathbf{U}_{f}(\mathbf{p},:)^{-1}\|$ is the magnification factor
and $\|(\mathbf{I}-\mathbf{U}_{f}^{T}\mathbf{U}_{f})\mathbf{f}\|$ is the
optimal error of approximating $\mathbf{f}$ in the span of $\mathbf{U}_{f}$,
which is obtained via the orthogonal projection of $\mathbf{f}$ onto
$\mathbf{U}_{f}$. The DEIM algorithm is designed to minimize $\eta$ using a
greedy approach.

###  2.4 DEIM for low-rank approximation of matrix differential equations

The DEIM algorithm, developed for the low-rank approximation of MDEs in [19],
deals with an MDE expressed as $d\mathbf{V}/dt=F(\mathbf{V})$, where
$\mathbf{V}\in\mathbb{R}^{N_{1}\times N_{2}}$ and
$F(\mathbf{V}):\mathbb{R}^{N_{1}\times
N_{2}}\rightarrow\mathbb{R}^{N_{1}\times N_{2}}$. The DLRA aims to find
solutions for the above MDE constrained to the manifold of rank $r$ matrices,
where $r\ll N_{1}$ and $r\ll N_{2}$. DLRA for MDE can be derived as a specific
case of DLRA for TDEs (refer to Eq. 6a-6b) by setting $d=2$ and
$r=r_{1}=r_{2}$. When the exact rank of $F(\mathbf{V})$ is large, the
computational cost of DLRA increases. For example, when $F(\sim)$ exhibits
non-polynomial nonlinearity, $F(\mathbf{V})$ is full rank, and the explicit
formation of $F(\mathbf{V})$ is needed. This results in the loss of
computational savings provided by DLRA.

The algorithm proposed in [19, Algorithm 1] constructs a low-rank
approximation of the matrix $\mathbf{F}=F(\mathbf{V})$ by selectively sampling
$r$ columns and $r$ rows from $\mathbf{F}$. This involves:

  1. 1.

Sampling the columns of the RHS matrix:
$\mathbf{F}(:,\mathbf{p}_{2})\in\mathbb{R}^{N_{1}\times r}$, where
$\mathbf{p}_{2}\in\mathbb{I}^{r}$ denotes the integer vector containing the
column indices.

  2. 2.

Conducting the QR decomposition of matrix
$\mathbf{F}(:,\mathbf{p}_{2})=\mathbf{Q}\mathbf{R}$ to construct an
orthonormal basis $\mathbf{Q}$ for the columns of matrix $\mathbf{F}$.

  3. 3.

Computing the rows of the RHS matrix
$\mathbf{F}(\mathbf{p}_{1},:)\in\mathbb{R}^{N_{1}\times r}$, where
$\mathbf{p}_{1}\in\mathbb{I}^{r}$ represents the integer row indices.

  4. 4.

Interpolating each column of matrix $\mathbf{F}$ onto the basis $\mathbf{Q}$
using the computed values at rows indexed by $\mathbf{p}_{1}$:
$\mathbf{F}\approx\hat{\mathbf{F}}=\mathbf{Q}\mathbf{Q}(\mathbf{p}_{1},:)^{-1}\mathbf{F}(\mathbf{p}_{1},:)$.

In the above algorithm, the row and column indices $\mathbf{p}_{1}$ and
$\mathbf{p}_{2}$ are determined using the DEIM algorithm:
$\mathbf{p}_{1}=\texttt{DEIM}(\mathbf{U}^{(1)}_{F})$ and
$\mathbf{p}_{2}=\texttt{DEIM}(\mathbf{U}^{(2)}_{F})$, where
$\mathbf{U}^{(1)}_{F}$ and $\mathbf{U}^{(2)}_{F}$ are the left and right
singular vectors of matrix $\hat{\mathbf{F}}$.

###  2.5 DEIM for Tucker tensor cross approximation

Cross tensor approximations are extensions of matrix CUR approximation
techniques, offering a practical approach for efficiently estimating low-rank
tensors [35]. In this paper, a novel cross tensor approximation technique is
proposed. The main steps of the proposed methodology are summarized in
Algorithm 1. We refer to our algorithm as DEIM fiber sampling (DEIM-FS). This
algorithm is presented for a three-dimensional tensor for simplicity, but it
can easily be generalized to higher-dimensional tensors. The algorithm
constructs a low-rank Tucker tensor approximation of $\mathcal{F}$ by only
sampling a few fibers along each mode of the tensor.

To this end, consider a Tucker low-rank approximation of $\mathcal{F}$ given
by:

$\mathcal{F}\approx\mathcal{S}_{\mathcal{F}}\times_{1}{\mathbf{U}^{(1)}_{\mathcal{F}}}\times_{2}{\mathbf{U}^{(2)}_{\mathcal{F}}}\times_{3}{\mathbf{U}^{(3)}_{\mathcal{F}}},$
(10)

where $\mathcal{F}\in\mathbb{R}^{N_{1}\times N_{2}\times N_{3}}$,
$\mathcal{S}_{\mathcal{F}}\in\mathbb{R}^{{r_{\mathcal{F}}}_{1}\times{r_{\mathcal{F}}}_{2}\times{r_{\mathcal{F}}}_{3}}$,
${\mathbf{U}_{\mathcal{F}}}^{(1)}\in\mathbb{R}^{N_{1}\times{r_{\mathcal{F}}}_{1}}$,
${\mathbf{U}^{(2)}_{\mathcal{F}}}\in\mathbb{R}^{N_{2}\times{r_{\mathcal{F}}}_{2}}$,
${\mathbf{U}^{(3)}_{\mathcal{F}}}\in\mathbb{R}^{N_{3}\times{r_{\mathcal{F}}}_{3}}$,
and
$\mathbf{r}_{\mathcal{F}}=(r_{\mathcal{F}_{1}},r_{\mathcal{F}_{2}},r_{\mathcal{F}_{3}})$
is the multi-rank of the Tucker tensor decomposition. The number of selected
fibers along the first, second, and third modes are denoted by
$r^{\prime}_{\mathcal{F}_{1}},r^{\prime}_{\mathcal{F}_{2}},r^{\prime}_{\mathcal{F}_{3}}$,
respectively. As we will explain below (See Remark 1), the number of selected
fibers must be greater than or equal to the target Tucker low-rank, i.e.,
$r^{\prime}_{\mathcal{F}_{i}}\geq r_{\mathcal{F}_{i}}$. According to the
numerical results of Section 3.1, we demonstrate that
$r^{\prime}_{\mathcal{F}_{i}}=r_{{\mathcal{F}}_{i}}+2$ is a reasonable choice
and this choice is used in all demonstrations in this paper. The indices of
the selected fibers are shown by vectors
$\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3}$, where $\mathbf{p}_{i}$ is an
integer vector containing $r^{\prime}_{\mathcal{F}_{i}}$.

Figure 1: Schematic of the DEIM-FS cross algorithm for a 3D tensor. For
simplicity, we assume that all selected fibers are adjacent to each other.

The schematic of the algorithm is shown in Figure 1, where
$\mathcal{C}_{1}=\mathcal{F}(:,\mathbf{p}_{2},\mathbf{p}_{3})\in\mathbb{R}^{N_{1}\times
r^{\prime}_{\mathcal{F}_{2}}\times
r^{\prime}_{\mathcal{F}_{3}}},~{}\mathcal{C}_{2}=\mathcal{F}(\mathbf{p}_{1},:,\mathbf{p}_{3})\in\mathbb{R}^{N_{2}\times
r^{\prime}_{\mathcal{F}_{1}}\times
r^{\prime}_{\mathcal{F}_{3}}},~{}\mathcal{C}_{3}=\mathcal{F}(\mathbf{p}_{1},\mathbf{p}_{2},:)\in\mathbb{R}^{N_{3}\times
r^{\prime}_{\mathcal{F}_{1}}\times r^{\prime}_{\mathcal{F}_{2}}}$ are the sub-
tensors formed by clustering the selected fibers of $\mathcal{F}$ along each
direction. For DEIM-FS, we assume that the first
$r^{\prime}_{\mathcal{F}_{i}}$ left singular vectors of matrix
$\mathcal{F}_{(i)}$ or some close approximation of these vectors are known. We
denote these singular vectors with
$\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}\in\mathbb{R}^{N_{i}\times
r^{\prime}_{\mathcal{F}_{i}}}$, where $i=1,2,3$ denotes different unfolding.
In Section 2.6, we present DEIM-FS (iterative), where the matrices of left
singular vectors are not needed. In the following, we explain all the steps of
DEIM-FS.

####  2.5.1 Computing the factor matrices

The factor matrices are computed using the DEIM-selected fibers of tensor
$\mathcal{F}$. To this end, the DEIM algorithm is applied to the left singular
vectors of $\mathcal{F}_{(i)}$ matrices:

$\mathbf{p}_{i}=\texttt{DEIM}(\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}),\quad
i=1,2,3.$

This generates the fiber indices for each mode of tensor $\mathcal{F}$.

The DEIM-selected fibers are extracted from tensor $\mathcal{F}$ and the
resulting subtensors are unfolded along the corresponding modes to obtain the
following matrices:

$\mathbf{C}_{1}=\Bigl{(}\mathcal{F}(:,\mathbf{p}_{2},\mathbf{p}_{3})\Bigr{)}_{(1)},\quad\mathbf{C}_{2}=\Bigl{(}\mathcal{F}(\mathbf{p}_{1},:,\mathbf{p}_{3})\Bigr{)}_{(2)},\quad\mbox{and}\quad\mathbf{C}_{3}=\Bigl{(}\mathcal{F}(\mathbf{p}_{1},\mathbf{p}_{2},:)\Bigr{)}_{(3)},$

where $\mathbf{C}_{1}\in\mathbb{R}^{N_{1}\times
r^{\prime}_{\mathcal{F}_{2}}r^{\prime}_{\mathcal{F}_{3}}}$,
$\mathbf{C}_{2}\in\mathbb{R}^{N_{2}\times
r^{\prime}_{\mathcal{F}_{1}}r^{\prime}_{\mathcal{F}_{3}}}$, and
$\mathbf{C}_{3}\in\mathbb{R}^{N_{3}\times
r^{\prime}_{\mathcal{F}_{1}}r^{\prime}_{\mathcal{F}_{2}}}$. In the next step,
the factor matrices are computed along each mode of the Tucker tensor
decomposition. To this end, we perform the SVD of the $\mathbf{C}_{i}$
matrices and truncate at rank $r_{\mathcal{F}_{i}}$:

$[\mathbf{U}^{(i)}_{\mathcal{F}},\sim,\sim]=\texttt{SVD}(\mathbf{C}_{i},r_{\mathcal{F}_{i}}),\quad
i=1,2,3,$ (11)

where $\mathbf{U}^{(i)}_{\mathcal{F}}\in\mathbb{R}^{N_{i}\times
r_{\mathcal{F}_{i}}}$ is the matrix of left singular vectors. The rank
$r_{\mathcal{F}_{i}}$ can be either fixed a priori or determined adaptively
based on accuracy requirements as explained in Section 2.5.3. The
$\mathbf{U}^{(i)}_{\mathcal{F}}$ matrices computed in this step are the factor
matrices of the cross Tucker tensor decomposition.

####  2.5.2 Computing the core tensor

Since the factor matrices are already computed, the optimal core tensor may be
obtained via the _orthogonal projection_ of $\mathcal{F}$ onto the factor
matrices. However, the orthogonal projection would require all entries of
$\mathcal{F}$, which is undesirable. The intersection tensor is denoted with
$\mathcal{W}$ (see Figure 1) and it is equal to:

$\mathcal{W}=\mathcal{F}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})\in\mathbb{R}^{r^{\prime}_{\mathcal{F}_{1}}\times
r^{\prime}_{\mathcal{F}_{2}}\times r^{\prime}_{\mathcal{F}_{3}}}.$

Note that the entries of the intersection tensor are a subset of the already-
computed fibers in the previous step.

The core tensor is calculated such that the difference between the cross
Tucker tensor approximation and the actual values of the tensor is minimized
at the DEIM intersection tensor ($\mathcal{W}$). More specifically, we seek to
find a core tensor such that

$E_{\mathcal{W}}(\mathcal{S}_{\mathcal{F}})=\big{\|}\mathcal{W}-\mathcal{S}_{\mathcal{F}}\times_{1}~{}\mathbf{U}^{(1)}_{\mathcal{F}}(\mathbf{p}_{1},:)\times_{2}~{}\mathbf{U}^{(2)}_{\mathcal{F}}(\mathbf{p}_{2},:)\times_{3}~{}\mathbf{U}^{(3)}_{\mathcal{F}}(\mathbf{p}_{3},:)\big{\|}^{2}_{F},$

is minimized. The above minimization problem has a closed-form solution and it
is computed via the least squares solution as shown below:

$\mathcal{S}_{\mathcal{F}}=\mathcal{W}\times_{1}~{}\mathbf{U}^{(1)}_{\mathcal{F}}(\mathbf{p}_{1},:)^{\dagger}\times_{2}~{}\mathbf{U}^{(2)}_{\mathcal{F}}(\mathbf{p}_{2},:)^{\dagger}\times_{3}~{}\mathbf{U}^{(3)}_{\mathcal{F}}(\mathbf{p}_{3},:)^{\dagger}.$
(12)

The above procedure can be also viewed as an _oblique projection_ of
$\mathcal{F}$ onto the factor matrices similar to matrix CUR decompositions
[31, 20].

We denote the resulting Tucker tensor decomposition obtained via Algorithm 1
by:

$\hat{\mathcal{F}}=\mathcal{S}_{\mathcal{F}}\times_{1}{\mathbf{U}^{(1)}_{\mathcal{F}}}\times_{2}{\mathbf{U}^{(2)}_{\mathcal{F}}}\times_{3}{\mathbf{U}^{(3)}_{\mathcal{F}}}.$
(13)

If rank adaptivity is not desired, no further computation is required and the
output of the algorithm is the Tucker tensor decomposition with the user-
specified multi-rank of
$\mathbf{r}_{\mathcal{F}}=(r_{\mathcal{F}_{1}},r_{\mathcal{F}_{2}},r_{\mathcal{F}_{3}})$.

####  2.5.3 Rank adaptivity

In certain applications, there could be a specific requirement for a low-rank
approximation error. Consequently, adjusting the rank becomes necessary to
meet a user-specified criterion. To address this requirement, we demonstrate
that the cross tensor approximation algorithm can be made rank-adaptive with
minor modifications. The error criterion is application-dependent. We consider
$\epsilon_{i}$

$\epsilon_{i}=\dfrac{\min({\mathbf{\Sigma}_{\mathcal{F}}}_{i})}{\|{\mathbf{\Sigma}_{\mathcal{F}}}_{i}\|_{F}}\quad
i=1,2,3~{},$ (14)

as the error proxy, where $\mathbf{\Sigma}_{\mathcal{F}_{i}}$ is the matrix of
singular values of the unfolded core tensor
$(\mathcal{S}_{\mathcal{F}})_{(i)}$ along each mode, i.e.,
$[\sim,\mathbf{\Sigma}_{\mathcal{F}_{i}},\sim]=\texttt{SVD}((\mathcal{S}_{\mathcal{F}})_{(i)},r^{\prime}_{\mathcal{F}_{i}})$,
$i=1,2,3$. This quantity measures the relative contribution of the $r$th rank.
The rank ($r_{\mathcal{F}_{i}}$) is adjusted or remains unchanged to maintain
$\epsilon$ within a desired range of
$\epsilon_{l}\leq\epsilon_{i}\leq\epsilon_{u}$, where $\epsilon_{u}$ and
$\epsilon_{l}$ are user-specified upper and lower thresholds. Rank increase or
decrease only requires truncating the SVD of matrix $\mathbf{C}_{i}$ (Eq. 11)
at $r_{\mathcal{F}_{i}}+1$ or $r_{\mathcal{F}_{i}}-1$, respectively. If
$r_{\mathcal{F}_{i}}$ is updated, then
$r^{\prime}_{\mathcal{F}_{i}}=r_{\mathcal{F}_{i}}+2$ must be updated
accordingly.

###### Remark 1

The size of the intersection tensor along each mode must be greater than or
equal to the corresponding Tucker rank, i.e.,
$r^{\prime}_{\mathcal{F}_{i}}\geq r_{\mathcal{F}_{i}}$. This constraint can be
explained by inspecting Eq. 12. If
$r^{\prime}_{\mathcal{F}_{i}}=r_{\mathcal{F}_{i}}$, it is easy to verify that
the above least squares problem becomes an interpolation problem, i.e., the
$\hat{\mathcal{F}}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})=\mathcal{F}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})$.
When $r^{\prime}_{\mathcal{F}_{i}}>r_{\mathcal{F}_{i}}$, Eq. 12 amounts to a
regression solution, i.e., an overdetermined system of equations. However, if
$r^{\prime}_{\mathcal{F}_{i}}<r_{\mathcal{F}_{i}}$, Eq. 12 becomes an
underdetermined system of equations which can result in poor or unstable
solutions.

1

Input:

$\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$: matrix of exact or approximate left
singular vectors of $\mathcal{F}_{(i)}$.

$r_{\mathcal{F}_{i}}$: target Tucker rank.

$\mathcal{F}$: function handle to compute fibers of the target $\mathcal{F}$

2

Output:
$\mathcal{S}_{\mathcal{F}},~{}\mathbf{U}^{(1)}_{\mathcal{F}},~{}\mathbf{U}^{(2)}_{\mathcal{F}},~{}\mathbf{U}^{(3)}_{\mathcal{F}}$

3

4$\mathbf{p}_{i}=$ DEIM($\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$), $\rhd$
Determine $\mathbf{p}_{1}\in\mathbb{I}^{r^{\prime}_{\mathcal{F}_{1}}}$,
$\mathbf{p}_{2}\in\mathbb{I}^{r^{\prime}_{\mathcal{F}_{2}}}$,
$\mathbf{p}_{3}\in\mathbb{I}^{r^{\prime}_{\mathcal{F}_{3}}}$ where
$r^{\prime}_{\mathcal{F}_{i}}={r_{\mathcal{F}}}_{i}+2$

5$\mathbf{C}_{1}=\Bigl{(}\mathcal{F}(:,\mathbf{p}_{2},\mathbf{p}_{3})\Bigr{)}_{(1)}$,
$\rhd$ Calculate $\mathbf{C}_{1}\in\mathbb{R}^{N_{1}\times
r^{\prime}_{\mathcal{F}_{2}}r^{\prime}_{\mathcal{F}_{3}}}$,
$\mathbf{C}_{2}\in\mathbb{R}^{N_{2}\times
r^{\prime}_{\mathcal{F}_{1}}r^{\prime}_{\mathcal{F}_{3}}}$,
$\mathbf{C}_{3}\in\mathbb{R}^{N_{3}\times
r^{\prime}_{\mathcal{F}_{1}}r^{\prime}_{\mathcal{F}_{2}}}$

$\mathbf{C}_{2}=\Bigl{(}\mathcal{F}(\mathbf{p}_{1},:,\mathbf{p}_{3})\Bigr{)}_{(2)}$,

$\mathbf{C}_{3}=\Bigl{(}\mathcal{F}(\mathbf{p}_{1},\mathbf{p}_{2},:)\Bigr{)}_{(3)}$

6$[\mathbf{U}^{(i)}_{\mathcal{F}},\sim,\sim]=\texttt{SVD}(\mathbf{C}_{i},r_{\mathcal{F}_{i}})$
$\rhd$ Calculate the left singular vectors of $\mathbf{C}_{i}$ and truncate at
rank $r_{\mathcal{F}_{i}}$

7$\mathcal{W}=\mathcal{F}(\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})$
$\rhd$ Form $\mathcal{W}\in\mathbb{R}^{r^{\prime}_{\mathcal{F}_{1}}\times
r^{\prime}_{\mathcal{F}_{2}}\times r^{\prime}_{\mathcal{F}_{3}}}$

8$\mathcal{S}_{\mathcal{F}}=\mathcal{W}\times_{1}~{}\mathbf{U}^{(1)}_{\mathcal{F}}(\mathbf{p}_{1},:)^{\dagger}\times_{2}~{}\mathbf{U}^{(2)}_{\mathcal{F}}(\mathbf{p}_{2},:)^{\dagger}\times_{3}~{}\mathbf{U}^{(3)}_{\mathcal{F}}(\mathbf{p}_{3},:)^{\dagger}$
$\rhd$ Calculate the core tensor ($\mathcal{S}_{\mathcal{F}}$)

Algorithm 1 DEIM-FS Tucker tensor low-rank approximation

###  2.6 Iterative DEIM-FS Tucker tensor approximation

Step 1 of Algorithm 1 requires access to the HOSVD factor matrices or some
close approximations to $\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$. As we will
discuss later, DLRA serves as one such example wherein the factor matrices
from the _previous_ time step are utilized in the DEIM algorithm to determine
which fibers should be sampled at the _current_ time step. However, for
problems in which these factor matrices are unknown, Algorithm 1 can be
applied iteratively with minor modifications, as explained below.

In the first iteration, Step 1 is skipped and instead, the fiber indices
($\mathbf{p}_{i}$) are chosen randomly. Then Step 2 is executed. In Step 3,
$\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$ are stored as the first
$r^{\prime}_{\mathcal{F}_{i}}$ left singular vectors of matrix
$\mathbf{C}_{i}$. Therefore, the first $r_{\mathcal{F}_{i}}$ columns of
$\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$ are identical to matrix
$\mathbf{U}_{\mathcal{F}}^{(i)}$. However, since
$r^{\prime}_{\mathcal{F}_{i}}\geq r_{\mathcal{F}_{i}}$ , more SVD columns are
stored in matrix $\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$. Then Steps 4 and 5
are executed. In the second iteration and the iterations after that Steps 1-5
of Algorithm 1 are executed with the only modification that
$\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$ are updated in Step 3. The iterations
can continue until the singular values of unfolded core tensor
(${\mathcal{S}_{\mathcal{F}}}_{(i)}$) converge up to a threshold value. In
particular, the convergence criterion is defined as:

$\dfrac{\Bigl{|}\|\boldsymbol{\Sigma}_{\mathcal{F}_{i}}^{k}\|_{F}-\|\boldsymbol{\Sigma}_{\mathcal{F}_{i}}^{k-1}\|_{F}\Bigr{|}}{\|\boldsymbol{\Sigma}_{\mathcal{F}_{i}}^{k}\|_{F}}<\epsilon,$
(15)

where $\boldsymbol{\Sigma}_{\mathcal{F}_{i}}^{K}$ is the matrix of singular
values of the $i^{\text{th}}$ unfolding of the core tensor in the
$k^{\text{th}}$ iteration and $\epsilon$ is the threshold value. We also note
that similar iterative approaches have been used in the past for tensor train
cross approximation [33].

###  2.7 Computational complexity

In this section, we present the computational complexity of DEIM-FS. In the
following analysis, we use the fact that $r_{\mathcal{F}_{i}}$ and
$r^{\prime}_{\mathcal{F}_{i}}$ are of $\mathcal{O}(r)$. We also assume
$N_{i}\sim\mathcal{O}(N)$.

In Step 1, the computational cost of finding the DEIM indices for a
$d$-dimensional tensor, scale with $\mathcal{O}(rdN)$. Computing and or
storing the fibers in Step 2 scale with $\mathcal{O}(dNr^{d-1})$.

The computational cost of performing SVD is
$\mathcal{O}(Nr^{2(d-1)})+\mathcal{O}(r^{3(d-1)})$ or
$\mathcal{O}(N^{2}r^{d-1})+\mathcal{O}(N^{3})$, whichever is smaller. One can
also use randomized SVD to obtain $\mathbf{U}^{(i)}_{\mathcal{F}}$ in order to
reduce the computational cost of performing SVD. Randomized SVD algorithms can
be particularly effective since the $\mathbf{C}_{i}$ matrices have a large
number of columns ($\mathcal{O}(r^{d-1})$) and only a small ($r$) left
singular vectors need to be computed accurately. We have not used randomized
SVD algorithms in any of the test cases in this paper.

Forming $\mathcal{W}$ does not require extra calculation as the values of the
intersection tensor are already calculated in Step 2 and can be extracted from
any of the $\mathbf{C}_{i}$ matrices. The computational cost of computing the
core tensor is $\mathcal{O}(dr^{d+1})$.

In summary, Steps 2 and 3 are the two most computationally expensive parts of
the DEIM-FS algorithm. As mentioned above, in Step 3, randomized SVD
algorithms can be utilized to further reduce the cost. On the other hand, Step
2 requires accessing elements of tensor $\mathcal{F}$. For problems where
accessing any element of tensor $\mathcal{F}$ requires additional computation,
the computational cost of computing the fibers can dominate the overall cost.
An example of this type of problem is DLRA, where computing any element of the
right-hand side tensor $\mathcal{F}=\mathscr{F}(\hat{\mathcal{V}})$ requires
applying a nonlinear discrete differential operator on tensor
$\hat{\mathcal{V}}$.

###  2.8 DEIM fiber sampling for dynamical low-rank approximation

The DEIM-FS algorithm can be employed to construct a low-rank Tucker tensor
approximation of $\mathscr{F}(\hat{\mathcal{V}})$ involved in the DLRA
evolution equations. As demonstrated, this approach leads to significant
computational savings when solving Eqs. (6a) and (6b). Step 1 of Algorithm 1
requires $\tilde{\mathbf{U}}^{(i)}_{\mathcal{F}}$, representing the exact or
approximate left singular vectors of matrices $\mathcal{F}_{(i)}$. To obtain
these vectors, we utilize $\mathbf{U}^{(i)}_{\mathcal{F}}$ from the previous
time step, which has already been calculated using the DEIM-FS algorithm. It
is important to note that the $\tilde{\mathbf{U}}^{(i)}_{\mathcal{F}}$
matrices are necessary solely for the DEIM algorithm to compute the fiber
indices, while the actual computation of the fibers is carried out for
$\mathscr{F}(\hat{\mathcal{V}})$ at the current time step.

In practice, utilizing $\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$ from the
previous time step yields excellent performance. The difference in accuracy
compared to cases where Algorithm 1 is used iteratively, as detailed in
Section 2.6, is negligible. The same approach was used in previous studies for
low-rank approximation of matrix differential equations [19, 20]. At $t=0$, if
no good approximation for $\tilde{\mathbf{U}}^{(i)}_{\mathcal{F}}$ exists,
DEIM-FS (iterative) may be used.

In the DLRA equations, the fibers of tensor $\mathcal{F}$ are calculated by
evaluating the function $\mathcal{F}=\mathscr{F}(\hat{\mathcal{V}})$. For the
TDEs obtained from discretizing a PDE, the function $\mathscr{F}(\sim)$
involves discrete differential operators. Since the derivative calculation
requires adjacent points of the DEIM-selected points, the adjacent points must
be determined in Step 2 of Algorithm 1. This step depends on the numerical
scheme used for the spatial discretization. For instance, when employing the
spectral element method, determining the derivative at a specific spatial
point requires access to the values of other points within the same element
[19]. As a result, calculating any fiber of $\mathcal{F}$ requires access to
the values to additional (adjacent) fibers of $\hat{\mathcal{V}}$. We denote
the additional fiber indices with $\mathbf{p}_{a_{i}}$. For example, for a
three-dimensional tensor, $\mathbf{C}_{1}$ is calculated as:

$\mathcal{F}(:,\mathbf{p}_{2},\mathbf{p}_{3})=\mathscr{F}\Bigl{(}\mathcal{S}\times_{1}~{}\mathbf{U}^{(1)}\times_{2}~{}\mathbf{U}^{(2)}([\mathbf{p}_{2},\mathbf{p}_{a_{2}}],:)\times_{3}~{}\mathbf{U}^{(3)}([\mathbf{p}_{3},\mathbf{p}_{a_{3}}],:)\Bigr{)}.$
(16)

Applying this methodology to DLRA enables a significant reduction in the
computational cost and memory. In these equations, $\mathcal{F}$ is $N^{d}$
tensor, while the memory advantage of DEIM-FS anbles storing $\mathcal{F}$ in
the Tucker compressed form. Accordingly, Eqs. (6a) and (6b) can be rewritten
as:

$\displaystyle\dot{\mathcal{S}}=\mathcal{S}_{\mathcal{F}}~{}\bigtimes_{i=1}^{d}~{}\bigl{(}{\mathbf{U}^{(i)}}^{T}\mathbf{U}^{(i)}_{\mathcal{F}}\bigr{)},$
(17a)
$\displaystyle\dot{\mathbf{U}}^{(i)}=\Bigl{(}\mathbf{I}-\mathbf{U}^{(i)}\mathbf{U}^{(i)^{T}}\Bigr{)}\mathbf{U}^{(i)}_{\mathcal{F}}~{}\Bigl{[}\mathcal{S}_{\mathcal{F}}\bigtimes_{k\neq
i}\bigl{(}{\mathbf{U}^{(k)}}^{T}\mathbf{U}^{(k)}_{\mathcal{F}}\bigr{)}\Bigr{]}_{(i)}~{}\mathcal{S}_{(i)}^{\dagger},$
(17b)

where $\mathcal{S}_{\mathcal{F}}$ and $\mathbf{U}^{(i)}_{\mathcal{F}}$ are the
outputs of the Algorithm 1. Another advantage of employing DEIM-FS for DLRA is
its simplification of the implementation of DLRA equations. This is due to the
tensor $\mathcal{F}$ being approximated in a black-box fashion and Eqs. (17a)
and (17b) are agnostic to the type of TDE being solved by DLRA. Hereinafter,
Eqs. (17a) and (17b) are referred to as DLRA-DEIM-FS.

###  2.9 Comparison to an existing fiber sampling algorithm

We compare our proposed algorithm with the FSTD algorithm [32]. While FSTD has
some similarities to DEIM-FS it also has some key differences with the
presented algorithm. We briefly review the FSTD algorithm here. In FSTD, the
cross Tucker model of $\mathcal{F}$ is given by:

$\mathcal{F}\approx\mathcal{W}\times_{1}\mathbf{C}_{1}\mathcal{W}_{(1)}^{\dagger}\times_{2}\mathbf{C}_{2}\mathcal{W}^{\dagger}_{(2)}\times_{3}\mathbf{C}_{3}\mathcal{W}^{\dagger}_{(3)},$
(18)

where $\mathbf{C}_{1}\in\mathbb{R}^{N_{1}\times
r^{\prime}_{\mathcal{F}_{2}}r^{\prime}_{\mathcal{F}_{3}}}$,
$\mathbf{C}_{2}\in\mathbb{R}^{N_{2}\times
r^{\prime}_{\mathcal{F}_{1}}r^{\prime}_{\mathcal{F}_{3}}}$,
$\mathbf{C}_{3}\in\mathbb{R}^{N_{3}\times
r^{\prime}_{\mathcal{F}_{1}}r^{\prime}_{\mathcal{F}_{2}}}$ are the unfolded
$\mathcal{C}_{1},\mathcal{C}_{2},\mathcal{C}_{3}$ fiber sub-tensors shown in
Figure 1, $\mathcal{W}\in\mathbb{R}^{r^{\prime}_{\mathcal{F}_{1}}\times
r^{\prime}_{\mathcal{F}_{2}}\times r^{\prime}_{\mathcal{F}_{3}}}$ is the
intersection tensor. According to the FSTD algorithm presented in [32], at the
first step, the index of the first fiber is initialized. Then, the indices of
the other fibers are determined based on a deterministic greedy algorithm.

We compare FSTD to DEIM-FS in terms of the information these algorithms
require about the tensors, their accuracy, and the computational cost. The
FSTD algorithm does not require any information about $\mathcal{F}$, whereas
DEIM-FS requires the exact or approximate left singular vectors of all
unfoldings of $\mathcal{F}$. The DEIM-FS (iterative) algorithm, on the other
hand, does not require left singular vectors of $\mathcal{F}$ unfolding.
However, computing DEIM-FS (iterative) requires iterations and it is more
expensive to compute in comparison to FSTD.

The DEIM-FS Tucker tensor models appear to be more accurate than FSTD Tucker
tensor models of the same rank. This is due to several factors: (i) The FSTD
becomes ill-conditioned as the rank increases. The same issue also exists for
matrix CUR decompositions, where the intersection matrix becomes singular as
rank increases [33]. This issue is mitigated in matrix CUR by performing QR
decomposition of selected columns or rows. The issue of ill-conditioning is
also mitigated in DEIM-FS, since the factor matrices are obtained by
performing SVD of the selected fibers. The SVD, similar to QR, results in a
set of orthonormal modes, and the matrices
$\mathbf{U}^{(i)}_{\mathcal{F}}(\mathbf{p}_{i},:)$ are well-conditioned
matrices. In fact, the DEIM algorithm is designed to maintain
$\|\mathbf{U}^{(i)}_{\mathcal{F}}(\mathbf{p}_{i},:)^{-1}\|_{2}$ as small as
possible. (ii) The choice of initial fibers in FSTD can have a significant
impact on the accuracy of the resulting Tucker model. We show this effect in
our numerical examples. We also show that DEIM-FS (iterative) is much less
sensitive to the initial choice of fibers.

From the computational cost point of view, both FSTD and DEIM-FS require
access to $\mathcal{O}(dr^{d-1}N)$ number of elements of $\mathcal{F}$. They
also have the same memory requirements. However, DEIM-FS requires more flops
due to the computation of SVD of $\mathbf{C}_{i}$ matrices. The FSTD algorithm
offers the advantage of storing the actual fibers of the target tensor,
thereby inheriting the structure of the target tensor. For instance, in the
case of sparse tensors, storing FSTD factor matrices in a sparse form can
potentially achieve a very high compression ratio.

##  3 Demonstration

###  3.1 Toy Examples

As our first example, we consider two three-dimensional functions as shown
below:

$\displaystyle\mathscr{F}_{1}(x_{1},x_{2},x_{3})=e^{-(x_{1}~{}x_{2}~{}x_{3})^{2}}\hskip
96.73918ptx_{1},x_{2},x_{3}\in[-1,1],$ (19a)
$\displaystyle\mathscr{F}_{2}(x_{1},x_{2},x_{3})=\frac{1}{\bigl{(}x_{1}^{b}+x_{2}^{b}+x_{3}^{b}\bigr{)}^{1/b}}\hskip
28.45274ptx_{1},x_{3}\in[1,300]~{}~{}~{}~{}x_{2}\in[1,400].$ (19b)

The tensor $\mathcal{F}_{1}$ is obtained by evaluation
$\mathscr{F}_{1}(x_{1},x_{2},x_{3})$ at 100 equally spaced elements of
$x_{1}$, $x_{2}$, and $x_{3}$ in their respective domains. Therefore,
$\mathcal{F}_{1}\in\mathbb{R}^{100\times 100\times 100}$. Similarly,
$\mathcal{F}_{2}\in\mathbb{R}^{300\times 400\times 300}$ is obtained by
evaluating $\mathscr{F}_{2}(x_{1},x_{2},x_{3})$ on a uniform grid in each
direction. Two choices of $b=3$ and $b=5$ are considered for
$\mathcal{F}_{2}$. In all demonstrations of DEIM-FS, the exact left singular
vectors of the unfolding of the target tensors are used in the DEIM algorithm.
To compare the accuracy and efficiency of DEIM-FS against HOSVD and FSTD,
[32], $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$ are approximated using the three
mentioned algorithms. Then the error between the approximated tensors and the
actual tensors is calculated. Denoting $\hat{\mathcal{F}}_{1}$ and
$\hat{\mathcal{F}}_{2}$ as the low-rank Tucker approximation of
$\mathcal{F}_{1}$ and $\mathcal{F}_{2}$, respectively, the error is defined
as:

$\mathcal{E}=\|\hat{\mathcal{F}}-\mathcal{F}\|_{F}.$ (20)

For DEIM-FS, we consider
${r_{\mathcal{F}}}_{1}={r_{\mathcal{F}}}_{2}={r_{\mathcal{F}}}_{3}={r_{\mathcal{F}}}$.
Hence,
$r^{\prime}_{\mathcal{F}_{1}}=r^{\prime}_{\mathcal{F}_{2}}=r^{\prime}_{\mathcal{F}_{3}}=r^{\prime}_{\mathcal{F}}$.
As mentioned in Remark 1, $r^{\prime}_{\mathcal{F}}\geq r_{\mathcal{F}}$. In
Figure 2a, we increase $r^{\prime}_{\mathcal{F}}$ for a fixed
${r_{\mathcal{F}}}$ to study the effect of increasing
$r^{\prime}_{\mathcal{F}}$. The results of Figure 2a indicate that increasing
$r^{\prime}_{\mathcal{F}}$ beyond $r_{\mathcal{F}}+2$ results in a negligible
reduction in error. We use $r^{\prime}_{\mathcal{F}}=r_{\mathcal{F}}+2$ for
rest of examples in this paper.

Figures 2b and 2c show the low-rank error versus rank for $\mathcal{F}_{1}$
and $\mathcal{F}_{2}$, respectively. As it can be seen, for both tensors, the
error of DEIM-FS method closely follows the HOSVD error. However, the FSTD
error is always greater than the DEIM-FS error. Moreover, the FSTD error
either does not decrease (see Figure 2c) or even increases (see Figure 2b) as
rank increases. This is due to the issue of ill-conditioning as discussed in
Section 2.9.

Figure 2: Toy Example: Comparison of the FSTD, DEIM-FS, DEIM-FS (iterative),
and HOSVD: (a) approximation errors of $\mathcal{F}_{1}$ and $\mathcal{F}_{2}$
versus $r^{\prime}_{\mathcal{F}}$ for a fixed rank of $r_{\mathcal{F}}$; (b)
approximation error of $\mathcal{F}_{1}$ versus rank; (c) approximation error
of $\mathcal{F}_{2},(b=3,5)$ versus rank; (d) approximation error of
$\mathcal{F}_{2},(b=5)$ versus rank; (e) effect of fiber initialization on the
approximation error of $\mathcal{F}_{2},(b=5)$ for a fixed rank of
$r_{\mathcal{F}}=20$ for 100 different initializations; (f) convergence of the
singular values of the core tensor ($\mathcal{S}_{\mathcal{F}}$) versus
iterations for DEIM-FS (iterative) algorithm.

As discussed in Section 2.6, in cases where
$\tilde{\mathbf{U}}_{\mathcal{F}}^{(i)}$ is not available, the DEIM-FS
(iterative) approach can be employed. Figure 2d illustrates that the error
obtained by the DEIM-FS (iterative) algorithm closely follows the errors of
the HOSVD and DEIM-FS methods. In the DEIM-FS (iterative) algorithm, the
initial fibers are randomly selected, while in DEIM-FS, the exact left
singular vectors are utilized.

In both the FSTD and DEIM-FS (iterative) algorithms, the initial fibers are
selected randomly. Figure 2e examines the effect of different initializations
on accuracy. The errors of Tucker models, for $\mathcal{F}_{2}$ resulting from
100 random initializations of both FSTD and DEIM-FS (iterative), are depicted
in Figure 2e. It is noticeable that the variance of errors obtained by FSTD is
significantly larger than those obtained by DEIM-FS (iterative).
Interestingly, for some fiber initializations in FSTD, the fiber indices
remain fixed on the initial fibers and do not update, leading to very large
errors ($\mathcal{E}\approx 30$). These cases are not displayed in Figure 2e
because we believe a minor fix could resolve this issue. In summary, our
observation reveals that the error in FSTD is highly sensitive to fiber
initializations, while DEIM-FS (iterative) displays a significantly smaller
variance in error.

For the 100 initializations presented in Figure 2e, the average number of
iterations required for convergence by DEIM-FS (iterative) is 5.48, with a
standard deviation of 1.24. Figure 2f illustrates the convergence of the
singular values of the core tensor ($\mathcal{S}_{\mathcal{F}}$) determined by
the DEIM-FS (iterative) algorithm.

###  3.2 DLRA for four-dimensional Fokker Planck equation

For the second test case, a four-dimensional Fokker Planck (FP) equation is
considered. The physics under consideration is governed by the Langevin form
of stochastic differential equations (SDE).

$d\mathbf{x}=\mathbf{f}(\mathbf{x},t)dt+\mathbf{S}(\mathbf{x},t)d\mathbf{w},$
(21)

where $\mathbf{x}=[x_{1},x_{2},x_{3},x_{4}]$ is the stochastic transport
variable, $\mathbf{f}(\mathbf{x},t)$ is the drift function,
$\mathbf{S}(\mathbf{x},t)$ is the diffusion matrix and $\mathbf{w}$ denotes
the Wiener-Levy process [40]. In this equation, $\mathbf{x}$ and
$\mathbf{f}(\mathbf{x},t)$ are $d$-dimensional vectors,
$\mathbf{S}(\mathbf{x},t)$ is a $d\times m$ matrix, and $\mathbf{w}$ is an
$m$-dimensional vector. The transitional probability density function (PDF) of
the stochastic variable is governed by its FP equation:

$\frac{\partial p(\mathbf{x},t)}{\partial
t}=-{\sum}_{i=1}^{d}\frac{\partial}{\partial
x_{i}}\bigl{[}f_{i}(\mathbf{x},t)p(\mathbf{x},t)\bigr{]}+\frac{1}{2}{\sum}_{i=1}^{d}{\sum}_{j=1}^{d}\frac{\partial^{2}}{\partial
x_{i}\partial x_{j}}\bigl{[}D_{ij}(x,t)p(x,t)\bigr{]},$ (22)

where $\mathbf{D}=[D_{ij}]$ is a $d\times d$ matrix and
$\mathbf{D}=\mathbf{S}\mathbf{S}^{T}$. We consider a 4-dimensional super-
symmetric PDF and we choose $r=r_{1}=r_{2}=r_{3}=r_{4}$ and
${r_{\mathcal{F}}}_{1}={r_{\mathcal{F}}}_{2}={r_{\mathcal{F}}}_{3}={r_{\mathcal{F}}}_{4}={r_{\mathcal{F}}}$.
We consider homogeneous Dirichlet boundary condition and a normal distribution
as the initial condition as:
$p_{0}(\mathbf{x})=\mathcal{N}\bigl{(}\boldsymbol{\mu}_{0},~{}\mathbf{C}_{0}\bigr{)}$,
where $\mathcal{N}$ is the normal distribution and $\boldsymbol{\mu}_{0}$ is
the mean set to
$\boldsymbol{\mu}_{0}=\begin{bmatrix}1.5&0.6&-0.3&-1.2\end{bmatrix}^{T}$, and
$\mathbf{C}_{0}$ is the covariance matrix, where $\mathbf{C}_{0_{ii}}=1$,
$i=1,\dots,4$ and $\mathbf{C}_{0_{ij}}=0.5$, $i,j=1,\dots,4$ and $i\neq j$. We
also consider $f_{i}=-\alpha x_{i}~{},~{}i=1,2,3,4$ and $\alpha=0.75$. The
solution $p(x,t)$ remains Gaussian when $D_{ij}$ is a constant. The analytical
expressions are obtained for the moments over time as:

$\boldsymbol{\mu}(t)=\boldsymbol{\mu}_{0}\exp(-\alpha t),$ (23)
$\mathbf{C}(t)=\frac{\mathbf{D}}{2\alpha}+\bigl{(}\mathbf{C}_{0}-\frac{\mathbf{D}}{2\alpha}\bigr{)}\exp(-2\alpha
t).$ (24)

For spatial discretization, the third-order spectral method is used. Each
domain is discretized via $N=N_{1}=N_{2}=N_{3}=N_{4}=61$ points, which
includes 20 elements with a second-order order polynomial approximation within
each of the elements. For time advancement, the fourth-order Runge-Kutta (RK4)
method is used. The domain under consideration is $x\in[-6,6]^{4}$ within the
time interval $t\in[0,8]$. The time advancement of $\Delta t=2\times 10^{-3}$
is chosen. The resulting TDE is solved using DLRA-DEIM-FS (Eqs. 17a and 17b)
with $r=5$. Although the algorithm can be adaptive, we solve this example with
fixed ${r_{\mathcal{F}}}=5$ and $r^{\prime}_{\mathcal{F}}=7$. Figure 3a shows
the temporal evolution of the singular values obtained from DLRA-DEIM-FS and
analytical solution. For the rest of the paper, the error is defined as the
relative error as:

$\mathcal{E}(t)=\dfrac{\|\hat{\mathcal{V}}(t)-\mathcal{V}(t)\|_{F}}{\|\mathcal{V}(t)\|_{F}}.$
(25)

Figure 3: Fokker Planck Equation: (a) first 5 singular values of analytical
solution and DLRA-DEIM-FS; (b) relative error evolution for different $r$ and
$r_{\mathcal{F}}$.

We conduct a convergence study by solving the DLRA equations with various $r$
and $r_{\mathcal{F}}$ values. The temporal evolution of the error is depicted
in Figure 3b. It is observed that the model with $r=5$ and $r_{\mathcal{F}}=5$
yields the best error. The statistical moments obtained by solving the FP
equation with DLRA-DEIM-FS can be compared to those obtained analytically
(refer to Eqs. (23) and (24)). Table 1 presents the comparison of the mean of
the variables and the covariance matrices at the final time step $(t=8)$.

Table 1: Comparison of the moments of the PDF at $t=8$ with $r=5$ and $r_{\mathcal{F}}=5$ | Mean | Covariance
---|---|---
Analytical | $\begin{bmatrix}0.0000&0.0000&0.0000&0.0000\end{bmatrix}$ | $\begin{bmatrix}0.6675&0.0272&0.0272&0.0272\\\ 0.0272&0.6675&0.0272&0.0272\\\ 0.0272&0.0272&0.6675&0.0272\\\ 0.0272&0.0272&0.0272&0.6675\\\ \end{bmatrix}$
DLRA-DEIM-FS | $\begin{bmatrix}0.0017&0.0041&-0.0007&-0.0028\end{bmatrix}$ | $\begin{bmatrix}0.6674&0.0272&0.0270&0.0271\\\ 0.0272&0.6673&0.0271&0.0272\\\ 0.0270&0.0271&0.6675&0.0270\\\ 0.0271&0.0272&0.0270&0.6674\\\ \end{bmatrix}$

Although this example is not a nonlinear TDE, there are still advantages to
using DLRA-DEIM-FS to solve this TDE rather than employing standard DLRA
techniques. To illustrate this, consider the right-hand side of the FP
equation, which involves the summation of 20 terms. Implementing standard DLRA
requires a highly intrusive and meticulous treatment of these terms to avoid
storing the full-dimensional right-hand side tensor. However, the DLRA-DEIM-FS
algorithm remains agnostic to the nature of the TDE, involving the evaluation
of sparse fibers of the right-hand side in a black-box fashion. Therefore,
DLRA-DEIM-FS is significantly easier to implement compared to the standard
practice for solving DLRA equations, even for linear TDEs.

Solving the FOM for this example requires storing $N^{4}=13,845,841$ floating-
point numbers for the right-hand side of the TDE. Using DLRA-DEIM-FS with
$r_{\mathcal{F}}=5$ requires storing $d\times
N\times{r^{\prime}_{\mathcal{F}}}^{d-1}+{r^{\prime}_{\mathcal{F}}}^{4}=86,093$
floating-point numbers, resulting in a memory compression ratio of
$13,845,841/86,093=160.8$.

###  3.3 DLRA for four-dimensional nonlinear advection equation

For the final demonstration, we consider a four-dimensional nonlinear
advection given below:

$\frac{\partial v(\mathbf{x},t)}{\partial t}=-\mathbf{b}\cdot\nabla
v(\mathbf{x},t)+s(v(\mathbf{x},t)),\quad\mathbf{x}\in[-5,5]^{4},$ (26)

where $\mathbf{b}=[-\sin(t),~{}\cos(t),~{}-\sin(\pi+t),~{}\cos(\pi+t)]$ is the
advection velocity and $s(v)=-0.1v/(1+v^{2})$ is the nonlinear source term.
The boundary condition is periodic and the initial condition is considered as
$v_{0}(\mathbf{x})=f(\mathbf{x})+g(\mathbf{x})$, where

$f(\mathbf{x})=e^{-(x_{1}-\frac{1}{2})^{2}}e^{-(x_{1}+\frac{x_{2}}{2}-\frac{1}{2})^{2}}e^{-(x_{3}-\frac{1}{2})^{2}}e^{-(x_{3}+\frac{x_{4}}{2}-\frac{1}{2})^{2}},$
$g(\mathbf{x})=e^{-(x_{1}+\frac{1}{2})^{2}}e^{-(x_{2}+\frac{1}{2})^{2}}e^{-(x_{3}+\frac{1}{2})^{2}}e^{-(x_{4}+\frac{1}{2})^{2}}.$

The second-order spectral method is used for spatial discretization
($\mathbf{x}$), and the fourth-order Runge-Kutta (RK4) method is employed for
temporal integration within the time interval $t\in[0,4]$. The time step is
set to $\Delta t=2\times 10^{-3}$. Each domain is discretized with
$N=N_{1}=N_{2}=N_{3}=N_{4}=85$ points, comprising 42 elements with a second-
order polynomial approximation within each element. The rank-adaptive DLRA-
DEIM-FS is utilized for approximating the right-hand side of the resulting
TDE. It is important to note that the rank ($r$) of the solution
$\hat{\mathcal{V}}(t)$ is not adaptive. The relative error with respect to the
norm (Eq. (25)) is used for error calculation. As the analytical solution is
unavailable in this case, the FOM is solved using $N=85$ and $\Delta t=2\times
10^{-3}$. The solution obtained from FOM is considered the ground truth.

Three distinct errors contribute to the overall error while using DLRA-DEIM-
FS: the DEIM-FS error for low-rank approximation of
$\mathscr{F}(\hat{\mathcal{V}})$, controlled by varying threshold values; the
DLRA error for approximating $\hat{\mathcal{V}}(t)$, controlled by $r$; and
the temporal integration error, influenced by adjusting $\Delta t$.

First, the TDE is solved with a fixed $r=6$ and various rank-adaptivity
thresholds: $\epsilon_{u}={10^{-2},10^{-4},10^{-6}}$ and
$\epsilon_{l}={10^{-3},10^{-5},10^{-7}}$. Figure 4a indicates that for
$(\epsilon_{l}=10^{-6},\epsilon_{u}=10^{-7})$, the DLRA error does not
decrease. This suggests that the error in approximating the right-hand side
tensor is not dominant, and the error reaches a saturation point due to the
DLRA low-rank approximation error in $\hat{\mathcal{V}}(t)$. This also
demonstrates that in DLRA it is not necessary to compute
$\mathscr{F}(\hat{\mathcal{V}})$ exactly.

Figure 4b illustrates the evolution of $r_{\mathcal{F}}$ values associated
with different $\epsilon_{l}$ and $\epsilon_{u}$ threshold values. As
expected, lower thresholds correspond to higher ranks. In Figure 4c, we
examine the effect of the rank $\hat{\mathcal{V}}(t)$ on the total error. In
this case, the threshold values for approximating
$\mathscr{F}(\hat{\mathcal{V}})$ are fixed at
$(\epsilon_{l}=10^{-5},\epsilon_{u}=10^{-4})$. The results demonstrate that
increasing $r$ enhances accuracy. Another study is also conducted to consider
the effect of the time step $(\Delta t)$ on the error of the DLRA-DEIM-FS
method. To this end, the DLRA rank and the thresholds are fixed: $r=11$ and
$(\epsilon_{l}=10^{-5},\epsilon_{u}=10^{-4})$. The TDE is then solved for
three different values of $\Delta t=\\{2\times 10^{-3},4\times 10^{-3},6\times
10^{-3}\\}$. Note that the $\Delta t$ of the FOM (the ground truth) was not
changed and it was fixed to be $2\times 10^{-3}$. As shown in Figure 4d the
effect of different time steps is considerable until $t\approx 2$. Afterward,
the low-rank approximation errors of $\hat{\mathcal{V}}(t)$ and
$\mathscr{F}(\hat{\mathcal{V}})$ dominates, and decreasing $\Delta t$ does not
change the error considerably.

Figure 5a presents a comparison between the first 11 singular values obtained
using DLRA-DEIM-FS and those of FOM, while Figure 5b shows the corresponding
evolution of $r_{\mathcal{F}}$ for this problem.

Figure 4: Four-dimensional nonlinear advection equation (Example 3): (a)
Relative error evolution with fixed $r$ for different $\epsilon_{l}$ and
$\epsilon_{u}$; (b) Evolution of $r_{\mathcal{F}}$ associated with the plots
of Figure 4a; (c) Relative error evolution for different $r$ with fixed
$\epsilon_{l}=10^{-5}$ and $\epsilon_{u}=10^{-4}$; (d) Relative error
evolution for different $\Delta t$ with fixed $r$, $\epsilon_{l}$, and
$\epsilon_{u}$. Figure 5: Four-dimensional nonlinear advection equation
(Example 3): (a) First 11 singular values of FOM and TDB-DEIM-FS; (b)
Evolution of $r_{\mathcal{F}}$ associated with the solved system in Figure 5a.

We also compute the following marginal function:

$\overline{v}(x_{3},x_{4},t)=\int_{-5}^{5}\int_{-5}^{5}{v(\mathbf{x},t)}^{2}~{}dx_{1}dx_{2},$
(27)

using both DLRA-DEIM-FS and FOM.

Figure 6 shows contour plots of $\overline{v}(x_{3},x_{4})$ obtained from FOM
(top row) and DLRA-DEIM-FS (bottom row) at two different time instances. The
first two columns of the $x_{3}$ and $x_{4}$ factor matrices are shown:
$\\{\mathbf{u}_{1}^{(3)},\mathbf{u}_{2}^{(3)}\\}$ are shown on the top side of
each panel and $\\{\mathbf{u}_{1}^{(4)},\mathbf{u}_{2}^{(4)}\\}$ are shown on
the right side of each panel. The first basis is displayed using solid blue
color, while the second basis is shown with red dashed lines. The factor
matrices of the FOM are obtained by performing HOSVD on the FOM solution
tensor. The selected fibers along the third and fourth modes of the tensor are
also shown with dashed lines. A good agreement between the DLRA-DEIM-FS and
FOM contour plots is observed. Also, it can be seen that the first two
$\mathbf{x}_{3}$ and $\mathbf{x}_{4}$ bases are localized and advect with the
solution. Moreover, the selected fibers are concentrated around the localized
solution.

Figure 6: Four-dimensional advection-reaction equation (Example 3): Top row:
FOM; Bottom row: DLRA-DEIM-FS. The DEIM sampled fibers along each mode are
shown with white dashed lines. On the top of each panel the first two
$x_{3}$-bases (the first two columns of $\mathbf{U}^{(3)}$ ) are shown. The
first basis is shown in solid blue line and the second basis is shown in
dashed red line. On the right side of each panel.

The PDE defined by Eq. 26 contains non-polynomial nonlinearity. Consequently,
in standard DLRA implementation, the explicit formation of the right-hand-side
tensor is inevitable because $\mathcal{F}=\mathscr{F}(\hat{\mathcal{V}})$ is
full rank despite $\hat{\mathcal{V}}$ having a low rank. The explicit
formation of $\mathcal{F}$ incurs substantial memory and computational costs,
scaling at least as $\mathcal{O}(N^{4})$, both in terms of memory allocation
and floating-point operations. The DLRA-DEIM-FS algorithm avoids this cost
because the full-rank tensor is never formed. Instead, a cross Tucker tensor
approximation of $\mathcal{F}$ is constructed on the fly. The computational
cost of computing $\mathcal{F}$ using DLRA-DEIM-FS scales linearly with $N$.
In Figure 7a the computational cost of solving DLRA evolution equations versus
$N$ using standard DLRA and DLRA-DEIM-FS are shown.

Figure 7b displays the required memory relative to FOM for storing
$\mathcal{F}=\mathscr{F}(\hat{\mathcal{V}})$ in Examples 2 and 3. The results
demonstrate a significant reduction in memory for both of these examples.

Figure 7: (a) Computational time vs $N$ (number of discretization points) for
4D nonlinear advection equation Four-dimensional advection-reaction equation
(Example 3); (b) Comparison of the relative memory for storing the RHS for 4D
Fokker Planck and 4D nonlinear advection equations.

##  4 Conclusions

Tucker tensor low-rank approximation is a building block of many important
tensor low-rank approximation algorithms. Determining the optimal (lowest
error) Tucker tensor model lacks a known closed solution. Instead, HOSVD is
commonly used, which yields near-optimal Tucker tensor models. However, HOSVD
requires access to all elements of the tensor, and its computational cost
scales at least with the size of the tensor, i.e., $\mathcal{O}(N^{d})$. In
this work, we present DEIM-FS — a cross algorithm that builds a rank-$r$
Tucker tensor model by sampling $\mathcal{O}(dr^{d-1})$ strategically selected
fibers. The fibers are selected using the DEIM algorithm. The computational
cost of DEIM-FS scales linearly with $N$. Our numerical results demonstrate
that DEIM-FS remains well-conditioned as $r$ increases and the error of the
algorithm closely follows the same-rank HOSVD errors.

We augment the DEIM-FS algorithm with two capabilities: (i) the introduction
of the DEIM-FS (iterative) algorithm, which, unlike DEIM-FS, does not require
access to the left singular vectors of the tensor unfolding. Consequently,
DEIM-FS (iterative) can be viewed as a black-box Tucker tensor model
algorithm. (ii) Additionally, we introduce the rank-adaptive DEIM-FS, where
the Tucker rank is adjusted to meet a specified accuracy threshold.

Finally, we introduce DLRA-DEIM-FS to tackle computational cost challenges
encountered when solving DLRA equations in nonlinear TDEs. The computational
savings of DLRA diminish in cases where the exact rank of the right-hand side
tensor ($\mathcal{F}(\hat{\mathcal{V}})$) is exceptionally high, even when
$\hat{\mathcal{V}}$ is low-rank, such as in nonlinear TDEs. We demonstrate
that DLRA-DEIM-FS mitigates these cost issues by constructing a low-rank
Tucker tensor approximation of the right-hand side tensor in the TDE.
Utilizing DEIM-FS, we substantially reduce the required memory and
computational cost involved in solving DLRA evolution equations.

## Acknowledgments

We thank Professor Peyman Givi for many stimulating discussions that led to a
number of improvements in this paper. This work is sponsored by the National
Science Foundation (NSF), USA under Grant CBET-2152803 and by the Air Force
Office of Scientific Research award no. FA9550-21-1-0247.

## Appendix 1

The DEIM pseudocode is presented via Algorithm 2. This algorithm is adopted
from [39].

Input:
$\mathbf{U}_{p}=\left[\begin{array}[]{llll}\mathbf{u}_{1}&\mathbf{u}_{2}&\cdots&\mathbf{u}_{p}\end{array}\right]$

Output: $\mathbf{I}_{p}$

1 $\left[\rho,\mathbf{I}_{1}\right]=\max\left|\mathbf{u}_{1}\right|$ $\rhd$
choose the first index;

2 $\mathbf{P}_{1}=\left[\mathbf{e}_{\mathbf{I}_{1}}\right]$ $\rhd$ construct
first measurement matrix;

3 for  _$i=2$ to $p$_ do

4
$\mathbf{P}_{i}^{T}\mathbf{U}_{i}\mathbf{c}_{i}=\mathbf{P}_{i}^{T}\mathbf{u}_{i+1}$
$\rhd$ calculate $c_{i}$;

5 $\mathbf{R}_{i+1}=\mathbf{u}_{i+1}-\mathbf{U}_{i}\mathbf{c}_{i}$ $\rhd$
compute residual;

6 $\left[\rho,\mathbf{I}_{i}\right]=\max\left|\mathbf{R}_{i+1}\right|$ $\rhd$
find index of maximum residual;

7
$\mathbf{P}_{i+1}=\left[\begin{array}[]{ll}\mathbf{P}_{i}&\mathbf{e}_{\mathbf{I}_{i}}\end{array}\right]$
$\rhd$ add new column to measurement matrix;

8 end for

Algorithm 2 DEIM Algorithm [39]

## References

  * [1] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.
  * [2] L. Grasedyck, D. Kressner, and C. Tobler, “A literature survey of low-rank tensor approximation techniques,” GAMM-Mitteilungen, vol. 36, pp. 53–78, 2020/07/06 2013.
  * [3] L. R. Tucker, “Some mathematical notes on three-mode factor analysis,” vol. 31, no. 3, pp. 279–311, 1966.
  * [4] R. A. Harshman, “Foundations of the PARAFAC procedure: Models and conditions for an ”explanatory” multi-modal factor analysis,” UCLA Working Papers in Phonetics, vol. 16, pp. 1–84, 1970.
  * [5] L. Grasedyck, “Hierarchical singular value decomposition of tensors,” SIAM Journal on Matrix Analysis and Applications, vol. 31, pp. 2029–2054, 2023/12/11 2010.
  * [6] I. V. Oseledets, “Tensor-train decomposition,” SIAM Journal on Scientific Computing, vol. 33, pp. 2295–2317, 2020/03/23 2011.
  * [7] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” SIAM Review, vol. 51, no. 3, pp. 455–500, 2009.
  * [8] O. Koch and C. Lubich, “Dynamical tensor approximation,” SIAM Journal on Matrix Analysis and Applications, vol. 31, no. 5, pp. 2360–2375, 2010.
  * [9] M. H. Beck, A. Jäckle, G. A. Worth, and H. D. Meyer, “The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets,” Physics Reports, vol. 324, pp. 1–105, 1 2000.
  * [10] H. Risken, The Fokker-Planck-Equation. Methods of Solution and Applications.  Springer Berlin, Heidelberg, 09 1996.
  * [11] A. M. Boelens, D. Venturi, and D. M. Tartakovsky, “Tensor methods for the Boltzmann-BGK equation,” Journal of Computational Physics, vol. 421, p. 109744, 2020.
  * [12] S. Dolgov, D. Kalise, and K. K. Kunisch, “Tensor decomposition methods for high-dimensional Hamilton-Jacobi-Bellman equations,” SIAM Journal on Scientific Computing, vol. 43, no. 3, pp. A1625–A1650, 2021.
  * [13] I. Gavrilyuk and B. N. Khoromskij, “Tensor numerical methods: Actual theory and recent applications,” Computational Methods in Applied Mathematics, vol. 19, no. 1, pp. 1–4, 2019.
  * [14] O. Koch and C. Lubich, “Dynamical low‐rank approximation,” SIAM Journal on Matrix Analysis and Applications, vol. 29, pp. 434–454, 2017/04/02 2007.
  * [15] L. Einkemmer and C. Lubich, “A low-rank projector-splitting integrator for the vlasov–poisson equation,” SIAM Journal on Scientific Computing, vol. 40, pp. B1330–B1360, 2023/08/15 2018.
  * [16] J. Hu and Y. Wang, “An adaptive dynamical low rank method for the nonlinear boltzmann equation,” Journal of Scientific Computing, vol. 92, no. 2, p. 75, 2022.
  * [17] M. Donello, M. H. Carpenter, and H. Babaee, “Computing sensitivities in evolutionary systems: A real-time reduced order modeling strategy,” SIAM Journal on Scientific Computing, pp. A128–A149, 2022/01/19 2022.
  * [18] D. Ramezanian, A. G. Nouri, and H. Babaee, “On-the-fly reduced order modeling of passive and reactive species via time-dependent manifolds,” Computer Methods in Applied Mechanics and Engineering, vol. 382, p. 113882, 2021.
  * [19] M. H. Naderi and H. Babaee, “Adaptive sparse interpolation for accelerating nonlinear stochastic reduced-order modeling with time-dependent bases,” Computer Methods in Applied Mechanics and Engineering, vol. 405, p. 115813, 2023.
  * [20] M. Donello, G. Palkar, M. H. Naderi, D. C. Del Rey Fernández, and H. Babaee, “Oblique projection for scalable rank-adaptive reduced-order modelling of nonlinear stochastic partial differential equations with time-dependent bases,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 479, p. 20230320, 2023/10/19 2023.
  * [21] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM Journal on Matrix Analysis and Applications, vol. 21, pp. 1253–1278, 2020/07/11 2000.
  * [22] S. Ahmadi-Asl, S. Abukhovich, M. G. Asante-Mensah, A. Cichocki, A. H. Phan, T. Tanaka, and I. Oseledets, “Randomized algorithms for computation of Tucker decomposition and higher order SVD (HOSVD),” IEEE Access, vol. 9, pp. 28684–28706, 2021.
  * [23] A. K. Saibaba, “HOID: Higher order interpolatory decomposition for tensors based on Tucker representation,”
  * [24] N. Vannieuwenhoven, R. Vandebril, and K. Meerbergen, “A new truncation strategy for the higher-order singular value decomposition,” SIAM Journal on Scientific Computing, vol. 34, pp. A1027–A1052, 2020/07/06 2012.
  * [25] S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, “A theory of pseudoskeleton approximations,” Linear Algebra and its Applications, vol. 261, no. 1, pp. 1–21, 1997.
  * [26] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, How to Find a Good Submatrix, pp. 247–256.
  * [27] S. Goreinov and E. Tyrtyshnikov, The maximal-volume concept in approximation by low-rank matrices, pp. 47–51.  2001\.
  * [28] D. Savostyanov, Polilinear approximation of matrices and integral equations.  PhD thesis, Dept. Math., INM RAS, Moscow, Russia, 2006.
  * [29] E. Tyrtyshnikov, “Incomplete cross approximation in the mosaic-skeleton method,” Computing, vol. 64, no. 4, pp. 367–380, 2000.
  * [30] M. W. Mahoney, “Randomized algorithms for matrices and data,” Foundations and Trends® in Machine Learning, vol. 3, no. 2, pp. 123–224, 2011.
  * [31] D. C. Sorensen and M. Embree, “A DEIM induced CUR factorization,” SIAM Journal on Scientific Computing, vol. 38, no. 3, pp. A1454–A1482, 2016.
  * [32] C. F. Caiafa and A. Cichocki, “Generalizing the column–row matrix decomposition to multi-way arrays,” Linear Algebra and its Applications, vol. 433, no. 3, pp. 557–573, 2010.
  * [33] I. Oseledets and E. Tyrtyshnikov, “TT-cross approximation for multidimensional arrays,” Linear Algebra and its Applications, vol. 432, no. 1, pp. 70–88, 2010.
  * [34] J. Ballani, L. Grasedyck, and M. Kluge, “Black box approximation of tensors in hierarchical Tucker format,” Linear Algebra and its Applications, vol. 438, no. 2, pp. 639–657, 2013.
  * [35] S. Ahmadi-Asl, C. F. Caiafa, A. Cichocki, A. H. Phan, T. Tanaka, I. Oseledets, and J. Wang, “Cross tensor approximation methods for compression and dimensionality reduction,” IEEE Access, vol. 9, pp. 150809–150838, 2021.
  * [36] G. Ceruti and C. Lubich, “Time integration of symmetric and anti-symmetric low-rank matrices and Tucker tensors,” BIT Numerical Mathematics, vol. 60, no. 3, pp. 591–614, 2020.
  * [37] G. Ceruti, J. Kusch, and C. Lubich, “A rank-adaptive robust integrator for dynamical low-rank approximation,” BIT Numerical Mathematics, vol. 62, no. 4, pp. 1149–1174, 2022.
  * [38] K. Manohar, B. W. Brunton, J. N. Kutz, and S. L. Brunton, “Data-driven sparse sensor placement for reconstruction: Demonstrating the benefits of exploiting known patterns,” IEEE Control Systems Magazine, vol. 38, no. 3, pp. 63–86, 2018.
  * [39] S. Chaturantabut and D. C. Sorensen, “Nonlinear model reduction via discrete empirical interpolation,” SIAM Journal on Scientific Computing, vol. 32, no. 5, pp. 2737–2764, 2010.
  * [40] S. Karlin and H. E. Taylor, A Second Course in Stochastic Processes.  New York, NY: Academic Press, 1981.