# Interpolatory tensorial reduced order models for parametric dynamical
systems

Alexander V. Mamonov Department of Mathematics, University of Houston,
Houston, Texas 77204 (avmamonov@uh.edu). Maxim A. Olshanskii Department of
Mathematics, University of Houston, Houston, Texas 77204
(maolshanskiy@uh.edu).

###### Abstract

The paper introduces a reduced order model (ROM) for numerical integration of
a dynamical system which depends on multiple parameters. The ROM is a
projection of the dynamical system on a low dimensional space that is both
problem-dependent and _parameter-specific_. The ROM exploits compressed tensor
formats to find a low rank representation for a sample of high-fidelity
snapshots of the system state. This tensorial representation provides ROM with
an orthogonal basis in a universal space of all snapshots _and_ encodes
information about the state variation in parameter domain. During the online
phase and for any incoming parameter, this information is used to find a
reduced basis that spans a parameter-specific subspace in the universal space.
The computational cost of the online phase then depends only on tensor
compression ranks, but not on space or time resolution of high-fidelity
computations. Moreover, certain compressed tensor formats enable to avoid the
adverse effect of parameter space dimension on the online costs (known as the
curse of dimension). The analysis of the approach includes an estimate for the
representation power of the acquired ROM basis. We illustrate the performance
and prediction properties of the ROM with several numerical experiments, where
tensorial ROM’s complexity and accuracy is compared to those of conventional
POD-ROM.

###### keywords:

Model order reduction, parametric PDEs, low-rank tensors, dynamical systems,
proper orthogonal decomposition

##  1 Introduction

In numerical optimal control, inverse modeling or uncertainty quantification,
one commonly needs to integrate a parameter-dependent dynamical system for
various values of the parameter vector. For example, inverse modeling may
require repeated solutions of the forward problem represented by a dynamical
system, along the search path in a high-dimensional parameter space. This may
lead to extreme-scale computations that, if implemented straightforwardly,
often result in overwhelming computational costs. Reduced order models (ROMs)
offer a possibility to alleviate these costs by replacing a high-fidelity
model with a low-dimensional surrogate model [3, 33]. Thanks to this practical
value and apparent success, ROMs for parametric dynamical systems have already
attracted considerable attention; see, e.g., [10, 39, 14, 5, 9, 13].

In this paper, we are interested in projection based ROMs that build the
surrogate model by projecting a high-fidelity model onto a low-dimensional
problem-dependent vector space [10]. Projection-based ROMs for dynamical
systems include such well-known model order reduction techniques as proper
orthogonal decomposition (POD) ROMs [54, 64] (and its variants such as POD-
DEIM [17] and balanced POD [61]) and PGD-ROMs [18, 19]. In these approaches
the basis for the projection space is computed by building on the information
about the dynamical system provided through high-fidelity solutions sampled
for certain time instances and/or parameter values, the so-called solution
snapshots. However, building a general low-dimensional space for all times and
parameters of interest might be challenging, if possible at all, for wide
parameter ranges and long times. Several studies aimed to address this
challenge: In [26, 25, 2] the authors considered partitioning strategies which
introduce a subdivision of the parameter domain and assign an individual local
reduced-order basis to each subdomain offline. Another idea [1, 65] is to
adapt precomputed equal-dimension reduced order spaces by interpolating them
for out-of-sample parameters along geodesics on the Grassmann manifold. The
present paper introduces a different approach that builds on recent
developments in tensor decompositions and low-rank approximations to quickly
compute parameter-specific reduced bases for projection based ROMs.

For parametrized systems of time-dependent differential equations, the
generated data (the input of ROM) naturally takes a form of a multi-
dimensional tensor of solution snapshots, with dimensionality $D+2$, where $D$
is a the dimension of the parameter space and 2 accounts for the spatial and
time-wise distributions. Modern ROMs often proceed by unfolding such tensors
into a matrix to perform standard POD based on truncated SVD. This leads to
the loss of information about the dependency of solutions on parameters. We
propose to overcome these issues by working directly with tensor data, and by
exploiting low-rank tensor approximations based on the canonical polyadic
(CP), high order SVD (HOSVD), and Tensor Train (TT) decompositions. The
approach consists of two stages. First, at the offline stage, the compressed
snapshot tensor is computed using one of the three tensor decompositions, thus
preserving the essential information about variation of the solution with
respect to parameters. Up to the compression accuracy, each of these
decompositions provides a (global) basis for the universal space spanned by
all observed snapshots. The so-called core of the compressed representation is
then transmitted to the second stage, referred to as the online stage. At the
online stage the transmitted part of compressed tensor allows for a fast
computation of a _parameter-specific_ reduced basis for any incoming out-of-
sample parameter vector through an interpolation and fast linear algebra
routines. The reduced order basis is then given in terms of its coordinates in
the global basis that can be stored offline. For CP and TT formats, the cost
of these computations is free of exponential growth with respect to the
parameter space dimension. On analysis side of this work, we prove an estimate
for prediction power of the parameter-specific reduced order basis. The
estimate explicitly depends on the approximation accuracy of the original
tensor by the compressed one, parameter interpolation error, and singular
values of a small-size parameter-specific matrix.

Despite an outstanding recent progress in numerical multi-linear algebra and,
in particular, in understanding tensor decompositions (see, e.g., review
articles [46, 29, 63]), the application of tensor methods in reduced order
modeling of dynamical systems is still rather scarce. We mention two reports
by Nouy [55, 56], who reviewed tensor compressed formats and discussed their
possible use for sparse function representation and reduced order modeling, as
well as a series of publications on the treatment in compressed tensor formats
of algebraic systems resulting from the stochastic and parametric Galerkin
finite element method, see e.g. [11, 7, 8, 49, 48]. A POD-ROM was combined
with a low-rank tensor representation of a mapping from a parameter space onto
an output domain in [42]. The authors of survey [10] observe that “The
combination of tensor calculus $\ldots$ and parametric model reduction
techniques for time dependent problems is still in its infancy, but offers a
promising research direction”. We believe the statement holds true, and the
present study contributes to this largely open research field.

The remainder of the paper is organized as follows. In Section 2 we set up a
parameter-dependent Cauchy problem and recall the basics of POD-ROM approach
that is needed for reference purpose later in the text. Section 3 introduces a
general idea of the interpolatory tensorial ROM and considers its realization
using three popular tensor compression formats. Details are worked out for a
Cartesian grid-based sampling of the parameter domain, and then the approach
is extended to a more general parameter sampling scheme. A separate subsection
discusses online–offline complexity and storage requirements of the method. An
estimate on the prediction power of the reduced order basis is proved in
Section 4. Numerical examples in Section 5 illustrate the analysis and
performance of the method. In particular, we compare the delivered accuracy
with standard POD-ROM that employs a global low-dimensional basis.

##  2 Parameterized Cauchy problem and the conventional POD-ROM

To fix ideas, consider the following multi-parameter initial value problem.
For a vector of parameters
$\mbox{\boldmath$\alpha$\unboldmath}=(\alpha_{1},\dots,\alpha_{D})$ from the
parameter domain $\mathcal{A}\subset\mathbb{R}^{D}$ find the trajectory
$\mathbf{u}=\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath}):[0,T)\to\mathbb{R}^{M}$
solving

$\mathbf{u}_{t}=F(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath}),\quad
t\in(0,T),\quad\text{and}~{}\mathbf{u}|_{t=0}=\mathbf{u}_{0},$ (1)

with a given continuous flow field
$F:(0,T)\times\mathbb{R}^{M}\times\mathcal{A}\to\mathbb{R}^{M}$. Hereafter we
denote all vector quantities by bold lowercase letters. We assume that the
unique solution exists on $(0,T)$ for all
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$. Examples considered in
this paper include parameter-dependent parabolic equations, in which case one
can think of (1) as a system of ODEs for nodal values of the finite volume or
finite element solution to the PDE problem, where material coefficients,
boundary conditions, or the computational domain (via a mapping into a
reference domain) are parameterized by $\alpha$.

We are interested in projection based ROMs, where for an arbitrary but fixed
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ an approximation to
$\mathbf{u}$ is sought as a solution to equations projected onto a reduced
space. Projection based approaches aim at retaining the structure of the model
and thus at preserving the physics present in the high-fidelity model [10].
Among the projection based approaches to model reduction for time-dependent
differential equations, Proper Orthogonal Decomposition (POD) and its variants
are likely the most widely used ROM technique, which provides tools to
represent trajectories of a dynamical system in a low-dimensional, problem-
dependent basis [43, 60, 50, 52]. We summarize the POD-ROM below for further
reference and for the purpose of comparison to our approach in Section 5.

Assume for a moment that $\alpha$ is _fixed_. The POD-ROM computes a
representative collection of states
$\boldsymbol{\phi}_{k}(\mbox{\boldmath$\alpha$\unboldmath})=\mathbf{u}(t_{k},\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M}$
at times $0\leq t_{1},\dots,t_{N}<T$, referred to as snapshots, through high-
fidelity numerical simulations. Next, one finds a parameter-specific low-
dimensional basis $\\{\mathbf{z}_{i}^{\rm
pod}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}\subset\mathbb{R}^{M}$,
$n\ll N$, referred to hereafter as the _reduced basis_ , such that the
projection subspace $\mbox{span}\big{\\{}\mathbf{z}_{1}^{\rm
pod}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\mathbf{z}_{n}^{\rm
pod}(\mbox{\boldmath$\alpha$\unboldmath})\big{\\}}$ approximates the snapshot
space
$\mbox{span}\\{\boldsymbol{\phi}_{1}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\boldsymbol{\phi}_{N}(\mbox{\boldmath$\alpha$\unboldmath})\\}$
in the best possible way.

To determine the reduced basis, form a matrix of snapshots

$\Phi_{\text{pod}}(\mbox{\boldmath$\alpha$\unboldmath})=[\boldsymbol{\phi}_{1}(\mbox{\boldmath$\alpha$\unboldmath}),\ldots,\boldsymbol{\phi}_{N}(\mbox{\boldmath$\alpha$\unboldmath})]\in\mathbb{R}^{M\times
N},$ (2)

compute its SVD

$\Phi_{\text{pod}}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{U}\Sigma\mathrm{V}^{T},$
(3)

and define $\mathbf{z}_{i}^{\rm pod}(\mbox{\boldmath$\alpha$\unboldmath})$,
$i=1,\dots,n$, to be the first $n$ left singular vectors of
$\Phi_{\text{pod}}(\mbox{\boldmath$\alpha$\unboldmath})$, i.e., the first $n$
columns of $\mathrm{U}$. Hereafter we denote all matrices with upright capital
letters. The singular values in $\Sigma$ provide information about the
approximation power of $\mbox{span}\big{\\{}\mathbf{z}_{1}^{\rm
pod}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\mathbf{z}_{n}^{\rm
pod}(\mbox{\boldmath$\alpha$\unboldmath})\big{\\}}$. We refer to [51, 43] and
references therein for a discussion about algebraically different ways to
define POD and their equivalence.

For parameters $\alpha$ varying in $\mathcal{A}$, a parametric POD-ROM builds
a global reduced basis by sampling the parameter domain, generating snapshots
for each sampled parameter value and proceeding with SVD (3) for a cumulative
matrix of all snapshots. Possible sampling strategies include using a
Cartesian grid in $\mathcal{A}$, Monte–Carlo methods, and greedy algorithms
based on a posteriori error estimates; see, e.g., [10, 39]. Regardless of the
sampling procedure, the resulting basis can accurately reproduce only the data
from which it originated. Without parameter-specificity, the basis may lack
robustness for out-of-sample parameters, i.e., away from the reference
simulations. This is a serious limitation for using POD based ROMs in inverse
modeling. We plan to address this limitation by introducing tensorial
techniques for finding reduced bases that are both problem- and parameter-
specific.

##  3 Tensorial ROMs

We first consider in Section 3.1 a Cartesian grid-based sampling of the
parameter domain $\mathcal{A}$ in the case when $\mathcal{A}$ is the
$D$-dimensional box

$\mathcal{A}=\bigotimes_{i=1}^{D}[\alpha_{i}^{\min},\alpha_{i}^{\max}],$ (4)

and the sampling points are placed at the nodes of a Cartesian grid. Next, we
describe three tensorial ROMs (TROMs) based on three different tensor
decompositions, canonical polyadic (CP, Section 3.5), high order SVD (HOSVD,
Section 3.6) and tensor train (TT, Section 3.7).

###  3.1 Cartesian grid-based parameter sampling

To generate the sampling set ${\widehat{\mathcal{A}}}$, we distribute $n_{i}$
nodes $\\{\widehat{\alpha}_{i}^{j}\\}_{j=1,\dots,n_{i}}$ within each of the
intervals $[\alpha_{i}^{\min},\alpha_{i}^{\max}]$ in (4) for $i=1,\dots,D$,
and define

${\widehat{\mathcal{A}}}=\left\\{\widehat{\boldsymbol{\alpha}}=(\widehat{\alpha}_{1},\dots,\widehat{\alpha}_{D})^{T}\,:\,\widehat{\alpha}_{i}\in\\{\widehat{\alpha}_{i}^{j}\\}_{j=1,\dots,n_{i}},~{}i=1,\dots,D\right\\}.$
(5)

Hereafter we use hats to denote parameters from the sampling set
${\widehat{\mathcal{A}}}$, and the cardinality of ${\widehat{\mathcal{A}}}$ is
denoted by

$K=\prod_{i=1}^{D}n_{i}.$ (6)

The corresponding snapshots
$\boldsymbol{\phi}_{k}(\widehat{\boldsymbol{\alpha}})=\mathbf{u}(t_{k},\widehat{\boldsymbol{\alpha}})$,
$\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}$, are organized in a
multi-dimensional array

$(\boldsymbol{\Phi})_{:,i_{1},\dots,i_{D},k}=\boldsymbol{\phi}_{k}(\widehat{\alpha}_{1}^{i_{1}},\dots,\widehat{\alpha}_{D}^{i_{D}}),$
(7)

which is a tensor of order $D+2$ and size $M\times n_{1}\times\dots\times
n_{D}\times N$. We reserve the first and the last indices of
$\boldsymbol{\Phi}$ for the spatial and temporal distributions, respectively.
All tensors hereafter are denoted with bold uppercase letters.

Unfolding $\boldsymbol{\Phi}$ along the first index in an $M\times(n_{1}\cdots
n_{D})N$ matrix and applying (truncated) SVD to determine the first $n$ left
singular vectors is equivalent to the POD with grid-based parameter sampling.
The disadvantage of this approach for ROM construction is that it neglects any
information about the dependence of snapshots on parameter variation reflected
in the tensor structure of $\boldsymbol{\Phi}$. To preserve this information,
we proceed with a compressed approximation $\widetilde{\boldsymbol{\Phi}}$ of
$\boldsymbol{\Phi}$ rather than with the low rank approximation of the
unfolded matrix.

###  3.2 Tensor compression and universal space

The notion of a tensor rank and low-rank tensor approximation is ambiguous and
later in this section we consider three popular compressed tensor formats. For
now we only assume that $\widetilde{\boldsymbol{\Phi}}$ satisfies

$\big{\|}\widetilde{\boldsymbol{\Phi}}-\boldsymbol{\Phi}\big{\|}_{F}\leq\widetilde{\varepsilon}\big{\|}\boldsymbol{\Phi}\big{\|}_{F}$
(8)

for some small $\widetilde{\varepsilon}>0$, where tensor Frobenius norm is
simply

$\|\boldsymbol{\Phi}\|_{F}:=\Big{(}\sum_{j=1}^{M}\sum_{i_{1}=1}^{n_{1}}\dots\sum_{i_{D}=1}^{n_{D}}\sum_{k=1}^{N}\boldsymbol{\Phi}_{j,i_{1},\dots,i_{D},k}^{2}\Big{)}^{1/2}$
(9)

The ”low-rank” (compressed) tensor $\widetilde{\boldsymbol{\Phi}}$ is computed
during the first, offline stage of TROM construction and _a part of_
$\widetilde{\boldsymbol{\Phi}}$ is passed on to the second, online stage which
uses this information about variation of snapshots with respect to changes in
parameters to compute a parameter-specific TROM.

We call _universal space_ the space $\widetilde{V}$ spanned by the first-mode
fibers of $\widetilde{\boldsymbol{\Phi}}$, i.e., $\widetilde{V}$ is the column
space of the mode-1 unfolding matrix. For the exact decomposition (i.e., for
$\widetilde{\varepsilon}=0$), $\widetilde{V}$ is the space of all observed
system states. In general, $\widetilde{V}$ depends on a compression format,
dimension of $\widetilde{V}$ does not exceed $M$ and depends on
$\widetilde{\varepsilon}$ and snapshot variation. We shall see that
$\widetilde{V}$ does approximate the full space of high-fidelity snapshots,
while the CP, HOSVD and TT formats all deliver an orthogonal basis for
$\widetilde{V}$. In the online stage of TROM we find a local (parameter-
specific) ROM basis by specifying its coordinates in $\widetilde{V}$.

###  3.3 In-sample prediction

During the online stage we wish to be able to approximately solve (1) for an
arbitrary parameter $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ in an
_$\alpha$ -specific_ reduced basis. For this step we need to introduce the
notion of $k$-mode tensor-vector product
$\boldsymbol{\Psi}\times_{k}\mathbf{a}$ of a tensor
$\boldsymbol{\Psi}\in\mathbb{R}^{N_{1}\times\dots\times N_{m}}$ of order $m$
and a vector $\mathbf{a}\in\mathbb{R}^{N_{k}}$: the resulting tensor
$\boldsymbol{\Psi}\times_{k}\mathbf{a}$ has order $m-1$ and size
$N_{1}\times\dots\times N_{k-1}\times N_{k+1}\times\dots\times N_{m}$.
Specifically, elementwise

$(\boldsymbol{\Psi}\times_{k}\mathbf{a})_{j_{1},\dots,j_{k-1},j_{k+1},\dots,j_{m}}=\sum_{j_{k}=1}^{N_{k}}\boldsymbol{\Psi}_{j_{1},\dots,j_{m}}a_{j_{k}}.$
(10)

For a moment, consider some
$\widehat{\boldsymbol{\alpha}}=\left(\widehat{\alpha}_{1},\dots,\widehat{\alpha}_{D}\right)^{T}$
from the sampling set ${\widehat{\mathcal{A}}}$ and define $D$ vectors,
$\mathbf{e}^{i}(\widehat{\boldsymbol{\alpha}})=\left(e_{1}^{i}(\widehat{\boldsymbol{\alpha}}),\dots,e_{n_{i}}^{i}(\widehat{\boldsymbol{\alpha}})\right)^{T}\in\mathbb{R}^{n_{i}}$,
$i=1,\dots,D$, as

