<figure><img src="image_1.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_2.png" /></figure>

Check for
updates

# A generalized Fourier transform by means of change of variables within multilinear approximation

Mathilde Chevreuil<fnref n="1" />* and Myriam Slama²

**Abstract**

The paper deals with approximations of periodic functions that play a significant role in harmonic analysis. The approach revisits the trigonometric polynomials, seen as combinations of functions, and proposes to extend the class of models of the combined functions to a wider class of functions. The key here is to use structured functions, that have low complexity, with suitable functional representation and adapted parametrizations for the approximation. Such representation enables to approximate multivariate functions with few eventually random samples. The new parametrization is determined automatically with a greedy procedure, and a low rank format is used for the approximation associated with each new parametrization. A supervised learning algorithm is used for the approximation of a function of multiple random variables in tree-based tensor format, here the particular Tensor Train format. Adaptive strategies using statistical error estimates are proposed for the selection of the underlying tensor bases and the ranks for the Tensor-Train format. The method is applied for the estimation of the wall pressure for a flow over a cylinder for a range of low to medium Reynolds numbers for which we observe two flow regimes: a laminar flow with periodic vortex shedding and a laminar boundary layer with a turbulent wake (sub-critic regime). The automatic re-parametrization enables here to take into account the specific periodic feature of the pressure.

**Keywords:** Adapted parametrization, Statistical learning, High-dimensional approximation, Tree tensor networks, Hierarchical tensor format Tensor train format, Adaptive algorithms

**Introduction**

The approximation of periodic functions plays a significant role in harmonic analysis. In the case of the dynamical response of structures, these responses can notably be highly perturbed by low variability on the model and it then becomes necessary to develop reliable and efficient tools for the prediction of the dynamical random response. We are here interested in constructing an approximation of a multivariate function with periodicity in one or more dimensions based on observations. This is of special interest for instance in uncertainty quantification for vibroacoustic problems where the structure is excited with a harmonic wall pressure field. The wall pressure field is a multivariate function which

<figure><img src="image_3.png" /></figure>

<footnote>*Correspondence:
mathilde.chevreuil@univ-
nantes.fr</footnote>

<footnote n="1">¹GeM UMR CNRS 6183,
Université de Nantes, Nantes,
France</footnote>

<footnote>Full list of author information is
available at the end of the article</footnote>
depends on the time and on a set of variables, such as the Reynolds number. In practice, the wall pressure is computed for different instances of the variables and consequently on different discrete time grids depending on the instances. When fine discrete models are involved, the evaluations of the model are costly and we may have access only to sparse information in terms of instances of the variables and of observation time interval.

The set of trigonometric polynomials is well adapted for representing periodic functions. Indeed real trigonometric functions of degree m are written in the following form:

$$v(t) = a_0 + \sum_{n=1}^{m} (a_n \cos(nt) + b_n \sin(nt))$$

where the terms are periodic functions with period $2\pi$. This class of functions has nice properties for approximation use, especially if a function f is continuous and given a tolerance $\epsilon$, there exists a trigonometric polynomial v such that $|f(t) - v(t)| < \epsilon$ for all t. As mentioned above, in many applications we have access to samples of the polynomial from which we want to determine its coefficients. The trigonometric polynomials are therefore linked to discrete-time signal processing, e.g. the Discrete-Time Fourier Transform (DTFT) converts a sequence of length N on an equally spaced time grid into a trigonometric polynomial of degree $N - 1$. The DTFT is extended to the d-dimensional case in the same manner. Given a sample of a multivariate function, the construction of an approximation of the function in the class of trigonometric functions has been widely addressed and the methods for constructing such representation generally depend on the discretization (see [1–3] and the references herein).

In the present paper an alternative approach is proposed in order to tackle such problems using a sample that is not necessarily structured. It is based on statistical learning methods [4] for multi-dimensional problems with s variables where the multivariate output function $u(x_1, \dots, x_s)$ of the model, identified with an order-s tensor, is approximated in a parametrized subset of functions

$$\mathcal{M} = \{v = \Psi(\mathbf{a}); \mathbf{a} \in A\} \quad (1)$$

where the parameter a belongs to some set of parameters A and $\Psi$ is a multilinear function with respect to the parameters a. The key idea is to propose an adapted parametrization with m new variables $z_i = g_i(x_1, \dots, x_s)$, $i = 1, \dots, m$, for the computation of the response so to obtain structured approximations with low complexity by exploiting the periodicity of the function in some dimensions. In the last decade, active subspaces [6] and basis adaptation methods [5] have been proposed to find low dimensional structures using adapted parametrizations with reduced dimension. The first class of methods consist in detecting the directions of strongest variability of a function using gradient evaluations, and then construct an approximation of the function exploiting the defined low-dimension subspace. In [7], active subspaces have been advantageously used for quantifying uncertainty of hypersonic flows around a cylinder. The second class of methods, namely basis adaptation methods, identifies the dominant effects in the form of linear combinations of the input variables and the adapted reduction of the representation is performed through a projection technique. In the current work the change of variables is extended to a wider class of functions g and is done with a method inspired from the projection pursuit method [8] which defines automatically and sequentially the new variables to add. The approximation with this possibly high-dimensional new set of variables is created by exploiting
specific low-dimensional structures of the function such as sparsity [9] and low rank [10] structures of the function to be approximated that enable the construction of an approximation using few samples as introduced in [11,12]. In the latter reference, the output is approximated in suitable low-rank tensor subsets of the form

$$ \mathcal{M} = \{ \nu = \Psi(\mathbf{a}^1, \dots, \mathbf{a}^L); \mathbf{a}^l \in \mathbb{R}^{n_l} \} \qquad (2) $$

