Approximation of high-dimensional periodic functions with Fourier-based methods
Daniel Potts Michael Schmischke
math.NA cs.NA 65T, 42B05


Robust On-Line ADP-based Solution of a Class of Hierarchical Nonlinear Differential Game
    Mohammad Reza Satouri, 
             Hamed Kebriaei, Senior Member, IEEE,
        Abolhassan Razminia 
        and Mohammad Javad Yazdanpanah
        
M. R. Satouri and H. Kebriaei are with the School of Electrical and Computer Engineering, College of Engineering, University of Tehran, P.O. Box 4395-515, Tehran, Iran (e-mail: satouri.mr@ut.ac.ir ; kebriaei@ut.ac.ir ).


A. Razminia is with the Dynamical Systems & Control (DSC) Research Lab.,
Department of Electrical Engineering, School of Engineering, Persian Gulf University, P.O. Box 75169, Boushehr, Iran (e-mail: razminia@pgu.ac.ir).  


M. J. Yazdanpanah is with the School of Electrical and Computer Engineering, College of Engineering, University of Tehran, P.O. Box 4395-515, Tehran, Iran, and also with the Control & Intelligent Processing Center of Excellence (e-mail: yazdan@ut.ac.ir).
     ^⋆  Laboratoire d'Informatique et des Systèmes, Université de Toulon, Aix Marseille Univ, CNRS UMR 7020, LIS, Marseille, France

boizot@univ-tln.fr; 

^†  Department of Mathematics, Lehigh University, Bethlehem, PA, USA. 

lus219@lehigh.edu; 

 

This research has been partially supported by the ANR SRGI (reference ANR-15-CE40-0018). 


This is a post-peer-review, pre-copyedit version of an article published in Journal of Mathematical Imaging and Vision (JIMIV). The final authenticated version is available online at: <http://dx.doi.org/10.1007/s10851-020-00990-5>

   
      
   

 
 

=======================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





	In this paper we propose an approximation method for high-dimensional 1-periodic functions based on the multivariate ANOVA decomposition. We provide an analysis on the classical ANOVA decomposition on the torus and prove some important properties such as the inheritance of smoothness for Sobolev type spaces and the weighted Wiener algebra. We exploit special kinds of sparsity in the ANOVA decomposition with the aim to approximate a function in a scattered data or black-box approximation scenario. This method allows us to simultaneously achieve an importance ranking on dimensions and dimension interactions which is referred to as attribute ranking in some applications. In scattered data approximation we rely on a special algorithm based on the non-equispaced fast Fourier transform (or NFFT) for fast multiplication with arising Fourier matrices. For black-box approximation we choose the well-known rank-1 lattices as sampling schemes and show properties of the appearing special lattices.




  ANOVA decomposition, high-dimensional approximation, Fourier approximation




	65T, 42B05




§ INTRODUCTION


The approximation of high-dimensional functions is an important and current topic with great interest in many applications. We consider a setting of periodic functions f, d ∈, over the torus  where certain data about the function is known. Here, we distinguish between a black-box setting, i.e., f can be evaluated at points ∈^d at a certain cost, and a scattered data setting, i.e., sampling points X ⊂^d and function values (f())_∈ X are given. Besides the natural question of wanting to find an approximation for f, we want to consider the question of interpretability, i.e., analyzing the importance of the dimensions and dimension interactions of the function. In applications this is sometimes referred to as an attribute ranking.

The main tool to achieve our goals is the analysis of variance (ANOVA) decomposition <cit.> which is an important model in the analysis of dimension interactions of multivariate, high-dimensional functions. It has proved useful in understanding the reason behind the success of certain quadrature methods for high-dimensional integration <cit.> and also infinite-dimensional integration <cit.>. The ANOVA decomposition decomposes a d-variate function in 2^d ANOVA terms where each term belongs to a subset of {1,2,…,d}. The single term depends only on the variables in the corresponding subset and the number of these variables is the order of the ANOVA term. In this paper we study the classical ANOVA decomposition for periodic functions and how it acts on the frequency domain. The decomposition is referred to as classical since it is based on an integral projection operator. In this setting we find relationships between ANOVA terms and the support of the frequencies as subsets of ^d. Moreover, we prove formulas for the representation of ANOVA terms and projections.     

Classical approximation methods cannot be applied for high-dimensional functions in general since the data required increases exponentially because of the curse of dimensionality. However, the observation has been made that in practical applications with multivariate, high-dimensional functions, often only the ANOVA terms of a low order are enough to describe a function, see e.g. <cit.>. This leads to the notion of a superposition dimension d_s ∈ that limits the order of the ANOVA terms involved. Using this as a sparsity assumption to circumvent the curse of dimensionality, we consider functions where the ANOVA decomposition is mostly supported on terms of a low order, i.e., the norm of the remaining decomposition weighed by the norm of f is small. This leads to a truncation of the decomposition through a superposition threshold. We consider how the previously described error can be related to the decay of Fourier coefficients and specifically the smoothness of f.

We present and analyze an approximation method that uses sensitivity analysis, cf. <cit.>, on the truncated ANOVA decomposition which is able to identify important ANOVA terms and incorporate this information in finding an approximation. The goal is to simplify the approximation model which yields benefits in reducing the influence of overfitting regarding the amount of data. We determine approximations of the Fourier coefficients of the function (or learn them) by solving least-squares problems. This is done through exploiting the special structure of the system matrix by identifying submatrices with the corresponding ANOVA terms. In the case of black-box approximations we are using rank-1 lattices as a spatial discretization, see e.g. <cit.>, and for scattered data approximation we rely on the iterative LSQR method <cit.> and the fast matrix-vector multiplications for Fourier matrices provided by the non-equispaced fast Fourier transform (or NFFT) introduced in <cit.>.

The paper is organized as follows. In <Ref> we introduce the classical ANOVA decomposition and study its behavior for periodic functions with regard to the Fourier system. We prove new formulas for the Fourier coefficients of projections in <ref> and ANOVA terms in <ref>. Moreover, we prove that functions in Sobolev type spaces and the weighted Wiener algebra inherit their smoothness to the ANOVA terms, see <ref> and <ref>. In <Ref> we consider the truncated ANOVA decomposition and prove formulas for their Fourier coefficients, see <Ref> and <ref>. We also give direct formulas for the truncated decomposition using the projections in <Ref> and <Ref>. Furthermore, we relate Sobolev type spaces and the weighted Wiener algebra to the previously introduced functions of low-dimensional structure and compute the errors in <ref> and <ref>. Specifically, we consider a class of product and order-dependent weights, see <cit.>, of functions with isotropic and dominating-mixed smoothness, cf. <cit.>, to obtain specific error bounds, see <ref> and <ref>. In <Ref> we present an approximation method for functions that are of an low-dimensional structure, cf. <ref>. We start by considering a black-box approximation scenario with rank-1 lattices as sampling schemes and show properties of the arising lattices in <ref>, <ref>, and <ref>. Furthermore, we discuss scattered data approximation in <ref>. The arising approximation errors are considered in <ref> with main results being <Ref>, <Ref>, <ref>, and <ref>. In <ref> we perform numerical experiments with a specific benchmark function.



§ PREREQUISITES AND NOTATION


We consider multivariate 1-periodic functions f^d with spatial dimension d ∈, which are square-integrable, i.e., elements of (). Those functions have a unique representation with regard to the Fourier system {^2πı·̨}_ k ∈^d as Fourier series

    f( x) = ∑_ k ∈^df ^2πıkx,

where f∫_ f( x) ^-2πıkxd x ∈, k ∈^d, are the Fourier coefficients of f. Given a finite index set I ⊂^d, we call the trigonometric polynomial

    S_I f () = ∑_∈̨If ^2πıkx

the Fourier partial sum of f with respect to the index set I.

In this paper we make use of indexing with sets. First, for a given spatial dimension d we denote with = { 1,2,…,d } the set of coordinate indices and subsets as bold small letters, e.g., ⊂̆. The complement of those subsets are always with respect to , i.e., = ∖$̆. For a vector∈^dwe define_ = ( x_i )_i ∈ ∈^. There remains a small ambiguity regarding the order of the components of_which can be clarified if one chooses a consistent ordering, e.g., ascending order which would be a natural choice.

Furthermore, we use thep-norm which is defined as

    p = { i ∈ x_i ≠ 0 }      p = 0 
    ( ∑_i=1^d x_i^p )^1/p    1 ≤ p < ∞
    max_ı∈x_i    p = ∞

for∈^d. Note that the case1 ≤p < ∞can be expanded to0 < p < 1, but then·pwould only be a quasi-norm. In the casep = 0,·pis not a norm at all.



 §.§ Rank-1 lattice


In the case of black-box approximation we are going to rely on rank-1 lattice as sampling schemes, see e.g. <cit.>. For a given lattice sizeM ∈and a generating vectorẕ ∈^dwe define a rank-1 lattice 
    Λ( ẕ, M ) {_j j/Mẕ1̱ j = 0,1,…,M-1 }.
 These lattices are useful in the evaluation of trigonometric polynomialsp ∈Π_I span{^2πı·̨∘∈̨I }over a finite index setI ⊂^dfor given Fourier coefficientsp. As discussed in <cit.>, we have 
    p(_j) = p(j/Mẕ1̱) = ∑_l=0^M-1( ∑_∈̨I 
    ·̨ẕ≡ l  Mp) ^2πıj l/M.
 The computation of the sum overlcan be realized trough a one-dimensional FFT and therefore the evaluation ofpat all lattice nodes can be done using only this single FFT. The arithmetic complexity of this evaluation is in𝒪(MlogM + d I).