$e^{i}_{j}(\widehat{\boldsymbol{\alpha}})=\begin{cases}1&\text{if}~{}\widehat{\alpha}_{i}=\widehat{\alpha}_{i}^{j}\\\
0&\text{otherwise}\end{cases},\qquad j=1,\dots,n_{i}.$ (11)

In other words, $\mathbf{e}^{i}(\widehat{\boldsymbol{\alpha}})$ encodes the
position of $\widehat{\alpha}_{i}$ among the grid nodes on
$[\alpha_{i}^{\min},\alpha_{i}^{\max}]$, $i=1,\ldots,D$.

Vectors $\mathbf{e}^{i}(\widehat{\boldsymbol{\alpha}})$ defined above allow us
to extract the snapshots corresponding to a particular
$\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}$. Specifically, we
introduce the following _extraction_ operation

$\Phi_{e}(\widehat{\boldsymbol{\alpha}})=\boldsymbol{\Phi}\times_{2}\mathbf{e}^{1}(\widehat{\boldsymbol{\alpha}})\times_{3}\mathbf{e}^{2}(\widehat{\boldsymbol{\alpha}})\dots\times_{D+1}\mathbf{e}^{D}(\widehat{\boldsymbol{\alpha}})\in\mathbb{R}^{M\times
N},$ (12)

which extracts from tensor of all snapshots $\boldsymbol{\Phi}$ the matrix of
snapshots (2) for the particular
$\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}$, i.e.,
$\Phi_{e}(\widehat{\boldsymbol{\alpha}})=\Phi_{\text{pod}}(\widehat{\boldsymbol{\alpha}})$.

Combining (12) with compressed approximation (8), we conclude that is should
be possible to extract from $\widetilde{\boldsymbol{\Phi}}$ the information
about the space spanned by the snapshots
$\\{\mathbf{u}(t_{i},\widehat{\boldsymbol{\alpha}})\\}_{i=1}^{N}$, for a
particular $\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}$ up to the
accuracy of approximation in (8). Indeed, let
$\boldsymbol{\phi}_{i}(\widehat{\boldsymbol{\alpha}})=\mathbf{u}(t_{i},\widehat{\boldsymbol{\alpha}})$,
$i=1,\dots,N$, and denote by
$\\{\mathbf{z}_{j}(\widehat{\boldsymbol{\alpha}})\\}_{j=1}^{\widetilde{N}}$,
$\widetilde{N}\leq N$, an orthonormal basis for the column space of

$\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})=\widetilde{\boldsymbol{\Phi}}\times_{2}\mathbf{e}^{1}(\widehat{\boldsymbol{\alpha}})\times_{3}\mathbf{e}^{2}(\widehat{\boldsymbol{\alpha}})\dots\times_{D+1}\mathbf{e}^{D}(\widehat{\boldsymbol{\alpha}})\in\mathbb{R}^{M\times
N},$ (13)

where
$\widetilde{N}=\mbox{rank}\big{(}\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})\big{)}$.
Then, it holds

$\sum_{i=1}^{N}\left\|\boldsymbol{\phi}_{i}-\sum_{j=1}^{\widetilde{N}}\langle\boldsymbol{\phi}_{i},\mathbf{z}_{j}\rangle\mathbf{z}_{j}\right\|^{2}_{\ell^{2}}\leq\widetilde{\varepsilon}^{2}\big{\|}\boldsymbol{\Phi}\big{\|}_{F}^{2},$
(14)

where we use the shortcut
$\boldsymbol{\phi}_{i}=\boldsymbol{\phi}_{i}(\widehat{\boldsymbol{\alpha}})$,
$\mathbf{z}_{i}=\mathbf{z}_{i}(\widehat{\boldsymbol{\alpha}})$. To establish
(14), consider (thin) SVD
$\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})=\widetilde{\mathrm{U}}\widetilde{\Sigma}\widetilde{\mathrm{V}}^{T}$
and compute

$\begin{split}\sum_{i=1}^{N}\left\|\boldsymbol{\phi}_{i}-\sum_{j=1}^{\widetilde{N}}\langle\boldsymbol{\phi}_{i},\mathbf{z}_{j}\rangle\mathbf{z}_{j}\right\|^{2}_{\ell^{2}}&=\left\|(\mathrm{I}-\widetilde{\mathrm{U}}\widetilde{\mathrm{U}}^{T}){\Phi}_{e}(\widehat{\boldsymbol{\alpha}})\right\|^{2}_{F}\\\
&=\left\|(\mathrm{I}-\widetilde{\mathrm{U}}\widetilde{\mathrm{U}}^{T})\left({\Phi}_{e}(\widehat{\boldsymbol{\alpha}})-\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})\right)\right\|^{2}_{F}\\\
&\leq\left\|\mathrm{I}-\widetilde{\mathrm{U}}\widetilde{\mathrm{U}}^{T}\right\|^{2}\left\|{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})-\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})\right\|^{2}_{F}\leq\left\|\Phi_{e}(\widehat{\boldsymbol{\alpha}})-\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})\right\|^{2}_{F}\\\
&=\left\|(\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}})\times_{2}\mathbf{e}^{1}(\widehat{\boldsymbol{\alpha}})\times_{3}\mathbf{e}^{2}(\widehat{\boldsymbol{\alpha}})\dots\times_{D+1}\mathbf{e}^{D}(\widehat{\boldsymbol{\alpha}})\right\|^{2}_{F}\\\
&\leq\left\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\right\|^{2}_{F}\leq\widetilde{\varepsilon}^{2}\big{\|}\boldsymbol{\Phi}\big{\|}_{F}^{2},\end{split}$

where we used linearity of (10) and the inequality
$\|\mathrm{A}\mathrm{B}\|_{F}\leq\|\mathrm{A}\|\|\mathrm{B}\|_{F}$ for
matrices $\mathrm{A}$, $\mathrm{B}$, and spectral matrix norm $\|\cdot\|$. We
also used $\|\mathrm{P}\|\leq 1$ for an orthogonal projection matrix
$\mathrm{P}=\mathrm{I}-\widetilde{\mathrm{U}}\widetilde{\mathrm{U}}^{T}$.

Since $\mathbf{z}_{i}(\widehat{\boldsymbol{\alpha}})\in\widetilde{V}$ for all
in-sample $\widehat{\boldsymbol{\alpha}}$, the bound in (14) provides an
estimate on how accurate the true snapshots can be approximated in the
universal space. Furthermore, from (14) we conclude that given
$\widetilde{\boldsymbol{\Phi}}$ and
$\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}$ we can obtain a
parameter-specific (quasi)-optimal reduced basis by taking the first $n$ left
singular vectors of $\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})$.
Representation power of this basis is determined by $\widetilde{\varepsilon}$
from (8) and $\sigma_{i}$, $i>n$, the singular values of
$\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})$. If
$\widetilde{\varepsilon}$ is sufficiently small, i.e., the snapshot tensor
admits an efficient low-rank representation, then the computed reduced basis
better represents the snapshot space for a given
$\widehat{\boldsymbol{\alpha}}$ than the first $n$ left singular vectors of
the unfolded snapshot matrix (i.e., better than the POD basis); see numerical
results in Section 5 which show up to several orders of accuracy gain for some
of the examples.

The arguments above apply only to parameter values
$\widehat{\boldsymbol{\alpha}}$ from the sampling set
${\widehat{\mathcal{A}}}\subset\mathcal{A}$. Next, we consider ROM basis
computation for an arbitrary
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ that may not necessarily
come from the training set ${\widehat{\mathcal{A}}}$, the so-called _out-of-
sample_ $\alpha$. Below we explore the option of building ROM basis for an
out-of-sample $\alpha$ using interpolation in the parameter space. This
approach is based on an assumption of smooth dependence of the solution
$\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath})$ of (1) on $\alpha$. We
refer to the corresponding tensorial ROMs as interpolatory TROMs.

###  3.4 Interpolatory TROM

To construct parameter-specific ROM basis for an arbitrary
$\mbox{\boldmath$\alpha$\unboldmath}=(\alpha_{1},\dots,\alpha_{D})^{T}\in\mathcal{A}$
we introduce the interpolation procedure defined by

$\mathbf{e}^{i}\,:\,\mbox{\boldmath$\alpha$\unboldmath}\to\mathbb{R}^{n_{i}},\quad
i=1,\dots,D.$ (15)

Entrywise, we write
$\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})=\big{(}e_{1}^{i}(\mbox{\boldmath$\alpha$\unboldmath}),\ldots,e_{n_{i}}^{i}(\mbox{\boldmath$\alpha$\unboldmath})\big{)}^{T}\in\mathbb{R}^{n_{i}}$.
The interpolation procedure should satisfy the following property. For a
smooth function ${g}:[\alpha_{i}^{\min},\alpha_{i}^{\max}]\to\mathbb{R}$, it
holds

${g}(\alpha_{i})\approx\sum_{j=1}^{n_{i}}e_{j}^{i}(\mbox{\boldmath$\alpha$\unboldmath}){g}(\widehat{\alpha}_{i}^{j}),\quad
i=1,\ldots,D,$ (16)

where $\widehat{\alpha}_{i}^{j}$, $j=1,\ldots,n_{i}$, are the grid nodes on
$[\alpha_{i}^{\min},\alpha_{i}^{\max}]$.

We further consider Lagrangian interpolation of order $p$: for a given
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ let
$\widehat{\alpha}_{i}^{i_{1}},\ldots,\widehat{\alpha}_{i}^{i_{p}}$ be the $p$
closest grid nodes to $\alpha_{i}$ on $[\alpha_{i}^{\min},\alpha_{i}^{\max}]$,
for $i=1,\ldots,D$. Then

$e_{j}^{i}(\mbox{\boldmath$\alpha$\unboldmath})=\begin{cases}\prod\limits_{\begin{subarray}{c}m=1,\\\
m\neq
k\end{subarray}}^{p}(\widehat{\alpha}_{i}^{i_{m}}-\alpha_{i})\Big{/}\prod\limits_{\begin{subarray}{c}m=1,\\\
m\neq
k\end{subarray}}^{p}(\widehat{\alpha}_{i}^{i_{m}}-\widehat{\alpha}_{i}^{j}),&\text{if
}j=i_{k}\in\\{i_{1},\ldots,i_{p}\\},\\\
\qquad\qquad\qquad\qquad\qquad\qquad\qquad 0,&\text{otherwise},\end{cases}$
(17)

are the entries of $\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})$ for
$j=1,\dots,n_{i}$. For the numerical experiments in Section 5 we use $p=2$ or
$3$, i.e., linear or quadratic interpolation. The Lagrangian interpolation is
not the only option and depending on parameter sampling and solution
smoothness other fitting procedures can be more suitable.

Vectors $\mathbf{e}^{i}$ extend the notion of position vectors
$\mathbf{e}^{i}$ defined in (11) for out-of-sample vectors. Indeed, from (17)
it is easy to see that both vectors coincide if
$\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}$ and so we use the
same notation hereafter. Therefore, we can define a snapshot matrix
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ through the
extraction–interpolation procedure:

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\widetilde{\boldsymbol{\Phi}}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\times_{3}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\dots\times_{D+1}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M\times
N},$ (18)

to generalize (12) for any
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$. Note that (18) and (12)
are the same for
$\mbox{\boldmath$\alpha$\unboldmath}=\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}\subset\mathcal{A}$,
while (18) defines $\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$
also for out-of-sample parameter vectors. If the low-rank representation of
the snapshot tensor is exact, i.e.,
$\widetilde{\boldsymbol{\Phi}}=\boldsymbol{\Phi}$, then
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ is the
interpolation of the snapshot matrices ${\Phi}_{\rm
pod}(\widehat{}\mbox{\boldmath$\alpha$\unboldmath})$.

Once $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ is fixed and
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ is given by (18),
our _parameter-specific reduced basis_
$\\{\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}$ is
defined as the first $n$ left singular vectors of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$. Later we
demonstrate that the coordinates of this basis in the universal space can be
calculated quickly (i.e., using only low-dimensional calculations) online
without actually computing
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$.

In a non-interpolatory TROM, a parameter-specific reduced basis can be
constructed as follows. Choose $p\geq 2$ and fix
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, then let
$\widehat{\alpha}_{i}^{i_{1}},\ldots,\widehat{\alpha}_{i}^{i_{p}}$ be the $p$
closest grid nodes to $\alpha_{i}$ on $[\alpha_{i}^{\min},\alpha_{i}^{\max}]$,
for $i=1,\ldots,D$, similarly to the interpolatory construction above. Define
the set

${\widehat{\mathcal{A}}}_{p}:=\left\\{\widehat{\boldsymbol{\alpha}}=(\widehat{\alpha}_{1},\dots,\widehat{\alpha}_{D})^{T}\,:\,\widehat{\alpha}_{i}\in\\{\widehat{\alpha}_{i}^{j}\\}_{j\in\\{i_{1},\ldots,i_{p}\\}},~{}i=1,\dots,D\right\\}\subset{\widehat{\mathcal{A}}}.$
(19)

Then, assemble a large matrix by concatenating
$\widetilde{\Phi}_{e}(\widehat{\boldsymbol{\alpha}})$ for all
$\widehat{\boldsymbol{\alpha}}\in{\widehat{\mathcal{A}}}_{p}$ and take
$\\{\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}$ to be
its first $n$ left singular vectors. Of course, hybrid strategies (e.g.,
interpolation only in some parameter directions) are also possible. For non-
interpolatory or hybrid TROMs it is also possible to compute local basis
online with only low-dimensional calculations following same steps as
considered below.

In the rest of paper we focus on the interpolatory TROM and consider three
well-known compressed formats for low rank tensor approximation
$\widetilde{\boldsymbol{\Phi}}\approx\boldsymbol{\Phi}$: canonical polyadic
(CP), Tucker, a.k.a. higher order SVD (HOSVD), and tensor train (TT)
decomposition formats. Note that the notion of tensor rank(s) differs among
these formats. When applied to TROM computation, these formats lead to
different offline computational costs to build
$\widetilde{\boldsymbol{\Phi}}$, different amounts of information transmitted
from the offline stage to the online stage (measured by the compression rate,
as explained in Section 3.9), and slightly varying amounts of online
computations for finding the reduced basis given an incoming
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$.

###  3.5 CP-TROM

The first tensor decomposition that we consider is the canonical polyadic
decomposition of a tensor into the sum of rank one tensors [40, 16, 45, 46].
In CP-TROM we approximate $\boldsymbol{\Phi}$ by the sum of $R$ (where $R$ is
the so-called canonical tensor rank) direct products of $D+2$ vectors
$\mathbf{u}^{r}\in\mathbb{R}^{M}$,
$\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i}\in\mathbb{R}^{n_{i}}$,
$i=1,\dots,D$, and $\mathbf{v}^{r}\in\mathbb{R}^{N}$,

$\boldsymbol{\Phi}\approx\widetilde{\boldsymbol{\Phi}}=\sum_{r=1}^{R}\mathbf{u}^{r}\circ\mbox{\boldmath$\sigma$\unboldmath}^{r}_{1}\circ\dots\circ\mbox{\boldmath$\sigma$\unboldmath}^{r}_{D}\circ\mathbf{v}^{r},$
(20)

or entry-wise

$(\widetilde{\boldsymbol{\Phi}})_{j,i_{1},\dots,i_{D},k}=\sum_{r=1}^{R}u_{j}^{r}\sigma_{1,i_{1}}^{r}\dots\sigma_{D,i_{D}}^{r}v_{k}^{r}.$

CP decomposition often delivers excellent compression. However, there are
well-known difficulties in determining the accurate canonical rank $R$ and
working with the CP format, see, e.g., [38, 22]. Since we are interested in
the approximation $\widetilde{\boldsymbol{\Phi}}$ to $\boldsymbol{\Phi}$, the
alternating least squares (ALS) algorithm [36, 46] can be used to minimize
$\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\|_{F}$ for a specified
target canonical rank $R$ to find the approximate factors
$\mathbf{u}^{r}\in\mathbb{R}^{M}$,
$\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i}\in\mathbb{R}^{n_{i}}$, and
$\mathbf{v}^{r}\in\mathbb{R}^{N}$, $r=1,\ldots,R$. We further assume $R\leq
M$, where $M$ is the dimension of high-fidelity snapshots.

Note that the second-mode product of a $D+2$-dimensional rank-one tensor
$\mathbf{u}^{r}\circ\mbox{\boldmath$\sigma$\unboldmath}^{r}_{1}\circ\dots\circ\mbox{\boldmath$\sigma$\unboldmath}^{r}_{D}\circ\mathbf{v}^{r}$
and a vector
$\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{n_{1}}$ is
the $D+1$-dimensional rank-one tensor
$\left\langle\mbox{\boldmath$\sigma$\unboldmath}^{r}_{1},\mathbf{e}\right\rangle(\mathbf{u}^{r}\circ\mbox{\boldmath$\sigma$\unboldmath}^{r}_{2}\circ\dots\circ\mbox{\boldmath$\sigma$\unboldmath}^{r}_{D}\circ\mathbf{v}^{r})$.
Proceeding with this computation for other modes, in the decomposition (20) we
find that the definition (18) yields representation of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ as the sum of rank
one matrices for any $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$:

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\sum_{r=1}^{R}s_{r}\mathbf{u}^{r}\circ\mathbf{v}^{r}\in\mathbb{R}^{M\times
N},~{}~{}\text{with}~{}~{}s_{r}=\prod_{i=1}^{D}\left\langle\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i},\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})\right\rangle\in\mathbb{R}.$
(21)

However, (21) is _not_ the SVD of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$, since vectors
$\mathbf{u}^{r}$ (and $\mathbf{v}^{r}$) are not necessarily orthogonal. To
avoid computing $\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$
and its SVD online, the following preparatory offline step is required.
Organize the vectors $\mathbf{u}^{r}$ and $\mathbf{v}^{r}$ from (20) into
matrices

$\widehat{\mathrm{U}}=[\mathbf{u}^{1},\dots,\mathbf{u}^{R}]\in\mathbb{R}^{M\times
R},\quad\widehat{\mathrm{V}}=[\mathbf{v}^{1},\dots,\mathbf{v}^{R}]\in\mathbb{R}^{N\times
R},$ (22)

and compute the thin QR factorizations

