Interpolatory tensorial reduced order models for parametric dynamical systems

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

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

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

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.

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

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

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}$.

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

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

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

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

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

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

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.

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

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

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

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:

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

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}$.

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}$,

or entry-wise

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}$:

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

and compute the thin QR factorizations

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}}$,

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

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

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.

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.

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

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

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,

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

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

Consider the thin SVD of the core matrix

Combining this with the representation above we get

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

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.

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

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

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

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

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

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

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