However, using a special kind of rank-1 lattice, it is possible to reconstruct the Fourier coefficientspby samplingpat the nodes inΛ(ẕ,M)in an exact and stable way. For an index setI ⊂^dwe define a reconstructing rank-1 latticeΛ(ẕ,M,I)as a rank-1 latticeΛ(ẕ,M)such that the condition 
    m̱·ẕ≢0  M   ∀m̱∈(I)∖{0̱}
 is fulfilled with 
    (I) {-̨ẖ,̨ẖ∈ I }
 being the difference set forI. Using the nodes of a reconstructing rank-1 latticeΛ(ẕ,M,I), the Fourier coefficients can be calculated as 
    p = 1/M∑_j=0^M p(j/Mẕ1̱) ^-2πı j·̨ẕ/M.
 The calculation of all Fourier coefficientsp,∈̨I, can then be realized trough a one-dimensional FFT and the computation of the products·̨ẕ. Consequently, the arithmetic complexity of this evaluation is again in𝒪(MlogM + d I).

The principle of this reconstruction can be generalized to functionsf ∈(^d) { f ∈L_1() f(^d) ∑_k ∈^d w()̨ f < ∞} ,w^d[1,∞), by taking the Fourier partial sumS_I ffor a suitable index setI ⊂^dand treating the evaluations offas the evaluations of the trigonometric polynomialS_I f. Using the same idea as before, we get 
    f≈f̂_1/M∑_j=0^M f(j/Mẕ1̱) ^-2πı j·̨ẕ/M
 with a a reconstructing rank-1 latticeΛ(ẕ, M, I). The error for each coefficient is 
    f̂_ = f + ∑_∈Λ^(ẕ, M)∖{0̱}+̨f
 with the integer dual lattice 
    Λ^(ẕ, M) {∈̨^d ·̨ẕ≡ 0  M }.
 Consequently, if∑_∈Λ^(ẕ, M)∖{0̱} +̨ fis small, the approximationsf̂_are close to the Fourier coefficientsf. For further details on this topic we refer to <cit.> and <cit.>.



§ THE CLASSICAL ANOVA DECOMPOSITION OF 1-PERIODIC FUNCTIONS


In this section we introduce the ANOVA decomposition, see e.g. <cit.>, and derive new results for the periodic setting specifically with regard to the decomposition acting on the frequency domain.

We start by defining the projection operator

    P_ f ( x_ ) ∫_^ f( x) x_

that integrates over the variablesx_. For> 0this operator maps a function from()toL_2(^)by the Cauchy-Schwarz inequality and the imageP_ fdepends only on the variablesx_ ∈^. In the case of=̆ ∅, the projection gives us the integral off. We define the index set

    { k ∈^d k_ =  0 }

which can be identified with^using the mappingk ↦k_. Note that we use the convention^∅ = { 0 }. We now prove a relationship between the Fourier coefficients ofP_ fandf.

	Let f ∈(^d) and ℓ∈^. Then
	
    ℓ̱P_ f = f

	for k ∈^d with k_ = ℓ and k_ =  0.

	We consolidate the two integrals and derive 
	
    ℓ̱P_ f   = ∫_^∫_^ f( x) x_ ^-2πıℓx_ ux_
       = ∫_ f( x) ^-2πıℓx_ ux
       = ∫_ f( x) ^-2πı kxx = f.
 Using <ref>, we are able to writeP_ fas both, ad-dimensional Fourier seriesP_ f (x) = ∑_k ∈ f   ^2πıkxand a-dimensional Fourier seriesP_ f (x_) = ∑_ℓ∈^ ℓP_ f   ^2πıℓx_.

Now, we recursively define the ANOVA term for⊆̆f_P_ f - ∑_ v ⊊ u  f_.
 There exists a direct formula for the ANOVA termsf_defined in (<ref>). 


	Let a ∈_0 and b ∈ with b > a. Then
	
    ∑_n=a^b-1 (-1)^n-a+1b-an-a = (-1)^b-a.

	We prove an equivalent form obtained through multiplication with (-1)^a and an index shift
	
    ∑_n=0^b-a-1 (-1)^n+a+1b-an = (-1)^b.

	Splitting the sum and applying the Binomial theorem yields
	
    ∑_n=0^b-a-1 (-1)^n+a+1b-an   = ∑_n=0^b-a (-1)^n+a+1b-an - (-1)^b+1
       = (-1)^a+1∑_n=0^b-a(-1)^nb-an_=(-1+1)^b-a=0 + (-1)^b
       = (-1)^b.

	Let f ∈(^d) with ⊆̆. Then
	
    f_ = ∑_⊂̌ (-1)^-P_ f.

	A proof based on properties of projection operators was given in <cit.> while we use combinatorial arguments. We prove this statement through structural induction over the cardinality of u̱. For u̱ = 0, i.e., u̱ = ∅, we have
	
    (-1)^0-0(P_∅)(x̱) = (P_∅)(x̱) = (P_∅)(x̱) - ∑_v̱⊊∅ f_v̱(x̱).
 
	Now, let (<ref>) be true for ⊂̌, = 0,1,…,m-1, m ∈{1,2,…,d}, and take a subset u̱⊂ with u̱ = m. We use the notation
	
    δ_w̱⊂v̱ = 
    			1   w̱⊂v̱
    
    			0   otherwise.

	and start from the recursive expression in (<ref>) to obtain
	
    f_u̱(x̱)    = (P_u̱f)(x̱) - ∑_v̱⊊u̱ f_v̱(x̱) 
    		= (P_u̱f)(x̱) - ∑_v̱⊊u̱∑_w̱⊂v̱ (-1)^v̱-w̱(P_w̱f)(x̱) 
       = (P_u̱f)(x̱) - ∑_v̱⊊u̱∑_w̱⊊u̱ (-1)^v̱-w̱(P_w̱f)(x̱)δ_w̱⊂v̱.

	We exchange the two sums and sum over the order of the ANOVA terms 
	
    ∑_v̱⊊u̱∑_w̱⊊u̱ (-1)^v̱-w̱(P_w̱f)(x̱)δ_w̱⊂v̱   = ∑_w̱⊊u̱(P_w̱f)(x̱)∑_v̱⊊u̱ (-1)^v̱-w̱δ_w̱⊂v̱
       = ∑_w̱⊊u̱(P_w̱f)(x̱)∑_n=w̱^m-1∑_v̱⊂u̱
    v̱ = n (-1)^v̱-w̱δ_w̱⊂v̱
       = ∑_w̱⊊u̱(P_w̱f)(x̱)∑_n=w̱^m-1(-1)^n-w̱∑_v̱⊂u̱
    v̱ = nδ_w̱⊂v̱.

	Applying <ref> yields the formula.


We proceed to present a relationship between the Fourier coefficients off_andf. Furthermore, we provef_ ∈(^). Therefore, we define the index set

    { k ∈^d k_ =  0, k_j ≠ 0  ∀ j ∈ u }

which can be identified with( ∖{ 0 } )^using the mappingk ↦k_. Here, we use the convention( ∖{ 0 } )^∅ = { 0 }. 

	Let ,̆⊆̌ with ≠̆$̌. Then∩𝔽_^(d) = ∅. Moreover, we have
	
    ^d = ⋃_⊆̆.

	Let ,̆⊆̌, ≠̆$̌, and w.l.o.g.≥. We assume there exists ak̃∈∩𝔽_^(d)and first consider the case∩̆=̌∅. Sincek̃∈we havek̃_ =  0and thereforek̃_ =  0. This contradictsk̃∈𝔽_^(d). In the case of∩̆≠̌∅there exists aj ∈∩$̆. Then k̃∈𝔽_^(d) implies that k̃_j = 0 which contradicts k̃∈.
	
	The inclusion ⋃_⊆̆⊆^d is trivial since ⊆^d for every ⊆̆. To prove ^d ⊆⋃_⊆̆ we take a k ∈^d and define =̆{ j ∈ k_j ≠ 0 }. Then k ∈ and therefore k ∈⋃_⊆̆.

	Let f ∈(^d) with ⊆̆ and ℓ∈^. Then
	
    ℓf_ = f     ℓ̱∈(∖{0})^
    δ_,̆∅·0f     ℓ̱= 0̱
    
    	0    otherwise

	for k ∈^d with k_ = ℓ and k_ =  0. Furthermore, f_∈(^).

	We begin by employing the direct formula (<ref>) to obtain
	
    ℓf_   = ∫_^ f_(_)  ^-2πıℓ̱__
       = ∫_^[∑_⊂̌ (-1)^-P_ f (_)] ^-2πıℓ̱__
       = ∑_⊂̌ (-1)^-∫_^P_ f (_) ^-2πıℓ̱__
       = ∑_⊂̌ (-1)^-c__̨(P_f) δ__̨∖̆,0̱.

	We go on to prove 0f_ = δ_,̆∅·0f. In this case, _̨ = 0̱ and δ__̨∖̆,0̱ = 1 for every ⊂̌$̆. By the Binomial Theorem, we have
	
    ℓf_   = ∑_⊂̌ (-1)^-c__̨(P_f) δ__̨∖̆,0̱ = 0f∑_⊂̌ (-1)^-
       = 0f∑_n=0^n (-1)^-n = 0f·δ_,̆∅.

	For the second case, we consider anℓ̱and with a set⊂$̆ such that ∅≠{ i ∈k_i = 0 }≠$̆. Thenδ__̨∖̆,0̱ = 1 ⟺^c∖̆⊂$̌ and with the Binomial Theorem we get
	
    ℓf_   = ∑_⊂̌ (-1)^-c__̨(P_f) δ__̨∖̆,0̱ = ∑_^c⊂⊂̌ (-1)^-c__̨(P_f) 
       = f∑_^c⊂⊂̌ (-1)^- = f∑_n=^c^-^cn-^c (-1)^-n
       = f∑_m=0^-^c-^cm (-1)^-^c-m = 0.

	For the case were the entries of ℓ̱ are all nonzero, only the addend where =̌$̆ is nonzero, i.e.,ℓf_ = f.