$\widehat{\mathrm{U}}={\mathrm{U}}\mathrm{R}_{U},\quad\widehat{\mathrm{V}}={\mathrm{V}}\mathrm{R}_{V},$
(23)

if $R>N$ let further $\mathrm{R}_{V}=\widehat{\mathrm{V}}$ and
${\mathrm{V}}=\mathrm{I}$. _The columns of ${\mathrm{U}}$ form an orthogonal
basis in the universal space $\widetilde{V}$_. Matrix ${\mathrm{U}}$ is stored
offline (${\mathrm{V}}$ is not used and can be dropped), while low-dimensional
matrices $\mathrm{R}_{U}$ and $\mathrm{R}_{V}$ together with vectors
$\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i}$ form the _online_ part of
$\widetilde{\boldsymbol{\Phi}}$,

$\mbox{online}(\widetilde{\boldsymbol{\Phi}})=\left\\{\mathrm{R}_{U},\mathrm{R}_{V}\in\mathbb{R}^{R\times
R},~{}\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i}\in\mathbb{R}^{n_{i}},~{}{\small
i=1,\dots,D}\right\\},$ (24)

which is transmitted to the online stage.

At the online stage and for any incoming
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, we compute the SVD of the
$R\times R$ _core_ matrix

$\mathrm{C}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{R}_{U}\mathrm{S}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{R}_{V}^{T},$
(25)

where
$\mathrm{S}(\mbox{\boldmath$\alpha$\unboldmath})=\mbox{diag}(s_{1},\dots,s_{R})$,
with $s_{r}$ from (21):
$\mathrm{C}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{U}_{c}\Sigma_{c}\mathrm{V}_{c}^{T}$.
Since

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})={\mathrm{U}}\mathrm{C}(\mbox{\boldmath$\alpha$\unboldmath}){\mathrm{V}}^{T}=\left({\mathrm{U}}\mathrm{U}_{c}\right)\Sigma_{c}\left({\mathrm{V}}\mathrm{V}_{c}\right)^{T}$
(26)

is the SVD of $\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$,
_the first $n$ columns of $\mathrm{U}_{c}$, denoted by
$\left\\{\mbox{\boldmath$\beta$\unboldmath}_{1}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\mbox{\boldmath$\beta$\unboldmath}_{n}(\mbox{\boldmath$\alpha$\unboldmath})\right\\}$,
are the coordinates of the local reduced basis in the universal space
$\widetilde{V}$._ The parameter-specific basis in the physical space is then
$\\{\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}$, with
$\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})={\mathrm{U}}\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})$.
Note that $\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})$ are not
actually computed.

Under certain assumptions on $F$, the dynamical system (1) is projected
offline onto $\widetilde{V}$ and passed to the online stage, where for any
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ it is further projected
onto the local basis
$\left\\{\mbox{\boldmath$\beta$\unboldmath}_{1}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\mbox{\boldmath$\beta$\unboldmath}_{n}(\mbox{\boldmath$\alpha$\unboldmath})\right\\}$.
This avoids any online computations with high-dimensional objects used in
high-fidelity simulations; see further discussion in Section 3.10.

We summarize the above in the following algorithm.

######  Algorithm 1 (CP-TROM).

  * •

Offline stage.
_Input:_ snapshot tensor $\boldsymbol{\Phi}\in\mathbb{R}^{M\times
n_{1}\times\ldots\times n_{D}\times N}$, target canonical rank $R$;
_Output:_ CP decomposition factors, universal basis matrix
${\mathrm{U}}\in\mathbb{R}^{M\times R}$, and upper triangular matrices
$\mathrm{R}_{U}$, $\mathrm{R}_{V}\in\mathbb{R}^{R\times R}$;
_Compute:_

    1. 1.

Use ALS algorithm to minimize
$\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\|_{F}$ to find CP
decomposition factors $\mathbf{u}^{r}$,
$\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i}$, and $\mathbf{v}^{r}$ of
$\widetilde{\boldsymbol{\Phi}}$ satisfying (20);

    2. 2.

Assemble matrices $\widehat{\mathrm{V}},\widehat{\mathrm{U}}$ as in (22) and
compute their thin QR factorizations (23) to obtain $\mathrm{R}_{U}$,
$\mathrm{R}_{V}$ and ${\mathrm{U}}$.

  * •

Online stage.
_Input:_ $\mbox{\rm online}(\widetilde{\Phi})$ as defined in (24), reduced
space dimension $n\leq R$, and parameter vector
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$;
_Output:_ Coordinates of the reduced basis in $\widetilde{V}$:
$\\{\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}\subset\mathbb{R}^{R}$;
_Compute:_

    1. 1.

Use (25) to assemble the core matrix
$\mathrm{C}(\mbox{\boldmath$\alpha$\unboldmath})$;

    2. 2.

Compute the SVD of the core matrix
$\mathrm{C}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{U}_{c}\Sigma_{c}\mathrm{V}_{c}^{T}$,
with
$\mathrm{U}_{c}=[\widetilde{\mathbf{u}}_{1},\widetilde{\mathbf{u}}_{2},\ldots,\widetilde{\mathbf{u}}_{R}]$;

    3. 3.

Set
$\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})=\widetilde{\mathbf{u}}_{i}$,
$i=1,\dots,n$.

Note that we do not have direct control over ALS algorithm to enforce a priori
CP decomposition accuracy
$\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\|_{F}<\widetilde{\varepsilon}\|\boldsymbol{\Phi}\|_{F}$.
One option is to rerun the offline stage for different trial values of $R$.
Given that the offline stage is computationally expensive, this may become
prohibitive in cases where a desired accuracy $\widetilde{\varepsilon}$ must
be strictly enforced. The other two variants of TROM presented below are free
from this limitation.

###  3.6 HOSVD-TROM

As we already mentioned, truncated variant of CP decomposition is not known to
satisfy any simple minimization property (unlike the SVD decomposition for
matrices). A classical tensor decomposition, known to deliver a
(quasi)-minimization property, is the so-called higher order SVD (HOSVD) [21].
In HOSVD-TROM variant we approximate the snapshot tensor with a Tucker tensor
[68, 46] $\widetilde{\Phi}$ of the form

$\boldsymbol{\Phi}\approx\widetilde{\boldsymbol{\Phi}}=\sum_{j=1}^{\widetilde{M}}\sum_{q_{1}=1}^{\widetilde{n}_{1}}\dots\sum_{q_{D}=1}^{\widetilde{n}_{D}}\sum_{k=1}^{\widetilde{N}}(\mathbf{C})_{j,q_{1},\dots,q_{D},k}\mathbf{u}^{j}\circ\mbox{\boldmath$\sigma$\unboldmath}^{q_{1}}_{1}\circ\dots\circ\mbox{\boldmath$\sigma$\unboldmath}^{q_{D}}_{D}\circ\mathbf{v}^{k},$
(27)

with $\mathbf{u}^{j}\in\mathbb{R}^{M}$,
$\mbox{\boldmath$\sigma$\unboldmath}_{i}^{q_{i}}\in\mathbb{R}^{n_{i}}$, and
$\mathbf{v}^{k}\in\mathbb{R}^{N}$. The numbers $\widetilde{M}$,
$\widetilde{n}_{1}$, $\ldots$, $\widetilde{n}_{D}$ and $\widetilde{N}$ are
referred to as Tucker ranks of $\widetilde{\boldsymbol{\Phi}}$. The HOSVD
delivers an efficient compression of the snapshot tensor, provided the size of
the _core tensor_
$\mathbf{C}\in\mathbb{R}^{\widetilde{M}\times\widetilde{n}_{1}\times\dots\times\widetilde{n}_{D}\times\widetilde{N}}$
is (much) smaller than the size of $\boldsymbol{\Phi}$.

In what follows, it is helpful to organize the column vectors from (27) into
matrices

$\begin{split}\mathrm{U}&=[\mathbf{u}^{1},\dots,\mathbf{u}^{\widetilde{M}}]\in\mathbb{R}^{M\times\widetilde{M}},\quad\mathrm{V}=[\mathbf{v}^{1},\dots,\mathbf{v}^{\widetilde{N}}]\in\mathbb{R}^{N\times\widetilde{N}},\\\
\mathrm{S}_{i}&=[\mbox{\boldmath$\sigma$\unboldmath}^{1}_{i},\dots,\mbox{\boldmath$\sigma$\unboldmath}^{\widetilde{n}_{i}}_{i}]^{T}\in\mathbb{R}^{\widetilde{n}_{i}\times{n}_{i}},\quad
i=1,\ldots,D.\end{split}$ (28)

In contrast with CP decomposition, HOSVD computes vectors $\mathbf{u}^{j}$,
$j=1,\ldots,\widetilde{M}$, and $\mathbf{v}^{k}$, $k=1,\ldots,\widetilde{N}$,
that are orthonormal. Therefore, _the columns of $\mathrm{U}$ form an
orthogonal basis in the universal reduced space $\widetilde{V}$._ The
dimension of this space is defined by the first Tucker rank, $\mbox{\rm
dim}(\widetilde{V})=\widetilde{M}$. The information about
$\widetilde{\boldsymbol{\Phi}}$ to be transmitted to the online stage includes
matrices $\mathrm{S}_{i}$ and the core tensor $\mathbf{C}$. Explicitly,

$\mbox{online}(\widetilde{\boldsymbol{\Phi}})=\left\\{\mathbf{C}\in\mathbb{R}^{\widetilde{M}\times\widetilde{n}_{1}\times\dots\times\widetilde{n}_{D}\times\widetilde{N}},~{}\mathrm{S}_{i}\in\mathbb{R}^{n_{i}\times\widetilde{n}_{i}},~{}{\small
i=1,\dots,D}\right\\}.$ (29)

To find coordinates of the local basis for
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, define the
$\alpha$-specific _core matrix_
$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ as

$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\mathbf{C}\times_{2}\left(\mathrm{S}_{1}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\right)\times_{3}\left(\mathrm{S}_{2}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\right)\dots\times_{D+1}\left(\mathrm{S}_{D}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\right)\in\mathbb{R}^{\widetilde{M}\times\widetilde{N}}.$
(30)

Using the definition of $k$-mode product, (18) and (27), one computes

$\begin{split}\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})&=\sum_{j=1}^{\widetilde{M}}\sum_{q_{1}=1}^{\widetilde{n}_{1}}\dots\sum_{q_{D}=1}^{\widetilde{n}_{D}}\sum_{k=1}^{\widetilde{N}}(\mathbf{C})_{j,q_{1},\dots,q_{D},k}\langle\mbox{\boldmath$\sigma$\unboldmath}^{q_{1}}_{1},\mathbf{e}^{1}(\alpha)\rangle\cdot\dotsc\cdot\langle\mbox{\boldmath$\sigma$\unboldmath}^{q_{D}}_{D},\mathbf{e}^{D}(\alpha)\rangle\,\mathbf{u}^{j}\circ\mathbf{v}^{k}\\\
&=\sum_{j=1}^{\widetilde{M}}\sum_{q_{1}=1}^{\widetilde{n}_{1}}\dots\sum_{q_{D}=1}^{\widetilde{n}_{D}}\sum_{k=1}^{\widetilde{N}}(\mathbf{C})_{j,q_{1},\dots,q_{D},k}\left(\mathrm{S}_{1}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{q_{1}}\cdot\dotsc\cdot\left(\mathrm{S}_{D}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{q_{D}}\mathbf{u}^{j}\circ\mathbf{v}^{k}\\\
&=\sum_{j=1}^{\widetilde{M}}\sum_{k=1}^{\widetilde{N}}\left(\mathbf{C}\times_{2}\left(\mathrm{S}_{1}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\right)\times_{3}\left(\mathrm{S}_{2}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\right)\dots\times_{D+1}\left(\mathrm{S}_{D}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\right)\right)_{jk}\mathbf{u}^{j}\circ\mathbf{v}^{k}=\mathrm{U}\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{V}^{T}\end{split}$

Consider the thin SVD of the core matrix

$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{U}_{c}\Sigma_{c}\mathrm{V}_{c}^{T}.$
(31)

Combining this with the representation above we get

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\left({\mathrm{U}}\mathrm{U}_{c}\right)\Sigma_{c}\left({\mathrm{V}}\mathrm{V}_{c}\right)^{T},$
(32)

which is the thin SVD of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ since both
matrices $\mathrm{U}$ and $\mathrm{V}$ are orthogonal. We conclude that _the
coordinates
$\left\\{\mbox{\boldmath$\beta$\unboldmath}_{1}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\mbox{\boldmath$\beta$\unboldmath}_{n}(\mbox{\boldmath$\alpha$\unboldmath})\right\\}$
of the local reduced basis in the universal space $\widetilde{V}$ are the
first $n$ columns of $\mathrm{U}_{c}$ from (31)._ The parameter-specific basis
is then $\\{\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}$,
with
$\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})={\mathrm{U}}\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})$
(not actually computed at the online stage).

To compute the low-rank HOSVD approximation (27) we employ the standard
algorithm [21] based on repeated computations of truncated SVD for unfolded
matrices. In particular, one may compute $\widetilde{\boldsymbol{\Phi}}$ with
either prescribed Tucker ranks or prescribed accuracy
$\widetilde{\varepsilon}$. Moreover, for fixed Tucker ranks one can show that
the recovered $\widetilde{\boldsymbol{\Phi}}$ satisfies a quasi-minimization
property [21] of the form

$\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\|\leq\sqrt{D+2}\|\boldsymbol{\Phi}-\boldsymbol{\Phi}^{\rm
opt}\|\quad\text{and}\quad\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\|\leq\left(\sum_{i=1}^{D+1}{\Delta_{i}^{2}}\right)^{\frac{1}{2}},$
(33)

where $\boldsymbol{\Phi}^{\rm opt}$ is the best approximation to
$\boldsymbol{\Phi}$ among all Tucker tensors of the given rank (such
approximation always exists), and ${\Delta_{i}}$ measures truncated SVD error
on the $i$-th step of the HOSVD. We summarize the above in the following
algorithm.

######  Algorithm 2 (HOSVD-TROM).

  * •

Offline stage.
_Input:_ snapshot tensor $\boldsymbol{\Phi}\in\mathbb{R}^{M\times
n_{1}\times\ldots\times n_{D}\times N}$ and target accuracy
$\widetilde{\varepsilon}$;
_Output:_ Compressed tensor ranks, HOSVD decomposition matrices as in (28),
and core tensor $\mathbf{C}$;
_Compute:_ Use algorithm [21] with prescribed accuracy
$\widetilde{\varepsilon}$ to compute HOSVD decomposition matrices and the core
tensor.

  * •

Online stage.
_Input:_ $\mbox{\rm online}(\widetilde{\boldsymbol{\Phi}})$ as defined in
(29), reduced space dimension $n\leq\min\\{\widetilde{M},\widetilde{N}\\}$,
and parameter vector $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$;
_Output:_ Coordinates of the reduced basis in $\widetilde{V}$:
$\\{\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}\subset\mathbb{R}^{\widetilde{M}}$;
_Compute:_

    1. 1.

Use the core tensor $\mathbf{C}$ and matrices $\mathrm{S}_{i}$,
$i=1,\ldots,D$, to assemble the core matrix
$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{\widetilde{M}\times\widetilde{N}}$
as in (30);

    2. 2.

Compute the SVD of the core matrix
$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{U}_{c}\Sigma_{c}\mathrm{V}_{c}^{T}$
with
$\mathrm{U}_{c}=[\widetilde{\mathbf{u}}_{1},\ldots,\widetilde{\mathbf{u}}_{\widetilde{M}}]$;

    3. 3.

Set
$\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})=\widetilde{\mathbf{u}}_{i}$,
$i=1,\ldots,n$.

###  3.7 TT-TROM

A third low-rank tensor decomposition of interest is the _Tensor Train_ (TT)
decomposition [58]. In TT-TROM we seek to approximate the snapshot tensor with
$\widetilde{\boldsymbol{\Phi}}$ in the TT-format, namely

$\boldsymbol{\Phi}\approx\widetilde{\boldsymbol{\Phi}}=\sum_{j_{1}=1}^{\widetilde{r}_{1}}\dots\sum_{j_{D+1}=1}^{\widetilde{r}_{D+1}}\mathbf{u}^{j_{1}}\circ\mbox{\boldmath$\sigma$\unboldmath}^{j_{1},j_{2}}_{1}\circ\dots\circ\mbox{\boldmath$\sigma$\unboldmath}^{j_{D},j_{D+1}}_{D}\circ\mathbf{v}^{j_{D+1}},$
(34)

with $\mathbf{u}^{j_{1}}\in\mathbb{R}^{M}$,
$\mbox{\boldmath$\sigma$\unboldmath}^{j_{i},j_{i+1}}_{i}\in\mathbb{R}^{n_{i}}$,
and $\mathbf{v}^{j_{D+1}}\in\mathbb{R}^{N}$, where the positive integers
$\widetilde{r}_{i}$ are referred to as the _compression ranks_ (or TT-ranks)
of the decomposition. For higher order tensors the TT format is in general
more efficient compared to HOSVD. This may be beneficial for large $D$, the
dimension of parameter space. In [58, 57] a stable algorithm for finding
$\widetilde{\boldsymbol{\Phi}}$ based on truncated SVD for a sequence of
unfolding matrices was introduced and the optimality property similar to (33)
was proved.

Once an optimal TT approximation (34) is computed, we organize the vectors and
matrices from (34) into matrices

$\mathrm{U}=[\mathbf{u}^{1},\dots,\mathbf{u}^{\widetilde{r}_{1}}]\in\mathbb{R}^{M\times\widetilde{r}_{1}},\quad\mathrm{V}=[\mathbf{v}^{1},\dots,\mathbf{v}^{\widetilde{r}_{D+1}}]\in\mathbb{R}^{N\times\widetilde{r}_{D+1}},$
(35)

and third order tensors $\mathbf{S}_{i}\in\mathbb{R}^{\widetilde{r}_{i}\times
n_{i}\times\widetilde{r}_{i+1}}$, defined entry-wise as

$(\mathbf{S}_{i})_{jkq}=(\mbox{\boldmath$\sigma$\unboldmath}^{jq}_{i})_{k},\quad
j=1,\ldots,\widetilde{r}_{i},\quad k=1,\ldots,n_{i},\quad
q=1,\ldots,\widetilde{r}_{i+1},$ (36)

for all $i=1,\ldots,D$. Note that matrix $\mathrm{U}$ is orthogonal and so
_its columns provide an orthogonal basis in the universal space
$\widetilde{V}$_. The dimension of $\widetilde{V}$ is defined by the first TT-
rank, $\mbox{dim}(\widetilde{V})=\widetilde{r}_{1}$.

