A sparse FFT approach for ODE with random coefficients

The paper presents a general strategy to solve ordinary differential equations (ODE), where some coefficient depend on the spatial variable and on additional random variables. The approach is based on the application of a recently developed dimension-incremental sparse fast Fourier transform. Since such algorithms require periodic signals, we discuss periodization strategies and associated necessary deperiodization modifications within the occurring solution steps. The computed approximate solutions of the ODE depend on the spatial variable and on the random variables as well. Certainly, one of the crucial challenges of the high-dimensional approximation process is to rate the influence of each variable on the solution as well as the determination of the relations and couplings within the set of variables. The suggested approach meets these challenges in a full automatic manner with reasonable computational costs, i.e., in contrast to already existing approaches, one does not need to seriously restrict the used set of ansatz functions in advance.


Introduction
During the last years, the concept of random variables has become a very popular tool to model uncertain properties mathematically.For instance, diffusion characteristics of inhomogeneous materials can be distinctly more accurately described by functions that additionally depend on random variables.One common area of application of these mathematical designs are diffusion coefficients in differential equations.Certainly, the additional random variables affect the solvability and -if exist -the solutions of the differential equations under consideration.Besides investigations on existence, uniqueness and regularity of solutions for specific mathematical problems that involve randomness, cf.e.g.[5,10,2,1], numerical solution approaches need to be developed in order to compute approximations of the desired solutions.Accordingly, the established numerical solution approaches for differential equations without random coefficients need to be -at least -extended in order to meet the new challenges that are caused by the randomness of the diffusion coefficient.Commonly, discretizations of the domain of the stochastic variables lead to discretized solutions that are used to compute solutions in polynomial spaces or finite element representations, cf.e.g.[16,8,19,9,7].One essential task in this approach is the choice of suitable polynomial spaces and corresponding basis polynomials, which can be extremely challenging for higher numbers of random variables and occuring dependencies within the random variables.Furthermore, preferable choices of the used basis functions can improve the efficiency of the arising computations.
In the recent literature, many solution approaches deal with differential equations with random parameters in their coefficients.Most commonly, one suggests to choose several fixed instances of the random variables and applies known solvers for the considered differential equations without random coefficients for each of those instances.One achieves a set of solutions and computes the quantities of interest, which may be the coefficients of a specific expansion of the full solution of the differential equation or simply the expectation function, from these solutions by applying stochastic estimators.
In this paper, we present a closed approach that deals with the spatial variable as well as the random variables simultaneously in order to solve an ordinary differential equation with a diffusion coefficient affected by randomness.In more detail, we consider the differential equation with homogeneous boundary conditions, where a : D a → R, D a := × 1+d ξ j=1 [α j , β j ] ⊂ R × R d ξ , is the diffusion coefficient, that depends on the spatial variable η as well as the random variables that are the components of the vector ξ ∈ R d ξ , and the right hand side f : D f → R, D f := [α 1 , β 1 ] ⊂ R, is a function that depends only on the spatial variable η.Here we would like to point out, that the suggested solution approach is not restricted to the homogeneous boundary conditions or the independence on the random variables ξ of the right hand side f .Simple modifications of the presented approach leads to a solution strategy even for more general settings.However, the restrictions will simplify the notations and help to preserve clarity.
The essential restriction 0 for all (η, ξ) ∈ D a guarantees the existence of a unique solution u(•, ξ) of (1.1) for each fixed ξ.Hence, we suggest to approximately compute the unique solution u of (1.1) by means of a dimension-incremental sparse fast Fourier transform (FFT) approach, cf.[18,15], and based on a direct reversion of the occurring derivatives.On the one hand, the assumptions on the differential equation (1.1) do not guarantee for periodic signals that has to be treated.On the other hand, the sparse FFT approaches consider the input signals as periodic signals and is more successful -in the sense of approximation rates, number of needed samples to ensure a specific accuracy, etc. -when dealing with smooth periodic signals.Therefore, the periodization of the arising signals will be necessary in order to compute good approximate solutions of (1.1).At this point we would like to highlight that we approximately compute a complete solution of (1.1), which can be used to subsequently approximate several quantities of interest.The crucial advantage of this approach is that these complete approximate solutions reveals detailed characteristics of the random variables, i.e., the influence of each single variable on the solution as well as the interaction between different variables.We stress that the suggested strategy automatically detects these detailed characteristics with reasonable computational costs.
The paper is organized as follows: First we roughly outline the concept of the dimensionincremental sparse FFT approach and indicate the basic properties of suitable periodization mappings in Section 2. Section 3 presents the suggested strategy to treat the considered problem which leads to an approximation of the solution of the considered differential equation.As mentioned above, one may be interested in specific quantities of interest of this solution.Thus, we demonstrate how to compute nth moments of the computed solution in Section 4. Section 5 contains various numerical examples, shows the operability of the suggested approach, and discusses advantages and disadvantages of the applied sparse FFT approaches.