With <ref> we have two equivalent series representations for the ANOVA termf_, thed-dimensional Fourier seriesf_ ( x) = ∑_ k ∈f ^2πıkxand the-dimensional Fourier seriesf_ (x_) = ∑_ℓ∈^ℓf_ ^2πıℓx_withℓf_as in <ref>. The ANOVA terms have the following important property.

	Let f ∈(^d) and ,̆⊆̌ with ≠̆$̌. Then the ANOVA termsf_andf_are orthogonal, i.e.,
	
    ⟨ f_, f_⟩_(^d) = 0.

	We employ <ref> and <ref> to deduce
	
    ⟨ f_, f_⟩_(^d)   = ⟨∑_ k ∈f ^2πıkx, ∑_ℓ∈ℓf ^2πıℓx⟩_(^d)
       =  ∑_ k ∈∑_ℓ∈fℓf δ_,̨ℓ = 0.


Having defined the ANOVA terms, we now go on to the ANOVA decomposition, cf. <cit.>.

	Let f ∈(^d), the ANOVA terms f_ as in (<ref>) and the set of coordinate indices = {1,2,…,d}. Then f can be uniquely decomposed as
	
    f(x̱) = f_∅ + ∑_i=1^d f_{i} (x_i) + ∑_i = 1^d-1∑_j = i+1^d f_{i,j} (x̱_{i,j})+ … + f_() = ∑_⊆̆ f_u̱(x̱_)

	which we call analysis of variance (ANOVA) decomposition.

	We use that ^d is the disjoint union of the sets  for ⊂̆ and obtain
	
    ∑_⊂̆ f_( _ )    = ∑_⊂̆∑_ k ∈kf ^2πıkx = ∑_ k ∈⋃_⊆̆kf ^2πıkx
       = ∑_ k ∈^dkf ^2πıkx = f(x̱).

	Since the union is disjoint, the decomposition is unique.

	The ANOVA decomposition (<ref>) depends strongly on the projection operator P_ f, see (<ref>). The integral operator considered in this paper leads to the so called classical ANOVA decomposition. Another important variant is the anchored decomposition where one chooses an anchor point c̱∈^d and the projection operator is then defined as 
	
    P_ f (_) = f(y̱),  y̱_ = _, y̱_ = c̱_.

	This decomposition can for example be used in methods for the integration of high-dimensional functions such as the multivariate decomposition method, see e.g. <cit.>. However, the error analysis may again be based on the classical ANOVA decomposition, see e.g. <cit.>. 


In <Ref> we have visualized the different frequency index sets,⊂̆, for a3-dimensional example.



 §.§ Variance and Sensitivity


In order to get a notion of the importance of single terms compared to the entire function, we define the variance of a function
    f∫_^d( f() - 0̱f)^2

for real-valuedf. In this case, we have the equivalent formulation

    f = f(^d)^2 - 0̱f^2

which yields a sensible definition for complex-valued functionsf. For the ANOVA termsf_with∅≠⊂̆we have0̱f_ = 0and therefore

    f_ = f_(^)^2.

	Let f ∈(^d). Then we obtain for the variance
	
    f = ∑_∅≠⊂̆f_.

	We show that the right-hand side equals the left-hand side by employing <ref> and <ref>
    ∑_∅≠⊂̆f_   = ∑_∅≠⊂̆∑_∈̨f^2 = ∑_∈̨⋃_∅≠⊆̆f^2 
       = ∑_∈̨^df^2 - 0̱f^2 = f(^d)^2 - 0̱f^2.


The global sensitivity indices 

    ff_/f∈ [0,1]

for∅≠⊂̆provide a comparable score to rank the importance of ANOVA terms against each other, cf. <cit.>. Clearly, we have∑_∅≠⊂̆f = 1by <ref>.

We now introduce one notion of effective dimensions as proposed in <cit.>. Given a fixedδ∈ (0,1], the general notion of superposition dimension is defined as the minimum 
    min{ s ∈∑_∅≠⊂̆
    ≤ sf_≥δf}.
 If we consider a particular Hilbert spaceH ⊂Ł_2(^d)with norm·H, we modify the superposition dimension in the sense of this space, see e.g. <cit.>. Forf ∈ Handδ∈ (0,1]we define the modified superposition dimension as 
    min{ s ∈sup_fH≤ 1∑_ > sf_Ł_2(^d)^2 ≤ 1 - δ}.


Finally, we investigate how the smoothness offtranslates to projectionsP_ fand ANOVA termsf_. For a different setting this has been discussed in <cit.> and therein called inheritance of smoothness. In our setting, we express smoothness through special subspaces of()and howfbeing an element of those spaces translates to the projectionsP_ fand ANOVA termsf_. In particular, we look at Sobolev type spaces, cf. <cit.>, 

    (^d) { f ∈() f(^d)(∑_ k ∈^d w^2()̨f^2)^1/2 < ∞}

and the weighted Wiener algebra
    (^d) { f ∈L_1() f(^d)∑_ k ∈^d w()̨f < ∞}

with a weight functionw^d[1,∞)for both cases. 


	Let f ∈(^d) with weight function w^d[1,∞). Then for any weight function 
w_^[1,∞) with
	
    w_(k_) ≤ w()̨ ∀∈̨

	we have P_ f ∈H^w_(^) and f_∈H^w_(^).

	We show that the norm P_ fH^w_(^) is finite by using <ref>
    ∑_ℓ∈^ w_^2(ℓ̱) ℓP_ f^2    = ∑_∈̨ w_^2(_̨) f^2 ≤∑_∈̨ w^2()̨f^2 
       ≤∑_∈̨^d w^2()̨f^2 = f(^d)^2 < ∞.
 
	Analogously, we employ <ref> to prove f_∈H^w_(^)
    ∑_ℓ∈^ w_^2(ℓ̱) ℓf_^2    = ∑_∈̨ w_^2(_̨) f^2 ≤∑_∈̨ w^2()̨f^2 
       ≤∑_∈̨^d w^2()̨f^2 = f(^d)^2 < ∞.

	Let f ∈(^d) with weight function w^d[1,∞). Then for any weight function w_^[1,∞) with
	
    w_(k_) ≤ w()̨ ∀∈̨

	we have P_ f ∈w_(^) and f_∈w_(^).

	We use <ref> to show that P_ f ∈(^)
    ∑_ℓ∈^ w_(ℓ̱) ℓf_   = ∑_∈̨ w_(_̨) kf≤∑_∈̨ w()̨kf
       ≤∑_∈̨^d w()̨kf = f(^d) < ∞.
 
	We utilize <ref> to prove f_∈w_(^)
    ∑_ℓ∈^ w_(ℓ̱) ℓf_   = ∑_∈̨ w_(_̨) kf≤∑_∈̨ w()̨kf
       ≤∑_∈̨^d w()̨kf = f(^d) < ∞.


The inheritance of smoothness has special significance with regard to the numerical realization of the method presented in <Ref>. It ensures that the ANOVA termsf_are at least as smooth as the functionfin consideration which is relevant for the quality of the approximation produced by the method.



§ TRUNCATED ANOVA DECOMPOSITION


The number of ANOVA terms of a function is equal to the cardinality of𝒫()=2^dand therefore grows exponentially in the dimension. This reflects the curse of dimensionality in a certain way and poses a problem for the approximation of a function. In this section we consider truncating the ANOVA decomposition, i.e., removing certain termsf_, and therefore creating a certain form of sparsity. We define a subset of ANOVA terms as a subset of the power set of, i.e.,U ⊂𝒫(D), such that the inclusion condition

    ∈̆U ⟹∀⊂̌∈̌U

holds, cf. <cit.>. This is necessary due to the recursive definition of the ANOVA terms, see (<ref>). 

For any subset of ANOVA termsUwe then define the truncated ANOVA decomposition as 

    T_U f ∑_∈̆U f_.

A specific truncation idea can be obtained by relating to the superposition dimension, see (<ref>). For a chosen superposition thresholdd_s ∈(that may or may not be equal to the superposition dimension), we defineU_d_s{⊂̆≤ d_s }andT_d_sT_U_d_s. We subsequently prove properties of bothT_Uin general andT_d_sin particular.

	Let f ∈(^d) and U ⊂𝒫(D) be a subset of ANOVA terms. Then T_U f∈(^d) and for ∈̨^d the Fourier coefficient is 
	
    T_U f = f      ∃∈̆U ∈̨
    
    	0    otherwise.

	Clearly, we have T_U f∈(^d). Let now ∈̨^d. Then there exists a _̆0 ⊂ such that ∈̨𝔽__̆0^(d). We employ <ref> and obtain
	
    T_U f   = ∫_^d( ∑_∈̆U f_(_) ) ^2πı·̨= ∑_∈̆U∫_^d f_(_) ^2πı·̨
       = ∑_∈̆U∫_^ f_(_) ^2πı_̨·__ δ__̨,0̱ = ∑_∈̆Uf δ__̨,0̱ (1-δ__̨,0̱) 
       = f      _̆0 ∈ U 
    
    	0    otherwise.

	Let f ∈(^d) and d_s ∈ a superposition threshold. Then T_d_sf∈(^d) and only the Fourier coefficients corresponding to d_s-sparse frequencies are nonzero, i.e.,  
	
    T_d_sf = f      _0 ≤ d_s 
    
    	0    otherwise.

	Since U_d_s is a subset of ANOVA terms, T_d_sf∈(^d) follows directly from <ref>. Moreover, ∃∈̆U_d_s∈̨⟺_0 ≤ d_s.