where $\Psi$ is a multilinear map with parameters $\mathbf{a}^l$, $l = 1, \dots, L$. This is a special case of (1) with $\mathbf{a} = (\mathbf{a}^1, \dots, \mathbf{a}^L)$ and $A = \mathbb{R}^{n_1} \times \dots \times \mathbb{R}^{n_L}$. The dimension of the parametrization $\sum_{l=1}^L n_l$ grows linearly with $m$ and thus makes possible the approximation with high dimension $m$.

The first part of the paper presents an interpretation of the trigonometric functions as a composition of functions $h \circ g$ with specific structured representations for $h$, and proposes a generalization of the representation for $h$. The second part is dedicated to the algorithm for constructing the rank structured approximation combined with the change of variable $g(x)$ used to handle periodic functions. Finally the last part illustrates the method on the pressure of a flow around a cylinder.

## Periodic functions

Let $u \in \mathcal{H}$ be a multivariate function, with $\mathcal{H}$ a Hilbert space, which depends on a set of independent variables $X = (X_1, \dots, X_s)$. In the present paper, we consider the specific case where the function $u$ is periodic with respect to one variable denoted $\tau$ so that we have $X = (\Xi, \tau)$, with $\Xi = (X_1, \dots, X_d)$ and $\tau = X_{d+1}$, and $s = d + 1$.

A variable $X_v$ has values in $\mathcal{X}_v$ and an associated measure $dp_v$, $1 \le v \le s$. The variable $X$ has values in $\mathcal{X} = \mathcal{X}_1 \times \dots \times \mathcal{X}_s$ and an associated measure $dp = dp_1 \times \dots \times dp_s$. The variables $X_v$, $v = 1, \dots, d$ can be random variables, $dp_v$ being in that case a probability measure on $\mathcal{X}_v$. Let $\mathcal{H}$ be the Hilbert space defined on $\mathcal{X}$, it is a tensor space $\mathcal{H} = \mathcal{H}_1 \otimes \dots \otimes \mathcal{H}_s$ with $\mathcal{H}_v$ a Hilbert space defined on $\mathcal{X}_v$. We consider that $\mathcal{H}_v \subset L_{p_v}^2(\mathcal{X}_v)$ is a finite dimensional subspace of square-integrable functions equipped with the norm $\|u\|_v^2 = \int_{\mathcal{X}_v} u^2 dp_v$ and $\mathcal{H}$ is a subspace of $L_p^2(\mathcal{X})$ equipped with the canonical norm $\|u\|^2 = \int_{\mathcal{X}} u^2 dp$. Let $\{\psi_i^v\}_{i=1}^{P_v}$ be an orthonormal basis of $\mathcal{H}_v$ and $\{\psi_{i_1}^1 \otimes \dots \otimes \psi_{i_s}^s\}_{(i_1, \dots, i_s) \in [1, \dots, P_1] \times \dots \times [1, \dots, P_d]}$ be an orthonormal basis of $\mathcal{H}$.

A natural representation of the periodic function $u$ can be obtained using the Fourier series. In the following, $x = (\xi, t)$ will denote an observation of $X$ with $\xi$ an observation of $\Xi = (X_1, \dots, X_d)$ and $t$ an observation of $\tau$, i.e. a point in the periodic dimension.

## Trigonometric functions as a composition of functions

Let us consider a T-periodic and continuous real valued function $u : \mathbb{R} \to \mathbb{R}$. It can be represented by its truncated Fourier series which is a sum of harmonic functions:

$$ u(t) \approx \sum_{n=0}^{m-1} a_n \cos(\omega_n t) + b_n \sin(\omega_n t) \qquad (3) $$
where the circular frequencies are such that $\omega_n = n\omega_f$ (multiplicative constraint on $\omega_n$), with $\omega_f = \frac{2\pi}{T}$, and the coefficients $a_n$ and $b_n$ defined as follows:

$$a_n = \frac{\omega_f}{2\pi} \int_{-\frac{\pi}{\omega_f}}^{\frac{\pi}{\omega_f}} u(t) \cos(\omega_n t) \, dt \quad \text{and} \quad b_n = \frac{\omega_f}{2\pi} \int_{-\frac{\pi}{\omega_f}}^{\frac{\pi}{\omega_f}} u(t) \sin(\omega_n t) \, dt, \qquad (4)$$

for $n = 0, \dots, m-1$.

The truncated Fourier series can be seen as a composition of functions of the form:

$$v(t) = h(g_1(t), \dots, g_m(t))$$

where $h$ is an additive model $h(z_1, \dots, z_m) = \sum_{n=1}^{m} h_n(z_n)$ with $h_n(z_n) = a_n \cos(z_n) + b_n \sin(z_n) \in V$, where $V = \text{span}\{1, \cos(\cdot), \sin(\cdot)\}$, and $g_n(t) = \omega_{n-1}t$ for $n = 1, \dots, m$.

We propose here to extend the Fourier series to a more generalized framework where the function $v$ is a multivariate function depending both on $t$ and $\xi$, where the circular frequencies $\omega_n$ are chosen adaptively with no multiplicative constraint and where $h$ is chosen in a wider set of functions than additive models.

## Generalizing Fourier series for the representation of a multivariate function in tensor format

Let us focus on a function $u(x) = u(\xi, t)$, periodic with respect to $t$. A natural representation of $u$ is thus by means of a relevant change of variables:

$$u(x) \approx h(g_1(x), \dots, g_m(x), \xi) \qquad (5)$$

