Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Parametric integration by magic point empirical interpolation

Parametric integration by magic point empirical interpolation Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 IMA Journal of Numerical Analysis (2019) 39, 315–341 doi: 10.1093/imanum/drx072 Advance Access publication on December 25, 2017 Maximilian Gaß Center for Mathematics, Technische Universitat ¨ Munc ¨ hen, Munich, Germany Kathrin Glau Center for Mathematics, Technische Universitat ¨ Munc ¨ hen, Munich, Germany and School of Mathematical Sciences, Queen Mary University of London, London, UK Corresponding author: k.glau@qmul.ac.uk [Received on 30 November 2015; revised on 13 October 2017] We derive analyticity criteria for explicit error bounds and an exponential rate of convergence of the magic point empirical interpolation method introduced by Barrault et al. (2004, An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations. C. R. Math., 339, 667–672). Furthermore, we investigate its application to parametric integration. We find that the method is well-suited to Fourier transforms and has a wide range of applications in such diverse fields as probability and statistics, signal and image processing, physics, chemistry and mathematical finance. To illustrate the method, we apply it to the evaluation of probability densities by parametric Fourier inversion. Our numerical experiments display convergence of exponential order, even in cases where the theoretical results do not apply. The magic point integration performs considerably well, when compared with Clenshaw–Curtis quadrature and with Chebyshev interpolation in the parameter space. Keywords: parametric integration; Fourier transform; parametric Fourier inversion; magic point interpola- tion; empirical interpolation; online–offline decomposition. 1. Introduction At the basis of a large variety of mathematical applications lies the computation of parametric integrals of the form I(h ) := h (z)μ(dz) for all p ∈ P, (1.1) p p where Ω ⊂ C is a compact integration domain, P a compact parameter set, (Ω,A, μ) a measure space with μ(Ω) < ∞ and h : P × Ω → C a bounded function. A prominent class of such integrals are parametric Fourier transforms iz,x f (z) = e f (x) dx for all p = (q, z) ∈ P, (1.2) q q with a parametric family of complex functions f defined on the compact set Ω. Today, Fourier transforms lie at the heart of applications in optics, electrical engineering, chemistry, probability, partial differential equations, statistics and finance. The application of Fourier transforms even sparked many advancements in those disciplines, underlining its impressive power. The introduction of Fourier analysis to image © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 316 M. GASS AND K. GLAU formation theory in Duffieux (1946) marked a turning point in optical image processing, as outlined, for instance, in Stark (1982). The application of Fourier transform to nuclear magnetic resonance (NMR) in Ernst & Anderson (1966) was a major breakthrough in increasing the sensitivity of NMR, for which Richard R. Ernst was awarded the Nobel prize in Chemistry in 1991. Fourier analysis also plays a fundamental role in statistics, in particular in the spectral representation of time series of stochastic data (see Brockwell & Davis, 2002). For an elucidation of other applications impacted by Fourier theory, we recommend (Kammler, 2007, Appendix 1). For parametric Fourier integrals of form (1.2) with fixed value of q and z being the only varying parameter, efficient numerical methods have been developed based on discrete Fourier transform (DFT) and fast Fourier transform (FFT). The latter was developed by Cooley & Tukey (1965) to obtain the Fourier transform f (z) for a large set of values z simultaneously. With regard to (1.2), this is the special case of a parametric Fourier transform, where q is fixed and z varies over a specific set. The immense impact of FFT highlights the exceptional usefulness of efficient methods for parametric Fourier integration. Shifting the focus to other examples, parametric integrals arise as generalized moments of parametric distributions in probability, statistics, engineering and computational finance. These expressions often appear in optimization routines, where they have to be computed for a large set of different parameter constellations. Prominent applications are model calibration and statistical learning algorithms based on moment estimation, regression and expectation–maximization (see, e.g., Hastie et al., 2009). The efficient computation of parametric integrals is also a cornerstone of the quantification of parameter uncertainty in generalized moments of form (1.1), which leads to an integral of the form h (z)μ(dz)P(dp), (1.3) P Ω where P is a probability distribution on the parameter space P. In all of these cases, parametric integrals of the form (1.1) have to be determined for a large set of parameter values. Therefore, efficient algorithms for parametric integration have a wide range of applica- bility. In the context of reduced basis methods for parametric partial differential equations, a magic point empirical interpolation method has been developed by Barrault et al. (2004) to treat nonlinearities that are expressed by parametric integrals. The applicability of this method to the interpolation of parametric functions has been demonstrated by Maday et al. (2009). In this article, we focus on the approximation of parametric integrals by magic point empirical interpolation in general, a method that we call magic point integration and – provide sufficient conditions for an exponential rate of convergence of magic point interpolation and integration in its degrees of freedom, accompanied by explicit error bounds; – translate these conditions to the special case of parametric Fourier transforms, which shows the broad scope of the results; – empirically demonstrate efficiency of the method in a numerical case study; and – compare its performance with the Clenshaw–Curtis quadrature and with Chebyshev interpolation in the parameter space. An in-depth case study related to finance is presented in Gaß et al. (2017) and supports our empirical findings. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 317 1.1 Sketch of the main idea To sketch the idea, let us start with an observation that is as simple as it is crucial. If the interpolation M M M M I (f ) := β e (z) ≈ f (z) with coefficients β ∈ C and functions e : Ω → C is a good m=1 m m m m approximation of f in L (Ω) in terms of stability and order of convergence, then M M I (f )(z) := β e (z)μ(dz) (1.4) m m m=1 is a good approximation for f (z)μ(dz) in the sense that few number of summands M are sufficient to obtain a high accuracy. An approximation of type (1.4), however, is only numerically efficient, if this is M M the case for the computation of both the coefficients β and the weights w := e (z)μ(dz). m m In classical quadrature rules, the interpolation I (f ) is a polynomial interpolation, for instance, the Chebyshev interpolation of f , which results in the well-known Clenshaw–Curtis quadrature rule. In this case, the Chebyshev coefficients can be easily and efficiently computed, and, thirdly the basis functions e are polynomials whence the weights w are known explicitly. We notice that classical quadrature rules yield the same quadrature nodes and basis functions for all integrands f on a given domain Ω. Yet as we often have more insight into our integration problem, namely that the integrand stems from a particular set {h | p ∈ P}, it is natural to ask: How can we systematically tailor a quadrature rule to a specific set of integrands? The conceptual difficulty arises from the fact that we want to obtain a method that autonomously finds and exploits the specific features of the given set of integrands. Therefore, we propose a two-step procedure in which the first step is a learning algorithm. Here, we choose the magic point interpolation method of Barrault et al. (2004). In a recursive greedy procedure, it delivers a nodal interpolation of functions tailored a given set {h | p ∈ P}. The learning ∗ ∗ M M procedure after M ∈ N steps delivers both the magic points z , ... , z and the functions θ , ... , θ on 1 M 1 M the basis of which the magic point interpolation operator ∗ M I (h)(p, z) := h (z )θ (z) (1.5) M p m m m=1 is constructed to approximate all functions from the set {h | p ∈ P} particularly well in the L -norm. This allows us to define the magic point integration operator ∗ M I (h)(p) := h (z ) θ (z)μ(dz). (1.6) M p m m m=1 This integration routine inherits its favourable numerical properties from the respective properties of the magic point interpolation. We will discuss convergence results in detail. To efficiently employ operator (2.2) for numerical integration, we split the procedure in an offline and an online phase. ∗ ∗ In the offline phase, we iteratively construct the magic points z , ... , z along with the functions 1 M M M M M θ , ... , θ , thus the operator (1.5). Moreover, we precompute all weights w := θ (z) dz. 1 M m Ω m The online phase consist in evaluating the sum ∗ M h (z )w . (1.7) m m m=1 Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 318 M. GASS AND K. GLAU 1.2 Related literature Another algorithm for parametric integration has been proposed by Antil et al. (2013). Their two-step greedy algorithm for reduced order quadratures is similar to the algorithm we propose, since it also divides the problem into an offline and an online phase, where an essential part of the offline phase consists in an application of an empirical interpolation. In contrast to the method we propose here, the authors explicitly abstain from applying the original magic point interpolation method of Barrault et al. (2004), and the resulting algorithm of Antil et al. (2013) significantly differs from magic point integra- tion. The main reason is that they consider parametric scalar products of functions in a Hilbert space to obtain a reduced basis. Their starting point thus is the representation of the integrands in an orthonormal basis of a Hilbert space. Following the spirit of reduced order modelling, they first consider the pro- jection to a large but finite dimensional subspace. In the first step of the offline phase, a reduced basis of the approximation space is generated. The second step of the offline phase is a greedy procedure, where the discrete empirical interpolation (DEIM) extracts empirical integration points from the reduced basis. Philosophically, both approaches differ as well: Antil et al. (2013) employ a two-stage procedure. First, the integration problem is discretized. Secondly, the discrete problem is reduced by exploiting additional knowledge. In contrast, the magic point integration that we propose systematically exploits insight into the infinite dimensional integration problem directly. Similarly to the procedure of Antil et al. (2013), low-rank tensor decompositions for parameter- dependent integrals are constructed in two steps, first employing a discretization of choice, and second reducing the tenor structure (see Ballani, 2012; the survey article by Grasedyck et al., 2013 and the further references therein). The two approaches, empirical interpolation and low-rank tensor decompo- sitions, are closely interrelated as discussed in Bebendorf et al. (2014). In contrast to these approaches, magic point integration uses the actual function set of interest as an input into the learning procedure. No a priori discretization method is involved. Instead, it is the learning procedure itself that deliv- ers the discretization. Moreover, the article Bebendorf et al. (2014) discusses an interesting and close relation between magic point interpolation and the adaptive cross approximation. This relation can be directly transferred to our framework, i.e., when we compare our ansatz relying on the magic point inter- polation of the parametric integrands to the use of the adaptive cross approximation of the integrands instead. Monte Carlo methods for parametric integrals have been developed and investigated in Heinrich (1998), Heinrich & Sindambiwe (1999), by introducing the multilevel Monte Carlo method. The key idea of this approach is to combine interpolation of the integrand in the parameter space with Monte Carlo integration in a hierarchical way. In contrast to our approach, the interpolation is known and fixed beforehand and, thus, is not tailored to the parametric family of integrands, which is the key feature of magic point integration. The remainder of the article is organized as follows. In Section 2, we revisit the magic point empirical interpolation method from Barrault et al. (2004), related error bounds and present the related integral approximation, which can be perceived from two different perspectives. On the one hand, it delivers a quadrature rule for parametric integrals and, on the other, an interpolation method for parametric integrals in the parameter space. In Section 3, we provide analyticity conditions on the parametric integrals that imply exponential order of convergence of the method. We focus on the case of parametric Fourier transforms in the subsequent Section 4. In a case study, we discuss the numerical implementation and its results in Section 5. Here, we use the magic point integration method to evaluate densities of a parametric class of distributions that are defined through their Fourier transforms. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 319 2. Magic point empirical interpolation for integration We now introduce the magic point empirical interpolation method for parametric integration to approxi- mate parametric integrals of the form (1.1). Before we closely follow Maday et al. (2009) to describe the interpolation method, let us state our basic assumptions that ensure the well-definedness of the iterative procedure. Assumption 2.1 Let (Ω, ) and (P, ) be compact, P × Ω (p, z) → h (z) bounded and ∞ ∞ p p → h be continuous, i.e., for every sequence p → p we have h − h → 0. Moreover, there exists p i p p ∞ p ∈ P such that the function h is not constantly zero. ∗ ∗ M M For M ∈ N, the method delivers iteratively the magic points z , ... , z , functions θ , ... , θ and 1 M 1 M constructs the magic point interpolation operator ∗ M I (h)(p, z) := h (z )θ (z), (2.1) M p m m m=1 which allows us to define the magic point integration operator ∗ M I (h)(p) := h (z ) θ (z) dz. (2.2) M p m m m=1 ∗ ∗ In the first step, M = 1, we define the first magic parameter p , the first magic point z and the first 1 1 basis function q by p := arg max ∞ , (2.3) p L p∈P z := arg max|h (z)|, (2.4) z∈Ω h (·) q (·) := . (2.5) h (z ) Notice that, thanks to Assumption 2.1, these operations are well-defined. For M ≥ 1 and 1 ≤ m ≤ M we define M M −1 M ∗ θ (z) := (B ) q (z), B := q (z ), (2.6) j m m jm jm j j=1 −1 M M where we denote by (B ) the entry in the jth line and mth column of the inverse of matrix B .By jm construction, B is a lower triangular matrix with unity diagonal and is thus invertible (see Barrault et al., 2004). Then, recursively, as long as there are at least M linearly independent functions in {h | p ∈ P}, the algorithm chooses the next magic parameter p according to a greedy procedure. It selects the parameter M Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 320 M. GASS AND K. GLAU so that h is the function in the set {h | p ∈ P} which is worst represented by its approximation with p p the previously identified M − 1 magic points and basis functions. So the Mth magic parameter is p := arg max h − I (h)(p, ·) ∞ . (2.7) p M−1 L p∈P In the same spirit, select the Mth magic point as ∗ ∗ z := arg max h (z) − I (h)(p , z) . (2.8) p M−1 M M z∈Ω The Mth basis function q is the residual r , normed to 1 when evaluated at the new magic point, i.e., M M r (z) := h (z) − I (h)(p , z), (2.9) M p M−1 r (·) q (·) := . (2.10) r (z ) Notice the well-definedness of the operations in the iterative step, thanks to Assumption 2.1, and the fact that the denominator in (2.10) is zero only if all functions in {h | p ∈ P} are perfectly represented by the interpolation I , in which case they span a linear space of dimension M − 1 or less. M−1 Let us understand the approach from two different perspectives. On the one hand, it is an interpolation method: Remark 2.2 (Interpolation of parametric integrals in the parameters) For each m = 1, ... , M, the function M m θ is a linear combination of snapshots h for j = 1, ... , M with coefficients β , which may be iteratively m j computed. We thus may rewrite formula (2.2)as M M ∗ m I (h)(p) = h (z ) β h (z) dz. (2.11) M p p m j m=1 j=1 Thus, for an arbitrary parameter p ∈ P, the integral h (z) dz is approximated by a linear combination of integrals h (z) dz for the magic parameters p . In other words, I (h) is an interpolation of the p M function p → h (z) dz. Here, the interpolation nodes are the magic parameters p and the coefficients ∗ m are given by h (z )β . In contrast to classic interpolation methods, the ‘magic nodes’ are tailored m=1 m j to the parametric set of integrands. On the other hand, magic point integration is an approximate integration rule: Remark 2.3 (Empirical quadrature rule for integrating parametric functions) Magic point integration is an interpolation method for integrating a parametric family of integrands over a compact domain. From this point of view, the magic point empirical interpolation of Barrault et al. (2004) provides a quadrature rule for integrating parametric functions. The weights are θ (z) dz and the nodes are the magic points z for m = 1, ... , M. As in the adaptive quadrature rules, these quadrature nodes are not fixed in advance. Although the adaptive quadrature rules iteratively determine the next quadrature node in view Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 321 of integrating a single function, the empirical integration method (2.2) tailors itself to the integration of agiven parametrized family of integrands. A discrete approach to empirical interpolation has been introduced by Chaturantabut & Sorensen (2010). Although the empirical interpolation is designed for parametric functions, the input data for DEIM is discrete. A canonical way to use DEIM for integration is to first choose a fixed grid in Ω for a discrete integration of all parametric integrands. Then, DEIM can be performed on the integrands evaluated on this grid. In contrast, using magic point empirical interpolation for parametric integration separates the choice of the integration grid from the selection of nodes p . For each function, a different integration discretization may be used. Indeed, in our numerical study, we leave the discretization choices regarding integration to Matlab’s quadgk routine. Fixing a grid for discrete integration beforehand might, e.g., become advantageous when the domain of integration Ω is high dimensional. 3. Convergence analysis of magic point integration We therefore summarize the results the convergence results for the magic point interpolation in the L -norm, which clearly extend to the convergence results for the magic point integration method. Inter- estingly, these results relate the convergence of the magic point interpolation to the best linear n-term approximation that is formally expressed by the Kolmogorov n-width, consider the monograph of Pinkus (1985) for a general reference. For a real or complex normed linear space X , and U ⊂ X , the Kolmogorov n-width is given by d (U , X ) = inf sup inf g − f , (3.1) U ∈E(X ,n) f ∈U n n g∈U where E (X , n) is the set of all n-dimensional subspaces of X . This concept allows to express the convergence behaviour of the interpolation subject to the ‘true complexity’ of set of functions of interest, if one agrees that the Kolmogorov n-width is a good measure of the intrinsic complexity. The approach is especially worthwhile when tackling the curse of dimensionality in the parameter or in the integration space. In the following section, we present a posteriori and a priori error bounds, mostly from the litera- ture, which are relevant for the integration method proposed. To obtain precise convergence rates, the Kolmogorov n-width needs to be estimated. We devote Section 3.2 below to this task and obtain widely applicable sufficient conditions for an exponential convergence of the magic point integration method. For univariate parametric Fourier integrals, the conditions boil down to exponential moment conditions. 3.1 Error bounds for the magic point interpolation We first present an a posteriori error bound, which follows as a corollary from the results in Maday et al. (2016). Let U :={h | p ∈ P} and X , = L (Ω), . Moreover, the convergence statement p ∞ involves the Lebesgue constant Λ := sup θ (z) , (3.2) z∈Ω m=1 with θ from equation (2.6). m Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 322 M. GASS AND K. GLAU Proposition 3.1 Let Ω ⊂ C , P and h : P × Ω → C be such that Assumption 2.1 is satisfied and that the sequence (Λ ) is monotonically increasing. Then, the following assertions hold: M M≥1 ∞ −α (a) If there exist constants c > 0 and α> 1, such that d (U , L (Ω)) ≤ c n for all n ≥ 1, then for 0 n 0 all M ≥ 1 and all h ∈ U,wehave 3α 3 1−α h − I (h) ≤ 3 · 2 c (1 + Λ ) M . M ∞ 0 M ∞ −c n (b) If there exist constants c , c > 0 and α> 0, such that d (U , L (Ω)) ≤ c e for all n ≥ 1, then 0 1 n 0 for all M ≥ 1 and all h ∈ U,wehave √ √ 2 −c M h − I (h) ≤ 2c (1 + Λ ) Me , M ∞ 0 M −2α−1 where c = c 2 . 2 1 Proof. The generalized empirical interpolation method (GEIM) has been introduced in (Maday et al., 2016, Section 2) to generalize the concept of the magic point interpolation to general Banach spaces instead of L (Ω). Indeed, it is elementary to verify that the magic point interpolation can be interpreted as a GEIM and that the assumptions P1–P3 in Maday et al. (2016, Section 2) are satisfied in our case. Moreover, the inequalities (9) in Maday et al. (2016) are satisfied for η = 1 and F = F = span{q , ... , q } for n 1 n all n ∈ N. We therefore can apply in Maday et al. (2016, Theorem 13). We notice that the Lebesgue constant Λ appearing in the latter theorem is differently defined as the operator norm of I , namely by n n I (ϕ) n ∞ n := sup ∞ . It is, however, elementary to verify that ≤ Λ = sup θ (z) . n n n ϕ∈L (Ω) z∈Ω m m=1 Therefore, we can replace the constant Λ in the result of Maday et al. (2016, Theorem 13) with the Lebesgue constant from equation (3.2). Since Λ is monotonically increasing by assumption, we can apply Corollary 11 and Lemma 12 from Maday et al. (2016) and obtain the analogue of Corollary 14 in Maday et al. (2016). As the last step of our proof, we estimate the constants derived in the latter corollary to simplify the expression. This is an elementary step, namely writing n = 4 + k, we see that = 2( + ) ≤ n/2 + 2 ≤ 3/2n. Next, we present a variant of Maday et al. (2009, Theorem 2.4), which is an a priori error bound for the magic point empirical interpolation method: Proposition 3.2 Let Ω ⊂ C , P and h : P × Ω → C be such that Assumption 2.1 is satisfied. If there exist constants α> log(4) and c > 0 such that ∞ −αM d U , L (Ω) ≤ ce , c α then for arbitrary ε> 0 and C := e + ε we have for all u ∈ U that −(α−log(4))M u − I (u) ≤ CMe . (3.3) The proof follows analogous to the proof of Maday et al. (2009, Theorem 2.4) and is in detail provided in Gaß et al. (2015a). In contrast to Maday et al. (2009, Theorem 2.4), Proposition 3.2 expresses the result directly in terms of the Kolmogorov n-width, which we prefer. We emphazize, however, that this is only a minor difference and the scientific contribution originates from Maday et al. (2009). Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 323 3.2 A priori error bounds under analyticity conditions On the basis of in Maday et al. (2009, Theorem 2.4), we identify two generic cases for which this result implies exponential convergence of magic point interpolation and integration. In the first case, the set of parametric functions consists of univariate functions that are analytic on a specific Bernstein ellipse. In the second case, the functions may be multivariate and do not need to satisfy higher-order regularity. The parameter set now is a compact subset of R, and the regularity of the parameter dependence allows an analytic extension to a specific Bernstein ellipse in the parameter space. To formulate our analyticity assumptions, we define the Bernstein ellipse B([−1, 1], ) with parameter > 1 as the open region in the complex plane bounded by the ellipse with foci ±1 and semiminor and semimajor axis lengths summing up to . Moreover, we define for b < b ∈ R the generalized Bernstein ellipse by B([b, b], ) := τ ◦ B([−1, 1], ), (3.4) [b,b] where the transform τ : C → C is given by [b,b] b − b b − b τ (x) := b + (1 −(x)) + i (x) for all x ∈ C. (3.5) b,b 2 2 For Ω ⊂ R we set B(Ω, ) := B([inf Ω, sup Ω], ). (3.6) Definition 3.3 Let f : X × X → C with X ⊂ C for m = 1, 2. We say f has the analytic property 1 2 m with respect to X ⊂ C in the first argument,if X ⊂ X and f has an extension f : X × X → C such 1 1 1 1 1 2 that for all fixed x ∈ X the mapping x → f (x , x ) is analytic in the inner of X and 2 2 1 1 1 2 1 sup |f (x , x )| < ∞. (3.7) 1 1 2 (x ,x )∈X ×X 1 2 1 2 We define the analytic property of f with respect to X in the second argument analogously. We denote I(h , μ) := h (z)μ(dz) (3.8) p p and for all p ∈ P and M ∈ N, I (h, μ)(p) := I (h)(p, z)μ(dz), (3.9) M M Ω Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 324 M. GASS AND K. GLAU where I is the magic point interpolation operator given by the iterative procedure from Section 2. Hence, ∗ M I (h, μ)(p) = h (z ) θ (z)μ(dz), (3.10) M p m m m=1 ∗ M where z are the magic points and θ are given by (2.6) for every m = 1, ... , M. m m Theorem 3.4 Let Ω ⊂ C , P and h : P × Ω → C be such that Assumption 2.1 holds and one of the following conditions is satisfied for some > 4, (a) d = 1, i.e., Ω is a compact subset of R, and h has the analytic property with respect to B(Ω, ) in the second argument, (b) P is a compact subset of R, and h has the analytic property with respect to B(P, ) in the first argument. Then, there exists a constant C > 0 such that for all p ∈ P and M ∈ N, −M h − I (h) ≤ CM( /4) , (3.11) −M sup I(h , μ) − I (h, μ)(p) ≤ Cμ(Ω)M( /4) . (3.12) p M p∈P Proof. Assume (a). Thanks to an affine transformation, we may without loss of generality assume Ω =[−1, 1]. As assumed, (p, z) → h (z) has an extension f : P × B(Ω, ) → C that is bounded p 1 and for every fixed p ∈ P is analytic in the Bernstein ellipse B(Ω, ). We exploit the analyticity prop- erty to relate the approximation error to the best n-term approximation of the set U :={h | p ∈ P}. This can conveniently be achieved by inserting an example of an interpolation method that is equipped with exact error bounds, and we choose Chebyshev projection for this task. From Trefethen (2013, Cheby Theorem 8.2) we obtain for each N ∈ N and the interpolation I defined in the appendix as the explicit error bound, Cheby −N sup f (p, ·) − I f (p, ·) ≤ c(f ) , (3.13) 1 1 1 p∈P with constant c(f ) := sup f (p, z) . Next, we can apply the general convergence result 1 1 (p,z)∈P×B(Ω, ) −1 from Maday et al. (2009, Theorem 2.4). Consulting their proof, we realize that −M sup h − I (h)(p, ·) ≤ CM( /4) , p M p∈P c(f ) with C = . Equation (3.12) follows by integration. The proof follows analogously under assumption (b).  Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 325 Remark 3.5 The proof makes the constant C explicit. Under assumption (a) it is given by C = max f (p, z) . (3.14) (p,z)∈P×B(Ω, ) 2( − 1) Under assumption (b), an analogous constant can be derived. Remark 3.6 Consider a multivariate integral of the form I (h, μ ⊗ μ ) := h(p, z)μ (dz)μ (dp), (3.15) M 1 2 1 2 P Ω D d with compact sets P ⊂ R and Ω ⊂ R equipped with finite Borel-measures μ and μ . Then, the 1 2 application of the magic point empirical interpolation as presented in Section 2 yields ∗ M I (h, μ ⊗ μ ) := h (z )μ (dp) θ (z)μ (dz), (3.16) M 1 2 p 2 1 m m P Ω m=1 and, under the assumptions of Theorem 3.4, we obtain −M I(h, μ ⊗ μ ) − I (h, μ ⊗ μ ) ≤ Cμ (Ω)μ (P)M( /4) . (3.17) 1 2 M 1 2 1 2 Remark 3.7 Theorem 3.4 has an obvious extension to functions (p, z) → h (z) that are piecewise analytic. To be more precise, we can assume that Ω =∪ Ω with piecewise disjunct and compact sets i=1 Ω ⊂ C such that for each i = 1, ... , n the mapping P × Ω (p, z) → h (z) has the analytic property i i p w.r.t. B(Ω , ) in the second argument. The proof can be generalized to this setting, and the error bound i i needs to be adapted in the obvious way. If the domain Ω is one-dimensional and z → h (z) enjoys desirable analyticity properties, the approximation error of magic point integration decays exponentially in the number M of interpolation nodes, independently of the dimension of the parameter space P. Vice versa, if the parameter space is one-dimensional and p → h (z) enjoys desirable analyticity properties, the approximation error decays exponentially in the number M of interpolation nodes, independently of the dimension of the integration space Ω. If the function set {h , p ∈ P} has a small Kolmogorov n-width, magic point interpolation and integration do not suffer from the curse of dimensionality—in contrast to standard polynomial approximation. This is owed to the fact that classic interpolation methods are adapted to very large function classes, whereas magic point interpolation is adapted to a parametrized class of functions. In return, this highly appealing convergence comes with an additional cost: The offline phase involves a loop over optimizations to determine the magic parameters and the magic points. Although the online phase is independent of the dimensionality, the offline phase is thus affected. Let us further point out that Theorem 3.4 is a result on the theoretical iterative procedure as described in Section 2. As we will discuss in Section 5.1 below, the implementation inevitably involves additional problem simplifications and approximations to perform the necessary optimizations. In particular, instead of the whole parameter space a training set is fixed in advance. For this more realistic setting, rigorous a posteriori error bounds have been developed for the empirical interpolation method (see Eftang et al., 2010). These bounds Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 326 M. GASS AND K. GLAU rely on derivatives in the parameters and straightforwardly translate to an error bound for magic point integration. In our implementation, we obtain the discrete training set by simulating according to the uniform distribution on the parameter space. For high-dimensional parameter spaces, simulation of pseudo-random numbers with quasi-Monte Carlo methods is highly promising. These methods are well-established and often reduce the computational complexity considerably. Also, other discretization techniques such as sparse grids can be employed. Moreover, the optimization problem (2.7), (2.8) itself can be implemented using more advanced techniques. In particular, for the applications that we consider in Section 5 and in Gaß et al. (2015a), the derivatives of the parametric integrands are explicitly available, which enables the use of gradient-based optimization techniques. 4. Magic points for Fourier transforms In various fields of applied mathematics, Fourier transforms play a crucial role and parametric Fourier inverse transforms of the form −izx f (x) = e f (z) dz (4.1) q q 2π D+1 need to be computed for different values p = (q, x) ∈ P = Q × X ⊂ R , where the function izy z → f (z) := e f (y) dy is well-defined and integrable for all q ∈ Q. q q A prominent example arises in the context of signal processing, where signals are filtered with the help of short-time Fourier transform (STFT) that is of the form −izt Φ (f )(z) := f (t)w(t − b)e dt, (4.2) with a window function w and parameter b. One typical example for windows are the Gaußian win- 2 2 −t /(2σ ) dows, w (t) = √ e , where an additional parameter appears, namely the standard deviation σ 2πσ of the normal distribution. The STFT with a Gaußian window has been introduced by Gabor (1946). His pioneering approach has become indispensable for time–frequency signal analysis. We refer to Debnath & Bhatta (2014) for historical backgrounds. We consider the truncation of the integral in (4.1) to a compact integration domain Ω =[Ω, Ω]⊂ R and choose the same domain Ω for all p = (q, x) ∈ P = Q × X . Thus, we are left to approximate for all p = (q, x) ∈ P the integral −izx I(h ) := h (z) dz, where h (z) = h (z) = e f (z). (4.3) p p p (q,x) q 2π Then, according to Section 2, we consider the approximation of I(h ) by magic point integration, i.e., by ∗ M I (h)(p) := h (z ) θ (z) dz, (4.4) M p m m m=1 ∗ M where z are the magic points and θ are given by (2.6) for every m = 1, ... , M. In such cases where the m m analyticity properties of q → f (z) or of z → f (z) are directly accessible, Theorem 3.4 can be applied to q q Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 327 estimate the error sup |I(h ) − I (h)(q, x)|. The following corollary offers a set of conditions (q,x) M (q,x)∈Q×X in terms of existence of exponential moments of the functions f . (η+ε)|x| Corollary 4.1 For η> 15/8 (Ω − Ω) and every parameter q ∈ Q assume e f (x) dx < ∞ for some ε> 0 and assume further η|x| sup e f (x) dx < ∞. q∈Q Then, −M h − I (h) ≤ C(η)M (η)/4 , (4.5) −M sup I(h ) − I (h)(q, x) ≤ C(η)|Ω|M (η)/4 , (4.6) (q,x) M (q,x)∈Q×X with 2η 2η (η) = + + 1, |Ω| |Ω| (η) η(−X ∨ X ) η|x| C(η) = e sup e f (x) dx, 4π( (η) − 1) q∈Q where |Ω| := Ω − Ω, X = inf X and X = sup X . Proof. From the theorem of Fubini and the lemma of Morera, we obtain that the mappings z → f (z) have analytic extensions to the complex strip R + i[−η, η]. We determine the value of (η). It has to be chosen such that the associated semiminor b of the Bernstein ellipse B([−1, 1], ) is mapped by τ onto the semiminor of the generalized Bernstein ellipse B([Ω, Ω], ) that is of length η. Using the [Ω,Ω] relation b = (4.7) between the semiminor b of a Bernstein ellipse and its ellipse parameter we solve Ω − Ω τ (ib ) = b = η (4.8) [Ω,Ω] for and find that for 2η 2η (η) = + + 1 (4.9) Ω − Ω Ω − Ω the Bernstein ellipse B(Ω, (η)) is contained in the strip of analyticity of f and by the choice of η we have (η) > 4. Hence, assumption (a) of Theorem 3.4 is satisfied. With regard to Remark 3.5, we also obtain the explicit form of the constant C(η), which proves the claim.  Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 328 M. GASS AND K. GLAU Similarly, we consider the case where Q is a singleton and −izx h (z) = h (z) = e f (z). (4.10) p x 2π Under this assumption we obtain an interesting additional assertion if the integration domain Ω is rather small. This case occurs, e.g., for approximations of STFT. Corollary 4.2 For η> 15/8 (X − X ),wehave −M h − I (h) ≤ C(η)M (η)/4 , (4.11) −M sup I(h ) − I (h)(x) ≤ C(η)|Ω|M (η)/4 , (4.12) x M x∈X with 2η 2η (η) = + + 1, |X | |X | (η) η(−Ω ∨ Ω) C(η) = e sup f (z) , 4π( (η) − 1) z∈Ω where |X | := X − X , X = inf X and X = sup X . The proof of the corollary follows similarly to the proof of Corollary 4.1. 5. Case study 5.1 Implementation and complexity The implementation of the algorithm requires further discretizations. For the parameter selection, we disc replace the continuous parameter space P by a discrete parameter cloud P . Additionally, for the disc selection of magic points we replace Ω by a discrete set Ω , instead of considering the whole continuous disc domain. Each function h is then represented by its evaluation on this discrete Ω and is thus represented by a finite-dimensional vector, numerically. The crucial step during the offline phase is finding the solution to the optimization problem ∗ M (p , z ) := arg max h (z) − h (z )θ (z) . (5.1) M+1 M+1 p p m m disc disc (p,z)∈P ×Ω m=1 In a continuous setup, problem (5.1) is solved by optimization routines. In our discrete implementation, disc however, we are able to consider all magic parameter candidates in the discrete parameter space P disc and all magic point candidates in the discrete domain Ω to solve problem (5.1). This results in a disc disc complexity of O(M ·|P |·|Ω |) during the offline phase for identifying the basis functions and magic points of the algorithm. We comment on our choices for the dimensionality of the involved discrete sets in the numerical section, later. Before magic point integration can be performed, the quantities θ (z) dz need to be computed for m = 1, ... , M. The complexity of this final step during the offline phase depends on the number of integration points that the integration method uses. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 329 5.2 Tempered stable distribution We test the parametric integration approach on the evaluation of the density of a tempered stable distri- bution as an example. This is a parametric class of infinitely divisible probability distributions. A random variable X is called infinitely divisible, if for each n ∈ N there exist n independent and identically dis- tributed random variables whose sum coincides with X in distribution. The rich class of infinitely divisible distributions, containing, for instance, the Normal and the Poisson distribution, is characterized through their Fourier transforms by the Levy–Khintchine ´ representation. Using this link, a variety of infinitely divisible distributions is specified by Fourier transforms. The class of tempered stable distributions is an example of this sort, namely it is defined by a parametric class of functions, which by the theorem of Levy–Khintchine ´ are known to be Fourier transforms of infinitely divisible probability distributions. Its density function, however, is not explicitly available. To evaluate it, a Fourier inversion has to be performed numerically. As we show below, magic point integration can be an adequate way to efficiently compute this Fourier inversion for different parameter values and at different points on the density’s domain. Following Kuchler & Tappe (2014), parameters + + − − + − α , λ , α , λ ∈ (0, ∞), β , β ∈ (0, 1) (5.2) define the distribution η on (R, B(R)) that we call a tempered stable distribution and write + + + − − − η = TS(q) = TS(α , β , λ ; α , β , λ ), (5.3) izx if its characteristic function ϕ (z) := e η (dx) is given by q q izx ϕ (z) = exp (e − 1)F (dx) , z ∈ R, (5.4) q q with Levy ´ measure + − α + α − −λ x −λ x F (dx) = e 1 + e 1 dx. (5.5) q x∈(0,∞) x∈(−∞,0) + − 1+β 1+β x x + − The tempered stable distribution is also defined for β , β ∈ (1, 2), where the characteristic function is given by a similar expression as in (5.4). We consider the tempered stable distribution in a special case by introducing the parameter space Q ={(C, G, M, Y ) | C > 0, G > 0, M > 0, Y ∈ (1, 2)} (5.6) and setting + − − + + − α = α = C, λ = G, λ = M, β = β = Y . (5.7) The resulting distribution γ = CGMY(C, G, M, Y ) is also known as CGMY distribution and is well established in finance (see Carr et al., 2002). The density function of a tempered stable or a CGMY As usual, we denote the parameters of the CGMY distribution with C, M, G, Y as a reference to the authors of Carr et al. (2002). It will be clear from the context whether M denotes the number of magic points or the CGMY parameter. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 330 M. GASS AND K. GLAU Table 1 Parameter intervals for the numerical study. The parameter Y is kept constant CG M Y x interval [1, 5][1, 8][1, 8] 1.1 [−1, 1] disc disc Fig. 1. Decay of the error (2.8)on P and Ω during the offline phase of the algorithm. distribution, respectively, is not known in the closed form. With q = (C, G, M, Y ) ∈ Q of (5.6) and z ∈ R, its Fourier transform, however, is explicitly given by Y Y Y Y ϕ (z) = exp CΓ(−Y ) (M − iz) − M + (G + iz) − G , (5.8) wherein Γ denotes the gamma function. By Fourier inversion, the density f can be evaluated, 1 1 −izx −izx f (x) = e ϕ (z) dz =  e ϕ (z) dz for all x ∈ R. (5.9) q q q 2π π R 0 5.3 Numerical results We restrict the integration domain of (5.9)to Ω =[0, 65] and apply the parametric integration algorithm on a subset of P = Q × R, (5.10) specified in Table 1. We draw 4000 random samples from P respecting the intervals bounds of Table 1 ∞ −12 and run the offline phase until an L accuracy of 10 is reached. Here, we compute the benchmark −14 values by numerical integration using Matlab’s quadgk routine with a tolerance 10 . The error decay during the offline phase is displayed by Fig. 1. We observe exponential error decay reaching the prescribed accuracy threshold at M = 40. Let us have a closer look at some of the basis functions q that are generated during the offline phase. Figure 2 depicts five basis functions that have been created at an early stage of the offline phase (left) Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 331 Fig. 2. Some basis functions q constructed early (left) and late (right) during the offline phase of the algorithm. together with five basis functions that have been identified towards the end of it (right). Note that all plotted basis functions are supported on a relatively small subinterval of Ω. They are numerically zero outside of it. Interestingly, the functions in both sets have intersections with the Ω axis rather close to the origin. Such intersections mark the location of magic points. Areas on Ω where these intersections accumulate thus reveal regions where the approximation accuracy of the algorithm is low. In these regions, the shapes across all parametrized integrands seem to be most diverse. In the right picture, we observe in comparison that at a later stage in the offline phase, magic points z are picked further away from the origin, as well. We interpret this observation by assuming that the more magic points are chosen close to the origin, the better the algorithm is capable of approximating integrands more precisely there. Thus, its focus shifts towards the right for increasing values of M where now smaller variations attract its attention. We now test the algorithm on parameters not contained in the training set. To this extent, we randomly draw 1000 parameter sets p = (C , G , M , Y , x ), i = 1, ... , 1000, uniformly distributed from the i i i i i i intervals given by Table 1. For each p , i = 1, ... , 1000, we evaluate f (x ) by Fourier inversion i (C ,G ,M ,Y ) i i i i i using Matlab’s quadgk routine. Additionally, we approximate f (x ) using the interpolation (C ,G ,M ,Y ) i i i i i operator I for all values m = 1, ... , M. For each such m we compute the absolute and the relative error. The maximal absolute error over the 1000 samples, i.e., the L ((p ) )-error for each i i=1,...,1000 m = 1, ... , M, is illustrated by Fig. 3. To be precise, we compute the benchmark values by numeri- cal Fourier inversion, restrict the resulting integration domain to Ω =[0, 65] and integrate numerically −14 with Matlab’s quadgk routine with a tolerance of 10 . The set-up of the numerical experiment does not fall in the scope of our theoretical result in several respects. We discretized both the parameter space and the set of magic point candidates. While Table 1 ensures a joint strip of analyticity R + i(−1, 1) that all integrands h share, a Bernstein ellipse B(Ω, ) with > 4 on which those integrands are analytic does not exist. In spite of all these limitations, we still observe that the exponential error decay of the 3 disc disc offline phase is maintained. From this, we draw two conclusions. First, both Ω and P appear This observation is confirmed by our empirical studies in Gaß et al. (2017), where the existence of such a shared strip of analyticity was sufficient for empirical exponential convergence. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 332 M. GASS AND K. GLAU Fig. 3. Out of sample error decay. For 1000 randomly drawn parameters (C , G , M , Y , x ), i = 1, ... , 1000 from the set prescribed i i i i i in Table 1, the CGMY density is evaluated by the parametric integration method and by numerical Fourier inversion via Matlab’s quadgk routine. For each number M of magic points, the maximal absolute difference between both is displayed. Fig. 4. All out of sample errors in detail. Left: The absolute errors achieved for each of the 1000 randomly drawn parameter sets. −3 Right: The relative errors for density values above 10 . sufficiently rich to represent their continuous counterparts. Secondly, the practical use of the algorithm extends beyond the analyticity conditions imposed by Theorem 3.4. Figure 4 presents the absolute errors and the relative errors for each parameter set p , individually, for the maximal number M magic points. Note that only relative errors for CGMY density values larger −3 than 10 have been plotted to exclude the influence of numerical noise. While individual outliers are not present among the absolute deviations, they dominate the relative deviations in contrast. This result is in line with the objective function that the algorithm minimizes. Finally, Fig. 5 displays all magic parameters identified during the offline phase, together with those randomly drawn parameter constellations that resulted in the ten smallest absolute errors (crosses) and the 10 largest absolute errors (stars). We observe that orange parameter constellations are found in areas densely populated by magic parameters, whereas green parameter constellations often do not have magic parameters in their immediate neighbourhood, at all. This result may surprise at first. On second thought, however, we understand that areas where magic parameters accumulate mark precisely those locations in the parameter space where approximation is especially challenging for the algorithm. Integrands associated with parameter constellations from there exhibit the largest variation. Approximation for white Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 333 5 8 8 1 3 4.5 4.5 0 1 1 1 -1 135 135 135 135 8 8 8 1 4.5 4.5 4.5 0 1 1 1 -1 135 1 4.5 8 14.5 8 14.5 8 8 8 8 1 4.5 4.5 4.5 0 1 1 1 -1 135 1 4.5 8 14.5 8 14.5 8 1 1 1 1 0 0 0 0 -1 -1 -1 -1 135 1 4.5 8 14.5 8 -1 01 Fig. 5. An illustration of those randomly drawn parameter constellations that resulted in the 10 worst (stars) and the 10 best (circles) absolute errors. The crosses (“+” signs) the magic parameters selected during the offline phase. areas, on the other hand, is already covered by the previously chosen magic parameters. Consequently, crosses are very likely to be found there. 5.4 Comparison of magic point integration and Clenshaw–Curtis Since in our numerical experiments, we compute univariate analytic integrands over a finite domain, it is interesting to compare the performance of the new integration routine to the Clenshaw–Curtis method. The latter is well-known for its excellence performance, particularly in case of analyticity. To compare both methods, we use again the setting introduced in Section 5.2. We randomly sample 1000 parameter constellations according to the uniform distribution in the set specified in Table 1.For these, we compute the benchmark integrals with Matlab’s quadgk routine for an absolute tolerance −14 of 10 and the maximal difference with respect to the results obtained by both the Clenshaw–Curtis quadrature rules and the magic point integration, for varying numbers of quadrature nodes, respectively, magic points, M. Figure 6 displays the resulting error decays. Here, the magic point integration method clearly out- performs the Clenshaw–Curtis quadrature rule. For less than 35 magic points, the accuracy of magic −10 point integration reaches already 10 , while the error of the Clenshaw–Curtis quadrature is still in the −12 second digit. The range of machine precision, i.e., an error of 10 , is achieved with about 40 magic points while about 200 Chebyshev nodes are required in the Clenshaw–Curtis routine. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 334 M. GASS AND K. GLAU Fig. 6. Comparison of the error decay of magic point integration versus Clenshaw–Curtis. The maximal absolute errors achieved for each of the 1000 randomly drawn parameter sets for varying numbers M of quadrature nodes. Left: Performance of both magic point integration and Clenshaw–Curtis for M = 1, ... , 41 quadrature nodes. Right: Performance of Clenshaw–Curtis for M = 1, ... , 200 quadrature nodes. Here, the magic point integration clearly outperforms the Clenshaw–Curtis quadrature. In accordance with the theoretical results, we observe an exponential error decay for both methods. For magic point integration, the rate of convergence, however, is much higher. This is not reflected by the current theoretical results. Indeed, in the proof of Theorem 3.4, we estimate the Kolmogorov n-width by the Chebyshev interpolation, i.e., by the Clenshaw–Curtis quadrature. The theoretical bound that we obtained is therefore larger than the error bound for the Clenshaw–Curtis quadrature. The numerical experiments reveal that the magic point interpolation is actually able to reduce the computational complexity by adapting a quadrature rule to the parametric family of integrands in the offline phase. Moreover, this gives evidence that the Kolmogorov n-width of the parametric family of integrands is indeed significantly smaller than the error obtained by the interpolation with Chebyshev polynomials. Figure 7 displays the magic points versus the Chebyshev nodes, scaled to the unit interval. It is interesting to see that there is not even a qualitative similarity between the two different sets of points. 5.5 Comparison of magic point integration and Chebyshev interpolation in the parameter space As highlighted in Remark 2.3 and Remark 2.2, magic point integration can be seen both as a quadrature rule and as an interpolation method in the parameter space. Therefore, we also compare the performance of the magic point integration with the Chebyshev interpolation in the parameter space. In our numerical experiments reported up to this point, we considered a five-dimensional parameter space as laid in Table 1. Tensorized Chebyshev interpolation, however, is exposed to the curse of dimensionality, and therefore magic point integration will clearly outperform this method. To gain more insight, we reduce the parameter space to a univariate, a bivariate and a three-dimensional setting, below. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 335 Fig. 7. Magic points versus Clenshaw–Curtis points for M = 40. 5.5.1 Comparison for a single free parameter. We fix the parameters C = 1, M = 4, Y = 1.1, x =−1, and only vary G ∈[1, 8]. For a uniform test grid of 100 parameters in [1, 8], we compute the −14 benchmark integrals with Matlab’s quadgk routine for an absolute tolerance of 10 , and the maximal difference with respect to the results obtained by both Chebyshev interpolation in the parameter G and magic point integration for a varying number of magic points and Chebyshev points M, respectively. We refer to the appendix for a precise definition of the Chebyshev interpolation that we implemented. Figure 8 displays the resulting error decays. Both methods perform comparably well. As expected, the convergence is exponential in the number of nodes. The magic point integration shows a slightly higher rate of decay as the Chebyshev interpolation. Figure 9 displays the first 10 magic parameters versus the first 10 Chebyshev nodes in G, scaled to the unit interval. Remarkably, we recognize similarities in the choices of points, as both sets cluster at the boundaries. We also recognize significant differences. The magic points first of all are not located at the Chebyshev points. Secondly, there is an asymmetry in their distribution with more points on the left side than on the right. 5.5.2 Comparison for two free parameters. We fix the parameters C = 1, M = 4, Y = 1.1, and vary G ∈[1, 8] and x ∈[−1, 1]. For a uniform test grid of 100 × 100 parameters in [1, 8]×[−1, 1],we −14 compute the benchmark integrals with Matlab’s quadgk routine for an absolute tolerance of 10 and the maximal difference with respect to the results obtained by both tensorized Chebyshev interpolation in G, x, and magic point integration for a varying number of magic points and Chebyshev points M, respectively. For the precise definition of the tensorized Chebyshev interpolation, we again refer to the appendix. Since we are using the tensorized Chebyshev interpolation, we have a grid of Chebyshev nodes. We specify N number of nodes in both directions, G and x, and obtain a total number of M = N number of nodes. Figure 10 displays the error decay of both the magic point integration and the tensorized Chebyshev interpolation in G, x. The magic point integration outperforms the Chebyshev interpolation by far. For −12 an accuracy in the range of machine precision, i.e., 10 , less than 25 magic parameters versus over 700 Chebyshev nodes are required. This is an empirical indicator that magic point integration has the potential to tackle the curse of dimensionality. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 336 M. GASS AND K. GLAU Fig. 8. Comparison of the error decay of magic point integration versus Chebyshev interpolation in the parameter space. The maximal absolute errors achieved for each of the 100 equidistant grid points in the parameter set G ∈[1, 8], for varying numbers M of interpolation nodes. Fig. 9. Magic parameters versus Chebyshev interpolation nodes for M = 19 interpolation nodes. Next, we investigate the performance in relation to the dimensionality of the parameter space in some more depth. Figure 10 shows that the convergence rate of magic point integration is considerably higher. Yet, it is not clear how much the method is still exposed to the curse of dimensionality in the parameter space. From Fig. 8, we may conclude that both methods perform comparably in the univariate setting. Therefore, we next display the error decays of the Chebyshev interpolation in a different scale. While for magic point integration we show the error for the number of magic points M, we show the error of the Chebyshev interpolation for N, the number of nodes chosen in each of the interpolation axes. The result is provided in Fig. 11. The similarity of both curves is striking. For M = N = 15, we observe roughly the same error of about −8 10 . To be clear, this means that 15 summands have been used for the magic point integration, to obtain Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 337 Fig. 10. Comparison of the error decay of magic point integration versus Chebyshev interpolation in the parameter space. The maximal absolute errors achieved for each of the 100 × 100 equidistant grid points in the parameter set (G, x) ∈[1, 8]×[−1, 1], for varying numbers M of interpolation nodes. Fig. 11. Comparison of the error decay of magic point integration versus Chebyshev interpolation in the parameter space; specifi- cations as in Fig. 10. For the Chebyshev interpolation the number of summands equals (N + 1) , while the number of summands in the magic point interpolation is M. This figure visualizes that the magic point integration seen as interpolation in the parameter space is not exposed to the curse of dimensionality—in contrast to the tensorized Chebyshev interpolation. −8 4 2 a maximal error of 10 over the 10 test parameter. For the same precision, (N + 1) = 256 summands are required for the tensorized Chebyshev interpolation in G, x. In other words, in these experiments, the magic point interpolation exhibits indeed no curse of dimensionality in the parameter space. Figure 12 displays the grid of tensorized Chebyshev nodes and of magic parameters for which roughly −12 the same accuracy of 10 has been achieved. This visualizes the complexity reduction of magic points versus Chebyshev interpolation. In contrast to the example of the univariate parameter space, the distri- bution of the magic parameters is totally different from the Chebyshev points. It is, moreover, apparent that they only appear at the boundaries of the domain. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 338 M. GASS AND K. GLAU Fig. 12. M = 23 Magic parameters versus (N + 1) Chebyshev interpolation nodes for N = 28, for which both methods achieve roughly the same accuracy. 5.5.3 Comparison for three free parameters with low-rank tensor techniques. In this section, we consider three varying parameters. We have already seen for the bivariate parameter space that the tensorized Chebyshev interpolation is exposed to the curse of dimensionality, in contrast to the magic point empirical interpolation. Low-rank tensor techniques start from the tensorized Chebyshev interpolation and exploit its compressible structure by a greedy procedure, which is similar to the offline phase in the magic point interpolation. It is, therefore, very interesting to investigate the performance of low-rank tensor techniques in this example. To do so, we use the chebfun toolbox for Matlab, see http://www.chebfun.org (accessed 16 November 2017), developed by Prof. Trefethen and his working group, particularly we use chebfun3 written by Benham Hashemi (see also Hashemi & Trefethen, 2017). We fix the parameters C = 1, Y = 1.1 and vary (G, M, x) ∈[1, 8]×[1, 8]×[−1, 1]. We uniformly sample a set of 1000 test parameters from [1, 8]×[1, 8]×[−1, 1], compute the benchmark integrals with −14 Matlab’s quadgk routine for an absolute tolerance of 10 and the maximal difference with respect to the results obtained by chebfun3 from in G, M, x for a number of absolute tolerances, namely for −1−k/2 (10 ) , and the magic point integration for varying numbers of magic points M. To compare the k=0,...,22 complexity of both interpolation methods, we compute the number M of summands of each chebfun3 object. Moreover, chebfun3 allows to extract the number M of summands required for the same precision without the use low-rank compression techniques. Figure 13 displays the error achieved by chebfun3, magic point integration for the number of summands M each. The third line, which is dashed, shows the error for chebfun3 for the number of summands that would be required for the tensorized Chebyshev interpolation. The figure shows that the low-rank tensor technique performs considerably better compared to the −2 tensorized Chebyshev interpolation. For instance, chebfun3 achieves an accuracy of about 10 with less than 200 summands, while the interpolation based on the full tensor structure requires about 500 summands. The performance of magic point integration outperforms both methods by far. This can be Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 339 Fig. 13. Magic integration and Chebyshev interpolation three varying parameters G, M, x. explained by the fact that the accuracy of the tensorized Chebyshev interpolation is determined by the resolution along the axes, which carries over to the low-rank approximations. This a priori choice of axes, however, is not present for the magic point interpolation. Instead, the set of basis functions is iteratively produced for the set of functions of interest. Acknowledgements The authors thank Bernard Haasdonk, Laura Iapichino, Daniel Wirtz and Barbara Wohlmuth for fruitful discussions, and two anonymous referees for their valuable comments and suggestions. Maximilian Gaß thanks the KPMG Center of Excellence in Risk Management and Kathrin Glau acknowledges the TUM Junior Fellow Fund for financial support. References Antil, H., Field, S. E., Herrmann, F., Nochetto, R. H. & Tiglio, M. (2013) Two-step greedy algorithm for reduced order quadratures. J. Sci. Comput., 57, 604–637. Ballani, J. (2012) Fast evaluation of singular BEM integrals based on tensor approximations. Numer. Math., 121, 433–460. Barrault, M., Maday, Y., Nguyen, N. C. & Patera, A. T. (2004) An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations. C. R. Math., 339, 667–672. Bebendorf M., Maday Y. & Stamm B. (2014) Comparison of Some Reduced Representation Approximations. In: A. Quarteroni & G. Rozza (eds) Reduced Order Methods for Modeling and Computational Reduction. MS&A - Modeling, Simulation and Applications, vol 9. Cham: Springer, pp. 67–100. Brockwell, P. J. & Davis, R. A. (2002) Introduction to Time Series and Forecastung, 2nd edn. New York: Springer. Carr, P., Geman, H., Madan, D. B. & Yor, M. (2002) The fine structure of asset returns: An empirical investigation. J. Bus., 75, 305–333. Chaturantabut, S. & Sorensen, D. C. (2010) Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput., 32, 2737–2764. Cooley, J. W. & Tukey, J. W. (1965) An algorithm for the machine calculation of complex Fourier series. Math. Comput., 19, 297–301. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 340 M. GASS AND K. GLAU Debnath, L. & Bhatta, D. (2014) Integral Transforms and Their Applications, 3rd edn. Boca Raton, USA: Apple Academic Press Inc. Duffieux, P.-M. (1946) L’Integr ´ ale de Fourier et ses Applications a ` l’Optique. Rennes: Oberthur. Eftang, J. L., Grepl, M. A. & Patera, A. T. (2010) A posteriori error bounds for the empirical interpolation method. Ser. I, 348, 575–579. Ernst, R. R. & Anderson, W. A. (1966) Application of Fourier transform spectroscopy to magnetic resonance. Rev. Sci. Instrum., 37, 93–102. Gabor, D. (1946) Theory of communication. J. Inst. Electr. Eng. III Radio Commun. Eng., 93, 429–457. Gaß, M., Glau, K., Mahlstedt, M. & Mair, M. (2015a) Chebyshev interpolation for parametrized option pricing. Working Paper. Available at URL http://arxiv.org/abs/1505.04648. Accessed 16 November 2017. Gaß, M., Glau, K. & Maikr, M. (2017) Magic points in finance: empirical integration for parametrized option pricing, SIAM J. Finan. Math., 8, 766–803. Grasedyck, L., Kressner, D. & Tobler, C. (2013) A literature survey of low-rank tensor approximation technique. GAMM-Mitt., 36, 53–78. Hashemi, B. & Trefethen, L. N. (2017) Chebfun in three dimensions. SIAM J. Sci. Comput., 39, C341–C363. Hastie, T., Tibshirani, R. & Friedman, J. (2009) The Elements of Statistical Learning, 2nd edn. Heidelberg: Springer. Heinrich, S. (1998) Monte Carlo complexity of global solution of integral equations. J. Complex., 14, 151–175. Heinrich, S. & Sindambiwe, E. (1999) Monte Carlo complexity of parametric integration. J. Complex., 15, 317–341. Kammler, D. W. (2007) A First Course in Fourier Analysis, 2nd edn. New York: Cambridge University Press. K¨ uchler, U. & Tappe, S. (2014) Exponential stock models driven by tempered stable processes. J. Econometrics, 181, 53–63. Maday, Y., Mula, O. & Turinici, G. (2016) Convergence analysis of the generalized empirical interpolation method. SIAM J. Numer. Analy., 54, 1713–1731. Maday, Y., Nguyen, C. N., Patera, A. T. & Pau, G. S. H. (2009) A general multipurpose interpolation procedure: the magic points. Comm. Pure Appl. Anal., 8, 383–404. Pinkus, A. (1985) n-widths in approximation theory. Ergebnisse der Mathematik und ihrer Grenzgebiete [Results in Mathematics and Related Areas (3)], Vol. 7. Heidelberg: Springer. Stark, H. (1982) Theory and measurement of the optical Fourier transform. Applications of Optical Fourier Transforms (H. Stark ed). New York: Academic Press, pp. 2–40. Trefethen, L. N. (2013) Approximation Theory and Approximation Practice. New York: SIAM Books. Appendix Since different versions of the Chebyshev interpolation are used in the literature, we introduce the version that we use. For interpolation univariate functions f : [−1, 1]→ R, the interpolation with Chebyshev polynomials of degree N is of the form Cheby I (f )(x) := c T (x) (A.1) j j j=0 0<j<N 2 k with coefficients c := f (x ) cos jπ , for j ≤ N, basis functions T (x) := cos j arccos(x) j k j k=0 N N for x ∈[−1, 1] and j ≤ N, where indicates that the first and last summands are halved. The Chebyshev nodes x = cos π are displayed in Figs 7 and 9. For an arbitrary compact parameter interval [x, x], interpolation (A.1) needs is adjusted by the appropriate linear transformation. We consider the tensorized Chebyshev interpolation of functions f : [−1, 1] → R. For a more general hyperrectangular domain [x , x ]× ... ×[x , x ], appropriate linear transformations need to be 1 D 1 D applied. Let N := (N , ... , N ) with N ∈ N for i = 1, ... , D. The interpolation, which has (N +1) 1 D i 0 i i=1 Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 341 summands, is given by Cheby I (f )(x) := c T (x), (A.2) j j j∈J where the summation index j is a multi-index ranging over J :={(j , ... , j ) ∈ N : j ≤ N for 1 D i i i = 1, ... , D}.For j = (j , ... , j ) ∈ J the basis function is T (x , ... , x ) = T (x ). The coefficients 1 D j 1 D j i i=1 i c for j = (j , ... , j ) ∈ J are given by j 1 D N N D D 1 1 D {0<j <N } i i 2 k (k ,...,k ) 1 D c = ... f (x ) cos j π , (A.3) j i N N i i i=1 k =0 k =0 i=1 1 D k k where the Chebyshev nodes x for the multi-index k = (k , ... , k ) ∈ J are given by x = (x , ... , x ) 1 D k k 1 D with the univariate Chebyshev nodes x = cos π for k = 0, ... , N and i = 1, ... , D. k i i For the precise definition of the approximation with chebfun3 we refer to Hashemi & Trefethen (2017). http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