The following lemma shows that the number of terms inU_d_sis polynomial indfor a fixedd_sand therefore allows us to circumvent the curse of dimensionality in terms of the number of sets. 
	We estimate the cardinality of U_d_s as follows
	
    U_d_s < ( d/d_s)^d_s,

	i.e., the number of terms in U_d_s has polynomial growth in d for fixed d_s ∈∖{d}.

    We estimate the sum as follows
    
    U_d_s = ∑_n=0^d_sdn≤∑_n=0^d_sd^n d_s^n/n!  d_s^n = ∑_n=0^d_s(d/d_s)^n d_s^n/n!≤(d/d_s)^d_s∑_n=0^d_sd_s^n/n!.

    Estimating the sum by the Taylor series for ^d_s yields the statement.


In the following we show direct formulas for the truncated ANOVA decomposition based on the projections similarly as for the ANOVA terms, see (<ref>). 
	Let f ∈(^d) and U ⊂𝒫(D) a subset of ANOVA terms. Then we have the direct formula
	
    T_U f = ∑_∈̆U∑_∈̌U 
    ⊂̆ (-1)^-P_ f.

	We apply equation (<ref>) and obtain immediately
	
    T_U f    = ∑_∈̆U f_ = ∑_∈̆U∑_⊂̌ (-1)^-P_ f = ∑_∈̆U∑_∈̌U (-1)^-P_ f δ_⊂̌
       = ∑_∈̌U∑_∈̆U 
    ⊂̌ (-1)^-P_f.

	Let f ∈(^d) and d_s ∈ a superposition threshold. Then we have the direct formula
	
    T_d_s f = ∑_⊂̆
    ≤ d_s[ ∑_n=^d_s (-1)^n-d-n-] P_ f.

	Since the equality
	
    ∑_∈̌U_d_s
    ⊂̆ (-1)^- = ∑_n=^d_s (-1)^n-d-n-,

	holds, we employ <ref> and the formula is proven.	


The truncated ANOVA decomposition plays a major role in our approximation approach presented in <Ref>. Therefore we are interested in functions that can be approximated well by a truncated ANOVA decomposition. Specifically, we are looking to characterize functions such that the truncation operation byT_U ffor different setsUretains most of the function, i.e., we have a relative error 
    f-T_U fH_1/fH_2 < ε

withε > 0, andH_1, H_2certain subspaces of(^d). It is especially interesting to characterize these functions by properties like the smoothness. To this end, we start by proving general bounds for Sobolev type spaces(^d)and the weighted Wiener algebra(^d)to later relate this to weight functionswdefined by specific kinds of smoothness.
 

Moreover, this can be related to the superposition dimensionfor aδ∈ (0,1], see (<ref>). LetH_1 = Ł_2(^d)andH_2 ∈{w(^d), w(^d)}for a weight functionw. If we choose truncation by a superposition thresholdd_s ∈then the bound on the right-hand sideε(d_s) ∈ (0,1)depends ond_s. Moreover, we have 
    sup_f ≠ 0f-T_d_s fŁ_2(^d)^2/fH_2^2 = sup_fH_2≤ 1∑_ > d_sf_Ł_2(^d)^2 < ε(d_s)
 which follows fromf-T_d_s fŁ_2(^d)^2 = ∑_ > d_sf_Ł_2(^d)^2. The modified superposition dimensionwill now be smaller or equal tomin{ d_s ∈ε(d_s) ≤ 1 - δ}, i.e., truncation by this minimum as superposition threshold is guaranteed to be effective in relation toδ.


	Let f ∈(^d) with weight function w^d[1,∞). Then
	
    f-T_U f(^d)/f(^d)≤1/min_∈̨⋃_⊆̆
    ∉̆U w()̨ .

	We employ Parseval's identity and <ref> to derive
	
    f-T_U f(^d)^2    = ∑_∈̨^df - kT_U f^2 = ∑_∈̨⋃_⊆̆
    ∉̆Uf^2 
       =  ∑_∈̨⋃_⊆̆
    ∉̆Uw^2()̨/w^2()̨f^2 
       ≤1/min_∈̨⋃_⊆̆
    ∉̆U w^2()̨f(^d)^2.

	Let f ∈(^d) with weight function w^d[1,∞). Then
	
    f-T_U fL_∞(^d)/f(^d)≤1/min_∈̨⋃_⊆̆
    ∉̆U w()̨.

	For f ∈(^d) with a weight function w^d[1,∞) such that {1/w()̨}_∈̨^d∈ℓ_2 we have
	
    f-T_U fL_∞(^d)≤√(∑_∈̨⋃_⊆̆
    ∉̆U1/w^2()̨)  f(^d).

	We estimate the L_∞-norm by the sum of the absolute values of the Fourier coefficients and then use <ref>
    f-T_U fL_∞(^d)   ≤∑_∈̨^df - kT_U f = ∑_∈̨⋃_⊆̆
    ∉̆Uf
       =  ∑_∈̨⋃_⊆̆
    ∉̆Uw()̨/w()̨f
       ≤1/min_∈̨⋃_⊆̆
    ∉̆U w()̨f(^d).

	Employing the Cauchy-Schwarz inequality in (<ref>) instead of extracting the minimum yields 
	
    f-T_U fL_∞(^d)≤√(∑_∈̨⋃_⊆̆
    ∉̆U1/w^2()̨)  f(^d).

	The condition {1/w()̨}_∈̨^d∈ℓ_2 assures that the sum which appears in the bound is finite.


In the following, we relate the truncation offby the operatorT_d_swith the smoothness off. To this end, we introduce the weights

    w^α,β()̨γ_^-1  (1+1)^α∏_s ∈ (1+k_s)^β

with=̨{ i ∈ k_i ≠ 0 }and parametersβ≥ 0, andα > -β. The parametersα, β, and the weightγ_,⊂̆, regulate the decay of the Fourier coefficients. Specifically, the parameterαis regulating the isotropic smoothness andβthe dominating mixed smoothness, cf. <cit.>. Moreover,γ̱controls the influence of the different dimensions. We choose a POD (product and order-dependent) structure forγ_such that 
    γ_ = Γ_∏_s ∈γ_s,
whereΓ̱∈ (0,1]^dis nonincreasing andγ̱= (γ_i)_i=1^d∈ (0,1]^d. The POD structure is motivated by the application of quasi-Monte Carlo methods for PDEs with random coefficients, cf. <cit.>. Similar weights for isotropic and dominating mixed smoothness have been considered in <cit.>. Moreover, the Sobolev type spaces may also be referred to as weighted Korobov spaces, cf. <cit.> for product weights and <cit.> for general weights.

We now use the previously obtained bounds for general weight functionswand derive results for the weightsw^α,βfrom (<ref>). We focus on the subsets of ANOVA termsU_d_sdefined by a superposition thresholdd_s ∈.

	Let f ∈w^α,β(^d) with weight function from (<ref>) with POD structure (<ref>), β≥ 0, α > -β, Γ̱∈ (0,1]^d, and γ̱∈ (0,1]^d. Then
	
    f-T_d_s fL_∞(^d)/fw^α,β(^d)≤Γ_d_s+1  (2+d_s)^-α 2^-β (d_s+1)∏_s=1^d_s+1γ_s^∗

	where γ̱^∗ is the non-increasing rearrangement of γ̱.

	We use <ref> and calculate the bound for the weight function w^α,β,γ̱ by computing the minimum
	
    M min_∈̨^d 
    0 > d_sΓ__0^-1  (1+_1)^α∏_s=1^d (1+k_s)^β∏_s∈γ_s^-1.

	Since Γ̱ is non-increasing by definition, Γ_d_s+1^-1 has to be equal to the smallest value. The frequencies in  have exactly  nonzero entries, therefore we get
	
    M = Γ_d_s+1^-1 (1+d_s+1)^α (1+1)^β (d_s+1)min_∈̨^d 
    0 > d_s∏_s∈γ_s^-1.

	The remaining product becomes minimal for the product of the d_s+1 smallest entries in γ̱ which yields the statement.

	Let n ∈ and γ̱∈ (0,1]^d. Then
	
    ∑_⊆̆
     = n∏_s∈γ_s^2 ≤γ̱_2^2n.

	We rewrite the sum as follows
	
    ∑_⊆̆
     = n∏_s∈γ_s^2 = ∑_i_1=1^dγ_i_1^2 ∑_i_2=i_1+1^dγ_i_2^2 ⋯∑_i_n=i_n-1+1^dγ_i_n^2.

	Then every single sum can be estimated by γ̱_2^2, i.e., ∑_i_j = i_j-1+1^d γ_i_j^2 ≤∑_i_j = 1^d γ_i_j^2 = γ̱_2^2 for j ∈{2,3,…,d} with equality for j = 1.

	Let f ∈H^w^α,β(^d) with weight function from (<ref>) with POD structure (<ref>), β≥ 0, α > -β, Γ̱∈ (0,1]^d, and γ̱∈ (0,1]^d. Then
	
    f-T_d_s f(^d)/fH^w^α,β(^d)≤Γ_d_s+1 (2+d_s)^-α 2^-β (d_s+1)∏_s=1^d_s+1γ_s^∗

	where γ̱^∗ = (γ_s^∗)_s=1^d is the non-increasing rearrangement of γ̱. For functions with isotropic smoothness α = 0 and dominating mixed smoothness β > 1/2 we have
	
    f-T_d_s fL_∞(^d)/fH^w^α,β(^d)≤√(∑_n=d_s+1^d 2^nΓ_n^2 (ζ(2β)-1)^nγ̱_2^2n)
	 
	where ζ is the Riemann zeta function. Exponential decay for Γ_s, i.e., Γ_s = c^s, 0 < c ≤ 1, such that the condition 
	
    γ̱_2 < 1/c√(2ζ(2β)-2)

	holds, yields the bound
	
    f-T_d_s fL_∞(^d)/fH^w^α,β(^d)≤( cγ̱_2√(2ζ(2β)-2))^d_s+1/√(1-2 c^2γ̱_2^2(ζ(2β)-1) ).

	The bound from statement (<ref>) is a consequence of <ref> and can be calculated analogously to the proof of <ref>. For the second statement, we calculate the constant in the bound from <Ref>. We use <Ref> and the product structure of the weights w^α,β()̨ to obtain
	
    ∑_∈̨⋃_⊆̆
     > d_s1/w^2()̨   = ∑_⊆̆
     > d_s∑_∈̨1/Γ_^-2(1+ k_s)^2β∏_s∈γ_s^-2
       = ∑_⊆̆
     > d_sΓ_^2∑_∈̨(∖{0})^1/(∏_s∈γ_s^-2)(∏_s=1^ (1+k_s)^2β)
       = ∑_⊆̆
     > d_sΓ_^2∏_s∈γ_s^2 ∑_k∈∖{0}1/ (1+k)^2β.

	We find an explicit form by replacing the sums with the Riemann zeta function
	
    ∏_s∈γ_s^2  ∑_k∈∖{0}1/ (1+k)^2β = ∏_s∈ 2  γ_s^2  ∑_k∈1/ (1+k)^2β = 2^(ζ(2β)-1)^∏_s∈γ_s^2.

	Applying <ref> then gives us the upper bound
	
    ∑_n=d_s+1^d 2^nΓ_n^2(ζ(2β)-1)^n∑_⊆̆
     = n∏_s∈γ_s^2 ≤∑_n=d_s+1^d 2^nΓ_n^2(ζ(2β)-1)^nγ̱_2^2n.

	If we choose an exponential decay for Γ_n, i.e., Γ_n c^n, 0 < c ≤ 1, the explicit upper bound becomes
	
    ∑_n=d_s+1^d 2^n c^2n(ζ(2β)-1)^nγ̱_2^2n = q^d_s+1/1-q (1-q^d-d_s)

	where q  2 c^2(ζ(2β)-1) γ̱_2^2 with 0 < q < 1 because of the condition (<ref>).


