Fast Computation of Fourier Integral Operators

We introduce a general purpose algorithm for rapidly computing certain types of oscillatory integrals which frequently arise in problems connected to wave propagation and general hyperbolic equations. The problem is to evaluate numerically a so-called Fourier integral operator (FIO) of the form $\int e^{2\pi i \Phi(x,\xi)} a(x,\xi) \hat{f}(\xi) \mathrm{d}\xi$ at points given on a Cartesian grid. Here, $\xi$ is a frequency variable, $\hat f(\xi)$ is the Fourier transform of the input $f$, $a(x,\xi)$ is an amplitude and $\Phi(x,\xi)$ is a phase function, which is typically as large as $|\xi|$; hence the integral is highly oscillatory at high frequencies. Because an FIO is a dense matrix, a naive matrix vector product with an input given on a Cartesian grid of size $N$ by $N$ would require $O(N^4)$ operations. This paper develops a new numerical algorithm which requires $O(N^{2.5} \log N)$ operations, and as low as $O(\sqrt{N})$ in storage space. It operates by localizing the integral over polar wedges with small angular aperture in the frequency plane. On each wedge, the algorithm factorizes the kernel $e^{2 \pi i \Phi(x,\xi)} a(x,\xi)$ into two components: 1) a diffeomorphism which is handled by means of a nonuniform FFT and 2) a residual factor which is handled by numerical separation of the spatial and frequency variables. The key to the complexity and accuracy estimates is that the separation rank of the residual kernel is \emph{provably independent of the problem size}. Several numerical examples demonstrate the efficiency and accuracy of the proposed methodology. We also discuss the potential of our ideas for various applications such as reflection seismology.Comment: 31 pages, 3 figure

Compressive sensing Petrov-Galerkin approximation of high-dimensional parametric operator equations

We analyze the convergence of compressive sensing based sampling techniques for the efficient evaluation of functionals of solutions for a class of high-dimensional, affine-parametric, linear operator equations which depend on possibly infinitely many parameters. The proposed algorithms are based on so-called "non-intrusive" sampling of the high-dimensional parameter space, reminiscent of Monte-Carlo sampling. In contrast to Monte-Carlo, however, a functional of the parametric solution is then computed via compressive sensing methods from samples of functionals of the solution. A key ingredient in our analysis of independent interest consists in a generalization of recent results on the approximate sparsity of generalized polynomial chaos representations (gpc) of the parametric solution families, in terms of the gpc series with respect to tensorized Chebyshev polynomials. In particular, we establish sufficient conditions on the parametric inputs to the parametric operator equation such that the Chebyshev coefficients of the gpc expansion are contained in certain weighted $\ell_p$-spaces for $0<p\leq 1$. Based on this we show that reconstructions of the parametric solutions computed from the sampled problems converge, with high probability, at the $L_2$, resp. $L_\infty$ convergence rates afforded by best $s$-term approximations of the parametric solution up to logarithmic factors.Comment: revised version, 27 page