While $\mathrm{U}$ is an orthogonal matrix, the columns of $\mathrm{V}$ are
orthogonal, but not necessarily orthonormal. Thus, we introduce a diagonal
scaling matrix

$\mathrm{W}_{c}=\mbox{diag}\left(\|\mathbf{v}^{1}\|,\ldots,\|\mathbf{v}^{\widetilde{r}_{D+1}}\|\right)\in\mathbb{R}^{\widetilde{r}_{D+1}\times\widetilde{r}_{D+1}}.$
(37)

The essential information about $\widetilde{\boldsymbol{\Phi}}$ to be
transmitted to the online phase includes $\mathbf{S}_{i}$ tensors and the
scaling factors:

$\mbox{online}(\widetilde{\boldsymbol{\Phi}})=\left\\{\mathbf{S}_{i}\in\mathbb{R}^{\widetilde{r}_{i}\times
n_{i}\times\widetilde{r}_{i+1}},~{}{\small
i=1,\dots,D},~{}\mathrm{W}_{c}\in\mathbb{R}^{\widetilde{r}_{D+1}\times\widetilde{r}_{D+1}}\right\\}.$
(38)

To find the coordinates of the local basis, we define the parameter-specific
_core matrix_
$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{\widetilde{r}_{1}\times\widetilde{r}_{D+1}}$
as the product

$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\prod_{i=1}^{D}\left(\mathbf{S}_{i}\times_{2}\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})\right).$
(39)

Using the definition of $k$-mode product, (18) and (34), one computes

$\begin{split}\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})&=\sum_{j_{1}=1}^{\widetilde{r}_{1}}\dots\sum_{j_{D+1}=1}^{\widetilde{r}_{D+1}}\langle\mbox{\boldmath$\sigma$\unboldmath}^{j_{1},j_{2}}_{1},\mathbf{e}_{1}(\mbox{\boldmath$\alpha$\unboldmath})\rangle\cdot\dotsc\cdot\langle\mbox{\boldmath$\sigma$\unboldmath}^{j_{D},j_{D+1}}_{D},\mathbf{e}_{D}(\mbox{\boldmath$\alpha$\unboldmath})\rangle\,\mathbf{u}^{j_{1}}\circ\mathbf{v}^{j_{D+1}}\\\
&=\sum_{j_{1}=1}^{\widetilde{r}_{1}}\dots\sum_{j_{D+1}=1}^{\widetilde{r}_{D+1}}\left(\mathbf{S}_{1}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{j_{1},j_{2}}\cdot\dotsc\cdot\left(\mathbf{S}_{D}\times_{D+1}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{j_{D},j_{D+1}}\,\mathbf{u}^{j_{1}}\circ\mathbf{v}^{j_{D+1}}\\\
&=\sum_{j_{1}=1}^{\widetilde{r}_{1}}\sum_{j_{D+1}=1}^{\widetilde{r}_{D+1}}\Big{(}\prod_{i=1}^{D}\left(\mathbf{S}_{i}\times_{2}\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})\right)\Big{)}_{j_{1}j_{D+1}}\,\mathbf{u}^{j_{1}}\circ\mathbf{v}^{j_{D+1}}=\mathrm{U}\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{V}^{T}.\end{split}$

Consider the SVD of the rescaled core matrix:

$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{W}_{c}=\mathrm{U}_{c}\Sigma_{c}\mathrm{V}_{c}^{T}.$
(40)

Using this and the above representation of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ we compute

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{U}\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{W}_{c}\mathrm{W}_{c}^{-1}\mathrm{V}^{T}=\left({\mathrm{U}}\mathrm{U}_{c}\right)\Sigma_{c}\left({\mathrm{V}}\mathrm{W}_{c}^{-1}\mathrm{V}_{c}\right)^{T}.$
(41)

The right-hand side of (41) is the thin SVD of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$, since matrices
${\mathrm{U}}$, $\mathrm{U}_{c}$, ${\mathrm{V}}\mathrm{W}_{c}^{-1}$, and
$\mathrm{V}_{c}$ are all orthogonal. We conclude that _the coordinates
$\left\\{\mbox{\boldmath$\beta$\unboldmath}_{1}(\mbox{\boldmath$\alpha$\unboldmath}),\dots,\mbox{\boldmath$\beta$\unboldmath}_{n}(\mbox{\boldmath$\alpha$\unboldmath})\right\\}$
of the local reduced basis in the universal space $\widetilde{V}$ are the
first $n$ columns of $\mathrm{U}_{c}$._ The parameter-specific basis is then
$\\{\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}$, with
$\mathbf{z}_{i}(\mbox{\boldmath$\alpha$\unboldmath})={\mathrm{U}}\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})$
(not actually computed at the online stage).

We summarize the above in the following algorithm.

######  Algorithm 3 (TT-TROM).

  * •

Offline stage.
_Input:_ snapshot tensor $\boldsymbol{\Phi}\in\mathbb{R}^{M\times
n_{1}\times\ldots\times n_{D}\times N}$ and target accuracy
$\widetilde{\varepsilon}$;
_Output:_ Compression ranks, TT decomposition matrices and third order tensors
as in (35)–(36);
_Compute:_ Use algorithm from [58] with prescribed accuracy
$\widetilde{\varepsilon}$ to compute TT decomposition (34).

  * •

Online stage.
_Input:_ $\mbox{\rm online}(\widetilde{\boldsymbol{\Phi}})$ as defined in
(38), reduced space dimension
$n\leq\min\\{\widetilde{r}_{1},\widetilde{r}_{D+1}\\}$, and parameter vector
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$;
_Output:_ Coordinates of the reduced basis in $\widetilde{V}$:
$\\{\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})\\}_{i=1}^{n}\subset\mathbb{R}^{\widetilde{r}_{1}}$;
_Compute:_

    1. 1.

Use tensors $\mathbf{S}_{i}$ to assemble the core matrix
$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{\widetilde{r}_{1}\times\widetilde{r}_{D+1}}$
as in (39);

    2. 2.

Compute the SVD of the scaled core matrix
$\mathrm{C}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{W}_{c}=\mathrm{U}_{c}\Sigma_{c}\mathrm{V}_{c}^{T}$
with
$\mathrm{U}_{c}=[\widetilde{\mathbf{u}}_{1},\ldots,\widetilde{\mathbf{u}}_{\widetilde{r}_{1}}]$;

    3. 3.

Set
$\mbox{\boldmath$\beta$\unboldmath}_{i}(\mbox{\boldmath$\alpha$\unboldmath})=\widetilde{\mathbf{u}}_{i}$,
$i=1,\ldots,n$.

###  3.8 General parameter sampling

Grid-based sampling of parameter space can be computationally expensive or not
applicable if the set of admissible parameters $\mathcal{A}$ is not a box (or
an image of a box) in Euclidean space. However, interpolatory TROMs introduced
above can be extended to accommodate a more general sampling set
${\widehat{\mathcal{A}}}$. If $\mathcal{A}$ does have the Cartesian structure
(a box or an image of a box), then one way to reduce offline computational
costs is to compute the snapshots for only a few parameter values from a
Cartesian grid ${\widehat{\mathcal{A}}}\subset\mathcal{A}$. To recover the
missing entries of the full snapshot tensor $\boldsymbol{\Phi}$, one may use a
low-rank tensor completion method, e.g., one of those studied in [53, 28, 41,
69, 6]. The low-rank completion can be performed for any of the three
compressed tensor formats considered above. We shall investigate this option
elsewhere. In this paper, we consider another (more general) approach.

With a slight abuse of notation, let
${{\widehat{\mathcal{A}}}}=\\{\widehat{\boldsymbol{\alpha}}_{1},\dots,\widehat{\boldsymbol{\alpha}}_{K}\\}\subset\mathcal{A}$
be a set of sampled parameter values. We assume that ${\widehat{\mathcal{A}}}$
is a frame in $\mathbb{R}^{D}$ and so $K\geq D$. Note that $K$ does not obey
(6) for a general sampling. Given an out-of-sample vector of parameters
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, let

$\mathbf{e}\,:\,\mbox{\boldmath$\alpha$\unboldmath}\to\mathbb{R}^{K}$ (42)

be the representation of $\alpha$ in ${\widehat{\mathcal{A}}}$, i.e.,

$\mbox{\boldmath$\alpha$\unboldmath}=\sum_{j=1}^{K}a_{j}\widehat{\boldsymbol{\alpha}}_{j},\quad\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath})=(a_{1},\dots,a_{K})^{T},$
(43)

with an additional constraint enforcing uniqueness of the representation.

Similarly to (16) for the Cartesian grid case, we assume that for a smooth
function ${g}:\mathcal{A}\to\mathbb{R}$ it holds

${g}(\mbox{\boldmath$\alpha$\unboldmath})\approx\sum_{j=1}^{K}a_{j}{g}(\widehat{\boldsymbol{\alpha}}_{j}).$
(44)

In Section 5.1 we describe one particular choice of
$\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath})$ that is used in all
numerical experiments reported in Section 5.

To assemble the snapshot tensor, for each
$\widehat{\boldsymbol{\alpha}}_{j}\in{\widehat{\mathcal{A}}}$, $j=1,\ldots,K$,
collect the snapshot vectors
$\mathbf{u}(t_{k},\widehat{\boldsymbol{\alpha}}_{j})=\left(u_{1}(t_{k},\widehat{\boldsymbol{\alpha}}_{j}),\ldots,u_{M}(t_{k},\widehat{\boldsymbol{\alpha}}_{j})\right)^{T}$,
$k=1,\dots,N$, and arrange them in a third order tensor
$\boldsymbol{\Phi}\in\mathbb{R}^{M\times K\times N}$ with entries

$(\boldsymbol{\Phi})_{ijk}=u_{i}(t_{k},\widehat{\boldsymbol{\alpha}}_{j}),\quad
i=1,\ldots,M,\quad j=1,\ldots,K,\quad k=1,\ldots,N.$ (45)

Then, for any $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, the
parameter-specific reduced basis is defined as the first $n$ left singular
vectors of

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\widetilde{\boldsymbol{\Phi}}\times_{2}\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath}),$
(46)

where $\widetilde{\boldsymbol{\Phi}}$ is a low rank approximation of the
snapshot tensor $\boldsymbol{\Phi}$ with entries (45). The same three
compressed tensor formats considered above (CP, HOSVD and TT) can be used for
$\widetilde{\boldsymbol{\Phi}}$, with TT format being inferior to HOSVD (for
3D tensors both HOSVD- and TT-decomposition are Tucker tensors). An orthogonal
basis in the universal space of $\widetilde{\boldsymbol{\Phi}}$ and
coordinates of a local parameter-specific basis in it are computed similarly
to the Cartesian grid sampling cases considered previously in Sections
3.5–3.7. The only difference is that all the calculations therein are
performed setting $D=1$ and replacing
$\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})$ with
$\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath})$. This includes the
optimality result (33), where the factor $\sqrt{D+2}$ becomes $\sqrt{3}$.

###  3.9 Complexity and compression analysis

Projection-based parametric ROM framework consists in general of the following
steps.

  * (i)

High-fidelity simulations of (1) to generate the snapshot tensor
$\boldsymbol{\Phi}$;

  * (ii)

Offline stage: computing the compressed approximation
$\widetilde{\boldsymbol{\Phi}}$ to $\boldsymbol{\Phi}$ in one of low-rank
tensor formats;

  * (iii)

Passing the $\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ part of the
compressed tensor to the online stage;

  * (iv)

Online stage: using $\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ to compute
the coordinates of the parameter-specific reduced basis for an input $\alpha$;

  * (v)

Solving (1) projected onto the reduced space.

Since steps (i) and (v) are common for all projection-based ROM approaches, we
focus below on the computational and storage/transmission costs invoked in
steps (ii)–(iv). The necessary details on step (v) are included in Section
3.10.

First, we discuss briefly the computational costs at the more expensive
offline stage. For CP-TROM, the standard algorithm for finding
$\widetilde{\boldsymbol{\Phi}}$ in CP format (20) is the ALS method [36] which
for a _given_ CP rank $R$ iteratively fits a rank $R$ tensor
$\widetilde{\Phi}$ by solving on each iteration $D+2$ least squares problems
for the factors $\mathbf{u}^{r}$,
$\mbox{\boldmath$\sigma$\unboldmath}^{r}_{i}$, $i=1,\dots,D$, and
$\mathbf{v}^{r}$, $r=1,\ldots,R$. While straightforward to implement, the
method is sensitive to the choice of initial guess and may converge slowly. We
refer the reader to [46] for a guidance on the literature on improving the
efficiency of ALS and possible alternatives. On the other hand, computing
$\widetilde{\boldsymbol{\Phi}}$ in either HOSVD or TT formats relies on
finding truncated SVDs for matrix unfoldings of $\boldsymbol{\Phi}$ [21, 58].
Therefore, the computational complexity and cost of step (ii) for HOSVD- and
TT-TROM is essentially the same as that of standard POD-ROM.

Second, to measure the amount of information transmitted to the online stage
at step (iii), we introduce the _compression factor_ CF, defined as

$\text{CF}=\frac{\\#(\boldsymbol{\Phi})}{\\#(\mbox{online}(\widetilde{\boldsymbol{\Phi}}))}\,,$
(47)

where we denote by $\\#(\boldsymbol{\Psi})$ the number of floating point
numbers needed to store a tensor $\boldsymbol{\Psi}$. Specifically,
$\\#(\boldsymbol{\Phi})=MKN$ is simply the total number of entries in
$\boldsymbol{\Phi}$, while $\\#(\mbox{online}(\widetilde{\boldsymbol{\Phi}}))$
is the number of entries needed to store all the factors passed to the online
stage, as defined in (24), (29) and (38) for CP-, HOSVD- and TT- TROMs, which
we summarize in Table 1.

Table 1: Number of entries needed to store $\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ for CP, HOSVD and TT formats. | $\\#(\mbox{online}(\widetilde{\boldsymbol{\Phi}}))$
---|---
Format | Cartesian grid-based | General
CP | $R\big{(}\sum\limits_{i=1}^{D}n_{i}+R+1\big{)}$ | $R\big{(}K+R+1\big{)}$
HOSVD | $\widetilde{N}\widetilde{M}\prod\limits_{i=1}^{D}\widetilde{n}_{i}+\sum\limits_{i=1}^{D}\widetilde{n}_{i}n_{i}$ | $\widetilde{N}\widetilde{M}\widetilde{n}_{1}+\widetilde{n}_{1}K$
TT | $\widetilde{r}_{D+1}+\sum\limits_{i=1}^{D}\widetilde{r}_{i}n_{i}\widetilde{r}_{i+1}$ | $\widetilde{r}_{2}+\widetilde{r}_{1}K\widetilde{r}_{2}$

Table 1 shows that the compression factor is largely determined by the
compression ranks. In turn, the ranks depend on $\widetilde{\varepsilon}$ and
variability of observed states.

Third, the computational complexity of finding $\alpha$-specific reduced basis
in step (iv) is determined by the interpolation procedure and the computation
of first $n$ left singular vectors of the core matrix. Since vectors
$\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})$ contain very few non-
zero entries, e.g., $p=2$ or $3$ of them for the Cartesian sampling, the
number of operations for computing core matrices
$\mathrm{C}(\mbox{\boldmath$\alpha$\unboldmath})$ for CP-, HOSVD- and TT-TROM
is

$O\left(R^{2}\right),~{}~{}O\Big{(}\widetilde{M}\widetilde{N}\prod\limits_{i=1}^{D}\widetilde{n}_{i}\Big{)},~{}~{}\text{and}~{}~{}O\Big{(}\sum\limits_{i=2}^{D}\widetilde{r}_{i-1}\widetilde{r}_{i}\widetilde{r}_{i+1}\Big{)},$
(48)

respectively. CP-, HOSVD- and TT-TROM algorithms proceed to compute the SVD of
small core matrices of sizes $R\times R$, $\widetilde{M}\times\widetilde{N}$
and $\widetilde{r}_{1}\times\widetilde{r}_{D+1}$, respectively. If a reduced
basis in the physical space is desired, then one finds its vectors as linear
combinations of columns of $\mathrm{U}$, which requires $O(MRn)$,
$O(M\widetilde{M}n)$ or $O(M\widetilde{r}_{1}n)$ operations for CP-, HOSVD- or
TT-TROM, respectively. Section 3.10 below discusses how these costs can be
avoided at the online phase. We note that for a fixed compression accuracy
$\widetilde{\epsilon}$, it is often observed in practice that the
corresponding ranks of HOSVD and TT formats satisfy
$\widetilde{M}\simeq\widetilde{r}_{1}$,
$\widetilde{N}\simeq\widetilde{r}_{D+1}$.

In summary, the computational costs of the offline stage for TROMs are
comparable to those of POD-ROM for a multi-parameter problem. At the online
stage complexity of all preparatory steps depends only on compressed tensor
ranks rather than the size of the snapshot tensor $\boldsymbol{\Phi}$. The
amount of information transmitted from offline to online stages is determined
by the compressed tensor ranks, as should be clear from Table 1.

###  3.10 TROM evaluation

Besides finding a suitable reduced basis, a fast evaluation of the reduced
model for any incoming $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ is
required for a reduced modeling scheme to be effective. Efficient
implementation of a projected parametric model is a well-known challenge that
have been addressed in the literature with various approaches; see, e.g., [12,
30, 17, 24, 15, 10, 39, 47]. The tensorial approach presented here does not
directly contribute to resolving this issue, but it does not make it harder
either and so techniques known from the literature can be adapted in the TROM
framework.

For example, assume that $F(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath})$
from (1) has an affine dependence on parameters and linear dependence on
$\mathbf{u}$:

$F(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath})=\sum_{i=1}^{P}f_{i}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{A}_{i}\mathbf{u},$

with some $f_{i}:\mathcal{A}\to\mathbb{R}$ and parameter-independent
$\mathrm{A}_{i}\in\mathbb{R}^{M\times M}$. We assume that $P$ is not too
large, at least independent of other dimensions. Then the _offline_ stage of
model reduction consists of projecting matrices onto the universal space by
computing $\widehat{\mathrm{A}}_{i}=\mathrm{U}^{T}\mathrm{A}_{i}\mathrm{U}$,
where $\mathrm{U}$ is an orthogonal basis matrix for $\widetilde{V}$ provided
by the tensor decompositions. The new matrices $\widehat{\mathrm{A}}_{i}$ have
the reduced size $T_{r}\times T_{r}$, with
$T_{r}\in\\{R,\widetilde{M},\widetilde{r}_{1}\\}$ for CP-, HOSVD- and TT-
TROMs, respectively.