The bound in <Ref> and (<ref>) in <Ref> are independent of the spatial dimensionsdof the functionsfas long as they have the same superposition threshold and the norm stays the same. This allows us to circumvent the curse of dimensionality here and use the ANOVA terms inU_d_sfor a superposition thresholdd_s ∈. The bound (<ref>) can also be considered ford →∞. The dependence on the dimensiondis contained within the normγ̱_2^2. Choosing a square-summable sequence{γ_ℓ}_ℓ∈results in an upper bound forγ̱_2for anyd →∞. In this case the bound can be made independent ofdby the condition (<ref>). 

<Ref> shows the different bounds for weightsw^α,βwithγ̱= (1/s)_s=1^9andΓ̱= (π^-s√(3)^s )_s=1^9, see (<ref>). With regard to the superposition dimensionforw^α,β(^d), cf. (<ref>), one may interpret this as follows: Givenf ∈w^α,β(^d), the valueε(α,β) ∈ (0,1)of the bound in part (a) of <Ref> tells us that forδ = 1 - ε(α,β)^2the superposition dimensionis smaller or equal to the superposition thresholdd_s = 3, e.g.,ε(0,1) ≈ 0.0008and thereforeδ = 0.99999936.



§ ANOVA APPROXIMATION METHOD


We consider the general problem of approximating a periodic functionf^dgiven certain function evaluations off. Specifically, we distinguish two approximation scenarios – black-box approximation and scattered data approximation. In the case of black-box approximation, we are able to evaluatefat any given point∈^d. Since the evaluations come at a certain cost, we aim to keep them minimal or require a certain trade-off. For scattered data approximation we have a finite set of nodesX ⊂^dand know the function valuesy̱ = (f())_∈ X. Here, one cannot add more nodes toXor choose the locations of the nodes. Both scenarios have a high relevance for problems in various applications. 

In this section, we consider an approximation scheme for high-dimensional, periodic functions of a low-dimensional structure, i.e., functions with a small superposition dimension∈for aδ∈ (0,1]that is close to one, cf. (<ref>). In this case the truncation byT_d_swith a small superposition thresholdd_s ∈will be effective. It has been observed that functions in many practical applications belong to such a class, see e.g. <cit.>. In <Ref> we have considered errors for functions of dominating-mixed and isotropic smoothness defined trough the decay of the Fourier coefficients and therefore obtained an upper bound for the modified superposition dimensionfrom (<ref>). Considering <ref>, we know that e.g. POD weights lead to a decay such that the functions are of a low-dimensional structure.

The approximation scheme can be viewed in both approximation scenarios although the details are different. We work for now with the node setXas well as function evaluationsy̱and keep in mind thatXmay also be chosen if we are in the black-box case. The first step is to reduce the ANOVA decomposition to the terms inU_d_s, i.e., we approximate 
    f ≈ f = ∑_∈̆U_d_s f_.
 The Fourier coefficients fcan only be nonzero if the frequency$̨ is at most d_s-sparse, i.e., 0≤ d_s, see <ref>. Based on this, we aim to approximate f by a Fourier partial sum S_I f with a finite index set 
    I ⊂{∈̨^d ∈̨U_d_s}.
 The challenge is to determine an appropriate index set I. To this end, we employ a special scheme to determine frequency locations based on the ANOVA terms and an importance ranking on them.

We call the first step active set detection and its aim is to determine an importance ranking on the terms f_ with ∈̆U_d_s based on the global sensitivity indices f, cf. (<ref>). This information is also highly relevant to interpret relations in our data X and . 

Based on the sensitivity indices we build an active set of ANOVA terms U ⊂ U_d_s. This relates to the importance of frequencies and therefore information on how to choose the index set I from (<ref>). Reducing the number of ANOVA terms and in turn the number of frequencies leads to a reduction of the model complexity. The effects of overfitting are therefore lessened. In <Ref> we consider the details of the active set detection and in <Ref> the approximation with an active set as well as approximation errors.



 §.§ Active set detection


The method assumes that the underlying function f is of a low-dimensional structure, i.e., f ≈ f for some superposition threshold d_s ∈. The goal in the active set detection step is to determine an importance ranking for the ANOVA terms. In order to do this, we choose an appropriate search index set. Since we have no a-priori knowledge about the importance of the ANOVA terms or the smoothness of the function f, we work with order-dependent finite index sets I_0 = { 0 }, I_1 ⊂ (∖{0}), …, I_d_s⊂ (∖{0})^d_s. This achieves that two ANOVA terms f_ and f_ with = are supported on equivalent index sets. We then use the projection operator 
    P_ I {∈̨^d _̨∈ I, _̨ = 0̱}
 to project the index sets and obtain 
    I(U_d_s) = ⋃_∈̆U_d_sP_I_.
 This leads to the approximation by a Fourier partial sum 
    f() ≈ f () ≈ S_I(U_d_s) f () = I(U_d_s)f.


The Fourier coefficients f in (<ref>) are unknown and we aim to determine approximations for them from the data X and y̱. To this end, we consider the least-squares problem 
    _sol = _∈^I(U_d_s) - _I(U_d_s)2^2

with Fourier matrix _I(U_d_s) = ( ^2πı·̨)_∈ X,∈̨I(U_d_s). If the Fourier matrix has full rank, the elements of the solution vector _sol = (f̂_)_∈̨I(U_d_s) are the unique least-squares approximation to the Fourier coefficients, i.e., f̂_≈f, with respect to X and . Depending on the approximation scenario, there are different methods of solving least-squares problems of the type (<ref>). We refer to <Ref> for details.

We use the approximate Fourier coefficients f̂_ to build the approximate Fourier partial sum 
    S_I(U_d_s) f () ≈ S_I(U_d_s)^X f () = I(U_d_s)f̂_

which provides an initial approximation to the function f. In order to achieve a Fourier matrix _I(U_d_s) with full rank and combat the effects of overfitting, we may need to severely limit the number of frequencies in the order-dependent sets I_1, I_2, …, I_d_s. Details on this will be considered in the following subsections for the specific approximation scenarios.

In order to determine an importance ranking on the ANOVA terms, we assume that the global sensitivity indices of S_I(U_d_s)^X f and f behave similarly, i.e., it holds that 
    _̆1S_I(U_d_s)^X f≤_̆2S_I(U_d_s)^X f⟹_̆1f≤_̆2f
 for _̆1, _̆2 ∈ U_d_s. This allows us to use a threshold vector ε̱∈ [0,1]^d_s to define an active set of ANOVA terms that only contains the important terms with respect to ε̱
    U_X,^(){⊂̌∃∈̆U_d_s⊂̌and S_I(U_d_s)^X f > ε_}.

The inclusion condition (<ref>) is fulfilled by definition. We reduce the ANOVA decomposition to this set of terms to determine an approximation for f in <Ref>.



 §.§ Least-squares approximation


In this section, we discuss the solution of least-squares problems of the form 
    min_∈^I(U) - _I(U)2^2
 
with a Fourier matrix _I(U) = ( ^2πı·̨)_∈ X,∈̨I(U). Here, U is an arbitrary subset of ANOVA terms and for each term we have a given finite frequency index set I_⊂ (∖{0})^. The set 
    I(U) = ⋃_∈̆UP_I_
 is obtained through the projections (<ref>). 