where $g_n : \mathbb{R}^{d+1} \to \mathbb{R}$, $n = 1, \dots, m$ are new variables chosen under the form

$$g_n(\xi, t) = \omega_n(\xi)t \qquad (6)$$

and where $h \in V^{\otimes m} \otimes \mathcal{H}_1 \otimes \dots \otimes \mathcal{H}_d : \mathbb{R}^{[0,T]^m \times \Xi_1 \times \dots \times \Xi_d} \to \mathbb{R}$. Representing the periodic function under the form (5) can lead to the definition of a high number $m$ of new variables so that we will consider subsets of low rank tensor formats for the high dimensional function $h$ of dimension $M = m + d$.

A $M$-dimensional function $v$ in a subset of tensors can be written:

$$v(z) = \Psi(z)(\mathbf{a}_1, \dots, \mathbf{a}_L) \qquad (7)$$

where $\Psi(x)$ is a multilinear map with parameters $(\mathbf{a}_1, \dots, \mathbf{a}_L)$. We consider here a model class of rank-structured functions associated with a notion of rank. A well known rank is the canonical rank associated to the sum of multiplicative models. The canonical rank of a function $h$ is the minimal integer $\text{rank}_C(v) = r$ such that

$$v(z) = \sum_{k=1}^{r} v_k^1(z_1) \dots v_k^M(z_M)$$

and we define the subset of canonical tensors:

$$\mathcal{T}_r^C = \{v \in \mathcal{H} : \text{rank}_C(v) \le r\}. \qquad (8)$$
It can be associated to the parametrized representation (7) with $L = M$, where $\mathbf{a}_l \in \mathbb{R}^{r \times n_l}$ with $n_l$ the dimension of the functional basis on which the functions $v_k^l$ are represented, $k = 1, \dots, r$ and $l = 1, \dots, M$. We can consider other notions of ranks which provide different models with lower complexity. The $\alpha$-rank of $v$, denoted by $\text{rank}_\alpha(v)$, is the minimal integer $r_\alpha$ such that

$$v(z) = \sum_{k=1}^{r_\alpha} v_k^\alpha(z_\alpha) v_k^{\alpha^c}(z_{\alpha^c})$$

with $\alpha \subset \{1, \dots, M\}$, and $z_\alpha$ and $z_{\alpha^c}$ complementary groups of variables. The $T$-rank of $v$, denoted by $\text{rank}_T(v) = \{\text{rank}_\alpha(v) : \alpha \in T\}$, is the tuple $r = \{r_\alpha\}_{\alpha \in T}$ such that

$$v(z) = \sum_{k=1}^{r_\alpha} v_k^\alpha(z_\alpha) v_k^{\alpha^c}(z_{\alpha^c}), \quad \forall \alpha \in T$$

where $T$ is a collection of subsets of $\{1, \dots, M\}$. We define the subset of rank-structured functions:

$$\mathcal{T}_r^T = \{v \in \mathcal{H} : \text{rank}_T(v) \le r\}. \qquad (9)$$

The complexity of the associated parametrized representations of tensors is linear with the dimension $M$ and polynomial with the ranks.

## Statistical learning method for approximating a function in tensor format with a change of variables

### Supervised statistical learning

We consider a model that returns a real-valued variable $Y = u(X)$. An approximation $v$ of the function $u$, also referred to as metamodel, can be obtained by minimizing the risk

$$\mathcal{R}(v) = \int_X \ell(u(x), v(x)) dp(x)$$

over a model class $\mathcal{M}$. The loss function $\ell$ measures a distance between the observation $u(x)$ and the prediction $v(x)$. In the case of the least squares method, it is chosen as $\ell(y, v(x)) = (y - v(x))^2$.

Let $S = \{(x^k, y^k) : 1 \le k \le N\}$ be a sample of $n$ realizations of $(X, Y)$. In practice, the approximation is constructed by minimizing the empirical risk

$$\mathcal{R}_S(v) = \frac{1}{N} \sum_{k=1}^{N} \ell(y^k, v(x^k)), \qquad (10)$$

taken as an estimator of the risk. A regularization term $R$ can be used for stability reasons when the training sample is small. An approximation $\tilde{u}$ of $u$ is then solution of

$$\min_{v \in \mathcal{M}} \frac{1}{N} \sum_{k=1}^{N} (y^k - v(x^k))^2 + \lambda R(v), \qquad (11)$$

with the regularization parameter $\lambda \ge 0$, chosen or computed. The accuracy of the metamodel $\tilde{u}$ is estimated using an estimator of the $L^2$ error. In practice, the number of numerical experiments is too small to sacrifice part of it for the error estimation. The
error is thus estimated using a k-fold cross validation estimator and more specifically the leave-one-out cross validation estimator [4] which can be easily evaluated by constructing one single metamodel [13]. Cross validation estimators can be used for model selection.

In the following, we present a method to determine $\tilde{u}(x)$ in a sequence of model classes

$$ \mathcal{M}(m) = \{v(x) = h(g_1(x), \dots, g_m(x), \xi); h \in \mathcal{M}_{hv}, g_i \in \mathcal{M}_g \text{ for } i = 1, \dots, m\} \quad (12) $$

with $\mathcal{M}_h$ a linear or multilinear model class and $g_i \in \mathcal{M}_g$ a linear model class. We consider here multilinear models and more specifically the tensor subset $\mathcal{T}_r^T$ in order to handle the possibly high dimensional $m + d$ problem. We briefly recall the learning algorithm in a tensor subset [12] and then present the automatic computation of the new variables $z_i = g_i(x)$, $i = 1, \dots, m$.

## Learning with tensor formats