For each of TROMs, denote by $\mathrm{U}_{c}(n)$ the matrix of the first $n$
columns of $\mathrm{U}_{c}$, left singular vectors of $\alpha$-specific core
matrices. During the _online_ stage one solves the system projected further on
the parameter-specific local basis:

$\mathbf{v}_{t}=\sum_{i=1}^{P}f_{i}(\mbox{\boldmath$\alpha$\unboldmath}){\mathrm{U}^{T}_{c}(n)}\widehat{\mathrm{A}}_{i}{\mathrm{U}_{c}(n)}\mathbf{v},$

where $\mathbf{v}(t)$ is the trajectory in a space spanned by the columns of
${\mathrm{U}_{c}(n)}\in\mathbb{R}^{T_{r}\times n}$, i.e., the corresponding
physical states are given by
$\mathbf{u}(t)=\mathrm{U}{\mathrm{U}_{c}(n)}\mathbf{v}(t)$. We see that online
computations depend only on reduced dimensions (tensor ranks) and the small
dimension $n$ of parameter-specific basis. This observation can be extended to
the case when $F$ has a low order polynomial non-linearity with respect to
$\mathbf{u}$. For example, quadratic nonlinear terms, as in Burgers or Navier-
Stokes equations, can be evaluated in $O(T_{r}^{2})$ operations on each time
step given a vector $\mathbf{v}$ in the local reduced basis.

To evaluate more general nonlinear terms, one can use a hyper-reduction
technique such as the discrete empirical interpolation method (DEIM) [17]. In
this approach, the nonlinear term is approximated in a basis of its snapshots.
As an example, consider
$F(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{A}\mathbf{u}+f(t,\mathbf{u}(t),\mbox{\boldmath$\alpha$\unboldmath})$,
with $\mathrm{A}\mathbf{u}$ representing linear and
$f(t,\mathbf{u}(t),\mbox{\boldmath$\alpha$\unboldmath})$ representing the non-
linear part of $F$. Define the snapshots
$\mathbf{f}_{i}=f(t_{i},\mathbf{u}(t_{i}),\mbox{\boldmath$\alpha$\unboldmath}_{i})$,
$i=1,\dots,N_{\rm DEIM}$ for some “greedy” choice of parameters and time
instances and high-fidelity solution $\mathbf{u}$. Denote by $\mathrm{Q}$ an
orthogonal basis matrix for
$\mbox{span}\\{\mathbf{f}_{1},\dots,\mathbf{f}_{N_{\rm DEIM}}\\}$, then DEIM
approximates

$f(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath})\approx\mathrm{Q}(\mathrm{P}\mathrm{Q})^{-1}\mathrm{P}f(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath}),$

where $\mathrm{P}^{T}$ is an $N_{\rm DEIM}\times M$ “selection” matrix such
that for any $f\in\mathbb{R}^{M}$, $\mathrm{P}f$ contains $N_{\rm DEIM}$
selected entries of $f$. A particular $\mathrm{P}$ corresponds to the choice
of spatial interpolation nodes, cf. [17]. In TROM one may pre-compute
$\widehat{\mathrm{Q}}=\mathrm{U}^{T}\mathrm{Q}$ and
$\widehat{\mathrm{A}}=\mathrm{U}^{T}\mathrm{A}\mathrm{U}$ during the offline
stage, then solve at the online stage

$\mathbf{v}_{t}={\mathrm{U}_{c}^{T}(n)}\widehat{\mathrm{A}}\mathrm{U}_{c}(n)\mathbf{v}+{\mathrm{U}^{T}_{c}(n)}\widehat{\mathrm{Q}}(\mathrm{P}\mathrm{Q})^{-1}\mathrm{P}f(t,\mathrm{U}{\mathrm{U}_{c}(n)}\mathbf{v},\mbox{\boldmath$\alpha$\unboldmath}),$

with costs depending on compressed tensor ranks, $n$, and the dimension of
DEIM space, but not on the dimensions of high-fidelity simulations. It is an
interesting question, whether the tensor technique can be applied to make the
DEIM space parameter-specific for more efficient reduce online computations.
We plan to address this question elsewhere.

##  4 Prediction analysis

In this section we assess the prediction power of the reduced basis
$\mathcal{Z}_{n}(\mbox{\boldmath$\alpha$\unboldmath})=\\{\mathbf{z}_{1},\ldots,\mathbf{z}_{n}\\}$
consisting of the first $n$ left singular vectors of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ from (18), for a
parameter
$\mbox{\boldmath$\alpha$\unboldmath}=(\alpha_{1},\ldots,\alpha_{D})^{T}\in\mathcal{A}$,
not necessarily from a sampling set; i.e.,
$[\mathbf{z}_{1},\ldots,\mathbf{z}_{n}]=\mathrm{U}{\mathrm{U}_{c}(n)}$.

For the discussion below we also need the following notation. Given an
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, we denote by
$\boldsymbol{\psi}_{i}=\mathbf{u}(t_{i},\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M}$,
$i=1,\dots,N$, the snapshots of a high-fidelity solution to (1) and let
$\Psi(\mbox{\boldmath$\alpha$\unboldmath})=[\boldsymbol{\psi}_{1},\dots,\boldsymbol{\psi}_{N}]\in\mathbb{R}^{M\times
N}$ be the corresponding snapshot matrix. Note that in practice the snapshots
for out-of-sample parameters are not available, so the matrix
$\Psi(\mbox{\boldmath$\alpha$\unboldmath})$ should be treated as unknown.

We estimate the prediction power of
$\mathcal{Z}_{n}(\mbox{\boldmath$\alpha$\unboldmath})$ in terms of the
quantity

$E_{n}(\mbox{\boldmath$\alpha$\unboldmath})=\frac{1}{NM}\sum_{i=1}^{N}\left\|\boldsymbol{\psi}_{i}-\sum_{j=1}^{n}\langle\boldsymbol{\psi}_{i},\mathbf{z}_{j}\rangle\mathbf{z}_{j}\right\|^{2}_{\ell^{2}},$
(49)

which measures how accurate the solution
$\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath})$ at time instances $t_{i}$
can be represented in the reduced basis for the arbitrary but fixed
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$. The scaling $1/(NM)$
accounts for the variation of dimensions $N$ and $M$, which may correspond to
the number of temporal and spatial degrees of freedom, respectively, if (1)
comes from a discretization of a parabolic PDE defined in a spatial domain
$\Omega$. In this case and for uniform grids, the quantity in (49) is
consistent with the $L^{2}(0,T,L^{2}(\Omega))$ norm.

Below we prove an estimate for $E_{n}(\mbox{\boldmath$\alpha$\unboldmath})$ in
terms of $\widetilde{\varepsilon}$ from (8), the singular values of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$ and interpolation
properties of $\mathbf{e}^{i}(\mbox{\boldmath$\alpha$\unboldmath})$,
$i=1,\ldots,D$. To make use of the latter, we introduce the following
quantities related to the interpolation procedure. For Cartesian grid-based
sampling we define the maximum grid step

$\delta_{i}=\max\limits_{1\leq j\leq
n_{i}-1}\left|\widehat{\alpha}_{j}^{i}-\widehat{\alpha}_{j+1}^{i}\right|,\quad
i=1,\ldots,D.$ (50)

Relation (17) implies that the interpolation procedure (15)–(17) is of order
$p$, i.e., for any sufficiently smooth
$f:[\alpha_{i}^{\min},\alpha_{i}^{\max}]\to\mathbb{R}$ it holds

$\sup_{a\in[\alpha_{i}^{\min},\alpha_{i}^{\max}]}\Big{|}f(a)-\sum_{j=1}^{n_{i}}e^{i}_{j}\left(a\mathbf{e}_{i}\right)f(\widehat{\alpha}_{i}^{j})\Big{|}\leq
C_{a}\|f^{(p)}\|_{C([\alpha_{i}^{\min},\alpha_{i}^{\max}])}\delta_{i}^{p},$
(51)

for $i=1,\dots,D$, where $\mathbf{e}_{i}\in\mathbb{R}^{n_{i}}$ is the
$i^{\mbox{\scriptsize th}}$ column of an $n_{i}\times n_{i}$ identity matrix.
The constant $C_{a}$ does not depend on $f$. We let
$\delta^{p}=\sum_{i=1}^{D}\delta_{i}^{p}$ and also assume that the
interpolation procedure is stable in the sense that

$\left(\sum_{j=1}^{n_{i}}\left|e^{i}_{j}(a\mathbf{e}_{i})\right|^{2}\right)^{\frac{1}{2}}\leq
C_{e},$ (52)

with some $C_{e}$ independent of $a\in[\alpha_{i}^{\min},\alpha_{i}^{\max}]$
and $i=1,\ldots,D$. For the example of linear interpolation with $p=2$ and
$\alpha_{i}^{\min}$, $\alpha_{i}^{\max}$ included among the grid nodes, bounds
(51) and (52) hold with $C_{a}=\frac{1}{8}$ and $C_{e}=1$.

To estimate $E_{n}(\mbox{\boldmath$\alpha$\unboldmath})$, consider the SVD of
$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M\times
N}$ given by

$\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\widetilde{\mathrm{U}}\widetilde{\Sigma}\widetilde{\mathrm{V}}^{T},~{}~{}\text{with}~{}~{}\widetilde{\Sigma}=\text{diag}(\widetilde{\sigma}_{1},\dots,\widetilde{\sigma}_{N}).$
(53)

Then $\mathrm{Z}=[\mathbf{z}_{1},\ldots,\mathbf{z}_{n}]\in\mathbb{R}^{M\times
n}$ is build as the first $n$ columns of $\widetilde{\mathrm{U}}$, the reduced
basis vectors, i.e., $\mathrm{Z}=\mathrm{U}\mathrm{U}_{c}$. Then,

$\displaystyle E_{n}(\mbox{\boldmath$\alpha$\unboldmath})$
$\displaystyle=\frac{1}{NM}\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})\Psi(\mbox{\boldmath$\alpha$\unboldmath})\right\|^{2}_{F}$
$\displaystyle\leq\frac{1}{NM}\left(\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})(\Psi(\mbox{\boldmath$\alpha$\unboldmath})-\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath}))\right\|_{F}+\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}\right)^{2}$
$\displaystyle\leq\frac{1}{NM}\left(\left\|{\Psi}(\mbox{\boldmath$\alpha$\unboldmath})-\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}+\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}\right)^{2},$
(54)

where we used triangle inequality and
$\|\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T}\|\leq 1$ for the spectral norm of the
projector. For the last term in (54), we observe

$\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}=\left\|\widetilde{\mathrm{U}}\;\text{diag}(0,\dots,0,\widetilde{\sigma}_{n+1},\dots,\widetilde{\sigma}_{N})\;\widetilde{\mathrm{V}}^{T}\right\|_{F}=\left(\sum_{j=n+1}^{N}\widetilde{\sigma}_{j}^{2}\right)^{\frac{1}{2}}.$
(55)

To handle the first term of (54), consider the extraction

$\Phi_{e}(\mbox{\boldmath$\alpha$\unboldmath})=\boldsymbol{\Phi}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\times_{3}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\dots\times_{D+1}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})$
(56)

and proceed using the triangle inequality

$\big{\|}\Psi(\mbox{\boldmath$\alpha$\unboldmath})-\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\big{\|}_{F}\leq\big{\|}\Psi(\mbox{\boldmath$\alpha$\unboldmath})-\Phi_{e}(\mbox{\boldmath$\alpha$\unboldmath})\big{\|}_{F}+\big{\|}\Phi_{e}(\mbox{\boldmath$\alpha$\unboldmath})-\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\big{\|}_{F}.$
(57)

We use the stability of interpolation (52) and (8) to bound the second term of
(57). Specifically,

$\begin{split}\left\|\Phi_{e}(\mbox{\boldmath$\alpha$\unboldmath})-\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}&=\left\|(\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}})\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\times_{3}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\dots\times_{D+1}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}\\\
&\leq\left\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\right\|_{F}\|\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\|_{\ell^{2}}\|\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\|_{\ell^{2}}\dots\|\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\|_{\ell^{2}}\\\
&\leq(C_{e})^{D}\left\|\boldsymbol{\Phi}-\widetilde{\boldsymbol{\Phi}}\right\|_{F}\leq(C_{e})^{D}\widetilde{\varepsilon}\left\|\boldsymbol{\Phi}\right\|_{F}.\end{split}$
(58)

It remains to handle the first term in (57). At this point we need more
precise assumptions on the smoothness of
$\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath})$, the solution of (1). In
particular,

$\mathbf{u}\in
C([0,T]\times\overline{\mathcal{A}})^{M},\quad\frac{\partial^{\mathbf{j}}\mathbf{u}}{\partial\alpha^{j_{1}}_{1}\dots\alpha^{j_{D}}_{D}}\in
C([0,T]\times\overline{\mathcal{A}})^{M},\quad|\mathbf{j}|\leq p.$ (59)

We note that (59) is guaranteed to hold if the unique solution to (1) exists
on $(0,T_{1})$ for all
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}_{1}$, with $T<T_{1}$ and
$\overline{\mathcal{A}}\subset\mathcal{A}_{1}$ and $F$ is continuous with
continuous partial derivatives in components of $\mathbf{u}$ and $\alpha$ of
order up to $p$ [37].

Using interpolation property (51), we compute

$\begin{split}\left(\boldsymbol{\Phi}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{:,i_{2},\dots,i_{D},k}&=\sum_{j=1}^{n_{1}}e^{1}_{j}(\mbox{\boldmath$\alpha$\unboldmath})\mathbf{u}(t_{k},\widehat{\alpha}_{1}^{j},\widehat{\alpha}_{2}^{i_{2}},\dots,\widehat{\alpha}_{D}^{i_{D}})\\\
&=\mathbf{u}(t_{k},\alpha_{1},\widehat{\alpha}_{2}^{i_{2}}\dots,\widehat{\alpha}_{D}^{i_{D}})+\Delta^{1}_{:,i_{2},\dots,i_{D},k},\end{split}$
(60)

with the remainder term obeying a component-wise bound

$|\Delta^{1}_{:,i_{2},\dots,i_{D},k}|\leq
C_{a}\sup_{a\in[\alpha_{1}^{\min},\alpha_{1}^{\max}]}\left|\frac{\partial^{p}\mathbf{u}}{\partial\alpha^{p}_{1}}(t_{k},a,\widehat{\alpha}_{2}^{i_{2}},\dots,\widehat{\alpha}_{D}^{i_{D}})\right|\delta_{1}^{p},$
(61)

where the absolute value of vectors is understood entry-wise. Analogously, we
compute

$\begin{split}\big{(}\boldsymbol{\Phi}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})&\times_{3}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\big{)}_{:,i_{3},\dots,i_{D},k}=\left((\boldsymbol{\Phi}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath}))\times_{2}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{:,i_{3},\dots,i_{D},k}\\\
&=\sum_{j=1}^{n_{2}}e^{2}_{j}(\mbox{\boldmath$\alpha$\unboldmath})\left(\mathbf{u}(t_{k},\alpha_{1},\widehat{\alpha}_{2}^{j},\widehat{\alpha}_{3}^{i_{3}},\dots,\widehat{\alpha}_{D}^{i_{D}})+\Delta^{1}_{:,j,i_{3},\dots,i_{D},k}\right)\\\
&=\mathbf{u}(t_{k},\alpha_{1},\alpha_{2},\widehat{\alpha}_{3}^{i_{3}},\dots,\widehat{\alpha}_{D}^{i_{D}})+\Delta^{2}_{:,i_{3},\dots,i_{D},k}+\sum_{j=1}^{n_{2}}e^{2}_{j}(\mbox{\boldmath$\alpha$\unboldmath})\Delta^{1}_{:,j,i_{3},\dots,i_{D},k},\end{split}$

with a component-wise bound for the remainder

$\begin{split}\Big{|}\Delta^{2}_{:,i_{3},\dots,i_{D},k}+\sum_{j=1}^{n_{2}}e^{2}_{j}(\mbox{\boldmath$\alpha$\unboldmath})&\Delta^{1}_{:,j,i_{3},\dots,i_{D},k}\Big{|}\\\
\leq
C_{a}&\sup_{a\in[\alpha_{2}^{\min},\alpha_{2}^{\max}]}\left|\frac{\partial^{p}\mathbf{u}}{\partial\alpha^{p}_{2}}(t_{k},\alpha_{1},a,\widehat{\alpha}_{3}^{i_{3}},\dots,\widehat{\alpha}_{D}^{i_{D}})\right|\delta_{2}^{p}\\\
+\;C_{e}\;C_{a}&\sup_{a\in[\alpha_{1}^{\min},\alpha_{1}^{\max}]}\left|\frac{\partial^{p}\mathbf{u}}{\partial\alpha^{p}_{1}}(t_{k},a,\widehat{\alpha}_{2}^{i_{2}},\dots,\widehat{\alpha}_{D}^{i_{D}})\right|\delta_{1}^{p}.\end{split}$
(62)

Applying the same argument repeatedly, we obtain