Sparse FFT
As mentioned above, we suggest to approximate the solution of an ordinary differential equation with random coefficients using the dimension-incremental FFT approach presented in [15].In this section, we declare the necessary notation and indicate the basic idea of this algorithm.
The aim of the dimension-incremental approach is the reconstruction of the Fourier coefficients pk , k ∈ I, of an arbitrarily chosen trigonometric polynomial p(x) := k∈I pk e 2πik•x , where the frequencies k are supported on a frequency set I ⊂ Z d of finite cardinality, i.e. |I| < ∞.In contrast to usual FFT algorithms, the challenge of sparse FFT algorithms is the efficient determination of the unknown frequency set I in addition to the Fourier coefficients pk using only sampling values of p.
Appropriate thresholding strategies within dimension-incremental sparse FFT algorithms allow for the treatment of general functions f ∈ L 1 (T d )∩C(T d ), i.e., the sparse FFT determines an approximation of the frequency index set I as well as an approximation of the (roughly) largest Fourier coefficients of the function f .Accordingly, the algorithms can be used in order to compute approximations SI of the Fourier partial sum for sufficiently smooth functions f .In this context, the Fourier partial sum S I [f ] is the truncated Fourier series of f , which implies the formal definition of the Fourier coefficients In general, the coefficients fk are just approximations of the Fourier coefficients c k (f ), since they are computed using only function evaluations of f and thus are disturbed at least by aliasing.
In order to compute both, the set I of the most significant frequencies k as well as approximations fk of the corresponding Fourier coefficients, a dimension-incremental approach was developed in [18,15], where the fundamental concept arises from a dimension-incremental method for the reconstruction of anharmonic trigonometric polynomials based on Prony's method, cf.[17].An outline of this concept can be found in [15,Sec. 2.2].
In this paper, we restrict the discussion to the in-and output of the algorithm, cf.Algorithm 1.We require a restricted search space Γ ⊂ Z d in frequency domain, where the significant Fourier coefficients are assumed to be supported.For simplicity and without crucial influence on the runtime of the algorithm, we can choose a tensor product box of equal edge lengths, i.e., we fix Γ = [−N, N ] d for a suitable edge length 2N + 1, N ∈ N. Since the used sampling nodes are chosen adaptively, we assume the function f being given as a black box.The parameter θ ∈ R + is a thresholding for the minimal absolute values that should be accepted as significant Fourier coefficient fk and its projections in lower dimensions.Additional sparsity parameters s, s local ∈ N restrict the algorithm to deal with at most s or r s local frequencies in each dimension-incremental step.Here, the parameter r is the number of projections that are used in each dimension-incremental step.Multiple projections are necessary in order to avoid detection failures caused by cancellations.Since we will use only function evaluations of the function f in order to compute an approximation, we have to apply suitable sampling strategies.For the case where we use multiple rank-1 lattices, the adaptive construction of the sampling set is affected by a certain small default probability.Therefore, it may happen that one has to start this construction of the sampling set more than once.The parameter b ∈ N can be used in order to restrict the number of restarts of the construction in each dimension-incremental step in order to guarantee the termination of the algorithm, cf.[15].However, this parameter is not restrictive during the computation, since even the choice b = 5 is not reached in practice.
The output of Algorithm 1 is the frequency set I and the corresponding approximated Fourier coefficients fk , k ∈ I, where SI f , cf. (2.1), is a good approximation of f when all significant frequencies are collected in I.
One crucial point of the dimension-incremental approach is the construction of spatial discretizations for trigonometric polynomials with frequencies in a certain, adaptively determined candidate set.Additional preferable properties of these spatial discretizations are • fast discrete Fourier transform algorithms, • fast construction methods for the spatial discretizations, and • low oversampling factors, i.e., the ratio of the number of sampling values to the cardinality of the candidate set should be low.
For high dimensional sparse trigonometric polynomials the concept of multiple rank-1 lattices, cf.[14,13] combines all these advantages, that are particularly beneficial to our targeted application.
For the sake of completeness, we give further details on the used sampling schemes.For a given generating vector z ∈ Z d and a lattice size M ∈ N, we define the rank-1 lattice where the modulo operation is applied componentwise.For a given frequency set I ⊂ Z d , |I| < ∞, the corresponding Fourier matrix is given by .
The dimension-incremental sparse FFT deals with different candidate sets of frequencies and asks for spatial discretizations for trigonometric polynomials with frequencies supported on these frequency sets.Additional requirements on Λ(z, M ) guarantees the spatial discretization property, i.e., the full column rank of the matrix A(Λ(z, M ), I).Due to the structure of Λ(z, M ), the computations of the matrix vector products involving A and its pseudo inverse can be performed by fast Fourier transform algorithms, cf.[11].These fast algorithms as well as the component-by-component construction algorithms for the used spatial discretizations, cf.[12], are the essential building blocks for the dimension-incremental sparse FFT based on single rank-1 lattices as spatial discretizations, which we denote by R1LsFFT, cf.[18] for details on that approach.A very similar approach is considered in [15], where the authors replaced the used sampling schemes by multiple rank-1 lattices, i.e., the spatial discretizations are constructed by the union of more than one rank-1 lattice, which provides -at least with high probabilityasymptotically lower oversampling factors as well as much faster construction approaches for spatial discretizations.Furthermore, fast Fourier transform algorithms for the evaluation and the reconstruction of trigonometric polynomials were developed, cf.[13].The corresponding dimension-incremental sparse FFT that uses these algorithms, i.e., the FFT algorithms as well as the construction algorithms for the spatial discretizations, is denoted by MR1LsFFT in the following.Recently a dimension-incremental sparse FFT based on random sampling were suggested, cf.[4].
The aforementioned dimension-incremental sparse FFT algorithms can be applied to periodic functions.Higher order smoothness of the treated functions often leads to smaller and thus preferable frequency sets I. Hence, we consider reasonable approaches to (smoothly) periodize non-periodic functions.