The following remark shows that the Fourier matrix can be structured with respect to the ANOVA terms. Moreover, we can decompose the matrix-vector multiplications with both, _I(U) and its adjoint _I(U)^∗. 
	Let _I(U) be a Fourier matrix with respect to a node set X and an index set I(U) with a subset of ANOVA terms U ⊂𝒫(D) and index sets I_⊂ (∖{0})^, ∈̆U. Then 
    F̱ = ( F̱__̆1   F̱__̆2   ⋯   F̱__̆n)
 
	where _̆1, _̆2, …, _̆n with n = U is a numbering of the subsets of coordinate indices in U such that f̱̂̂ = ( f̱̂̂__̆1   f̱̂̂__̆2   ⋯   f̱̂̂__̆n)^⊤. The Fourier matrices are F̱_ = ( ^2πıℓ̱·_)_∈ X,ℓ̱∈ I_. The matrix-vector product with F̱ can therefore be decomposed as
	
    F̱ = ∑_∈̆UF̱__
 with vector components _. For the adjoint product F̱^∗f̱ with a vector f̱∈^X we obtain the result â̱̂∈^I(U) by computing the products 
    â̱̂_ = F̱_^∗f̱,  ∀∈̆U.
 
	Then we have the result vector â̱̂ = ( â̱̂__̆1   â̱̂__̆2   ⋯   â̱̂__̆n)^⊤.


  §.§.§ Black-box scenario
 

In the case of black-box approximation, i.e., the set X can be chosen, we have to determine an appropriate special discretization for index sets of the type I(U). Here, we have different possibilities. One might think of rank-1 lattices that have been used for integration before, see e.g. <cit.>, and approximation, see e.g. <cit.>. For a general introduction to lattice rules, we refer to <Ref>. Sparse grid sampling related to the Smolyak algorithm is a further possibility, cf. <cit.>.

In the following, we focus on using reconstructing single rank-1 lattice for function approximation. If we have a reconstructing single rank-1 lattice Λ(ẕ, M, I(U)) ⊂^d for a generating vector ẕ∈^d and size M ∈ with respect to an index set I(U), then 
    _I(U)^∗_I(U) = M ·I
 with I the identity matrix, see <cit.>. Then the solution to problem (<ref>) is unique and given by the multiplication of the Moore-Penrose inverse _I(U)^† with y̱, see e.g. <cit.>. Through the property (<ref>) the Moore-Penrose inverse is simplified to 
    _I(U)^† = 1/M_I(U)^∗,
 i.e., a multiplication with the adjoint matrix. This allows us to efficiently compute approximations for the Fourier coefficients of f if the nodes form a reconstructing rank-1 lattice.

It remains the issue of determining such a reconstructing rank-1 lattice given an index set of type I(U). In <cit.> it was shown that reconstructing lattices exist if the lattice size M is sufficiently large. Since the evaluations of f come at a certain cost, it is necessary to consider the lattice size for our special types of index sets which we do in the following.

An important quantity to get estimations on the lattice size is the difference set (I(U)) from (<ref>) since the result <cit.> tells us that there exists a reconstructing rank-1 lattice with prime cardinality 
    I(U)≤ M ≤(I(U)).
 In the following, we proof properties and show estimates on the cardinality of both I(U) and (I(U)).


	Let U ⊂𝒫() be a subset of ANOVA terms and I_⊂^, ∈̆U, finite symmetric frequency sets. Then we have 
    (I(U)) = ⋃_∈̆U 
    ⊂̌{-̨ẖ∈̨P_I_, ẖ∈P_I_}.

	It is easy to see that ⋃_∈̆U 
⊂̌{-̨ẖ∈̨P_I_, ẖ∈P_I_}⊂(I(U)) since P_I_⊂ I(U) for every ∈̆U and ∈̌U for all ⊂̌∈̆U due to (<ref>). In order to show the other inclusion we take an element ℓ̱∈(I(U)). By the uniqueness property of the ANOVA decomposition we know that there exists ,̆∈̌U such that ℓ̱= -̨ẖ with ḵ∈P_I_ and ẖ∈P_I_. Taking the symmetry of the index sets I_ into account, we have proven the statement.

The following lemma gives an estimate for the size of the difference set of index sets of type I(U) if there exists an upper bound on the cardinality of the term dependent sets I_.

	Let U be a subset of ANOVA terms and I_⊂ (∖{0})^, ∈̆U, symmetric frequency sets. Then the cardinality of the difference set of I(U) is bounded by 
    (I(U))≤∑_∈̆U∑_⊂̌I_I_≤ 2^max_∈̆UUmax_∈̆UI_^2.

	We estimate the cardinality of the difference set by applying <Ref>(I(U))≤∑_∈̆U∑_⊂̌I_I_. Here, we do not have equality since the union in <Ref> is not necessarily disjoint. Applying the upper bound on the cardinality of the sets I_, we arrive at 
    ∑_∈̆U∑_⊂̌I_I_   ≤∑_∈̆U∑_⊂̌max_∈̆UI_^2 ≤max_∈̆UI_^2   2^max_∈̆U∑_∈̆U 1 
       ≤ 2^max_∈̆UUmax_∈̆UI_^2.

    The cardinality of U_d_s is bounded by (· d/d_s)^d_s, see <ref>. Therefore the estimate in (<ref>) becomes 
    (I(U_d_s))≤(2· d/d_s)^d_smax_∈̆UI_^2.


In the following, we consider special term-dependent frequency index sets of the structure 
    I_{ℓ̱∈ (∖{0})^ w()̨≤ N_ for ∈̨^d  with _̨ = ℓ̱, _̨ = 0̱}
 with a subset of coordinate indices ∅≠⊂̆, a weight function w^d[1,∞) and cut-off N_∈. For a given subset of ANOVA terms U ⊂𝒫() we estimate the cardinalities of both, I(U) and the difference set (I(U)).


	Let U ⊂𝒫() be a subset of ANOVA terms, I_∅ = {0̱}, and I_, ∅≠∈̆U, finite frequency sets as in (<ref>) for a weight function w^d[1,∞) and N_∈. Moreover, let h_min[1,∞) and h_max[1,∞) be functions such that 
    c h_min(N__̆min) ≤min_∈̆U∖{∅}I_ and max_∈̆UI_≤ C  h_max(N__̆max)
 with _̆min = _∈̆U∖{∅}I_, _̆max = _∈̆UI_, and 0 < c ≤ C. Then we have for the asymptotic behavior of the cardinality of I(U)
    c h_min(N__̆min) ≤I(U)/U≤ C  h_max(N__̆max).
 The constants do not depend on the spatial dimension d.

	Since the projected sets P_ I_, ∈̆U, are disjoint, we have 
    I(U) = ∑_∈̆UI_.
 In order to show the upper bound, we estimate the cardinality of each index set by h_max
    ∑_∈̆UI_≤∑_∈̆U C  h_max(N__̆max) ≤ C  h_max(N__̆max) ∑_∈̆U 1 = U  C  h_max(N__̆max).
 The lower bound follows with similar arguments.

    Let U ⊂𝒫() be a subset of ANOVA terms, I_∅ = {0̱}, and I_, ∅≠∈̆U, finite symmetric frequency sets as in (<ref>) for a weight function w^d[1,∞) and N_∈. Moreover, let h_max be a function as in <Ref>. Then
    
    (I(U))≤ C^2   2^max_∈̆UU h_max^2(N__̆max).

    The corollary is a direct consequence of <ref> and <ref>.


We may apply <cit.> to construct the reconstructing rank-1 lattice Λ(ẕ, M, I(U)) via a component-by-component approach. 
Choosing the set X = Λ(ẕ, M, I(U)) as sampling nodes yields a Moore-Penrose inverse of type (<ref>) and we are able to compute the solution to (<ref>) by multiplying with the adjoint Fourier matrix. This computation can be done efficiently using a lattice fast Fourier transform or LFFT, see <cit.>.



  §.§.§ Scattered data scenario


In this section, we consider the scenario of scattered data approximation, i.e., we have a fixed set of nodes X ⊂^d. Here, we aim to solve the least-squares problem (<ref>) with the iterative LSQR method <cit.>. Specifically, we are interested in the matrix-free variant, i.e., we do not have to construct the system matrix _I(U)∈^X,I(U) explicitly. The curse of dimensionality would quickly lead to the size of the matrix becoming intractable. The matrix-free variant requires two algorithms, one which takes a vector a̱∈^I(U) as an input and returns the result of the matrix-vector multiplication _I(U)a̱ and one that takes â̱̂∈^X as an input and returns the result of ^∗_I(U)â̱̂. If we take <ref> into account, it is only necessary to provide algorithms for fast multiplication with Fourier matrices _I_∈^X,I_, ∈̆U.

The existence of such algorithms depends on the choice of the specific index sets I_. For full grids, i.e., frequency sets of the type  I_ = G_N^ = {∈̨^: -N_/2≤ k_i ≤N_/2 - 1, i = 1,2,…,}, N_∈ 2,  the non-equispaced fast Fourier transform (NFFT) was introduced in <cit.>. Moreover, for hyperbolic cross index sets of the form 
    I_ = H_n^ = ⋃_j̱∈_0^
    j̱1 = nĜ_j̱
 with Ĝ_ṉ = ×_s=1^Ĝ_n_s and Ĝ_n_s = (-2^n_s-1, 2^n_s-1]^∩, we have the non-equispaced hyperbolic cross fast Fourier transform (NHCFFT), cf. <cit.>.



 §.§ Approximation with active set


Now that we have obtained the active set U_X,^() from (<ref>), we aim to construct an approximation using only these ANOVA terms. The global sensitivity indices S_I(U_d_s)^X f calculated from the approximation S_I(U_d_s)^X f in (<ref>) provide us with a basis to choose term-dependent frequency index sets I_⊂ (∖{0})^, ∅≠∈̆U_X,^(). A higher sensitivity index suggests that the term is more important to the function and therefore a larger corresponding index set could be advisable. 