$\begin{split}\left(\Phi_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{:,k}&=\left(\boldsymbol{\Phi}\times_{2}\mathbf{e}^{1}(\mbox{\boldmath$\alpha$\unboldmath})\times_{3}\mathbf{e}^{2}(\mbox{\boldmath$\alpha$\unboldmath})\dots\times_{D+1}\mathbf{e}^{D}(\mbox{\boldmath$\alpha$\unboldmath})\right)_{:,k}\\\
&=\mathbf{u}(t_{k},\alpha_{1},\alpha_{2},\dots,\alpha_{D})+\Delta_{:,k}=\left(\Psi(\mbox{\boldmath$\alpha$\unboldmath})\right)_{:,k}+\Delta_{:,k},\end{split}$
(63)

with a component-wise bound for the remainder

$\displaystyle\left|\Delta_{:,k}\right|\leq C_{a}\Big{(}$
$\displaystyle\sup_{a\in[\alpha_{D}^{\min},\alpha_{D}^{\max}]}\left|\frac{\partial^{p}\mathbf{u}}{\partial\alpha^{p}_{D}}(t_{k},\alpha_{1},\dots,\alpha_{D-1},a)\right|\delta_{D}^{p}+\dots$
$\displaystyle+\,(C_{e})^{D-2}$
$\displaystyle\sup_{a\in[\alpha_{2}^{\min},\alpha_{2}^{\max}]}\left|\frac{\partial^{p}\mathbf{u}}{\partial\alpha^{p}_{2}}(t_{k},\alpha_{1},a,\widehat{\alpha}_{3}^{i_{3}},\dots,\widehat{\alpha}_{D}^{i_{D}})\right|\delta_{2}^{p}$
$\displaystyle+\,(C_{e})^{(D-1)}$
$\displaystyle\sup_{a\in[\alpha_{1}^{\min},\alpha_{1}^{\max}]}\left|\frac{\partial^{p}\mathbf{u}}{\partial\alpha^{p}_{1}}(t_{k},a,\widehat{\alpha}_{2}^{i_{2}},\dots,\widehat{\alpha}_{D}^{i_{D}})\right|\delta_{1}^{p}\Big{)}.$

Using the definition of the Frobenius norm, we arrive at

$\left\|\Psi(\mbox{\boldmath$\alpha$\unboldmath})-\Phi_{e}(\mbox{\boldmath$\alpha$\unboldmath})\right\|_{F}\leq\sqrt{NM}\;C_{a}C_{\mathbf{u}}\max\left\\{(C_{e})^{(D-1)},1\right\\}\;\delta^{p},$
(64)

where $C_{\mathbf{u}}$ depends only on the smoothness of $\mathbf{u}$ with
respect to the variations of parameters $\alpha$. More precisely, we can take
$C_{\mathbf{u}}=\|\mathbf{u}\|_{C(0,T;C^{p}(\mathcal{A}))}$, which is bounded
due to assumption (59).

Summarizing (54)–(64), we proved the following result.

######  Theorem 1.

Assume the solution $\mathbf{u}$ to (1) satisfies (59),
${\widehat{\mathcal{A}}}$ is a Cartesian grid in parameter domain
$\mathcal{A}$. Then for any
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ the interpolatory TROM
reduced basis
$\mathcal{Z}(\mbox{\boldmath$\alpha$\unboldmath})=\\{\mathbf{z}_{1},\ldots,\mathbf{z}_{n}\\}$
delivers the following representation estimate

$\frac{1}{3NM}\sum_{i=1}^{N}\left\|\mathbf{u}(t_{i},\mbox{\boldmath$\alpha$\unboldmath})-\sum_{j=1}^{n}\left\langle\mathbf{u}(t_{i},\mbox{\boldmath$\alpha$\unboldmath}),\mathbf{z}_{j}\right\rangle\mathbf{z}_{j}\right\|^{2}_{\ell^{2}}\\\
\leq\frac{1}{NM}\left((C_{e})^{2D}\widetilde{\varepsilon}^{2}\left\|\boldsymbol{\Phi}\right\|_{F}+\sum_{i=n+1}^{N}\widetilde{\sigma}_{i}^{2}\right)+C_{a}C_{\mathbf{u}}\max\left\\{(C_{e})^{2(D-1)},1\right\\}\delta^{2p},$
(65)

with $C_{\mathbf{u}}=\|\mathbf{u}\|_{C(0,T;C^{p}(\mathcal{A}))}$ independent
of $\alpha$, sampling grid and $n$.

We summarize here the definitions of quantities that appear in Theorem 1: $N$
is a number of time steps for snapshot collection, while $M$ is the spatial
dimension of snapshots, i.e.,
$\mathbf{u}(t_{i},\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M}$,
$i=1,\ldots,N$; $\widetilde{\varepsilon}$ is the relative accuracy of the
snapshot tensor compression from (8); $\widetilde{\sigma}_{i}$ are the
singular values of $\widetilde{\Phi}_{e}(\mbox{\boldmath$\alpha$\unboldmath})$
from (18) (note that in TROMs we have an access to $\widetilde{\sigma}_{i}$ as
the singular values of core matrices (25), (31) and (40) for CP-, HOSVD-, and
TT-TROM, respectively);
$\delta=\big{(}\sum_{i=1}^{D}\delta_{i}^{p}\big{)}^{\frac{1}{p}}$ is the grid
step parameter of the Cartesian grid in $\mathcal{A}$; $p$ is both the number
of nearest grid points and the order of interpolation of the interpolation
procedure (15)–(17); $C_{e}$ is the interpolation stability constant from
(52); and $D$ is the dimension of parameter space.

For the general parameter sampling, prediction power analysis follows the same
lines as above, simply setting $D=1$. However, the order of interpolation $p$
is slightly more difficult to formalize, so instead of (51)–(52) we rather
assume

$\displaystyle\sup_{\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}}\left|f(\mbox{\boldmath$\alpha$\unboldmath})-\sum_{j=1}^{K}(\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath}))_{j}f(\widehat{\boldsymbol{\alpha}}_{j})\right|\leq\|f\|_{C^{p}(\mathcal{A})}\delta,\quad\Big{(}\sum_{j=1}^{K}|(\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath}))_{j}|^{2}\Big{)}^{\frac{1}{2}}\leq
C_{e}$ (66)

with some $\delta$ depending on ${\widehat{\mathcal{A}}}$. The prediction
estimate then becomes

$\frac{1}{3NM}\sum_{i=1}^{N}\left\|\mathbf{u}(t_{i},\mbox{\boldmath$\alpha$\unboldmath})-\sum_{j=1}^{n}\left\langle\mathbf{u}(t_{i},\mbox{\boldmath$\alpha$\unboldmath}),\mathbf{z}_{j}\right\rangle\mathbf{z}_{j}\right\|^{2}_{\ell^{2}}\\\
\leq\frac{1}{NM}\left((C_{e})^{2}\widetilde{\varepsilon}^{2}\left\|\boldsymbol{\Phi}\right\|_{F}+\sum_{i=n+1}^{N}\widetilde{\sigma}_{i}^{2}\right)+C_{\mathbf{u}}\max\\{(C_{e})^{2},1\\}\delta^{2},$
(67)

with $C_{\mathbf{u}}=\|\mathbf{u}\|_{C(0,T;C^{p}(\mathcal{A}))}$.

We finally, note that the feasibility of a sufficiently accurate lower rank
representation of $\boldsymbol{\Phi}$ depends on the smoothness of
$\mathbf{u}$ as a function of $\mathbf{x},$ $t$ and $\alpha$. This question
can be addressed by considering tensor decompositions of multivariate
functions, e.g. [34, 55]. For these functional CP, HOSVD and hierarchical
Tucker (including TT) formats, the dependence of compression ranks on
$\widetilde{\varepsilon}$ from (8) and the regularity (smoothness) of
$\mathbf{u}$ was studied in [32, 67, 62, 35, 66]. This compression property
for multivariate functions was exploited to effectively represent solutions of
parametric elliptic PDEs using tensor formats in [44, 20, 4, 27, 23] among
other publications.

##  5 Numerical experiments

We perform several numerical experiments to assess the performance of the
three TROM approches and compare them to the conventional POD-ROM. The testing
in Section 5.2 is performed for a dynamical system originating from a
discretization of linear parameter-dependent heat equation. In Section 5.3 a
similar set of tests is carried out for a time-dependent parameterized
advection-diffusion system.

###  5.1 General parameter sampling interpolation scheme

For the numerical examples in the general parameter sampling setting we employ
the following interpolation scheme. Fix an integer $q\geq D+1$ and let
$\mbox{\boldmath$\alpha$\unboldmath}=(\alpha_{1},\ldots,\alpha_{D})\in\mathcal{A}$
be an out-of-sample parameter vector. The interpolation scheme is based on the
weighted minimum norm fit over $q$ nearest neighbors of $\alpha$ in the
sampling set. Thus, we denote by
$\widehat{\boldsymbol{\alpha}}_{i_{1}},\ldots,\widehat{\boldsymbol{\alpha}}_{i_{q}}$
the $q$ closest parameter samples in ${\widehat{\mathcal{A}}}$ to $\alpha$ and
set
$d_{k}=\|\widehat{\boldsymbol{\alpha}}_{i_{k}}-\mbox{\boldmath$\alpha$\unboldmath}\|>0$,
$k=1,\ldots,q$. Next, define the weighting matrix

${\mathrm{D}}=\mbox{diag}(d_{1}^{-1},\ldots,d_{q}^{-1})\in\mathbb{R}^{q\times
q}.$

Also, assemble the matrix

${\mathrm{X}}=\begin{bmatrix}(\widehat{\boldsymbol{\alpha}}_{i_{1}})_{1}&(\widehat{\boldsymbol{\alpha}}_{i_{2}})_{1}&\cdots&(\widehat{\boldsymbol{\alpha}}_{i_{q}})_{1}\\\
(\widehat{\boldsymbol{\alpha}}_{i_{1}})_{2}&(\widehat{\boldsymbol{\alpha}}_{i_{2}})_{2}&\cdots&(\widehat{\boldsymbol{\alpha}}_{i_{q}})_{2}\\\
\vdots&\vdots&\vdots&\vdots\\\
(\widehat{\boldsymbol{\alpha}}_{i_{1}})_{D}&(\widehat{\boldsymbol{\alpha}}_{i_{2}})_{D}&\cdots&(\widehat{\boldsymbol{\alpha}}_{i_{q}})_{D}\\\
1&1&\cdots&1\end{bmatrix}\in\mathbb{R}^{D+1\times q}.$

Solve the weighted minimum norm fitting problem to obtain

$\widehat{\mathbf{a}}={\mathrm{D}}({\mathrm{X}}{\mathrm{D}})^{\dagger}\begin{bmatrix}\mbox{\boldmath$\alpha$\unboldmath}\\\
1\end{bmatrix}\in\mathbb{R}^{q}.$ (68)

Note that the last row of ${\mathrm{X}}$ and the last entry of
$(\mbox{\boldmath$\alpha$\unboldmath},1)^{T}$ enforces the condition that the
entries of $\widehat{\mathbf{a}}$ sum to one. Meanwhile, the presence of the
weighting matrix puts more emphasis on the neighbors of $\alpha$ that are
closest to it.

Once $\widehat{\mathbf{a}}$ is obtained from (68), we define $a_{j}$,
$j=1,\ldots,K$, the entries of
$\mathbf{e}(\mbox{\boldmath$\alpha$\unboldmath})$ as

$a_{j}=\begin{cases}\widehat{a}_{k},&\text{if
}j=i_{k}\in\\{i_{1},\ldots,i_{q}\\}\\\ 0,&\text{otherwise}\end{cases}$

Clearly, such construction enforces representation (43).

###  5.2 Parameterized heat equation

We first assess performance of the three TROM approaches on a dynamical system
resulting from the discretization of a heat equation

$w_{t}=\Delta w,$ (69)

in a rectangular domain with three holes
$\Omega=\Omega_{r}\setminus(\Omega_{1}\cup\Omega_{2}\cup\Omega_{3})\subset\mathbb{R}^{2}$,
where $\Omega_{r}=[0,10]\times[0,4]$, and the holes are
$\Omega_{1}=[1,3]\times[1,3]$, $\Omega_{2}=[4,6]\times[1,3]$,
$\Omega_{3}=[7,9]\times[1,3]$. The PDE and geometry of $\Omega$ follow that of
[31], while the boundary conditions are modified from those used in [31], as
described below.

We parametrize the system with $D=4$ parameters that enter the boundary
conditions. Convection boundary conditions are enforced on the left side of
the rectangle $\Gamma_{o}=0\times[0,4]$ and on the boundaries of each hole
$\partial\Omega_{j}$, $j=1,2,3$. Explicitly,

$\left.(\mathbf{n}\cdot\nabla w+\alpha_{1}(w-1)\,)\right|_{\Gamma_{o}}=0,$
(70)

and

$\left.\left(\mathbf{n}\cdot\nabla
w+\frac{1}{2}w\right)\right|_{\partial\Omega_{j}}=\frac{1}{2}\alpha_{j+1},\quad
j=1,2,3,$ (71)

i.e., the first parameter in
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathbb{R}^{4}$ is Biot number at
$\Gamma_{o}$ with a fixed outside temperature $t_{o}=1$, while the other three
parameters are the temperatures at $\partial\Omega_{j}$, $j=1,2,3$,
respectively, with Biot numbers equal to $\frac{1}{2}$ on all three hole
boundaries. The rest of the boundary of $\Omega$ is assumed to be insulated

$\left.(\mathbf{n}\cdot\nabla
w)\right|_{\partial\Omega_{r}\setminus\Gamma_{o}}=0.$ (72)

In (70)–(72), $\mathbf{n}$ is the outer unit normal. Observe that the boundary
conditions (70)–(72) can be combined into

$\left.(\mathbf{n}\cdot\nabla
w+q(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})w)\right|_{\partial\Omega}=g(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath}),$

for the appropriate choices of
$q(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ and
$g(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ defined on
$\partial\Omega$ with $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$, a
parameter domain that we take to be the 4D box
$\mathcal{A}=[0.01,0.5]\times[0,0.9]^{3}$. The initial temperature is taken to
be zero throughout $\Omega$.

The system (69)–(72) is discretized with $P_{2}$ finite elements on a quasi-
uniform triangulation of $\Omega$ resulting in $M=3,562$ spatial degrees of
freedom. The choice of standard nodal basis functions
$\\{\theta_{j}(\mathbf{x})\\}_{j=1}^{M}$ defines the mass
$\mathrm{M}\in\mathbb{R}^{M\times M}$ and stiffness
$\mathrm{K}\in\mathbb{R}^{M\times M}$ matrices, as well as boundary terms
$\mathrm{Q}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M\times M}$,
$\mathbf{g}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M}$ with
entries given by

$(\mathrm{Q})_{ij}(\mbox{\boldmath$\alpha$\unboldmath})=\int_{\partial\Omega}q(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\theta_{j}(\mathbf{x})\theta_{i}(\mathbf{x})ds_{\mathbf{x}},\quad(\mathbf{g})_{j}(\mbox{\boldmath$\alpha$\unboldmath})=\int_{\partial\Omega}g(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\theta_{j}(\mathbf{x})ds_{\mathbf{x}},$

for $i,j=1,\ldots,M$.

The vector-valued function of nodal values
$\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath}):[0,T)\times\mathcal{A}\to\mathbb{R}^{M}$
solves

$\mathrm{M}\mathbf{u}_{t}+\left(\mathrm{K}+\mathrm{Q}(\mbox{\boldmath$\alpha$\unboldmath})\right)\mathbf{u}=\mathbf{g}(\mbox{\boldmath$\alpha$\unboldmath}),$
(73)

i.e., it satisfies the dynamical system of the form (1) with

$F(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath})=-\mathrm{M}^{-1}\left(\mathrm{K}+\mathrm{Q}(\mbox{\boldmath$\alpha$\unboldmath})\right)\mathbf{u}+\mathrm{M}^{-1}\mathbf{g}(\mbox{\boldmath$\alpha$\unboldmath})$

and the initial condition
$\mathbf{u}(0,\mbox{\boldmath$\alpha$\unboldmath})=\mathbf{u}_{0}=\boldsymbol{0}\in\mathbb{R}^{M}$
corresponding to zero initial temperature condition for $w$.

We compute the snapshots
$\boldsymbol{\phi}_{k}=\mathbf{u}(t_{k},\mbox{\boldmath$\alpha$\unboldmath})$
by time-stepping (73) at $t_{k}=0.2k$, $k=1,2,\ldots,N$, with $N=100$ time
steps and $T=20$ using Crank-Nicolson scheme. Setting
$\Theta(\mathbf{x})=[\theta_{1}(\mathbf{x}),\ldots,\theta_{M}(\mathbf{x})]$
allows to express the solution
$w(t,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ of (69)–(72) as
$w(t,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})=\Theta(\mathbf{x})\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath})$,
hence the solution snapshots are

$w(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})=\Theta(\mathbf{x})\mathbf{u}(t_{k},\mbox{\boldmath$\alpha$\unboldmath}).$
(74)

The setting is illustrated in Figure 1, where we display the domain $\Omega$
along with solution $w(T,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$
corresponding to parameter values
$\mbox{\boldmath$\alpha$\unboldmath}=(0.5,0,0,0.9)^{T}$.

Fig. 1: Domain $\Omega$ and the solution
$w(T,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ of the heat equation
(69)–(72) corresponding to
$\mbox{\boldmath$\alpha$\unboldmath}=(0.5,0,0,0.9)^{T}$.

For an arbitrary but fixed $\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$
let $\mathrm{Z}=[\mathbf{z}_{1},\ldots,\mathbf{z}_{n}]\in\mathbb{R}^{M\times
n}$ be a matrix with columns being vectors constituting the reduced basis,
i.e. $\mathrm{Z}=\mathrm{U}\mathrm{U}_{c}$ for TROM. Then, the projection ROM
of (73) is

$\widetilde{\mathrm{M}}\widetilde{\mathbf{u}}_{t}+\left(\widetilde{\mathrm{K}}+\widetilde{\mathrm{Q}}(\mbox{\boldmath$\alpha$\unboldmath})\right)\widetilde{\mathbf{u}}=\widetilde{\mathbf{g}}(\mbox{\boldmath$\alpha$\unboldmath}),$
(75)

where

$\displaystyle\widetilde{\mathrm{M}}$
$\displaystyle=\mathrm{Z}^{T}\mathrm{M}\mathrm{Z}\in\mathbb{R}^{n\times
n},\quad\widetilde{\mathrm{K}}=\mathrm{Z}^{T}\mathrm{K}\mathrm{Z}\in\mathbb{R}^{n\times
n},$
$\displaystyle\widetilde{\mathrm{Q}}(\mbox{\boldmath$\alpha$\unboldmath})$
$\displaystyle=\mathrm{Z}^{T}\mathrm{Q}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{Z}\in\mathbb{R}^{n\times
n},\quad\widetilde{\mathbf{g}}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{Z}^{T}\mathbf{g}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{n},$

and the initial condition is
$\widetilde{\mathbf{u}}(0,\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{Z}^{T}\mathbf{u}_{0}=\boldsymbol{0}\in\mathbb{R}^{n}$.
As discussed in Section 3.10, the evaluation of (75) can be effectively split
between the offline and online stages. Solving (75) for
$\widetilde{\mathbf{u}}(t,\mbox{\boldmath$\alpha$\unboldmath})$ allows to
recover the approximate solution at times $t_{k}$ as

$\widetilde{w}(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})=\Theta(\mathbf{x})\;\mathrm{Z}\;\widetilde{\mathbf{u}}(t_{k},\mbox{\boldmath$\alpha$\unboldmath})\approx
w(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath}).$ (76)

####  5.2.1 In-sample prediction and compression study

We begin TROM assessment with in-sample prediction and compression study for
the linear parabolic system described Section 5.2. To measure TROM predictive
power and to compare it to that of POD-ROM, we sample $\mathcal{A}$ uniformly
in each direction with $n_{1}\times n_{2}\times n_{3}\times n_{4}=9\times
5\times 5\times 5$ samples, for a total of $K=1,125$ samples in the set
${\widehat{\mathcal{A}}}=\\{\widehat{\mbox{\boldmath$\alpha$\unboldmath}}_{1},\dots,\widehat{\mbox{\boldmath$\alpha$\unboldmath}}_{K}\\}$.
For each of the three TROMs and for POD-ROM we compute the following in-sample
prediction error

$E_{L^{2}({\widehat{\mathcal{A}}})}=\left(\frac{1}{MNK}\sum_{j=1}^{K}\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})\Phi_{e}(\widehat{\mbox{\boldmath$\alpha$\unboldmath}}_{j})\right\|_{F}^{2}\right)^{1/2},$
(77)