Periodization
The goal of a periodization is the approximation of a non-periodic function f : [α, β] → C using trigonometric polynomials that are naturally periodic and corresponding fast Fourier transform algorithms.Accordingly, we transform f to a periodic function g using a variable transform that has the following features ϕ is continuous in [0,1/2] and strictly increasing in (0,1/2), ( i.e., f (ϕ(x)) = g(x).In more detail, we are interested in approximations of the antiderivative of the function f , which leads to In order to compute F (t), t ∈ [α, β], we are interested in suitable approximations of g(x)ϕ ′ (x) for x ∈ [0, 1/2], which we want to realize using trigonometric polynomials.There are two different approaches to realize the computation of F .One point of view is to approximate g(x) and assume that ϕ ′ is constant almost everywhere in [0, 1/2], which leads to the well known tent transform approach [6,20].A more general approach will require additional assumptions on ϕ in order to obtain periodic smoothness of gϕ ′ that allow for suitable periodic approximations.For our purposes it is enough to deal with periodizations in one dimension.In higher dimensional settings, i.e., periodizations applied to a vector of variables, we simply apply the one-dimensional periodizations to each component of the vector.

Tent transform
The so-called tent-transform [6,20] is often used for periodization due to its simplicity.From a geometric point of view, the tent transform appends a mirror of the non-periodic function to the original function and dilates the resulting function such that its support is of length one.In addition, the new function is shifted such that its support is exactly [0, 1].In formula, the mapping . Certainly, this mapping ϕ is not continuously differentiable.Nevertheless, the constant first derivative ϕ ′ within (0, 1/2) and (1/2, 1) provides advantages within the integrals that we would like to deal with.For t ∈ [α, β] and Consequently, we only need to find an approximation of the antiderivative of the periodic function f • ϕ in order to achieve an approximation of an antiderivative of the non-periodic function f .

More general periodizations
In addition to the basic assumptions on the peridization mapping, cf.(2.2), we may assume periodic differentiability in order to obtain smoother integrands in (2.3).Higher order smoothness of the periodization could have positive effects for the approximation of the integrand gϕ ′ using trigonometric polynomials.Roughly speaking, the smoother the function, the faster the decay of the Fourier coefficients, i.e, the smaller the cardinality of the frequency index set of suitable approximating trigonometric polynomials.In some cases it may be enough to construct periodizations of a specific fixed smoothness, since the function f does not allow for higher order smoothness of g.E.g., splines of higher order seems to be ideally suited in order to guarantee the desired properties, cf.Example 2.1.In cases of functions f of higher but unknown smoothness, infinitely differentiable mappings ϕ may be an option to ensure that the periodization does not cause lower order smoothness of the integrand gϕ ′ .One suitable option for such a mapping is given in Example 2.2.However, the usage of more complicated mappings may imply disadvantages in the computation of the inverse mapping of the periodization.
Example 2.1 A spline of order four can be used to construct a periodization that is two times continuously differentiable.The mapping ϕ, plotted in Figure 2

Example 2.2
The cosine function can be used to construct a infinitely differentiable periodization mapping A corresponding plot for [α, β] = [−1, 1] can be found in Figure 2.1c.

ODE solver
Since the differentiations within the ODE (1.1) acts on only one variable, we revert the differentiation by integration and thus obtain a formal solution However, for high-dimensional variables ξ the computation of such a solution is a particular challenge.
We denote by F the antiderivative of f and obtain c 2 (ξ) = 0 since u * (α 1 , ξ) = 0 for homogeneous boundary conditions.The solution u * (t, ξ) changes to where we will use the term c 1 (ξ) in order to satisfy the boundary condition u * (β 1 , ξ) = 0.In particular, requirement (1.2) implies u 2 (β 1 , ξ) > 0 and thus fixing yields homogeneous boundary conditions for u * .Accordingly, for given f and a we need to compute suitable approximations of u 1 and u 2 .To this end, we will apply a dimensionincremental sparse FFT approach as described in Section 2.1.Since these FFT algorithms handles periodic signals, we need to periodize the upcoming functions.

Integration of the right hand side f
First, we determine the term F (α 1 ) − F (t) in (3.1) from above by approximating and integrating f .To this end, we periodize f using a suitable periodization ϕ, cf.Section 2.2, Accordingly, we obtain and approximate the integrand on the right hand side by a trigonometric polynomial which yields Consequently, we denote the approximation of the term F (α 1 ) − F (t) by F (t) and obtain 3.2 Approximating u 1 , u 2 , and c 1 We denote the integrands in (3.1) that determine u 1 and u 2 by v 1 and v 2 , respectively, i.e. we have .
First we consider the function v 1 and plug in the approximation F (η) of F (α 1 ) − F (η) from (3.4).This yields an approximation of v 1 v1 (η, ξ) := F (η) a(η, ξ) .Now, our goal is to construct an antiderivative of v1 with respect to η, which is an approximation of the antiderivative of v 1 .To this end, we construct a periodization of v1 using mappings ϕ η and ϕ ξ .We take into account the influence of the periodization during integration with respect to the first variable, which leads to the periodic integrand ṽ1 (x, y) = v1 (ϕ η (x), ϕ ξ (y)) ϕ ′ η (x).We compute a corresponding approximation SI which can be done by sparse FFT approaches, as described in Section 2.1, similar to those described in [18,15].The antiderivative of SI 1 [ ṽ1 ] with respect to x is given by , and we roll the periodization back, which leads to the approximation ȗ1 (t, ξ) := of u 1 given in (3.1).The analogous approach, but without approximating v 2 , leads to an approximation of u 2 (t, ξ) (3.8)The construction of ȗj (t, ξ) yields ȗj (α 1 , ξ) = 0, j = 1, 2. Consequently, each linear combination of ȗ1 and ȗ2 satisfies the homogeneous boundary condition in t = α 1 .A suitable approximation of c 1 (ξ), cf.(3.2), will lead to a linear combination of ȗ1 and ȗ2 that also satisfies the homogeneous boundary condition in t = β 1 .To this end, we periodize ȗ1 as well as ȗ2 and construct the approximation which are well defined due to the requirements on the diffusion coefficient a.We stress on the fact that the periodizations of ȗj do not coincide to the terms in (3.6), since these are non-periodic in general due to the terms that are linear in x.We approximate c1 using sparse FFT approaches by Altogether, an approximation of u * (t, ξ), cf.(3.1), is then given by ȗ(t, ξ) := ȗ1 (t, ξ) + c1 (ξ) ȗ2 (t, ξ), (3.11) which actually is built of three Fourier series combined with inverse mappings of the periodizations ϕ η and ϕ ξ .Algorithm 2 summarizes the approach stated above.

Tent transform in spatial domain
For the specific choice of the tent transform, cf.Section 2.2.1, for the periodization ϕ η and componentwise for ϕ ξ , we obtain some simplifications in the calculations above.Moreover, uniformly distributed random variables ξ j lead to additional simplifications due to the constant probability density.We observe and, hence, it is enough to compute an approximation of f = f • ϕ η , due to the equality Algorithm 2 Basic procedure for computing an approximation of the solution of (1.1) using a dimension-incremental sparse FFT approach Input: The approximation of SN [ f ] is preferable, since f ϕ ′ η is not continuous in the case where f (β 1 ) = 0 and, thus, problematic to approximate using trigonometric polynomials.Subsequent to the computation of SN [ f ] := N k=−N âk e 2πikx , the calculations of the antiderivative and deperiodization leads to The approximations ȗ1 and ȗ2 , cf. (3.7) and (3.8), can be computed in the exact same manner.Altogether, using the tent transform in spatial domain requires slight modifications in Algorithm 2 in lines 1, 3, and 4. The used sampling values must be computed with a factor 2 instead of ϕ ′ η .

Computing moments of the solution
In Section 3, we discussed a strategy for computing an approximate solution ȗ, cf.(3.11), of the ODE in (1.1).The computation of quantity of interests needs some further investigations.
For simplicity, we demonstrate one approach to compute approximations of the nth moments of the solution u * of the ODE in (1.1) based on the approximation ȗ.To this end, we denote the domain of the random variables by j=2 [α j , β j ].The nth moment of the solution u * of (1.1) is given by where p is the probability density function of the random variable vector ξ.Periodization yields where J is the involved Jacobian matrix.Assuming ϕ ξ is a periodization that acts on each component of ξ separately, cf.Section 2.2, the determinant of the Jacobian matrix is a tensor product function and we continue We approximate the integrand w n using a sparse FFT approach and achieve a Fourier partial sum Integrating S I 4 [w n ] instead of w n in (4.1) leads to the approximation of u * E n since each monomial that depends on y integrates to zero.

Numerical results
For our numerical tests, we use an example from [3] The parameters ξ k , k = 1, . . ., d ξ can be interpreted as random variables.Here we choose them to be uniformly distributed and we fix a 0 = 4.3 and γ = 2.In Figure 5.1 we (partially) plotted an approximation of the solution of this differential equation, where we restricted the number of random variables to two.
In order to demonstrate the applicability of the presented approach, we specify the settings of the applied algorithmic components.On the one hand, we restrict the numerical tests to the tent transform as periodization mapping, cf.Section 2.2.1, since this seems to be the most unfavourable choice due to its relatively low smoothness.On the other hand, we have to specify the applied sparse FFT approaches and the corresponding parameters.We choose three different sparsity levels s and refinements N for our approximated solutions ȗρ , ρ = I, II, III, cf.Section 3.2.Furthermore, we apply two different algorithms for computing approximate solutions denoted by ȗρ r1l and ȗρ mr1l namely the sFFT-algorithms that use sampling schemes that are r ank-1 lattices and multiple r ank-1 lattices, respectively.We call the corresponding sFFT algorithms R1LsFFT and MR1LsFFT.The basic structure of both algorithms is Increasing the number d ξ of random variables yields approximation problems of higher dimensionality.Clearly for practical applications, the number of random variables needs to be suitably bounded.The used diffusion coefficient a is build in such a way, that the influence of the random variable ξ j decreases with growing index j.Our first crucial task is to estimate the index j for which we can truncate the series expansion of a without losing significant information of a.In other words we would like to estimate a suitable number d ξ .
Example 5.1 To this end, we computed the approximation of a solution of (5.1) by our approach with a fixed large number d ξ = 40 of random variables, i.e., we treat an 41-dimensional approximation problem.We end up with an approximation ȗIII r1l as represented in (3.11).In order to simplify the considerations on the influence and the interactions on the variables of ȗ we apply periodizations and the sparse FFT approach on ȗ which leads in essence to a single Fourier sum representation of ȗ.The associated frequency set of this approximate solution -together with the absolute values of the occuring (Fourier) coefficients of this solutionallow for rating the random variables to their importance.In particular, if the expansion h j − l j , (k, h), (k ′ , l) ∈ I of the frequency set in direction j is zero -or very small and the corresponding coefficients almost zero in relation to the largest occuring coefficients -the solution does not or not significantly depend on the variable ξ j .Accordingly, leaving out this variable should not cause significant errors.Figure 5.2 indicates the expansions in each coordinate direction of the frequency set of ȗ for d ξ = 40.Obviously, the last 18 random variables have a very small expansion.We stress that the variables ξ 21 as well as ξ 22 have a significant frequency support but can be neglected due to the low order of magnitude of its Fourier coefficients.For these reasons, we restrict the number of random variables to d ξ = 20 in the following experiments.
As mentioned in the last example, we fix d ξ = 20.We solved (5.1) by the means of the sparse FFT approaches that uses single or multiple rank-1 lattices as spatial discretizations.The applied parameter constellations are presented in Tables 5.1 and 5.2.Both tables contains the total amount of samples that were used for the approximation of the functions ṽ1 , ṽ2 , and S I 4 [w 1 ] in columns six to eight for the different parameter settings as well.Moreover, the last columns of both tables present the cardinality of the full grids Ĝ21 N , where the sFFT algorithms search for the frequencies of the sparse representations of the computed approximations.
Example 5.2 We consider the average error of the computed approximations ȗρ † , † ∈ {r1l, mr1l}, of the solution u * for fixed spatial nodes η k .To this end we calculate the solution of (5.1) for n test = 20000 fixed randomly chosen parameters ξ i ∈ [−1, 1] 20 , i = 1, . . ., 20000 as grid functions defined on the uniform grid via numerical integration and an error bound of 10 −6 .We denote the corresponding solution by ǔ, i.e., we assume that the values ǔ(η k , ξ i ) are suitable approximations of the true solution and we use these function values for comparison against our approximations.For a first comparison, we consider the pointwise difference with respect to our approximated solution and calculate the mean, i.e.

Figure 5 . 2 :
Figure 5.2: Directional expansion of the frequency set I ⊂ Z 1+d ξ of ȗ for the random variables.

Table 5 .
2: Parameter settings, the number of samples M used for the computation of the approximations of ṽ1 , ṽ2 , and S I 4 [w 1 ] for d ξ = 20, cf.Sections 3.2 and 4, and the total cardinality of the box of frequency candidates for the sFFT algorithm that uses multiple rank-1 lattices as spatial discretizations for d ξ = 20..2show the parameters we used in columns two to five for the two different sFFT algorithms.The impacts of these parameters are shortly described in Section 2.1.