We project the index sets as before to obtain I(U_X,^()), see (<ref>). Note that in general and depending on the threshold ϵ̱, we have reduced the number of frequencies significantly. This is a sensible measure to reduce the effects of overfitting. Now, we approximate f by the Fourier partial sum 
    f() ≈T_U_X,^() f () ≈ S_I(U_X,^()) f () = I(U_X,^())f.
 The Fourier coefficients f are again unknown and we determine them by least-squares approximation from X and . The unique solution is given by 
    _sol = _∈^I(U_X,^()) - _I(U_X,^())2^2
 if the Fourier matrix _I(U_X,^()) = ( ^2πı·̨)_∈ X,∈̨I(U_X,^()) has full rank. Details on how to solve this system for scattered data and black-box approximation can be found in <ref>. We use the elements of the solutions vector _sol = (f̂_)_∈̨I(U_X,^()) to form the approximate Fourier partial sum and our solution 
    f() ≈ S_I(U_X,^()) f () ≈ S_I(U_X,^())^X f () = I(U_X,^())f̂_.
 The following algorithm summarizes the proposed method.



§ ERROR ANALYSIS


The error of our approximation method measured in the norm of some space H⊂Ł_2(^d) can be decomposed into multiple components by the triangle inequality 
    f - S_I(U)^X fH≤f - T_U fH_ANOVA truncation error+ T_U f - S_I(U)^X fH_approximation error
 for an active set of ANOVA terms U ⊂ U_d_s with superposition threshold d_s ∈. We distinguish between the ANOVA truncation error and the approximation error. Here, the analysis of the ANOVA truncation error is independent of the concrete approximation problem (<ref>) and the scenario (scattered data or black-box).



 §.§ ANOVA truncation error


The ANOVA truncation error is related to the truncation of the ANOVA decomposition to the set U_d_s with superposition threshold d_s ∈ and the active set U ⊂ U_d_s. We can separate the ANOVA truncation error as follows 
    f - T_U fH≤f - T_d_s fH_truncation by d_s + T_d_s f - T_U fH_active set truncation.
 Here, we bring T_d_s in with the aim to relate the error to our function class of low-order interactions, see (<ref>). To control the second term, we require assumptions on the sensitivity indices of the ANOVA terms in U_d_s∖ U. Since the error is only related to the structure of the function it can be considered independently of any specific approximation scenario like black-box or scattered data approximation. We show bounds for this error in the case that f is an element of a Sobolev type space (^d) or a Wiener algebra (^d) and H is (^d) or L_∞(^d).


    Let f ∈(^d) with a weight function w^d[1,∞) and superposition dimension , see (<ref>), for a δ∈ (0,1). If there exists a subset of ANOVA terms U ⊂ U_ such that f = f_/f < ε,  ε > 0, for every ∈̆U_∖ U then 
    f - T_U f(^d)/f(^d)≤√(1-δ) + √(U_∖ Uε).

    The ANOVA truncation error can be separated as in (<ref>). We prove an upper bound for the active set truncation. With Parseval's equality and the assumption on the global sensitivity indices, we estimate 
    T_ f - T_U f(^d)^2 = ∑_∈̆U_∖ U∑_∈̨f^2 ≤fU_∖ Uε.

    Clearly, we have f≤f(^d)^2 ≤f(^d)^2.

    Let f ∈(^d) with a weight function w^d[1,∞). If there exsists a subset of ANOVA terms U ⊂ U_d_s, d_s ∈, such that 
    ∑_∈̨f/∑_∈̨^df< ε_1, ε_1 >0
 for every ∈̆U_d_s∖ U and we have f-T_d_s fL_∞(^d)/f(^d) < ε_2,  ε_2 > 0,    then 
    f - T_U fL_∞(^d)/f(^d)≤ε_2 + √(U_d_s∖ Uε_1).

    We split the ANOVA truncation error as in (<ref>) and prove an upper bound for the second part. To this end, we estimate the L_∞ norm of f by the absolute values of its Fourier coefficients and apply (<ref>) to obtain
    
    T_d_s f - T_U fL_∞(^d)≤∑_∈̆U_d_s∖ U∑_∈̨f≤U_d_s∖ Uε_1 ∑_∈̨^df.

    Naturally, it holds that ∑_∈̨^df≤f(^d) which leads to the desired estimate.

Note that in order to prove a bound for the error in L_∞, we formulated a condition on an ℓ_1 equivalent of the global sensitivity indices f in accordance with the Wiener algebra norm.



 §.§ Approximation error


In this section, we focus on the approximation error which we separate into two parts as well 
    T_U f - S_I(U)^X fH≤T_U f - S_I(U) fH_truncation error + S_I(U) f - S_I(U)^X fH_aliasing error
 with H ∈{(^d), Ł_∞(^d)}, a subset of ANOVA terms U ⊂𝒫(), and a finite frequency index set I(U) ⊂^d of structure (<ref>) with sets I_ as in (<ref>). The truncation error remains independent of the approximation scenario and can be estimated by the norms in  and . 
     Let f ∈(^d), w^d[1,∞) a weight function, and I(U) ⊂^d a finite frequency index set of type (<ref>) with U ⊂𝒫(). Then the relative truncation error can be estimated as 
    T_U f - S_I(U) f(^d)/f(^d)≤1/min_∈̆U N_.

    If in addition we have ∑_∈̨^d1/w^2()̨ < ∞, we can estimate 
    T_U f - S_I(U) fŁ_∞(^d)/f(^d)≤√(∑_∈̆U∑_∈̨∖P_I_1/w^2()̨).

    In order to prove (<ref>) we employ Parseval's identity and use the weight w()̨
    T_U f - S_I(U) f(^d)^2    = ∑_∈̆U∑_∈̨∖P_I_f^2 = ∑_∈̆U∑_∈̨∖P_I_w^2()̨/w^2()̨f^2 
       ≤∑_∈̆U1/N_^2∑_∈̨∖P_I_ w^2()̨f^2 ≤1/min_∈̆U N_^2f(^d)^2.

    For the bound (<ref>) we estimate the norm by the absolute sum of the Fourier coefficients and use the Cauchy-Schwarz inequality 
    T_U f - S_I(U) fŁ_∞(^d)   = ∑_∈̨⋃_∈̆U∖ I(U)f = ∑_∈̨⋃_∈̆U∖ I(U)w()̨/w()̨f
       ≤f(^d)√(∑_∈̆U∑_∈̨∖P_I_1/w^2()̨).

     Let f ∈(^d) with w^d[1,∞) a weight function such that ∑_∈̨^d1/w^2()̨ < ∞, and I(U) ⊂^d a finite frequency index set of type (<ref>) with U ⊂𝒫() and sets I_ as in (<ref>). Then the relative truncation error can be estimated as 
    T_U f - S_I(U) fŁ_∞(^d)/f(^d)≤min{1/min_∈̆U N_, max_∈̆U√(∑_∈̨∖P_I_1/w^2()̨)}.

	The proof requires similar steps to the proof of <ref>. 


For the aliasing error in (<ref>), we start by considering the black-box approximation case where we solve the least-squares problem as described in <ref>. 
                                                                                                 
    Let f ∈(^d) with a weight function w^d[1,∞) such that ∑_∈̨^d1/w^2()̨ < ∞ and I(U) ⊂^d a finite frequency index set of type (<ref>) with sets I_ as in (<ref>). Moreover, we have a reconstructing rank-1 lattice Λ(ẕ, M, I(U)) for a generating vector ẕ∈^d and lattice size M ∈. Then the aliasing error can be estimated as 
    S_I(U) f - S_I(U)^Λ(ẕ, M, I(U)) f(^d)/f(^d)≤√(∑_∈̨^d∖ I(U)1/w^2()̨).
 Furthermore, if f ∈(^d) we get for the Ł_∞-norm 
    S_I(U) f - S_I(U)^Λ(ẕ, M, I(U)) fŁ_∞(^d)/f(^d)≤1/min_∈̆U N_.

    We show the bound (<ref>) by first applying Parseval's identity and (<ref>)
    S_I(U) f - S_I(U)^Λ(ẕ, M, I(U)) f(^d)^2    = ∑_∈̨I(U)f̂_ - f^2 
       = ∑_∈̨I(U)∑_ẖ∈Λ^(ẕ, M)∖{0̱}+̨ẖf^2.
 We then incorporate the weight and utilize the Cauchy-Schwarz inequality to obtain 
    S_I(U) f - S_I(U)^Λ(ẕ, M, I(U)) f(^d)^2    = ∑_∈̨I(U)∑_ẖ∈Λ^(ẕ, M)∖{0̱}w(+̨)/w(+̨)+̨ẖf^2 
       ≤∑_∈̨I(U)( ∑_ẖ∈Λ^(ẕ, M)∖{0̱} w^2(+̨)+̨ẖf^2 ) ×
           ×( ∑_ẖ∈Λ^(ẕ, M)∖{0̱}1/w^2(+̨ẖ))

	From <cit.> we know that for fixed ∈̨I(U) we have disjoint sets M_{+̨∈Λ^(ẕ, M)∖{0̱}}⊂^d ∖ I(U). This means we are able to estimate 
    ∑_ẖ∈Λ^(ẕ, M)∖{0̱} w^2(+̨)+̨ẖf^2    = ∑_ℓ̱∈ M_ w^2(ℓ̱)ℓ̱f^2 
       ≤∑_ℓ̱∈^d w^2(ℓ̱)ℓ̱f^2 = f(^d)^2
 such that S_I(U) f - S_I(U)^Λ(ẕ, M, I(U)) f(^d)^2 ≤f(^d)^2 ∑_∈̨I(U)∑_ẖ∈Λ^(ẕ, M)∖{0̱}1/w^2(+̨). Using that the sets M_ are disjoint and ⋃_∈̨I(U)M_⊂^d∖ I(U) yields
 	
    ∑_∈̨I(U)∑_ẖ∈Λ^(ẕ, M)∖{0̱}1/w^2(+̨)   = ∑_∈̨I(U)∑_ℓ̱∈ M_1/w^2(ℓ̱)
       = ∑_ℓ̱∈⋃_∈̨I(U) M_1/w^2(ℓ̱)≤∑_ℓ̱∈^d ∖ I(U)1/w^2(ℓ̱).

 	The Ł_∞-bound (<ref>) is obtained similarly to the method used in the proof of <cit.>. We proceed as follows 
    S_I(U) f - S_I(U)^Λ(ẕ, M, I(U)) fŁ_∞(^d)   ≤∑_∈̆U∑_∈̨∖P_I_f
       = ∑_∈̆U∑_∈̨∖P_I_w()̨/w()̨f
       ≤1/min_∈̆U N_∑_∈̆U∑_∈̨∖P_I_ w()̨f.

	The result is obtained through estimating the sum by f(^d).