to quantify the ability of the CP-TROM, HOSVD-TROM, TT-TROM local bases and
POD-ROM basis to represent original snapshots for in-sample parameter values.
Note that $E_{L^{2}({\widehat{\mathcal{A}}})}^{2}$ is the quantity from (14)
averaged over ${\widehat{\mathcal{A}}}$ and scaled by $(MN)^{-1}$.

Table 2: Cartesian grid-based sampling in-sample prediction and compression study results reporting prediction errors $E_{L^{2}({\widehat{\mathcal{A}}})}$ for ${\rm HOSVD},~{}{\rm TT},~{}{\rm POD}$, the number of compressed tensors elements transmitted to the online stage, and the corresponding compression factors CF. $\widetilde{\boldsymbol{\Phi}}$ ranks comprise: Tucker ranks for HOSVD-TROM in format $[\widetilde{M},\widetilde{n}_{1},\widetilde{n}_{2},\widetilde{n}_{3},\widetilde{n}_{4},\widetilde{N}]$, and compression ranks for TT-TROM in format $[\widetilde{r}_{1},\widetilde{r}_{2},\widetilde{r}_{3},\widetilde{r}_{4},\widetilde{r}_{5}]$. $\widetilde{\varepsilon}$ | 1e-4 | 1e-5 | 1e-6 | 1e-7 | 1e-9
---|---|---|---|---|---
$n$ | 12 | 16 | 19 | 23 | 30
| HOSVD | 3.45e-05 | 3.85e-06 | 2.78e-07 | 5.61e-08 | –
$E_{L^{2}({\widehat{\mathcal{A}}})}$ | TT | 6.15e-05 | 5.58e-06 | 5.21e-07 | 5.03e-08 | 5.41e-10
| POD | 1.67e-03 | 9.06e-04 | 5.84e-04 | 3.05e-04 | 7.76e-05
| HOSVD | $[34,4,2,$ | $[46,5,2,$ | $[57,6,2,$ | $[66,7,2,$ | –
$\widetilde{\boldsymbol{\Phi}}$ ranks |  | $2,2,12]$ | $2,2,16]$ | $2,2,20]$ | $2,2,23]$ |
| TT | $[34,35,30,$ | $[46,48,41,$ | $[57,61,51,$ | $[69,74,60,$ | $[99,97,80,$
|  | $21,12]$ | $29,16]$ | $36,19]$ | $43,23]$ | $57,30]$
| HOSVD | 13122 | 29515 | 54804 | 85101 | –
$\\#\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ | TT | 20382 | 37993 | 53877 | 86022 | 156607
| HOSVD | 3.05e+4 | 1.36e+4 | 7.31e+3 | 4.71e+3 | –
CF | TT | 1.97e+4 | 1.05e+4 | 7.42e+3 | 4.66e+3 | 2.56e+3
Table 3: General parameter sampling in-sample prediction and compression study results reporting prediction errors $E_{L^{2}({\widehat{\mathcal{A}}})}$ for ${\rm HOSVD},~{}{\rm TT},~{}{\rm POD}$, the number of compressed tensors elements transmitted to the online stage, and the corresponding compression factors CF. $\widetilde{\boldsymbol{\Phi}}$ ranks comprise: Tucker ranks for HOSVD-TROM in format $[\widetilde{M},\widetilde{N}]$, and compression ranks for TT-TROM in format $[\widetilde{r}_{1},\widetilde{r}_{2}]$. $\widetilde{\varepsilon}$ | 1e-4 | 1e-5 | 1e-6 | 1e-7 | 1e-9
---|---|---|---|---|---
$n$ | 12 | 16 | 19 | 23 | 30
| HOSVD | 5.69e-05 | 5.29e-06 | 4.66e-07 | 7.93e-08 | –
$E_{L^{2}({\widehat{\mathcal{A}}})}$ | TT | 5.63e-05 | 5.95e-06 | 5.72e-07 | 5.56e-08 | 5.43e-09
| POD | 2.05e-03 | 1.12e-03 | 5.84e-04 | 3.05e-04 | 8.46e-05
$\widetilde{\boldsymbol{\Phi}}$ ranks | HOSVD | $[33,11]$ | $[45,16]$ | $[56,19]$ | $[65,23]$ | –
| TT | $[32,11]$ | $[43,15]$ | $[55,19]$ | $[66,23]$ | $[95,29]$
| HOSVD | 11904 | 18450 | 26268 | 36680 | –
$\\#\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ | TT | 396011 | 725640 | 1175644 | 1633522 | 3099404
| HOSVD | 3.37e+4 | 2.17e+4 | 1.53e+4 | 1.09e+4 | –
CF | TT | 1.01e+3 | 5.52e+2 | 3.41e+2 | 2.45e+2 | 1.29e+2

We report in Tables 2 and 3 the in-sample prediction errors and compression
factors CF defined in (47) for Cartesian grid-based and general samplings,
respectively. For general sampling we organize the sampling parameters from
the Cartesian grid in 1D array leading to snapshot tensors or order 3.
Experimenting with the same number of randomly sampled parameters showed very
similar compression rates and in-sample prediction errors and so those are not
reported here. The results are reported in Tables 2 and 3 for a number of
decreasing values of $\widetilde{\varepsilon}$ and correspondingly increasing
$n$, such that $n\leq\min(\widetilde{N},\widetilde{r}_{D+1})$, where
$\widetilde{N}$ and $\widetilde{r}_{D+1}$ are the last Tucker and compression
ranks, respectively, for HOSVD- and TT-TROM. Available compression algorithm
failed to deliver the accuracy of 1e-9 for HOSVD, so we report only TT
statistics for this extreme value of $\widetilde{\varepsilon}$. Note that we
leave $E_{L^{2}({\widehat{\mathcal{A}}})}$ for CP-ROM out since there is no
direct way to control its relative error $\widetilde{\varepsilon}$, as
discussed at the end of Section 3.5. Instead, compression factors and
canonical ranks for CP-TROM are illustrated in Figure 2.

We observe in Tables 2 and 3 that both HOSVD- and TT-TROM outperform POD-ROM
in terms of prediction error up to four orders of magnitude for larger $n$ and
correspondingly small $\widetilde{\varepsilon}$. This result is consistent
across both Cartesian grid-based and general parameter samplings. The
$\\#\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ values tell us that for
Cartesian based sampling, HOSVD and TT are comparable in terms of memory and
data transmission requirements with HOSVD doing somewhat better for lower
representation accuracy for the snapshot tensor of order 6. Compression
achieved varies in $\widetilde{\varepsilon}$ (as should expected) and gives
more than 3 orders of saving even for the finest available representation
accuracy. If the Cartesian structure of $\widehat{\mathcal{A}}$ is abandoned
and snapshots are organized in tensors of order 3, then HOSVD format has a
clear advantage over TT in terms of compression achieved; see
$\\#\mbox{online}(\widetilde{\boldsymbol{\Phi}})$ and CF statistics in Table
3.

Fig. 2: CP-TROM compression factor CF and canonical rank $R$ as functions of
$\widetilde{\varepsilon}$. Left: Cartesian grid-based sampling; Right: general
sampling.

For CP-TROM we display in Figure 2 compression factors CF and canonical ranks
$R$ for a number of values of $\widetilde{\varepsilon}$. For comparison we
also show on the same plots CF for HOSVD, which shows that HOSVD-TROM is on
par with CP-TROM if Cartesian sampling allows to organize snapshots in a
higher order tensor. However, we were able to compute
$\widetilde{\boldsymbol{\Phi}}$ in CP compressed format only up to moderate
values of $\widetilde{\varepsilon}$, since the corresponding CP rank was
growing fast as $\widetilde{\varepsilon}$ decreases (see the left plot in
Figure 2). Smaller $\widetilde{\varepsilon}$ become feasible with CP if the
snapshots are organized in tensors of order 3, but in this case HOSVD-TROM
achieves much better compression than CP-TROM (see the right plot in Figure
2). We conclude that for this example with a relatively small number of
parameters ($D=4$) HOSVD-TROM appears to be the best performing TROM. We
finally note that for the same compression accuracy $\widetilde{\varepsilon}$
(if it was achieved) CP-TROM demonstrated very similar in-sample prediction
error as HOSVD- and TT-ROMs. This observation largely carries over to out-of-
sample representation studied next.

####  5.2.2 Out-of-sample prediction study

To quantify the ability of the CP-, HOSVD- and TT-TROM local bases to
represent the solution of (69) for arbitrary out-of-sample parameter values,
we use $E_{L^{2}(\mathcal{A})}$ which is defined as in (77) but with $\alpha$
(in place of $\widehat{\alpha}$) running through a large number of random
points from $\mathcal{A}$. We also use

$E_{L^{\infty}(\mathcal{A})}=\sup_{\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}}\left(\frac{1}{MN}\left\|(\mathrm{I}-\mathrm{Z}\mathrm{Z}^{T})\Phi_{e}({\mbox{\boldmath$\alpha$\unboldmath}})\right\|_{F}^{2}\right)^{1/2},$

for the maximum of representation error over the parameter domain. An estimate
of $E_{L^{\infty}(\mathcal{A})}$ is given by Theorem 1. Regarding constants
appearing in (65) we note that for our choice of the uniform grid in
$\mathcal{A}$ and $p=2,3$ one computes $C_{e}=1$, $C_{a}=\frac{1}{8}$ ($p=2$)
and $C_{a}=\frac{1}{48}$ ($p=2$), while constant $C_{\mathbf{u}}$ is hard to
evaluate. Experimentally we found that the grid in the forth parameter should
be finer than for the first three parameters to balance the observed error
suggesting that the solution is less smooth as a function of $\alpha_{4}$. To
study the error dependence on the parameter mesh size $\delta$, we reduce the
number of parameters to two, letting $\alpha_{1}=\alpha_{2}=\alpha_{3}$ in
(71).

TROM and POD prediction error vs $n$ | TROM prediction error vs $\delta$ | TROM prediction error vs $\widetilde{\varepsilon}$
---|---|---
|  |
Fig. 3: Out-of-sample prediction error of TROM versus ROM basis dimension $n$
(left plot), parameter mesh size $\delta$ (middle plot) and tensor compression
accuracy $\widetilde{\varepsilon}$ (right plot).

In Figure 3 we plot $E_{L^{\infty}(\mathcal{A})}$ and $E_{L^{2}(\mathcal{A})}$
versus ROM basis dimension $n$, parameter mesh size $\delta$ and tensor
compression accuracy $\widetilde{\varepsilon}$. The results were computed
using 100 randomly distributed parameters from $\mathcal{A}$ to evaluate the
error quantities for HOSVD-TROM. For TT-TROM and CP-TROM the error dependence
on $n$, $\delta$ and $\widetilde{\varepsilon}$ was virtually the same and are
omitted. Variants of TROM of course may differ by complexity. While the cost
of offline phase of computing $\widetilde{\boldsymbol{\Phi}}$ varies
significantly depending on the format, the online distribution of costs was
persistent for all three formats showing between 30% and 40% of online time
spent on finding $U_{c}$ (coordinates of local basis), less then 5% on
projection to $\alpha$-specific coordinates, and about 60% of time on the
integration of the projected system (this last step is common with standard
POD approach).

The left plot in Figure 3 shows the $E_{L^{\infty}(\mathcal{A})}$ error for
different values of $n$ and compares it to the error of POD-ROM. We use
$\widetilde{\varepsilon}$=1e-7 and $65\times 33$ parameter grid. Such fine
compression accuracy and grid allows us to isolate the effect of the second
term on the right hand side of (65). Indeed, we see that the error curve
follows closely the graph of “SVD remainder” (the maximum over all $\alpha$ of
the second term on the right-hand side of (65)) until the interpolation error
starts dominating for larger $n$. As can be expected, the interpolation error
for $p=3$ starts affecting $E_{L^{\infty}(\mathcal{A})}$ for larger $n$ than
the interpolation error for $p=2$.

The middle plot in Figure 3 demonstrates that both
$E_{L^{\infty}(\mathcal{A})}$ and $E_{L^{2}(\mathcal{A})}$ for TROM decrease
as $O(\delta^{2})$ for $p=2$ (computed with $\widetilde{\varepsilon}$=1e-7 and
$n=\widetilde{N}$) just as predicted by (65). Likewise, the right plot in
Figure 3 gives evidence for $O(\widetilde{\varepsilon})$ decrease of
$E_{L^{\infty}(\mathcal{A})}$ and $E_{L^{2}(\mathcal{A})}$ (computed with
$65\times 33$ parameter grid and $n=\widetilde{N}$) in accordance to (65) as
long the first term on the right hand side of (65) dominates.

####  5.2.3 Out-of-sample TROMs vs POD-ROM performance: heat equation

Performance of TROM and POD-ROM may vary for different out-of-sample parameter
values. Therefore, assessment of TROM vs POD-ROM performance in this section
is conducted in a statistical setting. The quantities of interest are computed
for $N_{r}\gg 1$ out-of -sample realizations of
$\mbox{\boldmath$\alpha$\unboldmath}^{(r)}\in\mathcal{A}$,
$r=1,2,\ldots,N_{r}$, where we use $N_{r}=200$ for the numerical studies
below. Realizations
$\mbox{\boldmath$\alpha$\unboldmath}^{(r)}=\left(\alpha_{1}^{(r)},\alpha_{2}^{(r)},\alpha_{3}^{(r)},\alpha_{4}^{(r)}\right)^{T}$
are drawn at random from $\mathcal{A}=[0.01,0.5]\times[0,0.9]^{3}$ with each
$\alpha_{i}$ distributed uniformly on $[\alpha_{i}^{\min},\alpha_{i}^{\max}]$,
$i=1,2,3,4$.

For statistical tests we use the following quantities to measure performance
of TROM. First, we introduce the relative $L^{\infty}(0,T,L^{2}(\Omega))$ ROM
solution error

