Multiscale Discrete Approximation of Fourier Integral Operators

Abstract

We develop a discretization and computational procedures for approximation of the action of Fourier integral operators the canonical relations of which are graphs. Such operators appear, for instance, in the formulation of imaging and inverse scattering of seismic reflection data. Our discretization and algorithms are based on a multiscale low-rank expansion of the action of Fourier integral operators using the dyadic parabolic decomposition of phase space and on explicit constructions of low-rank separated representations using prolate spheroidal wave functions, which directly reflect the geometry of such operators. The discretization and computational procedures connect to the discrete almost symmetric wave packet transform. Numerical wave propagation and imaging examples illustrate our computational procedures. F. ANDERSSON, M. V. DE HOOP, AND H. WENDT 3d−1 2 log(N )), or O(D N d log(N )) if D is the number of significant tiles in the dyadic parabolic decomposition of u, valid in arbitrary finite dimension d. Our separated representation is expressed in terms of geometric attributes of the canonical relation of the FIO. We make use of prolate spheroidal wave functions (PSWFs) in connection with the dyadic parabolic decomposition, while the propagation of singularities or canonical transformation is accounted for via an unequally spaced FFT (USFFT). The use of PSWFs was motivated by the work of Beylkin and Sandberg [7] and the proposition of an efficient algorithm for their numerical evaluation by Xiao, Rokhlin, and Yarvin [62]. We note that it is also possible to obtain low-rank separated representations of the complex exponential in (1.1) purely numerically at the cost of losing the explicit relationship with the geometry. The algorithm presented here can be applied to computing parametrices of hyperbolic evolution and wave equations; we show that then our approximation corresponds to the solution of the paraxial wave equation in curvilinear coordinates, MULTISCALE DISCRETE APPROXIMATION OF FIOs 113 i.e., directionally developed relative to the central wave vector. However, it also forms the basis of a computational procedure following the construction of weak solutions of Cauchy initial value problems for the wave equation if the medium is C 2,1 , in which, in addition, a Volterra equation needs to be solved (de Hoop et al. [18]).We derive our discretization from the (inverse) transform based on discrete almost symmetric wave packets [25]. The connection of our algorithm to discrete almost symmetric wave packets is important in imaging and inverse scattering applications, where the FIOs act on data (u in the above). The wave packets can aid in regularizing the data from a finite set of samples through sparse optimization (instead of standard interpolation, for example) [14,15,17,58]. Moreover, the mentioned connection enables multiscale imaging and, in the context of directional pointwise regularity analysis [1,30,32,33,34,35,42], the numerical estimation and study of propagation of scaling exponents by the ...

Multiscale Discrete Approximations of Fourier Integral Operators Associated with Canonical Transformations and Caustics

“…It is based on the discretization and approximation, to accuracy O(2 k/2 ), of the action ofF ij on a wave packet ' j,⌫,k (x), Expansion of the cutoff functions. To numerically evaluate (3.6) in reasonable time, we will use separated (in y andξ) representations ofǍ ij (y,ξ) andŠ ij (y,ξ) [2,10,11]. Such representations can be obtained by restricting the integration overξ to domains following a dyadic parabolic decomposition.…”

“…The number R of expansion terms is controlled by the prescribed accuracy ε of the tensor product representation. For a detailed description of the box-algorithm and its implementation, we refer to [2]. Based on this tensor product representation, it is possible to group computations and to evaluate the action ofF ij in Step 3 for all data wave packets of the same frequency box Bν ,k at once instead of for each ϕ γ individually.…”

“…An alternative approach for obtaining solutions in the vicinity of caustics has been proposed previously [2,19,20] for the special case of Fourier integral operators corresponding to parametrices of evolution equations for isotropic media. It consist in a re-decomposition strategy following a multi-product representation of the propagator.…”

“…To arrive at such an algorithm, we construct a universal oscillatory integral representation of the kernels of these Fourier integral operators by introducing singularity resolving diffeomorphisms where caustics occur. The universal representation is of the form such that the algorithm based on the dyadic parabolic decomposition of phase space previously developed by the authors applies [2]. We refer to [7,8,10,11] for related computational methods aiming at the evaluation of the action of Fourier integral operators.…”

“…We continue this procedure until the caustic is covered with overlapping sets, associated with diffeomorphisms for the corresponding rank deficiencies. Then we repeat the steps for each caustic and arrive at a collection of open sets covering the canonical relation.In the special case of Fourier integral operators corresponding to parametrices of evolution equations, for isotropic media, an alternative approach for obtaining solutions in the vicinity of caustics, based on a re-decomposition strategy following a multi-product representation of the propagator, has been proposed previously [2,19,20]. Unlike multi-product representations, our construction does not involve a subdivision of the evolution parameter and yields a single-step computation.…”