Let $z \in \mathbb{R}^M$, an approximation of $u$ in a tensor subset (2) can be obtained by minimizing the empirical least squares risk:

$$ \min_{\mathbf{a}^1, \dots, \mathbf{a}^L} \frac{1}{N} \sum_{k=1}^{N} \left( y^k - \Psi(z^k)(\mathbf{a}_1, \dots, \mathbf{a}_L) \right)^2 + \sum_{i=1}^{L} \lambda_i R_i(\mathbf{a}_i) \quad (13) $$

where $\lambda_i R_i(\mathbf{a}_i)$ are regularization functions. Problem (13) is solved using an alternating minimization algorithm which consists in successively solving an optimization problem on $\mathbf{a}^j$

$$ \min_{\mathbf{a}_j} \frac{1}{N} \sum_{k=1}^{N} \left( y^k - \Psi(z^k)(\dots, \mathbf{a}_j, \dots) \right)^2 + \lambda_j R_j(\mathbf{a}_j) \quad (14) $$

for fixed parameters $\mathbf{a}_i$, $i \neq j$. Introducing the linear map $\Psi^j(z)(\mathbf{a}_j) = \Psi(z)(\mathbf{a}_1, \dots, \mathbf{a}_L)$, problem (14) yields the following learning problem with a linear model

$$ \min_{\mathbf{a}_j} \frac{1}{N} \sum_{k=1}^{N} \left( y^k - \Psi^j(z^k)(\mathbf{a}_j) \right)^2 + \lambda_j R_j(\mathbf{a}_j). \quad (15) $$