In the following we consider the approximation error for scattered data approximation with a fixed node set X ⊂^d. Previously, we assumed that the index set I(U) and the node set X are such that the Fourier matrix _I(U) has full rank. In this case the least-squares problem (<ref>) has a unique solution. Assuming that the nodes in X are i.i.d. random variables that are uniformly distributed in ^d, it is possible to achieve good bounds on the approximation error, see <cit.>. 


    Let f ∈(^d) with a weight function w^d[1,∞) such that ∑_∈̨^d1/w^2()̨ < ∞, X ⊂^d a finite set of i.i.d. uniformly distributed points, y̱ = (f())_∈ X, and I(U) ⊂^d a finite frequency index set of type (<ref>) with U ⊂𝒫() a subset of ANOVA terms and sets I_ as in (<ref>). If for the number of frequencies we have I(U)≤X/7rlogX, r > 0, then 
    sup_f(^d)≤ 1∑_∈ X(f-S_I(U) f)()^2/X≤ 5 max( θ_I(U)^2,  8rκ^2logX/X∑_∈̨^d∖ I(U)1/w^2()̨)
 with a probability of at least 1-3X^1-r for θ_I(U) = f - S_I(U)fŁ_2(^d) and κ = 1+√(5)/2.

    The setting of this lemma is a special case of <cit.>.

The following theorem deals with the actual approximation error by incorporating the previous lemma.

    Let f ∈(^d) with a weight function w^d[1,∞) such that ∑_∈̨^d1/w^2()̨ < ∞, I(U) ⊂^d a finite frequency index set of type (<ref>) with sets I_ as in (<ref>). Moreover, U ⊂𝒫(), and S_I(U)^X f are the corresponding approximate Fourier partial sum obtained through the scattered data approximation method described in <ref>. If the elements of X ⊂^d are i.i.d. random variables uniformly distributed on ^d and for the number of frequencies we have I(U)≤X/7rlogX, r > 0, then
    
    S_I(U) f - S_I(U)^X fŁ_2(^d)/f(^d)≤√( 8 max( θ_I(U)^2, κ^2 logX/X∑_∈̨^d∖ I(U)1/w^2()̨) )
 with a probability of at least 1-3X^1-r for θ_I(U) = f - S_I(U)fŁ_2(^d) and κ = 1+√(5)/2.

    We denote the Fourier coefficients with ĉ̱̂ = (f)_∈̨I(U) and the approximate Fourier coefficients computed by <ref> with f̱̂̂ = (f̂_)_∈̨I(U). With Parseval's identity as well as the Moore-Penrose inverse we obtain 
    S_I(U) f - S_I(U)^X fŁ_2(^d)   = √(∑_∈̨I(U)f̂_-f^2) = f̱̂̂ - ĉ̱̂2
       = ( _I(U)^∗_I(U))^-1_I(U)^∗y̱ - ĉ̱̂2
       = ( _I(U)^∗_I(U))^-1_I(U)^∗( y̱ - _I(U)ĉ̱̂)2.

    We use the properties of the spectral norm and estimate further 
    ≤( _I(U)^∗_I(U))^-1_I(U)^∗2y̱ - _I(U)ĉ̱̂2
       = ( _I(U)^∗_I(U))^-1_I(U)^∗2√(∑_∈ X(f-S_I(U) f)()^2. )
 
    Applying <cit.> yields 
    sup_f(^d)≤ 1S_I(U) f - S_I(U)^X fŁ_2(^d)≤sup_f(^d)≤ 1√(2/X∑_∈ X(f-S_I(U) f)()^2 ).
 Finally, we use <ref> to obtain our bound 
    sup_f(^d)≤ 1S_I(U) f - S_I(U)^X fŁ_2(^d)≤√( 8 max( θ_I(U)^2, κ^2 logX/X∑_∈̨^d∖ I(U)1/w^2()̨) )
 with a probability of at least 1-3X^1-r. 


This concludes the consideration of the error of the presented method in both approximation scenarios. We were able to achieve bounds for  and Ł_∞ for functions in weighted Wiener algebras and Sobolev type spaces.



§ NUMERICAL RESULTS


We present numerical results for the method described in <Ref> for a test function f[0,1)^9, 
    f()  B_2(x_1) B_4(x_5) + B_2(x_2) B_4(x_6) + B_2(x_3) B_4(x_7) + B_2(x_4) B_4(x_8) B_6(x_9) ,

where B_2, B_4 and B_6 are parts of univariate, shifted, scaled and dilated B-splines of order 2, 4, and 6, respectively, see <Ref> for illustration. Their Fourier series is given by 
    B_j(x)  c_j ∑_k ∈sinc^j(π· k/j) cos( π· k )  ^2πı k · x
 with sinc(x) sin(x)/x and the constants c_2 √(3/4), c_4 √(315/604), c_6 √(277200/655177) such that B_jŁ_2(^d) = 1. This allows the direct computation of the Fourier coefficients f and the norm fŁ_2(^d). The ANOVA terms f_ are only nonzero for ∈̆U^∗𝒫({ 1,5 }) ∪𝒫({ 2,6 }) ∪𝒫({ 3,7 })∪𝒫({ 4,8,9 }). The function f therefore has an exact low-dimensional structure for d_s = 3, i.e., T_3 f = f. This leads to d_s = 3 being the optimal choice for the superposition threshold with no error caused by ANOVA truncation since it corresponds to the superposition dimension  for δ = 1, see (<ref>). In an approximation scenario with an unknown function f this information is of course not known.

We consider two errors 
    ε_ℓ_2 = y̱ - (S_I(U_d_s)^X f ())_∈ X2/y̱2, and  ε_L_2 = f - S_I(U_d_s)^X fŁ_2(^9)/fŁ_2(^9).
 Here, the error ε_ℓ_2 can be regarded as a training error since it is taken at the given sampling set X and the error ε_L_2 as a type of generalization error since it measures the error in the Fourier coefficients. Since our goal is to find the important ANOVA terms, i.e., the terms in U^∗, we expect to have an interval (or gap) in which to choose the order-dependent threshold ε̱∈ [0,1]^d_s. Therefore, we define 
    I^(j) = ∅     assumption (<ref>) is not fulfilled
    
    (a^(j), b^(j))    assumption (<ref>) is fulfilled

with 1 ≤ j ≤ d_s and 
    a^(j)   max{S_I(U_d_s)^X f∈̆U_d_s∖ U^∗,  = j }, 
    
    b^(j)   min{S_I(U_d_s)^X f∈̆U^∗,  = j }.
 Here, the assumption (<ref>) is to be understood for every order of terms, i.e., for $̆ and$̌ with =  = j.


	The norm occurring in the error ε_L_2 can be calculated using Parseval's identity 
    f - S_I(U_d_s)^X fŁ_2^2 = fŁ_2^2 + ∑_∈̨I(U_d_s)f - f̂_^2  - ∑_∈̨I(U_d_s)f^2
 which is possible since we know the exact Fourier coefficients and the norm of the function f. In general, this error cannot be computed.


 §.§ Scattered Data Approximation


For our numerical experiments we use one sampling set X ⊂^9 of uniformly distributed nodes with M X = 2.5 · 10^6, and an evaluation vector y̱ = ( f() )_∈ X. We are going to start by choosing three as the superposition threshold d_s while later reducing it to two which allows us to see the effect of truncating an ANOVA term. Our primary aim for now is to detect the ANOVA terms in U^∗ which we achieve using the first step of our method, see <ref>. To this end, we choose a frequency index set I(U_d_s) ⊂^9, cf. (<ref>), through order-dependent sets I_0 = {0}, I_1 = {-N_1/2,…,N_1/2-1}, I_2 = {-N_2/2,…,N_2/2-1}^2, and I_3 = {-N_3/2,…,N_3/2-1}^3 with N_1,N_2,N_3 ∈ 2. The method gives us an approximation S_I(U_d_s)^X f.

Results of numerical experiments with the function f from (<ref>) and different choices for the bandwidths N_1, N_2, and N_3 are displayed in <Ref>. They show that it is indeed possible to detect the ANOVA terms in U^∗ using trigonometric polynomials of small degrees. Moreover, both errors are roughly of the same order. Since our number of samples M is fixed, we are looking for values Ṉ such that one balances the effects of underfitting and overfitting. The experiments suggest that the choice in examples 5 and 8 is close to optimal. In <Ref> we depicted the global sensitivity indices S_I(U_d_s)^X f, cf. <ref>, for example 8 from <Ref>. The one-dimensional sets {i}, i=1,…,9, all have large indices as they are all in U^∗ while the two dimensional sets {1,5}, {2,6}, {3,7}, {4,8}, {4,9}, {8,9}∈ U^∗ are clearly separated from the two dimensional sets in U_d_s∖ U^∗. The same holds for the one three-dimensional term {4,8,9}∈ U^∗. The size of the intervals I^(j) suitable to choose the parameters ε_j is especially relevant since it separates important from unimportant terms.