Parametric integration by magic point empirical interpolation

Loading next page...
1
 
/lp/ou_press/parametric-integration-by-magic-point-empirical-interpolation-ZZ386DqIf6

References (31)

Publisher
Oxford University Press
Copyright
Copyright © 2022 Institute of Mathematics and its Applications
ISSN
0272-4979
eISSN
1464-3642
DOI
10.1093/imanum/drx072
Publisher site
See Article on Publisher Site

Abstract

Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 IMA Journal of Numerical Analysis (2019) 39, 315–341 doi: 10.1093/imanum/drx072 Advance Access publication on December 25, 2017 Maximilian Gaß Center for Mathematics, Technische Universitat ¨ Munc ¨ hen, Munich, Germany Kathrin Glau Center for Mathematics, Technische Universitat ¨ Munc ¨ hen, Munich, Germany and School of Mathematical Sciences, Queen Mary University of London, London, UK Corresponding author: k.glau@qmul.ac.uk [Received on 30 November 2015; revised on 13 October 2017] We derive analyticity criteria for explicit error bounds and an exponential rate of convergence of the magic point empirical interpolation method introduced by Barrault et al. (2004, An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations. C. R. Math., 339, 667–672). Furthermore, we investigate its application to parametric integration. We find that the method is well-suited to Fourier transforms and has a wide range of applications in such diverse fields as probability and statistics, signal and image processing, physics, chemistry and mathematical finance. To illustrate the method, we apply it to the evaluation of probability densities by parametric Fourier inversion. Our numerical experiments display convergence of exponential order, even in cases where the theoretical results do not apply. The magic point integration performs considerably well, when compared with Clenshaw–Curtis quadrature and with Chebyshev interpolation in the parameter space. Keywords: parametric integration; Fourier transform; parametric Fourier inversion; magic point interpola- tion; empirical interpolation; online–offline decomposition. 1. Introduction At the basis of a large variety of mathematical applications lies the computation of parametric integrals of the form I(h ) := h (z)μ(dz) for all p ∈ P, (1.1) p p where Ω ⊂ C is a compact integration domain, P a compact parameter set, (Ω,A, μ) a measure space with μ(Ω) < ∞ and h : P × Ω → C a bounded function. A prominent class of such integrals are parametric Fourier transforms iz,x f (z) = e f (x) dx for all p = (q, z) ∈ P, (1.2) q q with a parametric family of complex functions f defined on the compact set Ω. Today, Fourier transforms lie at the heart of applications in optics, electrical engineering, chemistry, probability, partial differential equations, statistics and finance. The application of Fourier transforms even sparked many advancements in those disciplines, underlining its impressive power. The introduction of Fourier analysis to image © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 316 M. GASS AND K. GLAU formation theory in Duffieux (1946) marked a turning point in optical image processing, as outlined, for instance, in Stark (1982). The application of Fourier transform to nuclear magnetic resonance (NMR) in Ernst & Anderson (1966) was a major breakthrough in increasing the sensitivity of NMR, for which Richard R. Ernst was awarded the Nobel prize in Chemistry in 1991. Fourier analysis also plays a fundamental role in statistics, in particular in the spectral representation of time series of stochastic data (see Brockwell & Davis, 2002). For an elucidation of other applications impacted by Fourier theory, we recommend (Kammler, 2007, Appendix 1). For parametric Fourier integrals of form (1.2) with fixed value of q and z being the only varying parameter, efficient numerical methods have been developed based on discrete Fourier transform (DFT) and fast Fourier transform (FFT). The latter was developed by Cooley & Tukey (1965) to obtain the Fourier transform f (z) for a large set of values z simultaneously. With regard to (1.2), this is the special case of a parametric Fourier transform, where q is fixed and z varies over a specific set. The immense impact of FFT highlights the exceptional usefulness of efficient methods for parametric Fourier integration. Shifting the focus to other examples, parametric integrals arise as generalized moments of parametric distributions in probability, statistics, engineering and computational finance. These expressions often appear in optimization routines, where they have to be computed for a large set of different parameter constellations. Prominent applications are model calibration and statistical learning algorithms based on moment estimation, regression and expectation–maximization (see, e.g., Hastie et al., 2009). The efficient computation of parametric integrals is also a cornerstone of the quantification of parameter uncertainty in generalized moments of form (1.1), which leads to an integral of the form h (z)μ(dz)P(dp), (1.3) P Ω where P is a probability distribution on the parameter space P. In all of these cases, parametric integrals of the form (1.1) have to be determined for a large set of parameter values. Therefore, efficient algorithms for parametric integration have a wide range of applica- bility. In the context of reduced basis methods for parametric partial differential equations, a magic point empirical interpolation method has been developed by Barrault et al. (2004) to treat nonlinearities that are expressed by parametric integrals. The applicability of this method to the interpolation of parametric functions has been demonstrated by Maday et al. (2009). In this article, we focus on the approximation of parametric integrals by magic point empirical interpolation in general, a method that we call magic point integration and – provide sufficient conditions for an exponential rate of convergence of magic point interpolation and integration in its degrees of freedom, accompanied by explicit error bounds; – translate these conditions to the special case of parametric Fourier transforms, which shows the broad scope of the results; – empirically demonstrate efficiency of the method in a numerical case study; and – compare its performance with the Clenshaw–Curtis quadrature and with Chebyshev interpolation in the parameter space. An in-depth case study related to finance is presented in Gaß et al. (2017) and supports our empirical findings. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 317 1.1 Sketch of the main idea To sketch the idea, let us start with an observation that is as simple as it is crucial. If the interpolation M M M M I (f ) := β e (z) ≈ f (z) with coefficients β ∈ C and functions e : Ω → C is a good m=1 m m m m approximation of f in L (Ω) in terms of stability and order of convergence, then M M I (f )(z) := β e (z)μ(dz) (1.4) m m m=1 is a good approximation for f (z)μ(dz) in the sense that few number of summands M are sufficient to obtain a high accuracy. An approximation of type (1.4), however, is only numerically efficient, if this is M M the case for the computation of both the coefficients β and the weights w := e (z)μ(dz). m m In classical quadrature rules, the interpolation I (f ) is a polynomial interpolation, for instance, the Chebyshev interpolation of f , which results in the well-known Clenshaw–Curtis quadrature rule. In this case, the Chebyshev coefficients can be easily and efficiently computed, and, thirdly the basis functions e are polynomials whence the weights w are known explicitly. We notice that classical quadrature rules yield the same quadrature nodes and basis functions for all integrands f on a given domain Ω. Yet as we often have more insight into our integration problem, namely that the integrand stems from a particular set {h | p ∈ P}, it is natural to ask: How can we systematically tailor a quadrature rule to a specific set of integrands? The conceptual difficulty arises from the fact that we want to obtain a method that autonomously finds and exploits the specific features of the given set of integrands. Therefore, we propose a two-step procedure in which the first step is a learning algorithm. Here, we choose the magic point interpolation method of Barrault et al. (2004). In a recursive greedy procedure, it delivers a nodal interpolation of functions tailored a given set {h | p ∈ P}. The learning ∗ ∗ M M procedure after M ∈ N steps delivers both the magic points z , ... , z and the functions θ , ... , θ on 1 M 1 M the basis of which the magic point interpolation operator ∗ M I (h)(p, z) := h (z )θ (z) (1.5) M p m m m=1 is constructed to approximate all functions from the set {h | p ∈ P} particularly well in the L -norm. This allows us to define the magic point integration operator ∗ M I (h)(p) := h (z ) θ (z)μ(dz). (1.6) M p m m m=1 This integration routine inherits its favourable numerical properties from the respective properties of the magic point interpolation. We will discuss convergence results in detail. To efficiently employ operator (2.2) for numerical integration, we split the procedure in an offline and an online phase. ∗ ∗ In the offline phase, we iteratively construct the magic points z , ... , z along with the functions 1 M M M M M θ , ... , θ , thus the operator (1.5). Moreover, we precompute all weights w := θ (z) dz. 1 M m Ω m The online phase consist in evaluating the sum ∗ M h (z )w . (1.7) m m m=1 Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 318 M. GASS AND K. GLAU 1.2 Related literature Another algorithm for parametric integration has been proposed by Antil et al. (2013). Their two-step greedy algorithm for reduced order quadratures is similar to the algorithm we propose, since it also divides the problem into an offline and an online phase, where an essential part of the offline phase consists in an application of an empirical interpolation. In contrast to the method we propose here, the authors explicitly abstain from applying the original magic point interpolation method of Barrault et al. (2004), and the resulting algorithm of Antil et al. (2013) significantly differs from magic point integra- tion. The main reason is that they consider parametric scalar products of functions in a Hilbert space to obtain a reduced basis. Their starting point thus is the representation of the integrands in an orthonormal basis of a Hilbert space. Following the spirit of reduced order modelling, they first consider the pro- jection to a large but finite dimensional subspace. In the first step of the offline phase, a reduced basis of the approximation space is generated. The second step of the offline phase is a greedy procedure, where the discrete empirical interpolation (DEIM) extracts empirical integration points from the reduced basis. Philosophically, both approaches differ as well: Antil et al. (2013) employ a two-stage procedure. First, the integration problem is discretized. Secondly, the discrete problem is reduced by exploiting additional knowledge. In contrast, the magic point integration that we propose systematically exploits insight into the infinite dimensional integration problem directly. Similarly to the procedure of Antil et al. (2013), low-rank tensor decompositions for parameter- dependent integrals are constructed in two steps, first employing a discretization of choice, and second reducing the tenor structure (see Ballani, 2012; the survey article by Grasedyck et al., 2013 and the further references therein). The two approaches, empirical interpolation and low-rank tensor decompo- sitions, are closely interrelated as discussed in Bebendorf et al. (2014). In contrast to these approaches, magic point integration uses the actual function set of interest as an input into the learning procedure. No a priori discretization method is involved. Instead, it is the learning procedure itself that deliv- ers the discretization. Moreover, the article Bebendorf et al. (2014) discusses an interesting and close relation between magic point interpolation and the adaptive cross approximation. This relation can be directly transferred to our framework, i.e., when we compare our ansatz relying on the magic point inter- polation of the parametric integrands to the use of the adaptive cross approximation of the integrands instead. Monte Carlo methods for parametric integrals have been developed and investigated in Heinrich (1998), Heinrich & Sindambiwe (1999), by introducing the multilevel Monte Carlo method. The key idea of this approach is to combine interpolation of the integrand in the parameter space with Monte Carlo integration in a hierarchical way. In contrast to our approach, the interpolation is known and fixed beforehand and, thus, is not tailored to the parametric family of integrands, which is the key feature of magic point integration. The remainder of the article is organized as follows. In Section 2, we revisit the magic point empirical interpolation method from Barrault et al. (2004), related error bounds and present the related integral approximation, which can be perceived from two different perspectives. On the one hand, it delivers a quadrature rule for parametric integrals and, on the other, an interpolation method for parametric integrals in the parameter space. In Section 3, we provide analyticity conditions on the parametric integrals that imply exponential order of convergence of the method. We focus on the case of parametric Fourier transforms in the subsequent Section 4. In a case study, we discuss the numerical implementation and its results in Section 5. Here, we use the magic point integration method to evaluate densities of a parametric class of distributions that are defined through their Fourier transforms. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 319 2. Magic point empirical interpolation for integration We now introduce the magic point empirical interpolation method for parametric integration to approxi- mate parametric integrals of the form (1.1). Before we closely follow Maday et al. (2009) to describe the interpolation method, let us state our basic assumptions that ensure the well-definedness of the iterative procedure. Assumption 2.1 Let (Ω, ) and (P, ) be compact, P × Ω (p, z) → h (z) bounded and ∞ ∞ p p → h be continuous, i.e., for every sequence p → p we have h − h → 0. Moreover, there exists p i p p ∞ p ∈ P such that the function h is not constantly zero. ∗ ∗ M M For M ∈ N, the method delivers iteratively the magic points z , ... , z , functions θ , ... , θ and 1 M 1 M constructs the magic point interpolation operator ∗ M I (h)(p, z) := h (z )θ (z), (2.1) M p m m m=1 which allows us to define the magic point integration operator ∗ M I (h)(p) := h (z ) θ (z) dz. (2.2) M p m m m=1 ∗ ∗ In the first step, M = 1, we define the first magic parameter p , the first magic point z and the first 1 1 basis function q by p := arg max ∞ , (2.3) p L p∈P z := arg max|h (z)|, (2.4) z∈Ω h (·) q (·) := . (2.5) h (z ) Notice that, thanks to Assumption 2.1, these operations are well-defined. For M ≥ 1 and 1 ≤ m ≤ M we define M M −1 M ∗ θ (z) := (B ) q (z), B := q (z ), (2.6) j m m jm jm j j=1 −1 M M where we denote by (B ) the entry in the jth line and mth column of the inverse of matrix B .By jm construction, B is a lower triangular matrix with unity diagonal and is thus invertible (see Barrault et al., 2004). Then, recursively, as long as there are at least M linearly independent functions in {h | p ∈ P}, the algorithm chooses the next magic parameter p according to a greedy procedure. It selects the parameter M Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 320 M. GASS AND K. GLAU so that h is the function in the set {h | p ∈ P} which is worst represented by its approximation with p p the previously identified M − 1 magic points and basis functions. So the Mth magic parameter is p := arg max h − I (h)(p, ·) ∞ . (2.7) p M−1 L p∈P In the same spirit, select the Mth magic point as ∗ ∗ z := arg max h (z) − I (h)(p , z) . (2.8) p M−1 M M z∈Ω The Mth basis function q is the residual r , normed to 1 when evaluated at the new magic point, i.e., M M r (z) := h (z) − I (h)(p , z), (2.9) M p M−1 r (·) q (·) := . (2.10) r (z ) Notice the well-definedness of the operations in the iterative step, thanks to Assumption 2.1, and the fact that the denominator in (2.10) is zero only if all functions in {h | p ∈ P} are perfectly represented by the interpolation I , in which case they span a linear space of dimension M − 1 or less. M−1 Let us understand the approach from two different perspectives. On the one hand, it is an interpolation method: Remark 2.2 (Interpolation of parametric integrals in the parameters) For each m = 1, ... , M, the function M m θ is a linear combination of snapshots h for j = 1, ... , M with coefficients β , which may be iteratively m j computed. We thus may rewrite formula (2.2)as M M ∗ m I (h)(p) = h (z ) β h (z) dz. (2.11) M p p m j m=1 j=1 Thus, for an arbitrary parameter p ∈ P, the integral h (z) dz is approximated by a linear combination of integrals h (z) dz for the magic parameters p . In other words, I (h) is an interpolation of the p M function p → h (z) dz. Here, the interpolation nodes are the magic parameters p and the coefficients ∗ m are given by h (z )β . In contrast to classic interpolation methods, the ‘magic nodes’ are tailored m=1 m j to the parametric set of integrands. On the other hand, magic point integration is an approximate integration rule: Remark 2.3 (Empirical quadrature rule for integrating parametric functions) Magic point integration is an interpolation method for integrating a parametric family of integrands over a compact domain. From this point of view, the magic point empirical interpolation of Barrault et al. (2004) provides a quadrature rule for integrating parametric functions. The weights are θ (z) dz and the nodes are the magic points z for m = 1, ... , M. As in the adaptive quadrature rules, these quadrature nodes are not fixed in advance. Although the adaptive quadrature rules iteratively determine the next quadrature node in view Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 321 of integrating a single function, the empirical integration method (2.2) tailors itself to the integration of agiven parametrized family of integrands. A discrete approach to empirical interpolation has been introduced by Chaturantabut & Sorensen (2010). Although the empirical interpolation is designed for parametric functions, the input data for DEIM is discrete. A canonical way to use DEIM for integration is to first choose a fixed grid in Ω for a discrete integration of all parametric integrands. Then, DEIM can be performed on the integrands evaluated on this grid. In contrast, using magic point empirical interpolation for parametric integration separates the choice of the integration grid from the selection of nodes p . For each function, a different integration discretization may be used. Indeed, in our numerical study, we leave the discretization choices regarding integration to Matlab’s quadgk routine. Fixing a grid for discrete integration beforehand might, e.g., become advantageous when the domain of integration Ω is high dimensional. 3. Convergence analysis of magic point integration We therefore summarize the results the convergence results for the magic point interpolation in the L -norm, which clearly extend to the convergence results for the magic point integration method. Inter- estingly, these results relate the convergence of the magic point interpolation to the best linear n-term approximation that is formally expressed by the Kolmogorov n-width, consider the monograph of Pinkus (1985) for a general reference. For a real or complex normed linear space X , and U ⊂ X , the Kolmogorov n-width is given by d (U , X ) = inf sup inf g − f , (3.1) U ∈E(X ,n) f ∈U n n g∈U where E (X , n) is the set of all n-dimensional subspaces of X . This concept allows to express the convergence behaviour of the interpolation subject to the ‘true complexity’ of set of functions of interest, if one agrees that the Kolmogorov n-width is a good measure of the intrinsic complexity. The approach is especially worthwhile when tackling the curse of dimensionality in the parameter or in the integration space. In the following section, we present a posteriori and a priori error bounds, mostly from the litera- ture, which are relevant for the integration method proposed. To obtain precise convergence rates, the Kolmogorov n-width needs to be estimated. We devote Section 3.2 below to this task and obtain widely applicable sufficient conditions for an exponential convergence of the magic point integration method. For univariate parametric Fourier integrals, the conditions boil down to exponential moment conditions. 3.1 Error bounds for the magic point interpolation We first present an a posteriori error bound, which follows as a corollary from the results in Maday et al. (2016). Let U :={h | p ∈ P} and X , = L (Ω), . Moreover, the convergence statement p ∞ involves the Lebesgue constant Λ := sup θ (z) , (3.2) z∈Ω m=1 with θ from equation (2.6). m Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 322 M. GASS AND K. GLAU Proposition 3.1 Let Ω ⊂ C , P and h : P × Ω → C be such that Assumption 2.1 is satisfied and that the sequence (Λ ) is monotonically increasing. Then, the following assertions hold: M M≥1 ∞ −α (a) If there exist constants c > 0 and α> 1, such that d (U , L (Ω)) ≤ c n for all n ≥ 1, then for 0 n 0 all M ≥ 1 and all h ∈ U,wehave 3α 3 1−α h − I (h) ≤ 3 · 2 c (1 + Λ ) M . M ∞ 0 M ∞ −c n (b) If there exist constants c , c > 0 and α> 0, such that d (U , L (Ω)) ≤ c e for all n ≥ 1, then 0 1 n 0 for all M ≥ 1 and all h ∈ U,wehave √ √ 2 −c M h − I (h) ≤ 2c (1 + Λ ) Me , M ∞ 0 M −2α−1 where c = c 2 . 2 1 Proof. The generalized empirical interpolation method (GEIM) has been introduced in (Maday et al., 2016, Section 2) to generalize the concept of the magic point interpolation to general Banach spaces instead of L (Ω). Indeed, it is elementary to verify that the magic point interpolation can be interpreted as a GEIM and that the assumptions P1–P3 in Maday et al. (2016, Section 2) are satisfied in our case. Moreover, the inequalities (9) in Maday et al. (2016) are satisfied for η = 1 and F = F = span{q , ... , q } for n 1 n all n ∈ N. We therefore can apply in Maday et al. (2016, Theorem 13). We notice that the Lebesgue constant Λ appearing in the latter theorem is differently defined as the operator norm of I , namely by n n I (ϕ) n ∞ n := sup ∞ . It is, however, elementary to verify that ≤ Λ = sup θ (z) . n n n ϕ∈L (Ω) z∈Ω m m=1 Therefore, we can replace the constant Λ in the result of Maday et al. (2016, Theorem 13) with the Lebesgue constant from equation (3.2). Since Λ is monotonically increasing by assumption, we can apply Corollary 11 and Lemma 12 from Maday et al. (2016) and obtain the analogue of Corollary 14 in Maday et al. (2016). As the last step of our proof, we estimate the constants derived in the latter corollary to simplify the expression. This is an elementary step, namely writing n = 4 + k, we see that = 2( + ) ≤ n/2 + 2 ≤ 3/2n. Next, we present a variant of Maday et al. (2009, Theorem 2.4), which is an a priori error bound for the magic point empirical interpolation method: Proposition 3.2 Let Ω ⊂ C , P and h : P × Ω → C be such that Assumption 2.1 is satisfied. If there exist constants α> log(4) and c > 0 such that ∞ −αM d U , L (Ω) ≤ ce , c α then for arbitrary ε> 0 and C := e + ε we have for all u ∈ U that −(α−log(4))M u − I (u) ≤ CMe . (3.3) The proof follows analogous to the proof of Maday et al. (2009, Theorem 2.4) and is in detail provided in Gaß et al. (2015a). In contrast to Maday et al. (2009, Theorem 2.4), Proposition 3.2 expresses the result directly in terms of the Kolmogorov n-width, which we prefer. We emphazize, however, that this is only a minor difference and the scientific contribution originates from Maday et al. (2009). Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 323 3.2 A priori error bounds under analyticity conditions On the basis of in Maday et al. (2009, Theorem 2.4), we identify two generic cases for which this result implies exponential convergence of magic point interpolation and integration. In the first case, the set of parametric functions consists of univariate functions that are analytic on a specific Bernstein ellipse. In the second case, the functions may be multivariate and do not need to satisfy higher-order regularity. The parameter set now is a compact subset of R, and the regularity of the parameter dependence allows an analytic extension to a specific Bernstein ellipse in the parameter space. To formulate our analyticity assumptions, we define the Bernstein ellipse B([−1, 1], ) with parameter > 1 as the open region in the complex plane bounded by the ellipse with foci ±1 and semiminor and semimajor axis lengths summing up to . Moreover, we define for b < b ∈ R the generalized Bernstein ellipse by B([b, b], ) := τ ◦ B([−1, 1], ), (3.4) [b,b] where the transform τ : C → C is given by [b,b] b − b b − b τ (x) := b + (1 −(x)) + i (x) for all x ∈ C. (3.5) b,b 2 2 For Ω ⊂ R we set B(Ω, ) := B([inf Ω, sup Ω], ). (3.6) Definition 3.3 Let f : X × X → C with X ⊂ C for m = 1, 2. We say f has the analytic property 1 2 m with respect to X ⊂ C in the first argument,if X ⊂ X and f has an extension f : X × X → C such 1 1 1 1 1 2 that for all fixed x ∈ X the mapping x → f (x , x ) is analytic in the inner of X and 2 2 1 1 1 2 1 sup |f (x , x )| < ∞. (3.7) 1 1 2 (x ,x )∈X ×X 1 2 1 2 We define the analytic property of f with respect to X in the second argument analogously. We denote I(h , μ) := h (z)μ(dz) (3.8) p p and for all p ∈ P and M ∈ N, I (h, μ)(p) := I (h)(p, z)μ(dz), (3.9) M M Ω Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 324 M. GASS AND K. GLAU where I is the magic point interpolation operator given by the iterative procedure from Section 2. Hence, ∗ M I (h, μ)(p) = h (z ) θ (z)μ(dz), (3.10) M p m m m=1 ∗ M where z are the magic points and θ are given by (2.6) for every m = 1, ... , M. m m Theorem 3.4 Let Ω ⊂ C , P and h : P × Ω → C be such that Assumption 2.1 holds and one of the following conditions is satisfied for some > 4, (a) d = 1, i.e., Ω is a compact subset of R, and h has the analytic property with respect to B(Ω, ) in the second argument, (b) P is a compact subset of R, and h has the analytic property with respect to B(P, ) in the first argument. Then, there exists a constant C > 0 such that for all p ∈ P and M ∈ N, −M h − I (h) ≤ CM( /4) , (3.11) −M sup I(h , μ) − I (h, μ)(p) ≤ Cμ(Ω)M( /4) . (3.12) p M p∈P Proof. Assume (a). Thanks to an affine transformation, we may without loss of generality assume Ω =[−1, 1]. As assumed, (p, z) → h (z) has an extension f : P × B(Ω, ) → C that is bounded p 1 and for every fixed p ∈ P is analytic in the Bernstein ellipse B(Ω, ). We exploit the analyticity prop- erty to relate the approximation error to the best n-term approximation of the set U :={h | p ∈ P}. This can conveniently be achieved by inserting an example of an interpolation method that is equipped with exact error bounds, and we choose Chebyshev projection for this task. From Trefethen (2013, Cheby Theorem 8.2) we obtain for each N ∈ N and the interpolation I defined in the appendix as the explicit error bound, Cheby −N sup f (p, ·) − I f (p, ·) ≤ c(f ) , (3.13) 1 1 1 p∈P with constant c(f ) := sup f (p, z) . Next, we can apply the general convergence result 1 1 (p,z)∈P×B(Ω, ) −1 from Maday et al. (2009, Theorem 2.4). Consulting their proof, we realize that −M sup h − I (h)(p, ·) ≤ CM( /4) , p M p∈P c(f ) with C = . Equation (3.12) follows by integration. The proof follows analogously under assumption (b).  Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 325 Remark 3.5 The proof makes the constant C explicit. Under assumption (a) it is given by C = max f (p, z) . (3.14) (p,z)∈P×B(Ω, ) 2( − 1) Under assumption (b), an analogous constant can be derived. Remark 3.6 Consider a multivariate integral of the form I (h, μ ⊗ μ ) := h(p, z)μ (dz)μ (dp), (3.15) M 1 2 1 2 P Ω D d with compact sets P ⊂ R and Ω ⊂ R equipped with finite Borel-measures μ and μ . Then, the 1 2 application of the magic point empirical interpolation as presented in Section 2 yields ∗ M I (h, μ ⊗ μ ) := h (z )μ (dp) θ (z)μ (dz), (3.16) M 1 2 p 2 1 m m P Ω m=1 and, under the assumptions of Theorem 3.4, we obtain −M I(h, μ ⊗ μ ) − I (h, μ ⊗ μ ) ≤ Cμ (Ω)μ (P)M( /4) . (3.17) 1 2 M 1 2 1 2 Remark 3.7 Theorem 3.4 has an obvious extension to functions (p, z) → h (z) that are piecewise analytic. To be more precise, we can assume that Ω =∪ Ω with piecewise disjunct and compact sets i=1 Ω ⊂ C such that for each i = 1, ... , n the mapping P × Ω (p, z) → h (z) has the analytic property i i p w.r.t. B(Ω , ) in the second argument. The proof can be generalized to this setting, and the error bound i i needs to be adapted in the obvious way. If the domain Ω is one-dimensional and z → h (z) enjoys desirable analyticity properties, the approximation error of magic point integration decays exponentially in the number M of interpolation nodes, independently of the dimension of the parameter space P. Vice versa, if the parameter space is one-dimensional and p → h (z) enjoys desirable analyticity properties, the approximation error decays exponentially in the number M of interpolation nodes, independently of the dimension of the integration space Ω. If the function set {h , p ∈ P} has a small Kolmogorov n-width, magic point interpolation and integration do not suffer from the curse of dimensionality—in contrast to standard polynomial approximation. This is owed to the fact that classic interpolation methods are adapted to very large function classes, whereas magic point interpolation is adapted to a parametrized class of functions. In return, this highly appealing convergence comes with an additional cost: The offline phase involves a loop over optimizations to determine the magic parameters and the magic points. Although the online phase is independent of the dimensionality, the offline phase is thus affected. Let us further point out that Theorem 3.4 is a result on the theoretical iterative procedure as described in Section 2. As we will discuss in Section 5.1 below, the implementation inevitably involves additional problem simplifications and approximations to perform the necessary optimizations. In particular, instead of the whole parameter space a training set is fixed in advance. For this more realistic setting, rigorous a posteriori error bounds have been developed for the empirical interpolation method (see Eftang et al., 2010). These bounds Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 326 M. GASS AND K. GLAU rely on derivatives in the parameters and straightforwardly translate to an error bound for magic point integration. In our implementation, we obtain the discrete training set by simulating according to the uniform distribution on the parameter space. For high-dimensional parameter spaces, simulation of pseudo-random numbers with quasi-Monte Carlo methods is highly promising. These methods are well-established and often reduce the computational complexity considerably. Also, other discretization techniques such as sparse grids can be employed. Moreover, the optimization problem (2.7), (2.8) itself can be implemented using more advanced techniques. In particular, for the applications that we consider in Section 5 and in Gaß et al. (2015a), the derivatives of the parametric integrands are explicitly available, which enables the use of gradient-based optimization techniques. 4. Magic points for Fourier transforms In various fields of applied mathematics, Fourier transforms play a crucial role and parametric Fourier inverse transforms of the form −izx f (x) = e f (z) dz (4.1) q q 2π D+1 need to be computed for different values p = (q, x) ∈ P = Q × X ⊂ R , where the function izy z → f (z) := e f (y) dy is well-defined and integrable for all q ∈ Q. q q A prominent example arises in the context of signal processing, where signals are filtered with the help of short-time Fourier transform (STFT) that is of the form −izt Φ (f )(z) := f (t)w(t − b)e dt, (4.2) with a window function w and parameter b. One typical example for windows are the Gaußian win- 2 2 −t /(2σ ) dows, w (t) = √ e , where an additional parameter appears, namely the standard deviation σ 2πσ of the normal distribution. The STFT with a Gaußian window has been introduced by Gabor (1946). His pioneering approach has become indispensable for time–frequency signal analysis. We refer to Debnath & Bhatta (2014) for historical backgrounds. We consider the truncation of the integral in (4.1) to a compact integration domain Ω =[Ω, Ω]⊂ R and choose the same domain Ω for all p = (q, x) ∈ P = Q × X . Thus, we are left to approximate for all p = (q, x) ∈ P the integral −izx I(h ) := h (z) dz, where h (z) = h (z) = e f (z). (4.3) p p p (q,x) q 2π Then, according to Section 2, we consider the approximation of I(h ) by magic point integration, i.e., by ∗ M I (h)(p) := h (z ) θ (z) dz, (4.4) M p m m m=1 ∗ M where z are the magic points and θ are given by (2.6) for every m = 1, ... , M. In such cases where the m m analyticity properties of q → f (z) or of z → f (z) are directly accessible, Theorem 3.4 can be applied to q q Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 327 estimate the error sup |I(h ) − I (h)(q, x)|. The following corollary offers a set of conditions (q,x) M (q,x)∈Q×X in terms of existence of exponential moments of the functions f . (η+ε)|x| Corollary 4.1 For η> 15/8 (Ω − Ω) and every parameter q ∈ Q assume e f (x) dx < ∞ for some ε> 0 and assume further η|x| sup e f (x) dx < ∞. q∈Q Then, −M h − I (h) ≤ C(η)M (η)/4 , (4.5) −M sup I(h ) − I (h)(q, x) ≤ C(η)|Ω|M (η)/4 , (4.6) (q,x) M (q,x)∈Q×X with 2η 2η (η) = + + 1, |Ω| |Ω| (η) η(−X ∨ X ) η|x| C(η) = e sup e f (x) dx, 4π( (η) − 1) q∈Q where |Ω| := Ω − Ω, X = inf X and X = sup X . Proof. From the theorem of Fubini and the lemma of Morera, we obtain that the mappings z → f (z) have analytic extensions to the complex strip R + i[−η, η]. We determine the value of (η). It has to be chosen such that the associated semiminor b of the Bernstein ellipse B([−1, 1], ) is mapped by τ onto the semiminor of the generalized Bernstein ellipse B([Ω, Ω], ) that is of length η. Using the [Ω,Ω] relation b = (4.7) between the semiminor b of a Bernstein ellipse and its ellipse parameter we solve Ω − Ω τ (ib ) = b = η (4.8) [Ω,Ω] for and find that for 2η 2η (η) = + + 1 (4.9) Ω − Ω Ω − Ω the Bernstein ellipse B(Ω, (η)) is contained in the strip of analyticity of f and by the choice of η we have (η) > 4. Hence, assumption (a) of Theorem 3.4 is satisfied. With regard to Remark 3.5, we also obtain the explicit form of the constant C(η), which proves the claim.  Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 328 M. GASS AND K. GLAU Similarly, we consider the case where Q is a singleton and −izx h (z) = h (z) = e f (z). (4.10) p x 2π Under this assumption we obtain an interesting additional assertion if the integration domain Ω is rather small. This case occurs, e.g., for approximations of STFT. Corollary 4.2 For η> 15/8 (X − X ),wehave −M h − I (h) ≤ C(η)M (η)/4 , (4.11) −M sup I(h ) − I (h)(x) ≤ C(η)|Ω|M (η)/4 , (4.12) x M x∈X with 2η 2η (η) = + + 1, |X | |X | (η) η(−Ω ∨ Ω) C(η) = e sup f (z) , 4π( (η) − 1) z∈Ω where |X | := X − X , X = inf X and X = sup X . The proof of the corollary follows similarly to the proof of Corollary 4.1. 5. Case study 5.1 Implementation and complexity The implementation of the algorithm requires further discretizations. For the parameter selection, we disc replace the continuous parameter space P by a discrete parameter cloud P . Additionally, for the disc selection of magic points we replace Ω by a discrete set Ω , instead of considering the whole continuous disc domain. Each function h is then represented by its evaluation on this discrete Ω and is thus represented by a finite-dimensional vector, numerically. The crucial step during the offline phase is finding the solution to the optimization problem ∗ M (p , z ) := arg max h (z) − h (z )θ (z) . (5.1) M+1 M+1 p p m m disc disc (p,z)∈P ×Ω m=1 In a continuous setup, problem (5.1) is solved by optimization routines. In our discrete implementation, disc however, we are able to consider all magic parameter candidates in the discrete parameter space P disc and all magic point candidates in the discrete domain Ω to solve problem (5.1). This results in a disc disc complexity of O(M ·|P |·|Ω |) during the offline phase for identifying the basis functions and magic points of the algorithm. We comment on our choices for the dimensionality of the involved discrete sets in the numerical section, later. Before magic point integration can be performed, the quantities θ (z) dz need to be computed for m = 1, ... , M. The complexity of this final step during the offline phase depends on the number of integration points that the integration method uses. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 329 5.2 Tempered stable distribution We test the parametric integration approach on the evaluation of the density of a tempered stable distri- bution as an example. This is a parametric class of infinitely divisible probability distributions. A random variable X is called infinitely divisible, if for each n ∈ N there exist n independent and identically dis- tributed random variables whose sum coincides with X in distribution. The rich class of infinitely divisible distributions, containing, for instance, the Normal and the Poisson distribution, is characterized through their Fourier transforms by the Levy–Khintchine ´ representation. Using this link, a variety of infinitely divisible distributions is specified by Fourier transforms. The class of tempered stable distributions is an example of this sort, namely it is defined by a parametric class of functions, which by the theorem of Levy–Khintchine ´ are known to be Fourier transforms of infinitely divisible probability distributions. Its density function, however, is not explicitly available. To evaluate it, a Fourier inversion has to be performed numerically. As we show below, magic point integration can be an adequate way to efficiently compute this Fourier inversion for different parameter values and at different points on the density’s domain. Following Kuchler & Tappe (2014), parameters + + − − + − α , λ , α , λ ∈ (0, ∞), β , β ∈ (0, 1) (5.2) define the distribution η on (R, B(R)) that we call a tempered stable distribution and write + + + − − − η = TS(q) = TS(α , β , λ ; α , β , λ ), (5.3) izx if its characteristic function ϕ (z) := e η (dx) is given by q q izx ϕ (z) = exp (e − 1)F (dx) , z ∈ R, (5.4) q q with Levy ´ measure + − α + α − −λ x −λ x F (dx) = e 1 + e 1 dx. (5.5) q x∈(0,∞) x∈(−∞,0) + − 1+β 1+β x x + − The tempered stable distribution is also defined for β , β ∈ (1, 2), where the characteristic function is given by a similar expression as in (5.4). We consider the tempered stable distribution in a special case by introducing the parameter space Q ={(C, G, M, Y ) | C > 0, G > 0, M > 0, Y ∈ (1, 2)} (5.6) and setting + − − + + − α = α = C, λ = G, λ = M, β = β = Y . (5.7) The resulting distribution γ = CGMY(C, G, M, Y ) is also known as CGMY distribution and is well established in finance (see Carr et al., 2002). The density function of a tempered stable or a CGMY As usual, we denote the parameters of the CGMY distribution with C, M, G, Y as a reference to the authors of Carr et al. (2002). It will be clear from the context whether M denotes the number of magic points or the CGMY parameter. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 330 M. GASS AND K. GLAU Table 1 Parameter intervals for the numerical study. The parameter Y is kept constant CG M Y x interval [1, 5][1, 8][1, 8] 1.1 [−1, 1] disc disc Fig. 1. Decay of the error (2.8)on P and Ω during the offline phase of the algorithm. distribution, respectively, is not known in the closed form. With q = (C, G, M, Y ) ∈ Q of (5.6) and z ∈ R, its Fourier transform, however, is explicitly given by Y Y Y Y ϕ (z) = exp CΓ(−Y ) (M − iz) − M + (G + iz) − G , (5.8) wherein Γ denotes the gamma function. By Fourier inversion, the density f can be evaluated, 1 1 −izx −izx f (x) = e ϕ (z) dz =  e ϕ (z) dz for all x ∈ R. (5.9) q q q 2π π R 0 5.3 Numerical results We restrict the integration domain of (5.9)to Ω =[0, 65] and apply the parametric integration algorithm on a subset of P = Q × R, (5.10) specified in Table 1. We draw 4000 random samples from P respecting the intervals bounds of Table 1 ∞ −12 and run the offline phase until an L accuracy of 10 is reached. Here, we compute the benchmark −14 values by numerical integration using Matlab’s quadgk routine with a tolerance 10 . The error decay during the offline phase is displayed by Fig. 1. We observe exponential error decay reaching the prescribed accuracy threshold at M = 40. Let us have a closer look at some of the basis functions q that are generated during the offline phase. Figure 2 depicts five basis functions that have been created at an early stage of the offline phase (left) Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 331 Fig. 2. Some basis functions q constructed early (left) and late (right) during the offline phase of the algorithm. together with five basis functions that have been identified towards the end of it (right). Note that all plotted basis functions are supported on a relatively small subinterval of Ω. They are numerically zero outside of it. Interestingly, the functions in both sets have intersections with the Ω axis rather close to the origin. Such intersections mark the location of magic points. Areas on Ω where these intersections accumulate thus reveal regions where the approximation accuracy of the algorithm is low. In these regions, the shapes across all parametrized integrands seem to be most diverse. In the right picture, we observe in comparison that at a later stage in the offline phase, magic points z are picked further away from the origin, as well. We interpret this observation by assuming that the more magic points are chosen close to the origin, the better the algorithm is capable of approximating integrands more precisely there. Thus, its focus shifts towards the right for increasing values of M where now smaller variations attract its attention. We now test the algorithm on parameters not contained in the training set. To this extent, we randomly draw 1000 parameter sets p = (C , G , M , Y , x ), i = 1, ... , 1000, uniformly distributed from the i i i i i i intervals given by Table 1. For each p , i = 1, ... , 1000, we evaluate f (x ) by Fourier inversion i (C ,G ,M ,Y ) i i i i i using Matlab’s quadgk routine. Additionally, we approximate f (x ) using the interpolation (C ,G ,M ,Y ) i i i i i operator I for all values m = 1, ... , M. For each such m we compute the absolute and the relative error. The maximal absolute error over the 1000 samples, i.e., the L ((p ) )-error for each i i=1,...,1000 m = 1, ... , M, is illustrated by Fig. 3. To be precise, we compute the benchmark values by numeri- cal Fourier inversion, restrict the resulting integration domain to Ω =[0, 65] and integrate numerically −14 with Matlab’s quadgk routine with a tolerance of 10 . The set-up of the numerical experiment does not fall in the scope of our theoretical result in several respects. We discretized both the parameter space and the set of magic point candidates. While Table 1 ensures a joint strip of analyticity R + i(−1, 1) that all integrands h share, a Bernstein ellipse B(Ω, ) with > 4 on which those integrands are analytic does not exist. In spite of all these limitations, we still observe that the exponential error decay of the 3 disc disc offline phase is maintained. From this, we draw two conclusions. First, both Ω and P appear This observation is confirmed by our empirical studies in Gaß et al. (2017), where the existence of such a shared strip of analyticity was sufficient for empirical exponential convergence. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 332 M. GASS AND K. GLAU Fig. 3. Out of sample error decay. For 1000 randomly drawn parameters (C , G , M , Y , x ), i = 1, ... , 1000 from the set prescribed i i i i i in Table 1, the CGMY density is evaluated by the parametric integration method and by numerical Fourier inversion via Matlab’s quadgk routine. For each number M of magic points, the maximal absolute difference between both is displayed. Fig. 4. All out of sample errors in detail. Left: The absolute errors achieved for each of the 1000 randomly drawn parameter sets. −3 Right: The relative errors for density values above 10 . sufficiently rich to represent their continuous counterparts. Secondly, the practical use of the algorithm extends beyond the analyticity conditions imposed by Theorem 3.4. Figure 4 presents the absolute errors and the relative errors for each parameter set p , individually, for the maximal number M magic points. Note that only relative errors for CGMY density values larger −3 than 10 have been plotted to exclude the influence of numerical noise. While individual outliers are not present among the absolute deviations, they dominate the relative deviations in contrast. This result is in line with the objective function that the algorithm minimizes. Finally, Fig. 5 displays all magic parameters identified during the offline phase, together with those randomly drawn parameter constellations that resulted in the ten smallest absolute errors (crosses) and the 10 largest absolute errors (stars). We observe that orange parameter constellations are found in areas densely populated by magic parameters, whereas green parameter constellations often do not have magic parameters in their immediate neighbourhood, at all. This result may surprise at first. On second thought, however, we understand that areas where magic parameters accumulate mark precisely those locations in the parameter space where approximation is especially challenging for the algorithm. Integrands associated with parameter constellations from there exhibit the largest variation. Approximation for white Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 333 5 8 8 1 3 4.5 4.5 0 1 1 1 -1 135 135 135 135 8 8 8 1 4.5 4.5 4.5 0 1 1 1 -1 135 1 4.5 8 14.5 8 14.5 8 8 8 8 1 4.5 4.5 4.5 0 1 1 1 -1 135 1 4.5 8 14.5 8 14.5 8 1 1 1 1 0 0 0 0 -1 -1 -1 -1 135 1 4.5 8 14.5 8 -1 01 Fig. 5. An illustration of those randomly drawn parameter constellations that resulted in the 10 worst (stars) and the 10 best (circles) absolute errors. The crosses (“+” signs) the magic parameters selected during the offline phase. areas, on the other hand, is already covered by the previously chosen magic parameters. Consequently, crosses are very likely to be found there. 5.4 Comparison of magic point integration and Clenshaw–Curtis Since in our numerical experiments, we compute univariate analytic integrands over a finite domain, it is interesting to compare the performance of the new integration routine to the Clenshaw–Curtis method. The latter is well-known for its excellence performance, particularly in case of analyticity. To compare both methods, we use again the setting introduced in Section 5.2. We randomly sample 1000 parameter constellations according to the uniform distribution in the set specified in Table 1.For these, we compute the benchmark integrals with Matlab’s quadgk routine for an absolute tolerance −14 of 10 and the maximal difference with respect to the results obtained by both the Clenshaw–Curtis quadrature rules and the magic point integration, for varying numbers of quadrature nodes, respectively, magic points, M. Figure 6 displays the resulting error decays. Here, the magic point integration method clearly out- performs the Clenshaw–Curtis quadrature rule. For less than 35 magic points, the accuracy of magic −10 point integration reaches already 10 , while the error of the Clenshaw–Curtis quadrature is still in the −12 second digit. The range of machine precision, i.e., an error of 10 , is achieved with about 40 magic points while about 200 Chebyshev nodes are required in the Clenshaw–Curtis routine. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 334 M. GASS AND K. GLAU Fig. 6. Comparison of the error decay of magic point integration versus Clenshaw–Curtis. The maximal absolute errors achieved for each of the 1000 randomly drawn parameter sets for varying numbers M of quadrature nodes. Left: Performance of both magic point integration and Clenshaw–Curtis for M = 1, ... , 41 quadrature nodes. Right: Performance of Clenshaw–Curtis for M = 1, ... , 200 quadrature nodes. Here, the magic point integration clearly outperforms the Clenshaw–Curtis quadrature. In accordance with the theoretical results, we observe an exponential error decay for both methods. For magic point integration, the rate of convergence, however, is much higher. This is not reflected by the current theoretical results. Indeed, in the proof of Theorem 3.4, we estimate the Kolmogorov n-width by the Chebyshev interpolation, i.e., by the Clenshaw–Curtis quadrature. The theoretical bound that we obtained is therefore larger than the error bound for the Clenshaw–Curtis quadrature. The numerical experiments reveal that the magic point interpolation is actually able to reduce the computational complexity by adapting a quadrature rule to the parametric family of integrands in the offline phase. Moreover, this gives evidence that the Kolmogorov n-width of the parametric family of integrands is indeed significantly smaller than the error obtained by the interpolation with Chebyshev polynomials. Figure 7 displays the magic points versus the Chebyshev nodes, scaled to the unit interval. It is interesting to see that there is not even a qualitative similarity between the two different sets of points. 5.5 Comparison of magic point integration and Chebyshev interpolation in the parameter space As highlighted in Remark 2.3 and Remark 2.2, magic point integration can be seen both as a quadrature rule and as an interpolation method in the parameter space. Therefore, we also compare the performance of the magic point integration with the Chebyshev interpolation in the parameter space. In our numerical experiments reported up to this point, we considered a five-dimensional parameter space as laid in Table 1. Tensorized Chebyshev interpolation, however, is exposed to the curse of dimensionality, and therefore magic point integration will clearly outperform this method. To gain more insight, we reduce the parameter space to a univariate, a bivariate and a three-dimensional setting, below. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 335 Fig. 7. Magic points versus Clenshaw–Curtis points for M = 40. 5.5.1 Comparison for a single free parameter. We fix the parameters C = 1, M = 4, Y = 1.1, x =−1, and only vary G ∈[1, 8]. For a uniform test grid of 100 parameters in [1, 8], we compute the −14 benchmark integrals with Matlab’s quadgk routine for an absolute tolerance of 10 , and the maximal difference with respect to the results obtained by both Chebyshev interpolation in the parameter G and magic point integration for a varying number of magic points and Chebyshev points M, respectively. We refer to the appendix for a precise definition of the Chebyshev interpolation that we implemented. Figure 8 displays the resulting error decays. Both methods perform comparably well. As expected, the convergence is exponential in the number of nodes. The magic point integration shows a slightly higher rate of decay as the Chebyshev interpolation. Figure 9 displays the first 10 magic parameters versus the first 10 Chebyshev nodes in G, scaled to the unit interval. Remarkably, we recognize similarities in the choices of points, as both sets cluster at the boundaries. We also recognize significant differences. The magic points first of all are not located at the Chebyshev points. Secondly, there is an asymmetry in their distribution with more points on the left side than on the right. 5.5.2 Comparison for two free parameters. We fix the parameters C = 1, M = 4, Y = 1.1, and vary G ∈[1, 8] and x ∈[−1, 1]. For a uniform test grid of 100 × 100 parameters in [1, 8]×[−1, 1],we −14 compute the benchmark integrals with Matlab’s quadgk routine for an absolute tolerance of 10 and the maximal difference with respect to the results obtained by both tensorized Chebyshev interpolation in G, x, and magic point integration for a varying number of magic points and Chebyshev points M, respectively. For the precise definition of the tensorized Chebyshev interpolation, we again refer to the appendix. Since we are using the tensorized Chebyshev interpolation, we have a grid of Chebyshev nodes. We specify N number of nodes in both directions, G and x, and obtain a total number of M = N number of nodes. Figure 10 displays the error decay of both the magic point integration and the tensorized Chebyshev interpolation in G, x. The magic point integration outperforms the Chebyshev interpolation by far. For −12 an accuracy in the range of machine precision, i.e., 10 , less than 25 magic parameters versus over 700 Chebyshev nodes are required. This is an empirical indicator that magic point integration has the potential to tackle the curse of dimensionality. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 336 M. GASS AND K. GLAU Fig. 8. Comparison of the error decay of magic point integration versus Chebyshev interpolation in the parameter space. The maximal absolute errors achieved for each of the 100 equidistant grid points in the parameter set G ∈[1, 8], for varying numbers M of interpolation nodes. Fig. 9. Magic parameters versus Chebyshev interpolation nodes for M = 19 interpolation nodes. Next, we investigate the performance in relation to the dimensionality of the parameter space in some more depth. Figure 10 shows that the convergence rate of magic point integration is considerably higher. Yet, it is not clear how much the method is still exposed to the curse of dimensionality in the parameter space. From Fig. 8, we may conclude that both methods perform comparably in the univariate setting. Therefore, we next display the error decays of the Chebyshev interpolation in a different scale. While for magic point integration we show the error for the number of magic points M, we show the error of the Chebyshev interpolation for N, the number of nodes chosen in each of the interpolation axes. The result is provided in Fig. 11. The similarity of both curves is striking. For M = N = 15, we observe roughly the same error of about −8 10 . To be clear, this means that 15 summands have been used for the magic point integration, to obtain Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 337 Fig. 10. Comparison of the error decay of magic point integration versus Chebyshev interpolation in the parameter space. The maximal absolute errors achieved for each of the 100 × 100 equidistant grid points in the parameter set (G, x) ∈[1, 8]×[−1, 1], for varying numbers M of interpolation nodes. Fig. 11. Comparison of the error decay of magic point integration versus Chebyshev interpolation in the parameter space; specifi- cations as in Fig. 10. For the Chebyshev interpolation the number of summands equals (N + 1) , while the number of summands in the magic point interpolation is M. This figure visualizes that the magic point integration seen as interpolation in the parameter space is not exposed to the curse of dimensionality—in contrast to the tensorized Chebyshev interpolation. −8 4 2 a maximal error of 10 over the 10 test parameter. For the same precision, (N + 1) = 256 summands are required for the tensorized Chebyshev interpolation in G, x. In other words, in these experiments, the magic point interpolation exhibits indeed no curse of dimensionality in the parameter space. Figure 12 displays the grid of tensorized Chebyshev nodes and of magic parameters for which roughly −12 the same accuracy of 10 has been achieved. This visualizes the complexity reduction of magic points versus Chebyshev interpolation. In contrast to the example of the univariate parameter space, the distri- bution of the magic parameters is totally different from the Chebyshev points. It is, moreover, apparent that they only appear at the boundaries of the domain. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 338 M. GASS AND K. GLAU Fig. 12. M = 23 Magic parameters versus (N + 1) Chebyshev interpolation nodes for N = 28, for which both methods achieve roughly the same accuracy. 5.5.3 Comparison for three free parameters with low-rank tensor techniques. In this section, we consider three varying parameters. We have already seen for the bivariate parameter space that the tensorized Chebyshev interpolation is exposed to the curse of dimensionality, in contrast to the magic point empirical interpolation. Low-rank tensor techniques start from the tensorized Chebyshev interpolation and exploit its compressible structure by a greedy procedure, which is similar to the offline phase in the magic point interpolation. It is, therefore, very interesting to investigate the performance of low-rank tensor techniques in this example. To do so, we use the chebfun toolbox for Matlab, see http://www.chebfun.org (accessed 16 November 2017), developed by Prof. Trefethen and his working group, particularly we use chebfun3 written by Benham Hashemi (see also Hashemi & Trefethen, 2017). We fix the parameters C = 1, Y = 1.1 and vary (G, M, x) ∈[1, 8]×[1, 8]×[−1, 1]. We uniformly sample a set of 1000 test parameters from [1, 8]×[1, 8]×[−1, 1], compute the benchmark integrals with −14 Matlab’s quadgk routine for an absolute tolerance of 10 and the maximal difference with respect to the results obtained by chebfun3 from in G, M, x for a number of absolute tolerances, namely for −1−k/2 (10 ) , and the magic point integration for varying numbers of magic points M. To compare the k=0,...,22 complexity of both interpolation methods, we compute the number M of summands of each chebfun3 object. Moreover, chebfun3 allows to extract the number M of summands required for the same precision without the use low-rank compression techniques. Figure 13 displays the error achieved by chebfun3, magic point integration for the number of summands M each. The third line, which is dashed, shows the error for chebfun3 for the number of summands that would be required for the tensorized Chebyshev interpolation. The figure shows that the low-rank tensor technique performs considerably better compared to the −2 tensorized Chebyshev interpolation. For instance, chebfun3 achieves an accuracy of about 10 with less than 200 summands, while the interpolation based on the full tensor structure requires about 500 summands. The performance of magic point integration outperforms both methods by far. This can be Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 339 Fig. 13. Magic integration and Chebyshev interpolation three varying parameters G, M, x. explained by the fact that the accuracy of the tensorized Chebyshev interpolation is determined by the resolution along the axes, which carries over to the low-rank approximations. This a priori choice of axes, however, is not present for the magic point interpolation. Instead, the set of basis functions is iteratively produced for the set of functions of interest. Acknowledgements The authors thank Bernard Haasdonk, Laura Iapichino, Daniel Wirtz and Barbara Wohlmuth for fruitful discussions, and two anonymous referees for their valuable comments and suggestions. Maximilian Gaß thanks the KPMG Center of Excellence in Risk Management and Kathrin Glau acknowledges the TUM Junior Fellow Fund for financial support. References Antil, H., Field, S. E., Herrmann, F., Nochetto, R. H. & Tiglio, M. (2013) Two-step greedy algorithm for reduced order quadratures. J. Sci. Comput., 57, 604–637. Ballani, J. (2012) Fast evaluation of singular BEM integrals based on tensor approximations. Numer. Math., 121, 433–460. Barrault, M., Maday, Y., Nguyen, N. C. & Patera, A. T. (2004) An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations. C. R. Math., 339, 667–672. Bebendorf M., Maday Y. & Stamm B. (2014) Comparison of Some Reduced Representation Approximations. In: A. Quarteroni & G. Rozza (eds) Reduced Order Methods for Modeling and Computational Reduction. MS&A - Modeling, Simulation and Applications, vol 9. Cham: Springer, pp. 67–100. Brockwell, P. J. & Davis, R. A. (2002) Introduction to Time Series and Forecastung, 2nd edn. New York: Springer. Carr, P., Geman, H., Madan, D. B. & Yor, M. (2002) The fine structure of asset returns: An empirical investigation. J. Bus., 75, 305–333. Chaturantabut, S. & Sorensen, D. C. (2010) Nonlinear model reduction via discrete empirical interpolation. SIAM J. Sci. Comput., 32, 2737–2764. Cooley, J. W. & Tukey, J. W. (1965) An algorithm for the machine calculation of complex Fourier series. Math. Comput., 19, 297–301. Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 340 M. GASS AND K. GLAU Debnath, L. & Bhatta, D. (2014) Integral Transforms and Their Applications, 3rd edn. Boca Raton, USA: Apple Academic Press Inc. Duffieux, P.-M. (1946) L’Integr ´ ale de Fourier et ses Applications a ` l’Optique. Rennes: Oberthur. Eftang, J. L., Grepl, M. A. & Patera, A. T. (2010) A posteriori error bounds for the empirical interpolation method. Ser. I, 348, 575–579. Ernst, R. R. & Anderson, W. A. (1966) Application of Fourier transform spectroscopy to magnetic resonance. Rev. Sci. Instrum., 37, 93–102. Gabor, D. (1946) Theory of communication. J. Inst. Electr. Eng. III Radio Commun. Eng., 93, 429–457. Gaß, M., Glau, K., Mahlstedt, M. & Mair, M. (2015a) Chebyshev interpolation for parametrized option pricing. Working Paper. Available at URL http://arxiv.org/abs/1505.04648. Accessed 16 November 2017. Gaß, M., Glau, K. & Maikr, M. (2017) Magic points in finance: empirical integration for parametrized option pricing, SIAM J. Finan. Math., 8, 766–803. Grasedyck, L., Kressner, D. & Tobler, C. (2013) A literature survey of low-rank tensor approximation technique. GAMM-Mitt., 36, 53–78. Hashemi, B. & Trefethen, L. N. (2017) Chebfun in three dimensions. SIAM J. Sci. Comput., 39, C341–C363. Hastie, T., Tibshirani, R. & Friedman, J. (2009) The Elements of Statistical Learning, 2nd edn. Heidelberg: Springer. Heinrich, S. (1998) Monte Carlo complexity of global solution of integral equations. J. Complex., 14, 151–175. Heinrich, S. & Sindambiwe, E. (1999) Monte Carlo complexity of parametric integration. J. Complex., 15, 317–341. Kammler, D. W. (2007) A First Course in Fourier Analysis, 2nd edn. New York: Cambridge University Press. K¨ uchler, U. & Tappe, S. (2014) Exponential stock models driven by tempered stable processes. J. Econometrics, 181, 53–63. Maday, Y., Mula, O. & Turinici, G. (2016) Convergence analysis of the generalized empirical interpolation method. SIAM J. Numer. Analy., 54, 1713–1731. Maday, Y., Nguyen, C. N., Patera, A. T. & Pau, G. S. H. (2009) A general multipurpose interpolation procedure: the magic points. Comm. Pure Appl. Anal., 8, 383–404. Pinkus, A. (1985) n-widths in approximation theory. Ergebnisse der Mathematik und ihrer Grenzgebiete [Results in Mathematics and Related Areas (3)], Vol. 7. Heidelberg: Springer. Stark, H. (1982) Theory and measurement of the optical Fourier transform. Applications of Optical Fourier Transforms (H. Stark ed). New York: Academic Press, pp. 2–40. Trefethen, L. N. (2013) Approximation Theory and Approximation Practice. New York: SIAM Books. Appendix Since different versions of the Chebyshev interpolation are used in the literature, we introduce the version that we use. For interpolation univariate functions f : [−1, 1]→ R, the interpolation with Chebyshev polynomials of degree N is of the form Cheby I (f )(x) := c T (x) (A.1) j j j=0 0<j<N 2 k with coefficients c := f (x ) cos jπ , for j ≤ N, basis functions T (x) := cos j arccos(x) j k j k=0 N N for x ∈[−1, 1] and j ≤ N, where indicates that the first and last summands are halved. The Chebyshev nodes x = cos π are displayed in Figs 7 and 9. For an arbitrary compact parameter interval [x, x], interpolation (A.1) needs is adjusted by the appropriate linear transformation. We consider the tensorized Chebyshev interpolation of functions f : [−1, 1] → R. For a more general hyperrectangular domain [x , x ]× ... ×[x , x ], appropriate linear transformations need to be 1 D 1 D applied. Let N := (N , ... , N ) with N ∈ N for i = 1, ... , D. The interpolation, which has (N +1) 1 D i 0 i i=1 Downloaded from https://academic.oup.com/imajna/article/39/1/315/4774576 by DeepDyve user on 19 July 2022 PARAMETRIC INTEGRATION WITH MAGIC POINTS 341 summands, is given by Cheby I (f )(x) := c T (x), (A.2) j j j∈J where the summation index j is a multi-index ranging over J :={(j , ... , j ) ∈ N : j ≤ N for 1 D i i i = 1, ... , D}.For j = (j , ... , j ) ∈ J the basis function is T (x , ... , x ) = T (x ). The coefficients 1 D j 1 D j i i=1 i c for j = (j , ... , j ) ∈ J are given by j 1 D N N D D 1 1 D {0<j <N } i i 2 k (k ,...,k ) 1 D c = ... f (x ) cos j π , (A.3) j i N N i i i=1 k =0 k =0 i=1 1 D k k where the Chebyshev nodes x for the multi-index k = (k , ... , k ) ∈ J are given by x = (x , ... , x ) 1 D k k 1 D with the univariate Chebyshev nodes x = cos π for k = 0, ... , N and i = 1, ... , D. k i i For the precise definition of the approximation with chebfun3 we refer to Hashemi & Trefethen (2017).

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Jan 25, 2019

There are no references for this article.