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

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.


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 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 + m n=1 (a n cos(nt) + b n sin(nt)) where the terms are periodic functions with period 2π. This class of functions has nice properties for approximation use, especially if a function f is continuous and given a tolerance ε, there exists a trigonometric polynomial v such that |f (t) − v(t)| < ε 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][2][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 , · · · , x s ) of the model, identified with an order-s tensor, is approximated in a parametrized subset of functions where the parameter a belongs to some set of parameters A and 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 , . . . , x s ), i = 1, . . . , 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 subspace [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 where is a multilinear map with parameters a l , l = 1, . . . , L. This is a special case of (1) with a = (a 1 , . . . , a L ) and A = R n 1 × . . . × R n L . The dimension of the parametrization L l=1 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•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 ∈ H be a multivariate function, with H a Hilbert space, which depends on a set of independent variables X = (X 1 , . . . , X s ). In the present paper, we consider the specific case where the function u is periodic with respect to one variable denoted τ so that we have X = ( , τ ), with = (X 1 , . . . , X d ) and τ = X d+1 , and s = d + 1.
A variable X ν has values in X ν and an associated measure dp ν , 1 ≤ ν ≤ s. The variable X has values in X = X 1 × · · · × X s and an associated measure dp = dp 1 × · · · × dp s . The variables X ν , ν = 1, . . . , d can be random variables, dp ν being in that case a probability measure on X ν . Let H be the Hilbert space defined on X , it is a tensor space H = H 1 ⊗ · · · ⊗ H s with H ν a Hilbert space defined on X ν . We consider that H ν ⊂ L 2 p ν (X ν ) is a finite dimensional subspace of square-integrable functions equipped with the norm u 2 ν = X ν u 2 dp ν and H is a subspace of L 2 p (X ) equipped with the canonical norm u 2 = X u 2 dp. Let {ψ ν i } P ν i=1 be an orthonormal basis of H ν and {ψ 1 i 1 ⊗ · · · ⊗ ψ s i s } (i 1 ,...,i s )∈[1,...,P 1] ×···×[1,...,P d ] be an orthonormal basis of H. A natural representation of the periodic function u can be obtained using the Fourier series. In the following, x = (ξ , t) will denote an observation of X with ξ an observation of = (X 1 , . . . , X d ) and t an observation of τ , 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 : R → R. It can be represented by its truncated Fourier series which is a sum of harmonic functions: a n cos(ω n t) + b n sin(ω n t) where the circular frequencies are such that ω n = nω f (multiplicative constraint on ω n ), with ω f = 2π T , and the coefficients a n and b n defined as follows: for n = 0, . . . , m − 1.
The truncated Fourier series can be seen as a composition of functions of the form: where h is an additive model h(z 1 , . . . , z m ) = m n=1 h n (z n ) with h n (z n ) = a n cos(z n ) + b n sin(z n ) ∈ V , where V = span 1, cos(·), sin(·) , and g n (t) = ω n−1 t for n = 1, . . . , 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 ξ , where the circular frequencies ω 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(ξ , t), periodic with respect to t. A natural representation of u is thus by means of a relevant change of variables: where g n : R d+1 → R, n = 1, . . . , m are new variables chosen under the form 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: where (x) is a multilinear map with parameters (a 1 , . . . , 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 rank and we define the subset of canonical tensors: It can be associated to the parametrized representation (7) with L = M, where a l ∈ R r×n l with n l the dimension of the functional basis on which the functions v l k are represented, k = 1, . . . , r and l = 1, . . . , M. We can consider other notions of ranks which provide different models with lower complexity. The α-rank of v, denoted by rank α (v), is the minimal integer r α such that where T is a collection of subsets of {1, . . . , M}. We define the subset of rank-structured functions: The complexity of the associated parametrized representations of tensors is linear with the dimension M and polynomial with the ranks.

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 over a model class M. The loss function 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 In practice, the approximation is constructed by minimizing the empirical risk 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ũ of u is then solution of with the regularization parameter λ ≥ 0, chosen or computed. The accuracy of the metamodelũ 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ũ(x) in a sequence of model classes with M h a linear or multilinear model class and g i ∈ M g a linear model class. We consider here multilinear models and more specifically the tensor subset T T r 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

Learning with tensor formats
Let z ∈ R M , an approximation of u in a tensor subset (2) can be obtained by minimizing the empirical least squares risk: where λ i R i (a i ) are regularization functions. Problem (13) (15) is a convex optimization problem known as Lasso [14] or basis pursuit [15]. The 1 -norm is a sparsity inducing regularization function that yields a solution 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 i α ) α∈T , the strategy consists in estimating and studying the truncation error min rank α (v)≤r i α R(v)−R(u) for different α in T , and choosing to increase the ranks r i α associated with the indices α 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 with w n = (w n,1 , . . . , w n,p ) ∈ R p the vector of parameters of the representation of g n on an orthonormal functional basis {ϕ j } p j=1 of H, and h a M-dimensional function in the model class of rank structured formats T C r or T T r that can be represented with a parametrized multilinear map with parameters a l , l = 1, . . . , L. The new set of variables z = (z 1 , . . . , z m , ξ ) is such that z n = g n (x), n = 1, . . . , 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 ∈ T T r and g ∈ H we have h • g ∈ L 2 p (X ). The approximationũ of the form (5) is thus parametrized as follows: with x k = (ξ k , t k ). The solution of this problem is found by alternatively solving for fixed (w 1 , . . . , w m ) using a learning algorithm with rank adaptation [11,12] and for fixed (a 1 , . . . , 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.
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 T T r where T is a dimension partition tree over {1, . . . , 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 T T r (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 ∈ T T r (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 ϕ such that g n = w n ϕ in (16) can be identified to (6) by choosing The circular frequencies in (6) are then expressed as: On the other hand, we choose bases of trigonometric functions {ψ n i } P n i=1 for the representation of h in the dimensions associated with the new variables z n , n = 1, . . . , 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: This constraint is imposed for all values taken by ξ in S. Using expression (23), it is recast under the form where w = (w 1 , . . . , w m ) ∈ R p×m , A ∈ R N ×p is the array of evaluations of ϕ and B ∈ R N ×m is a full array with values −2π/T max . The optimization problem (21) for the computation of the parameters w is replaced with the constrained optimization problem min 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, , τ ) with τ the time, Re the Reynolds number and the angular steps, that take values in X = X Re × X × X τ and we have d = 2. The pressure is evaluated on a tensor grid with The accuracy of the metamodelũ 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: where R S test (v) is the empirical risk defined in (10). Considering the least squares estimator,

Low Reynolds numbers
We first consider the approximation of the wall pressure for low Reynolds numbers for which the data is given 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 looses 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.

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 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 θ and Re, • for h: polynomial bases of degree 20 for θ 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 (ξ ), . . . , g 5 (ξ ), ξ ) and TT-ranks r = 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

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 (ξ , t), . . . , g m (ξ , t), ξ ) chosen in appropriate classes of functions. The functions g i (ξ , 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.