$\begin{split}R_{X}(\mbox{\boldmath$\alpha$\unboldmath})&=\frac{\max\limits_{k=1,\ldots,N}\left\|\widetilde{w}(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})-w(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\right\|_{L^{2}(\Omega)}}{\max\limits_{k=1,\ldots,N}\left\|w(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\right\|_{L^{2}(\Omega)}}\\\
&\approx\frac{\sup\limits_{t\in[0,T]}\left\|\widetilde{w}(t,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})-w(t,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\right\|_{L^{2}(\Omega)}}{\sup\limits_{t\in[0,T]}\left\|w(t,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\right\|_{L^{2}(\Omega)}},\end{split}$
(78)

which we compute for each realization
$\mbox{\boldmath$\alpha$\unboldmath}^{(r)}$, $r=1,2,\ldots,N_{r}$, for both
POD-ROM and each of the three TROMs with X$\in\\{$POD, CP, HOSVD, TT$\\}$. The
true and reduced order snapshots for (78) are computed as in (74) and (76),
respectively. We report the mean, minimum, and standard deviation of the three
relative _gain_ distributions

$G_{\text{X}}^{(r)}=\frac{R_{\text{POD}}\big{(}\mbox{\boldmath$\alpha$\unboldmath}^{(r)}\big{)}}{R_{\text{X}}\big{(}\mbox{\boldmath$\alpha$\unboldmath}^{(r)}\big{)}},\quad
r=1,2,\ldots,N_{r},$ (79)

for X$\in\\{$CP, HOSVD, TT$\\}$, which quantify the error decrease of CP-TROM,
HOSVD-TROM and TT-TROM, respectively, relative to POD-ROM. We study the
dependency of (79) with respect to $K$, the number of sampled parameter values
in $\mathcal{A}$, and $n$, the dimension of the reduced space. The results are
reported for both Cartesian grid-based sampling and general parameter
sampling.

Table 4: Statistics of relative gain (79) for various values of $K$, the number of sampled parameter values in $\mathcal{A}$. The study is performed with $n=10$. | General parameter sampling | Cartesian grid-based sampling
---|---|---
$K$ | $G_{X}$ | CP | HOSVD | TT | CP | HOSVD | TT
| mean | 24.17 | 24.17 | 24.17 | 24.76 | 25.08 | 25.08
$135=$ | min | 0.49 | 0.49 | 0.49 | 0.56 | 0.56 | 0.56
$5\times 3^{3}$ | std | 17.33 | 17.33 | 17.33 | 16.88 | 17.31 | 17.32
| mean | 34.60 | 34.60 | 34.61 | 35.21 | 35.52 | 35.51
$1000=$ | min | 2.80 | 2.80 | 2.80 | 1.72 | 1.72 | 1.72
$8\times 5^{3}$ | std | 15.36 | 15.36 | 15.37 | 15.03 | 15.11 | 15.11
| mean | 38.14 | 38.15 | 38.15 | 37.80 | 38.80 | 38.80
$3430=$ | min | 3.81 | 3.81 | 3.81 | 4.20 | 4.45 | 4.43
$10\times 7^{3}$ | std | 14.20 | 14.21 | 14.22 | 12.96 | 13.61 | 13.62
Table 5: Statistics of relative gain (79) for $n=10$ and $20$. The study is performed for $K=10\times 7\times 7\times 7=3430$. | General parameter sampling | Cartesian grid-based sampling
---|---|---
$n$ | $G_{X}$ | CP | HOSVD | TT | CP | HOSVD | TT
| mean | 38.14 | 38.15 | 38.15 | 37.80 | 38.80 | 38.80
$10$ | min | 3.81 | 3.81 | 3.81 | 4.20 | 4.45 | 4.43
| std | 14.20 | 14.21 | 14.22 | 12.96 | 13.61 | 13.62
| mean | 155.00 | 158.97 | 161.75 | 49.80 | 155.65 | 154.03
$20$ | min | 4.59 | 4.54 | 4.55 | 1.51 | 5.26 | 5.23
| std | 513.33 | 530.50 | 557.86 | 39.48 | 551.92 | 541.88

We present in Table 4 the dependence of relative gain statistics on the values
of $K$ (statistics in the table were computed setting
$\widetilde{\varepsilon}=10^{-5}$ as the targeted accuracy of HOSVD-TROM, TT-
TROM, and $R=250$ as the targeted rank for CP-TROM). We observe that as $K$
increases, TROMs become both more accurate on average and more robust. The
robustness is observed in both the increase of the minimum relative gain and
the decrease of its standard deviation. On average, for $K=3,430$, all three
TROMs are almost $40$ times more accurate compared to POD-ROM. The performance
difference between the TROMs themselves is basically negligible for this
particular study.

The performance of TROMs in the example above is limited by the relatively
small value of $n=10$. The effect of increasing $n$ to $20$ while keeping
$K=3,430$ is shown in Table 5 (statistics in the table were computed setting
$\widetilde{\varepsilon}=10^{-7}$ as the targeted accuracy of HOSVD-TROM, TT-
TROM, and $R=250$ as the targeted rank for CP-TROM). While the worst case
scenario stays relatively unchanged, the average accuracy gain by TROMs is
over two orders of magnitude. An outlier here is CP-TROM that underperforms
HOSVD- and TT-TROM in case of Cartesian grid-based parameter sampling. Aside
from that, performance difference of TROMs for general and Cartesian samplings
is negligible. It is also interesting to see if the interpolation of
approximate snapshots from the lower-rank tensor
$\widetilde{\boldsymbol{\Phi}}$ alone, i.e., without finding a reduced local
basis and solving the projected problem, gives reasonable approximation to
high-fidelity solutions for out-of-sample parameters. Such interpolation-only
predicted solution for incoming
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ is given by columns of
$\widetilde{\Phi}(\mbox{\boldmath$\alpha$\unboldmath})$. Repeating the
experiment with Cartesian grid-based sampling and other parameters the same as
used for results in Tables 4 and 79, we find that the mean relative gain of
HOSVD-tROM compared to interpolation-only approach is
$\\{2.10,\,1.51,\,1.13\\}$ for $n=10$, $K=\\{135,1000,3430\\}$ and
$\\{7.64,\,7.60,\,7.54\\}$ for $n=20$ and same values of $K$. The numbers were
very close for other two TROMs. We see that TROM based on solving projected
problem in general gives more accurate results then pure interpolation,
especially if more vectors are included in the reduced basis. If $n$ is fixed,
then for sufficiently fine sampling the interpolation-only approach delivers
the same (or even better) accuracy.

Table 6: Statistics of relative gain (79) for POD-ROM with POD-Greedy reduced basis for $n=10$ and $20$. The study is performed for $K=10\times 7\times 7\times 7=3430$. | General parameter sampling | Cartesian grid-based sampling
---|---|---
$n$ | $G_{X}$ | CP | HOSVD | TT | CP | HOSVD | TT
| mean | 53.14 | 53.14 | 53.14 | 54.02 | 54.12 | 54.12
$10$ | min | 7.01 | 7.01 | 7.01 | 7.43 | 7.42 | 7.42
| std | 18.63 | 18.63 | 18.63 | 18.03 | 18.03 | 18.03
| mean | 188.57 | 195.57 | 198.77 | 62.18 | 193.66 | 191.86
$20$ | min | 7.39 | 7.35 | 7.36 | 1.87 | 8.51 | 8.47
| std | 571.50 | 607.44 | 634.94 | 49.07 | 649.81 | 639.61

We conclude the numerical study for the parameterized heat equation with a
comparison between the three TROMs and another variant of POD-ROM often used
in practice, the so-called greedy POD or POD-Greedy approach to computing the
reduced basis [59]. Replacing the conventional POD-ROM computation with POD-
Greedy algorithm, we perform the same out-of-sample performance study as
presented in Table 5. The resulting relative gain statistics are reported in
Table 6. Qualitatively, the results are very similar to those in Table 5.
However, quantitatively we observe $20\%$ to $40\%$ increased relative gain
for all TROMs. This is consistent with the fact that POD-Greedy reduced basis
is sub-optimal compared to the conventional POD-ROM reduced basis computed
from the snapshots corresponding to all parameter values in
${\widehat{\mathcal{A}}}$.

In terms of computational performance for this specific setting, POD-Greedy
algorithm was found to be significantly slower in the offline stage than all
three TROM approaches and the conventional POD-ROM even if one includes into
the offline cost of TROM and POD-ROM the computation of $KN$ snapshots
$\boldsymbol{\phi}_{k}(\widehat{\mbox{\boldmath$\alpha$\unboldmath}}_{j})$,
$k=1,\ldots,N$, $j=1,\ldots,K$. Indeed, the bulk of computational cost of POD-
Greedy approach is in the evaluation of error estimator that has to be
performed at each of its $n$ iterations for all
$\widehat{\mbox{\boldmath$\alpha$\unboldmath}}_{j}$, $j=1,\ldots,K,$ to
determine the sample with the largest error. In turn, each evaluation of the
estimator requires the computation of the residual of the chosen time-stepping
scheme (e.g., Crank-Nicolson) that needs $O(N)$ matrix-vector products of a
dense $M\times n$ matrix and a vector in $\mathbb{R}^{n}$. Thus, error
estimator evaluations alone account for $O(KNMn^{2})$ operations of POD-Greedy
offline stage. We note that in other situations, when computing the snapshopts
is much more expensive compared to the evaluation of the error estimator
(e.g., when the matrices of systems to be solved on each time step are dense),
both POD-ROM and TROM may benefit from a greedy approach to parameter
sampling.

###  5.3 Advection-diffusion PDE

In the second numerical example, for assessing performance of the three TROM
approaches we are interested in a case with a higher order of parameter space
compared to $D=4$ in Section 5.2. To that end we set up a dynamical system
resulting from the discretization of a linear advection-diffusion equation

$w_{t}=\nu\Delta
w-\boldsymbol{\eta}(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\cdot\nabla
w+f(\mathbf{x}),$ (80)

in a unit square domain $\Omega=[0,1]\times[0,1]\subset\mathbb{R}^{2}$,
$\mathbf{x}=(x_{1},x_{2})^{T}\in\Omega$. Here $\nu$ is a constant diffusion
coefficient, $\boldsymbol{\eta}:\Omega\times\mathcal{A}\to\mathbb{R}^{2}$ is
the advection field and $f(\mathbf{x})$ is a Gaussian source

$f(\mathbf{x})=\frac{1}{2\pi\sigma_{s}^{2}}\exp\left(-\cfrac{(x_{1}-x_{1}^{s})^{2}+(x_{2}-x_{2}^{s})^{2}}{2\sigma_{s}^{2}}\right),$
(81)

where we take $\sigma_{s}=0.05$, $x_{1}^{s}={x_{2}^{s}}=0.25$. We enforce
homogeneous Neumann boundary conditions and zero initial condition

$\left.\left(\mathbf{n}\cdot\nabla w\right)\right|_{\partial\Omega}=0,\quad
w(0,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})=0.$ (82)

The model is parametrized with $D=9$ parameters with only the advection field
$\boldsymbol{\eta}$ depending on
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathbb{R}^{9}$. The advection field is
given as follows

$\boldsymbol{\eta}(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})=\begin{pmatrix}\eta_{1}(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\\\
\eta_{2}(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\end{pmatrix}=\begin{pmatrix}\cos\alpha_{9}\\\
\sin\alpha_{9}\end{pmatrix}+\frac{1}{\pi}\begin{pmatrix}\partial_{x_{2}}h(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\\\
-\partial_{x_{1}}h(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})\end{pmatrix},$
(83)

where $h(\mathbf{x})$ is the cosine trigonometric polynomial

$\begin{split}h(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})=&\;\;\;\;\alpha_{1}\cos(\pi
x_{1})+\alpha_{2}\cos(\pi x_{2})+\alpha_{3}\cos(\pi x_{1})\cos(\pi x_{2})\\\
&+\alpha_{4}\cos(2\pi x_{1})+\alpha_{5}\cos(2\pi x_{2})+\alpha_{6}\cos(2\pi
x_{1})\cos(\pi x_{2})\\\ &+\alpha_{7}\cos(\pi x_{1})\cos(2\pi
x_{2})+\alpha_{8}\cos(2\pi x_{1})\cos(2\pi x_{2}).\end{split}$ (84)

Here $\alpha_{9}$ determines the angle of the dominant advection direction,
while parameters $\alpha_{i}$, $i=1,\ldots,8$, introduce perturbations into
the advection field. The parameter domain is a 9D box, see Section 5.3.1 for
details.

The system (80)–(82) is discretized using $P_{2}$ finite elements on a grid
with either $M=1,893$ or $M=4,797$ nodes (depending on a particular
experiment) using the standard nodal basis functions
$\\{\theta_{j}(\mathbf{x})\\}_{j=1}^{M}$ that define the mass
$\mathrm{M}\in\mathbb{R}^{M\times M}$, stiffness
$\mathrm{K}\in\mathbb{R}^{M\times M}$ and advection
$\mathrm{H}(\mbox{\boldmath$\alpha$\unboldmath})\in\mathbb{R}^{M\times M}$
matrices, and the source vector $\mathbf{f}\in\mathbb{R}^{M}$. The vector-
valued function of nodal values
$\mathbf{u}(t,\mbox{\boldmath$\alpha$\unboldmath}):[0,T)\times\mathcal{A}\to\mathbb{R}^{M}$
solves

$\mathrm{M}\mathbf{u}_{t}+\left(\mathrm{K}+\mathrm{H}(\mbox{\boldmath$\alpha$\unboldmath})\right)\mathbf{u}=\mathbf{f},$
(85)

i.e., it satisfies the dynamical system of the form (1) with

$F(t,\mathbf{u},\mbox{\boldmath$\alpha$\unboldmath})=-\mathrm{M}^{-1}\left(\mathrm{K}+\mathrm{H}(\mbox{\boldmath$\alpha$\unboldmath})\right)\mathbf{u}+\mathrm{M}^{-1}\mathbf{f},$
(86)

and the initial condition
$\mathbf{u}(0,\mbox{\boldmath$\alpha$\unboldmath})=\boldsymbol{0}\in\mathbb{R}^{M}$.

Similarly to the experiments in Section 5.2, we compute the snapshots
$\boldsymbol{\phi}_{k}=\mathbf{u}(t_{k},\mbox{\boldmath$\alpha$\unboldmath})$
by time-stepping (85) at $t_{k}=(1/30)k$, $k=1,2,\ldots,N$, with $N=30$ time
steps and $T=1$ using Crank-Nicolson scheme. Then, the physical solution
snapshots have the form (74). The setting is illustrated in Figure 4 where we
display the advection field for a random realization of
$\mbox{\boldmath$\alpha$\unboldmath}\in\mathcal{A}$ and the corresponding
solution $w(T,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$.

Advection field $\boldsymbol{\eta}(\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ | Solution $w(T,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$
---|---
|
Fig. 4: Advection field (83) (left) and the solution
$w(T,\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ of (80)–(82) (right)
corresponding to a random realization of $\alpha$ from $\mathcal{B}$ and
$\nu=0.01$.

Projection ROM of (85) is obtained similarly to that of (73) using the matrix
of reduced basis vectors
$\mathrm{Z}=[\mathbf{z}_{1},\ldots,\mathbf{z}_{n}]\in\mathbb{R}^{M\times n}$.
Specifically,

$\widetilde{\mathrm{M}}\widetilde{\mathbf{u}}_{t}+\left(\widetilde{\mathrm{K}}+\widetilde{\mathrm{H}}(\mbox{\boldmath$\alpha$\unboldmath})\right)\widetilde{\mathbf{u}}=\widetilde{\mathbf{f}},$
(87)

where $\widetilde{\mathrm{M}}$ and $\widetilde{\mathrm{K}}$ are defined as in
section 5.2, whereas
$\widetilde{\mathrm{H}}(\mbox{\boldmath$\alpha$\unboldmath})=\mathrm{Z}^{T}\mathrm{H}(\mbox{\boldmath$\alpha$\unboldmath})\mathrm{Z}\in\mathbb{R}^{n\times
n}$, $\widetilde{\mathbf{f}}=\mathrm{Z}^{T}\mathbf{f}\in\mathbb{R}^{n}$, and
the initial condition is
$\widetilde{\mathbf{u}}(0,\mbox{\boldmath$\alpha$\unboldmath})=\boldsymbol{0}\in\mathbb{R}^{n}$.
Efficient evaluation of (87) was discussed in Section 3.10. Solving (87) for
$\widetilde{\mathbf{u}}(t,\mbox{\boldmath$\alpha$\unboldmath})$ allows to
compute the approximate solution snapshots
$\widetilde{w}(t_{k},\mathbf{x},\mbox{\boldmath$\alpha$\unboldmath})$ exactly
as in (76).

####  5.3.1 Out-of-sample TROM performance: advection-diffusion equation

The testing of TROMs for out-of-sample parameters for the advection-diffusion
system (80)–(82) is performed similarly to that for the heat equation in
Section 5.2.3. In particular, we compute the average over $\mathcal{A}$ of the
relative $L^{\infty}(0,T;L^{2}(\Omega))$ and $L^{2}(0,T;H^{1}(\Omega))$ errors
for $\widetilde{w}$ and use the statistical behavior of the TROM vs POD gain
(79) to compare the three TROM variants to conventional POD-ROM.

First, we test TROM in the following setting. The flow behavior is balanced
between advection and diffusion, with a diffusion coefficient $\nu=0.1$ and
$M=1893$. The parameter domain is
$\mathcal{A}=[-0.05,0.05]^{8}\times[0.1\pi,0.3\pi]$ sampled at a Cartesian
grid with $K=3^{8}\times 9=59,049$ points to obtain the sampling set
${\widehat{\mathcal{A}}}$. We draw $N_{r}=200$ out-of-sample realizations
$\mbox{\boldmath$\alpha$\unboldmath}^{(r)}$ from $\mathcal{A}$ with each
$\alpha_{i}^{(r)}$, $i=1,\ldots,9$, distributed uniformly in its corresponding
interval. The averaged relative $L^{\infty}(0,T;L^{2}(\Omega))$ finite element
error of HOSVD-TROM is evaluated as $N_{r}^{-1}\sum\limits_{r=1}^{N_{r}}R_{\rm
HOSVD}(\mbox{\boldmath$\alpha$\unboldmath}^{(r)})$ with $R_{\rm
HOSVD}(\mbox{\boldmath$\alpha$\unboldmath})$ defined in (78). The average
relative $L^{2}(0,T;H^{1}(\Omega))$ finite element error is computed in the
same way after modifying $R_{\rm HOSVD}(\mbox{\boldmath$\alpha$\unboldmath})$
accordingly. The TROM vs POD gain statistic was defined in (79).

For such large values of $K$ and therefore large $\boldsymbol{\Phi}$,
$\\#(\boldsymbol{\Phi})=$ 3.3534e+09, the algorithm we use for finding CP
decomposition turns out to be the most memory-intensive and runs out of
memory. Thus, in what follows we only report the results for HOSVD- and TT-
TROM approaches.

Table 7: Tensor compression ranks, averaged relative error of the HOSVD-TROM finite element solutions, statistics of the gain (79) for $\nu=0.1$, $K=3^{8}\times 9=59,049$. |  $n=10$, $\widetilde{\varepsilon}=10^{-5}$ |  $n=12$, $\widetilde{\varepsilon}=10^{-6}$ |  $n=13$, $\widetilde{\varepsilon}=10^{-7}$
---|---|---|---
HOSVD ranks | [78,3,3,3,3,3, | [116,3,3,3,3,3, | [153,3,3,3,3,3,
| 3,3,3,6,11] | 3,3,3,7,12] | 3,3,3,9,13]
TT ranks | [77,99,117,121,121, | [116,162,204,219,218, | [179,251,330,364,362,
| 107,85,62,37,11] | 186,139,91,50,12] | 298,209,126,63,14]
TROM FE error |  |  |
$L^{\infty}(0,T;L^{2}(\Omega))$ | 1.15e-03 | 8.36e-04 | 7.95e-04
$L^{2}(0,T;H^{1}(\Omega))$ | 9.86e-04 | 9.46e-04 | 9.14e-04
$G_{X}$ | mean | std | min | mean | std | min | mean | std | min
HOSVD |  9.17 |  5.87 | 3.11 |  11.32 |  10.42 | 1.16 |  10.69 |  9.06 | 1.28
TT | 9.13 | 5.85 | 3.11 | 11.40 | 10.61 | 1.17 | 10.65 | 9.06 | 1.22

We present in Table 7 the averaged relative error of the HOSVD-TROM finite
element solutions (for TT-TROM the errors were very close and so are skipped)
and the behavior of the gain statistics when tensor compression error
$\widetilde{\varepsilon}$ decreases, while simultaneously increasing $n$ to be
slightly less or equal to Tucker rank $\widetilde{N}$ for HOSVD or compression
rank $\widetilde{r}_{D+1}$ for TT, respectively. We observe in Table 7 a
relatively weak dependence of the errors and TROM vs POD gain mean on the
choice of $\widetilde{\varepsilon}$ and $n$, hence, we conclude that higher
tensor compression error and smaller $n$ are more beneficial, since they
correspond to higher compression factors and possible faster run times for the
offline stage of TROM algorithms.

For the second example we choose an advection-dominated flow with a smaller
diffusion coefficient $\nu=0.01$ and $M=4,797$. The parameter domain is
$\mathcal{A}=[-0.01,0.01]^{8}\times[0.1\pi,0.5\pi]$ sampled on a Cartesian
grid with $K=20\times 2^{8}=5,120$ points to obtain ${\widehat{\mathcal{A}}}$.
This gives the snapshot tensor with $\\#({\boldsymbol{\Phi}})=$ 1.2280e+09
entries. We draw $N_{r}=100$ out-of-sample realizations
$\mbox{\boldmath$\alpha$\unboldmath}^{(r)}$ from $\mathcal{A}$ with each
$\alpha_{i}^{(r)}$, $i=1,\ldots,9$, distributed uniformly in its corresponding
interval.

Table 8: Tensor compression ranks, number of elements passed to the online