If $R_j(\mathbf{a}_j) = \|\mathbf{a}_j\|_1$, with $\|\mathbf{v}\|_1 = \sum_{i=1}^{\#v} |v_i|$ the $l_1$-norm, problem (15) is a convex optimization problem known as Lasso [14] or basis pursuit [15]. The $l_1$-norm is a sparsity inducing regularization function that yields a solution $\mathbf{a}_j$ of (15) that may have coefficients equal to zero. The Lasso is solved using the modified least angle regression algorithm (LARS in [16]).

The algorithm to solve Problem (13) is described in [11] for the canonical tensor format, and in [12] for the tree-based tensor format, which is a special case of rank-structured tensor where $T$ is a dimension partition tree. Adaptive algorithms are proposed to automatically select the tuple yielding a good convergence of the approximation with respect to its complexity. For tree-based tensor formats, at an iteration $i$, given an approximation $v^i$ of $u$ with $T$-rank $(r_\alpha^i)_{\alpha \in T}$, the strategy consists in estimating and studying the truncation error $\min_{\text{rank}_\alpha(v)} \|R(v) - \mathcal{R}(u)\|$ for different $\alpha$ in $T$, and choosing to increase the ranks $r_\alpha^i$ associated with the indices $\alpha$ yielding the highest errors. The algorithm and more details in the tree-based tensor format case can be found in [12, 17].
### Learning method with automatic definition of new variables

We now present the method used to automatically search an adapted parametrization of the problem by looking for favored directions in the space of the $d+1$ input variables. It consists in writing the approximation under the form (5) where $g_n(x)$ can be represented with a parametrized linear map

$$g_n(x) = \mathbf{w}_n^\top \varphi(x) = \sum_{j=1}^{p} w_{n,j} \varphi_j(x) \qquad (16)$$

with $\mathbf{w}_n = (w_{n,1}, \dots, w_{n,p})^\top \in \mathbb{R}^p$ the vector of parameters of the representation of $g_n$ on an orthonormal functional basis $\{\varphi_j\}_{j=1}^p$ of $\mathcal{H}$, and $h$ a $M$-dimensional function in the model class of rank structured formats $\mathcal{T}_r^C$ or $\mathcal{T}_r^T$ that can be represented with a parametrized multilinear map

$$h(z) = \Psi(z)(\mathbf{a}^1, \dots, \mathbf{a}^L) \qquad (17)$$

with parameters $\mathbf{a}^l$, $l = 1, \dots, L$. The new set of variables $z = (z_1, \dots, z_m, \xi)$ is such that $z_n = g_n(x)$, $n = 1, \dots, m$.

The method is a Projection Pursuit like method [8] that is generalized to a larger class of models for $h$ than just additive models. It is shown in [18] that under reasonable conditions, for $h \in \mathcal{T}_r^T$ and $g \in \mathcal{H}$ we have $h \circ g \in L_p^2(\mathcal{X})$. The approximation $\tilde{u}$ of the form (5) is thus parametrized as follows:

$$\tilde{u}(x) = \Psi(\mathbf{w}_1^\top \varphi(x), \dots, \mathbf{w}_m^\top \varphi(x), \xi)(\mathbf{a}^1, \dots, \mathbf{a}^L). \qquad (18)$$

Let $(z_1, \dots, z_{m-1})$ be an initial set of variables. A new variable $z_m = g_m(x)$ is introduced using Algorithm 1.

**Algorithm 1** Addition of a new variable $z_m$.

**Require:** orthonormal set of functions $\{\varphi_j\}_{j=1}^p$, parametrization $w = (w_1, \dots, w_{m-1})$ of the set of variables $(z_1, \dots, z_{m-1})$

**Ensure:** parametrization $w = (w_1, \dots, w_m)$ of $(z_1, \dots, z_m)$

1: set $w_m$ such that $w_m^\top \varphi(x) = 1 \quad \forall x \in \mathbb{R}^d$

2: set $w = (w_1, \dots, w_m)$

3: set $z_m = w_m^\top \varphi(\cdot)$

The parameters $\mathbf{a}_l$, $l = 1, \dots, L$, and $\mathbf{w}_n$, $n = 1, \dots, m$, solve the minimization problem (11) over the model class $\mathcal{M}(m)$:

$$\min_{\{\mathbf{a}_l\}_{l=1}^L, \{\mathbf{w}_n\}_{n=1}^m} \frac{1}{N} \sum_{k=1}^N \left( y^k - \Psi \left( \mathbf{w}_1^\top \varphi(x^k), \dots, \mathbf{w}_m^\top \varphi(x^k), \xi^k \right) (\mathbf{a}^1, \dots, \mathbf{a}^L) \right)^2 + \sum_{i=1}^L \lambda_i R_i(\mathbf{a}_i) \qquad (19)$$

with $x^k = (\xi^k, t^k)$. The solution of this problem is found by alternatively solving

$$\min_{\mathbf{a}_1, \dots, \mathbf{a}_L} \frac{1}{N} \sum_{k=1}^{N} \left( y^k - \Psi \left( \mathbf{w}_1^\top \varphi(x^k), \dots, \mathbf{w}_m^\top \varphi(x^k), \xi^k \right) (\mathbf{a}^1, \dots, \mathbf{a}^L) \right)^2 + \sum_{i=1}^{L} \lambda_i R_i(\mathbf{a}_i) \quad (20)$$
for fixed ($\mathbf{w}_1, \dots, \mathbf{w}_m$) using a learning algorithm with rank adaptation [11,12] and

$$ \min_{\mathbf{w}_1 \in \mathbb{R}^p, \dots, \mathbf{w}_m \in \mathbb{R}^p} \frac{1}{N} \sum_{k=1}^{N} \left( y^k - \Psi \left( \mathbf{w}_1^\top \varphi(x^k), \dots, \mathbf{w}_m^\top \varphi(x^k), \xi^k \right) (\mathbf{a}^1, \dots, \mathbf{a}^L) \right)^2 \quad (21) $$

for fixed ($\mathbf{a}^1, \dots, \mathbf{a}^L$). The optimization problem (21) is a nonlinear least squares problem that is solved with a Gauss-Newton algorithm. The overall algorithm is presented in Algorithm 2.

**Algorithm 2** Learning algorithm for an approximation with automatic definition of the new variables.

Require: sample $S = \{(x^k, y^k) : 1 \le k \le N\}$, collection of M parametrizations $\Psi_i(\cdot)(\mathbf{a}^1, \dots, \mathbf{a}^{L_i})$ of functions of $i + d$ variables, initial values for the parameters ($\mathbf{a}^1, \dots, \mathbf{a}^{L_i}$), basis for the representation of the $g_i$, $i = 1, \dots, m_{\max}$

Ensure: approximation $v(x) = \Psi(\mathbf{w}_1^\top \varphi(x), \dots, \mathbf{w}_m^\top \varphi(x), \xi)(\mathbf{a}^1, \dots, \mathbf{a}^L)$ with selected optimal $m \in \{1, \dots, m_{\max}\}$

1: **for** $i = 1, \dots, m_{\max}$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;add one variable $g_i(x) = \mathbf{w}_i^\top \varphi(x)$ using Algorithm 1

3: &nbsp;&nbsp;&nbsp;&nbsp;set $z = (\mathbf{w}_1^\top \varphi(x), \dots, \mathbf{w}_m^\top \varphi(x), \xi)$

4: &nbsp;&nbsp;&nbsp;&nbsp;**while** not converged **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute the approximation $\Psi_i(z)(\mathbf{a}^1, \dots, \mathbf{a}^{L_m})$ solution of (20)

6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;compute $\mathbf{w}_1, \dots, \mathbf{w}_m$, solutions of (21)

7: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;set $z = (\mathbf{w}_1^\top \varphi(x), \dots, \mathbf{w}_m^\top \varphi(x), \xi)$

8: &nbsp;&nbsp;&nbsp;&nbsp;**end while**

9: **end for**

10: select $m$ which minimizes $\mathcal{R}_V(v_m)$ with $V$ a validation set independent of $S$

11: set $\Psi \leftarrow \Psi_m, L \leftarrow L_m$

In step 5 of algorithm 2, the parametrization of the model is to be selected in a collection of parametrizations. In this paper, we consider the rank structured function $\mathcal{T}_r^T$ where $T$ is a dimension partition tree over $\{1, \dots, m\}$, which corresponds to the model class of functions in tree-based tensor format [19], a particular class of tensor networks [20]. The new node associated with the new variable $z_i = g_i(x)$ is added at the top of the tree. The representation of a function $v$ in $\mathcal{T}_r^T(\mathcal{H})$ requires the storage of a number of parameters that depends both on the collection $T$ and on the associated $T$-rank $r$, as well as on the dimensions of the functional spaces in each dimension. To reduce the number of coefficients that need to be computed during the learning process of a tensor $v$, one can then represent $v \in \mathcal{T}_r^T(\mathcal{H})$ for different collections $T$ and associated $T$-ranks $r$, choosing the ones yielding the smallest storage complexity. Furthermore, this adaptation can prove useful when dealing with changes of variables, as introduced in the previous subsection, because it can remove the difficulty of how to add a variable: no matter what the initial ordering of the variables is, an adaptation procedure may be able to find an optimal one yielding a smaller storage complexity. When $T$ is a dimension partition tree, a stochastic algorithm is presented in [12] for trying to reduce the storage complexity of a tree-based tensor format at a given accuracy. This adaptation is not considered in the paper.

## Change of variables for periodic functions

In this section we present a generalization of the Fourier series where one does not need a structured sample (e.g. a grid) in the variable $t$. Indeed when the function to approximate is known to be periodic with respect to $t$, the periodicity of the approximation can be forced. It is done on the one hand by introducing a functional basis $\varphi$ such that $g_n = \mathbf{w}_n^\top \varphi$
in (16) can be identified to (6) by choosing

$$ \varphi_j(x) = \varphi_j^{(\xi)}(\xi)t \qquad (22) $$

where {$\varphi_j^{(\xi)}$} is a d-dimensional tensorized orthogonal basis of $\mathcal{H}_1 \otimes \dots \otimes \mathcal{H}_d$. The circular frequencies in (6) are then expressed as:

$$ \omega_n(\xi) = \sum_{j=1}^{p} w_{n,j} \varphi_j^{(\xi)}(\xi) = \mathbf{w}_n^\top \varphi^{(\xi)}(\xi). \qquad (23) $$

On the other hand, we choose bases of trigonometric functions {$\psi_i^n$} for the representation of h in the dimensions associated with the new variables $z_n$, $n = 1, \dots, m$.

Let $T_{max}$ be the maximal width of the observation interval in the dimension $t$. Supposing this interval is large enough to include the largest periods of the periodic functions that can be learned, the guarantee for the approximation not to have larger periods is to constrain the circular frequencies in (6) such that:

$$ \omega_n \geq \frac{2\pi}{T_{max}}, \quad n = 1, \dots, m. \qquad (24) $$

This constraint is imposed for all values taken by $\xi$ in $S$. Using expression (23), it is recast under the form

$$ -Aw \leq -B \qquad (25) $$

where $w = (w_1, \dots, w_m) \in \mathbb{R}^{p \times m}$, $A \in \mathbb{R}^{N \times p}$ is the array of evaluations of $\varphi_j^{(\xi)}(\xi)$ for the values $(x_1^k, \dots, x_d^k)^N$ of $\xi$ in the training set $S$:

$$ A = \begin{bmatrix} \varphi_1^{(\xi)}(x_1^1, \dots, x_d^1) & \dots & \varphi_p^{(\xi)}(x_1^1, \dots, x_d^1) \\ \vdots & \ddots & \vdots \\ \varphi_1^{(\xi)}(x_1^N, \dots, x_d^N) & \dots & \varphi_p^{(\xi)}(x_1^N, \dots, x_d^N) \end{bmatrix}, \qquad (26) $$

and $B \in \mathbb{R}^{N \times m}$ is a full array with values $-2\pi/T_{max}$. The optimization problem (21) for the computation of the parameters $w$ is replaced with the constrained optimization problem

$$ \min_{w_1 \in \mathbb{R}^{p}, \dots, w_m \in \mathbb{R}^{p}} \frac{1}{N} \sum_{k=1}^{N} \left( y^k - \Psi \left( w_1 \varphi(x^k), \dots, w_m \varphi(x^k), \xi^k \right) (a^1, \dots, a^L) \right)^2 \qquad (27) $$

which is solved with a NonLinear Programming (NLP) method.

## Application

The method is applied to the prediction of the wall pressure $p$ for a flow over a cylinder for two ranges (low and medium) of Reynolds numbers for which we observe two flow regimes: a laminar flow with periodic vortex shedding and a laminar boundary layer with a turbulent wake (sub-critic regime). The automatic re-parametrization enables here to take into account the specific periodic feature of the pressure $p$ with the time. The variables of the problem are $X = (Re, \Theta, \tau)$ with $\tau$ the time, $Re$ the Reynolds number and $\Theta$ the angular steps, that take values in $\mathcal{X} = \mathcal{X}_{Re} \times \mathcal{X}_{\Theta} \times \mathcal{X}_{\tau}$ and we have $d = 2$. The pressure is evaluated on a tensor grid with
<figure><img src="image_4.png" /><figcaption>Fig. 1 Example of a linear tree with five variables yielding the tensor train format [21]:
T = {{1}, {1, 2}, ..., {1, 2, 3, 4}}</figcaption></figure>

* at low Reynolds numbers: 300 time steps in $\mathcal{X}_\tau = [27 \text{ s}, 29.4 \text{ s}]$, 50 simulations of Re chosen uniformly on $\mathcal{X}_{Re} = [70, 200]$ and 128 angular steps in $\mathcal{X}_\Theta = [0, 2\pi[$,

* at higher Reynolds numbers: 1500 time steps, 11 simulations of Re chosen on $\mathcal{X}_{Re} = [7000, 13000]$ and 320 angular steps in $\mathcal{X}_\Theta = [0, 2\pi[$.

A representation (5) with new variables $g_n(\xi, t)$, where $\xi = (Re, \theta)$, is computed with Algorithm 2 where the optimization problem (21) of step 6 is replaced with the constrained optimization problem (27). The new variables $g_n$ are represented on a basis of functions as in (22) where $\{\varphi_j^{(\xi)}\}_{j=1}^P$ is a polynomial basis with maximal partial degree 2 ($P = 9$). We choose for the model class of $h$ the tensor train (TT) format associated to the linear tree $T = \{\{1\}, \{1, 2\}, \dots, \{1, \dots, M-1\}\}$ (see Fig. 1) where $M = m + d$, the associated parametrized subset is

$$ \{v = \Psi(\mathbf{a}^1, \dots, \mathbf{a}^M); \mathbf{a}^l \in \mathbb{R}^{r_{l-1} \times n_l \times r_l}, l = 1, \dots, M\} $$

with $r = (r_0, r_1, \dots, r_M)$ the TT-rank where $r_0 = r_M = 1$, and $n_l$ is the dimension of the functional basis $\{\psi_l\}_{l=1}^{n_l}$ in dimension $l$ for the representation of $h$. Here we use trigonometric bases $\{1, \cos(x), \sin(x)\}$ in the dimensions of $z_i = g_i(t, \theta, Re)$ and polynomial bases in the dimensions of $\theta$ and $Re$. The sequential quadratic programming (SQL) method is used to solve (27).
The accuracy of the metamodel $\tilde{u}$ is estimated using the unbiased test error estimator based on a test set $S_{test}$ of $N_{test}$ realizations of $(X, Y)$ independent of the training set $S$:

$$e_{test}^2 = \frac{\mathcal{R}_{S_{test}}(\tilde{u})}{\mathcal{R}_{S_{test}}(0)} \qquad (28)$$

where $\mathcal{R}_{S_{test}}(\nu)$ is the empirical risk defined in (10). Considering the least squares estimator, we have $e_{test}^2 = \frac{\frac{1}{N_{test}} \sum_{k=1}^{N_{test}} (y^k - \tilde{u}(x^k))^2}{\frac{1}{N_{test}} \sum_{k=1}^{N_{test}} (y^k)^2}$.

## Low Reynolds numbers

We first consider the approximation of the wall pressure for low Reynolds numbers for which the data is given for 50 values of the Reynolds number in $\mathcal{X}_{Re} = [70, 200]$. The approximation is constructed with the following setting:

* for $g_i$: polynomial bases with a maximal degree of 1 for $t$ and 2 for $\theta$ and $Re$,

* for $h$: polynomial bases of degree 14 for $\theta$ and 3 for $Re$,

* training sample $S$ using 20 simulations of $Re$ and considering only the 181 first time steps.

The algorithm provided an approximation with dimension $M = 6: X = (g_1(\xi), \dots, g_4(\xi), \xi)$ and TT-ranks $r = [1 3 4 4 2 1]$. The model error was estimated using two different samples:

* estimation of the approximation error on the sample $S_{test}$ with 30 simulations of $Re$ on the train time range (consisting of the 181 first time steps): $e_{test}^2 = 1.14\%$,

* estimation of the extrapolation error on the sample $S_{extra}$ with 30 simulations of $Re$ on the whole time range (consisting of the 300 available time steps): $e_{extra}^2 = 1.19\%$.

On Fig. 2 are plotted the predictions with blue crosses and the observations with circles, the red ones were used for learning the approximation and the green ones for estimating the model error. We observe a very good match of the predictions with the observations even beyond the train time range.

The approximation constructed using the proposed change of variable is able to extrapolate the wall pressure beyond the time range used for training the approximation. This extrapolation was made possible by introducing the constraint (24). As an illustration, Fig. 3 shows the observations (in blue) versus the predictions (in red) on a longer time interval obtained without constraint. The approximation obviously loses its periodicity.

Table 1 summarizes the results obtained with the proposed approach and those obtained using the change of variable without using the tensor train format. One can observe that tensor formats enable to break the curse of dimensionality and thus to ease the learning of the approximation based on observations. Indeed the storage complexity of tensor train formats is $O(MnR^2)$ where $n$ is the order of the dimension of the representation space in each dimension and $R$ is the order of the rank. That is, it grows only linearly with the dimension $M$ and quadratically with the rank whereas the storage complexity without using tensor formats, i.e. on the polynomial chaos, grows factorially or exponentially with the dimension $M$. Exploiting low complexity representations as low rank structures is necessary to address the problem when the dimension $M$ increases with the definition of new variables.
<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 2** Low Reynolds numbers. Observations and predictions with respect to time for fixed values of the other parameters. The predictions are obtained using tensor formats combined with constrained changes of variables</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 3** Low Reynolds numbers. Observations and predictions with respect to time for fixed values of the other parameters. The predictions are obtained using tensor formats combined with changes of variables (unconstrained)</figcaption></figure>

<table><caption>**Table 1** Storage complexity and test error of the approximations obtained with the change of variable and with and without using Tensor Train (TT) format</caption><thead><tr><td></td><td>Storage complexity</td><td>Test error e<sup>2</sup><sub>test</sub> (%)</td></tr></thead><tbody><tr><td>TT format + change of variables</td><td>594</td><td>1.14</td></tr><tr><td>Polynomial chaos + change of variables</td><td>8482</td><td>5.42</td></tr></tbody></table>

## High Reynolds numbers

We now consider the approximation of the wall pressure for higher Reynolds numbers for which the data is given for 11 values of the Reynolds number in $\mathcal{X}_{Re} = [7000, 13000]$. The approximation is constructed with the following setting:

* for $g_i$: polynomial bases with a maximal degree of 1 for $t$ and 2 for $\theta$ and $Re$,

* for $h$: polynomial bases of degree 20 for $\theta$ and 6 for $Re$,

* training sample $S$ using 8 simulations of $Re$ and considering only the 1000 first time steps.

The algorithm provided an approximation with dimension $M = 7: X = (g_1(\xi), \dots, g_5(\xi), \xi)$ and TT-ranks $r = [1 \ 3 \ 5 \ 5 \ 5 \ 1 \ 1]$. The model error was estimated using two different samples:

* estimation of the approximation error on the sample $S_{test}$ with 3 simulations of $Re$ on the train time range (consisting of the 1000 first time steps): $e_{test}^2 = 2.97\%$,

* estimation of the extrapolation error on the sample $S_{extra}$ with 3 simulations of $Re$ on the whole time range (consisting of the 1500 available time steps): $e_{extra}^2 = 3.02\%$.

On Fig. 4 are plotted the predictions with blue crosses and the observations with circles, the red ones were used for learning the approximation and the green ones for estimating the
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 4** High Reynolds numbers. Observations and predictions with respect to time for fixed values of the other parameters. The predictions are obtained using tensor formats combined with constrained changes of variables</figcaption></figure>
model error. Again we observe a very good match of the predictions with the observations even beyond the train time range. The TT-rank of the approximation is low making possible the approximation using few samples of the Reynolds number $R_e$.

## Conclusion

This paper presents a new strategy to approximate multivariate functions with periodicity. It gives the principles of the method based on the combination of functions $h(g_1(\xi, t), \dots, g_m(\xi, t), \xi)$ chosen in appropriate classes of functions. The functions $g_i(\xi, t)$ define new variables of the multivariate functions $h$ which is here represented in the class of rank structured functions. Algorithms are proposed for constructing the approximation based on observations of the function, a constraint is added for the definition of the new parameters to promote periodicity of the representation. The numerical simulations yield good results. An analysis on the convergence of the approximation is to be studied.

## Acknowledgements

Authors' knowledge the support of the CNRS GdR 3587 AMORE.

## Authors' contributions

MC suggested the methodology, MS performed its implementation and run the numerical examples during her post-doctorate position at GeM. All authors contributed to the data analysis and discussed the content of the article. Both authors read and approved the final manuscript.

## Funding

Mathilde Chevreuil and Myriam Slama are grateful for the financial support provided by the French National Research Agency (ANR) and the Direction Générale de l'Armement (DGA) (MODUL'O π project, Grant ASTRID Number ANR-16-ASTR-0018).

## Declarations

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹GeM UMR CNRS 6183, Université de Nantes, Nantes, France, ²LOMC UMR CNRS 6294, Université Le Havre Normandie, Normandie Université, Le Havre, France.

Received: 29 December 2020 Accepted: 13 July 2021

Published online: 04 August 2021

## References

1. Bass RF, Grochenig K. Random sampling of multivariate trigonometric polynomials. SIAM J. Math. Anal. 2006;36(3):773-95. https://doi.org/10.1137/S0036141003432316.

2. Kammerer L, Potts D, Volkmer T. Approximation of multivariate periodic functions by trigonometric polynomials based on rank-1 lattice sampling. Journal of Complexity. 2015;31(4):543-76. https://doi.org/10.1016/j.jco.2015.02.004.

3. Briand T. Trigonometric polynomial interpolation of images. Image Processing On Line. 2019. https://doi.org/10.5201/ipol.2019.273.

4. Hastie T, Tibshirani R, Friedman J. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. : Springer; 2009.

5. Tipreddy R, Ghanem R. Basis adaptation in homogeneous chaos spaces. Journal of Computational Physics. 2014;259:304-17.

6. Constantine PG, Dow E, Wang Q. Active subspace methods in theory and practice: applications to kriging surfaces. SIAM Journal on Scientific Computing. 2014;36(4):1500-24.

7. Cortesi A, Constantine P, Magin T, Congedo PM. Forward and backward uncertainty quantification with active subspaces: application to hypersonic flows around a cylinder. INRIA: Technical report; 2017.

8. Friedman JH, Stuetzle W. Projection pursuit regression. Journal of the American statistical Association. 1981;76(376):817-23.

9. Hastie T, Tibshirani R, Wainwright M. Statistical Learning with Sparsity: the Lasso and Generalizations. : CRC Press; 2015.

10. Kolda TG, Bader BW. Tensor decompositions and applications. SIAM Review. 2009;51(3):455-500. https://doi.org/10.1137/07070111X.

11. Chevreuil M, Lebrun R, Nouy A, Rai P. A least-squares method for sparse low rank approximation of multivariate functions. SIAM/ASA Journal on Uncertainty Quantification. 2015;3(1):897-921. https://doi.org/10.1137/13091899X.

12. Grelier E, Nouy A, Chevreuil M. Learning with tree-based tensor formats; 2019. arXiv:1811.04455.
13. Cawley GC, Talbot NLC. Fast exact leave-one-out cross-validation of sparse least-squares support vector machines. Neural Networks. 2004;17:1467-75.

14. Tibshirani R. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. 1996;58(1):267-88.

15. Chen SS, Donoho DL, Saunders MA. Atomic decomposition by basis pursuit. SIAM Journal on Scientific Computing. 1999;20:33-61.

16. Efron B, Hastie T, Johnstone I, Tibshirani R. Least angle regression. The Annals of Statistics. 2004;32(2):407-99.

17. Nouy A, Grelier E, Giraldi L. Approximation toolbox. Zenodo. 2020. https://doi.org/10.5281/zenodo.3653970.

18. Grelier E. Learning with tree-based tensor formats: Application to uncertainty quantification in vibroacoustics. PhD thesis, Centrale Nantes; 2019.

19. Falcó A, Hackbusch W, Nouy A. Tree-based tensor formats. SeMA Journal; 2018.

20. Cichocki A, Lee N, Oseledets I, Phan A-H, Zhao Q, Mandic D. Tensor networks for dimensionality reduction and large-scale optimization: Part 1 low-rank tensor decompositions. Foundations and Trends in * Machine Learning. 2016;9(4-5):249-429.

21. Oseledets I, Tyrtyshnikov E. Recursive decomposition of multidimensional tensors. Doklady Math; 2009.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Submit your manuscript to a SpringerOpen® journal and benefit from:
▶ Convenient online submission
▶ Rigorous peer review
▶ Open access: articles freely available online
▶ High visibility within the field
▶ Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com