Lucas Monzón

Fast convolution with the free space Helmholtz Green's function

We construct an approximation of the free space Greens function for the Helmholtz equation that splits the application of this operator between the spatial and the Fourier domains, as in Ewalds method for evaluating lattice sums. In the spatial domain we convolve with a sum of decaying Gaussians with positive coefficients and, in the Fourier domain, we multiply by a band-limited kernel. As a part of our approach, we develop new quadratures appropriate for the singularity of Greens function in the Fourier domain. The approximation and quadratures yield a fast algorithm for computing volumetric convolutions with Greens function in dimensions two and three. The algorithmic complexity scales as , where is selected accuracy, κ is the number of wavelengths in the problem, d is the dimension, and C is a constant. The algorithm maintains its efficiency when applied to functions with singularities. In contrast to the Fast Multipole Method, as , our approximation makes a transition to that of the free space Greens function for the Poisson equation. We illustrate our approach with examples.

Fast algorithms for Helmholtz Green's functions

The formal representation of the quasi-periodic Helmholtz Greens function obtained by the method of images is only conditionally convergent and, thus, requires an appropriate summation convention for its evaluation. Instead of using this formal sum, we derive a candidate Greens function as a sum of two rapidly convergent series, one to be applied in the spatial domain and the other in the Fourier domain (as in Ewalds method). We prove that this representation of Greens function satisfies the Helmholtz equation with the quasi-periodic condition and, furthermore, leads to a fast algorithm for its application as an operator. We approximate the spatial series by a short sum of separable functions given by Gaussians in each variable. For the series in the Fourier domain, we exploit the exponential decay of its terms to truncate it. We use fast and accurate algorithms for convolving functions with this approximation of the quasi-periodic Greens function. The resulting method yields a fast solver for the Helmholtz equation with the quasi-periodic boundary condition. The algorithm is adaptive in the spatial domain and its performance does not significantly deteriorate when Greens function is applied to discontinuous functions or potentials with singularities. We also construct Helmholtz Greens functions with Dirichlet, Neumann or mixed boundary conditions on simple domains and use a modification of the fast algorithm for the quasi-periodic Greens function to apply them. The complexity, in dimension d≥2, of these algorithms is (κd log κ+C(log ε−1)d), where ε is the desired accuracy, κ is proportional to the number of wavelengths contained in the computational domain and C is a constant. We illustrate our approach with examples.

Grids and transforms for band-limited functions in a disk

We develop fast discrete Fourier transforms (and their adjoints) from a square in space to a disk in the Fourier domain. Since our new transforms are not unitary, we develop a fast inversion algorithm and derive corresponding estimates that allow us to avoid iterative methods typically used for inversion. We consider the eigenfunctions of the corresponding band-limiting and space-limiting operator to describe spaces on which these new transforms can be inverted and made useful. In the process, we construct polar grids which provide quadratures and interpolation with controlled accuracy for functions band-limited within a disk. For rapid computation of the involved trigonometric sums we use the unequally spaced fast Fourier transform, thus yielding fast algorithms for all new transforms. We also introduce polar grids motivated by linearized scattering problems which are obtained by discretizing a family of circles. These circles are generated by using a single circle passing through the origin and rotating this circle with the origin as a pivot. For such grids, we provide a fast algorithm for interpolation to a near optimal grid in the disk, yielding an accurate adjoint transform and inversion algorithm.

Trigonometric identities and sums of separable functions

Modern computers have made commonplace many calculations that were impossible to imagine a few years ago. Still, when you face a problem with a high physical dimension, you immediately encounter the Curse of Dimensionality [1, p.94]. This curse is that the amount of computing power that you need grows exponentially with the dimension. The simplest manifestation of this curse appears when you try to represent a function by its sample values on a grid. If a function of one variable requires N samples, then an analogous function of n variables will need a grid of N samples. Thus, even relatively small problems in high dimensions are still unreasonably expensive. A method has been proposed in [2] to address this problem, based on approximating a function by a sum of separable functions:

On the Design of Highly Accurate and Efficient IIR and FIR Filters

We describe a systematic method for designing highly accurate and efficient infinite impulse response (IIR) and finite impulse response (FIR) filters given their specifications. In our approach, we first meet the specifications by constructing an IIR filter with, possibly, a large number of poles. We then construct, for any given accuracy, an optimal IIR version of such filter (with a minimal number of poles). Finally, also for any given accuracy, we convert the IIR filter to an efficient FIR filter cascade (either serial or parallel). Since in this FIR approximation the non-causal part of the IIR filter only introduces an additional delay (as a function of the desired accuracy), our IIR construction does not have to enforce causality. Thus, we obtain a simple method for constructing linear phase filters if the specifications so require. All of these procedures are accomplished via robust, fast algorithms. We provide several illustrative examples of our method.

Rational approximations for tomographic reconstructions

We use optimal rational approximations of projection data collected in x-ray tomography to improve image resolution. Under the assumption that the object of interest is described by functions with jump discontinuities, for each projection we construct its rational approximation with a small (near optimal) number of terms for a given accuracy threshold. This allows us to augment the measured data, i.e., double the number of available samples in each projection or, equivalently, extend (double) the domain of their Fourier transform. We also develop a new, fast, polar coordinate Fourier domain algorithm which uses our nonlinear approximation of projection data in a natural way. Using augmented projections of the Shepp–Logan phantom, we provide a comparison between the new algorithm and the standard filtered back-projection algorithm. We demonstrate that the reconstructed image has improved resolution without additional artifacts near sharp transitions in the image.

A new inversion method for NMR signal processing

We present a new, semi-analytic inversion method for nuclear magnetic resonance (NMR) log measurements. Our method represents multiwait-time measurements via short sums of exponentials. The resulting sparse T2 distribution requires fewer T2 relaxation times than present in linearized inversion methods. The T1 relaxation times, and corresponding amplitudes are estimated via convex optimization and a semi-analytic algorithm. We obtain an efficient way to represent the NMR data that can be utilized to estimate petrophysical properties and for compression in logging-while-drilling applications.

Fast and accurate propagation of coherent light

We describe a fast algorithm to propagate, for any user-specified accuracy, a time-harmonic electromagnetic field between two parallel planes separated by a linear, isotropic and homogeneous medium. The analytical formulation of this problem (ca 1897) requires the evaluation of the so-called Rayleigh–Sommerfeld integral. If the distance between the planes is small, this integral can be accurately evaluated in the Fourier domain; if the distance is very large, it can be accurately approximated by asymptotic methods. In the large intermediate region of practical interest, where the oscillatory Rayleigh–Sommerfeld kernel must be applied directly, current numerical methods can be highly inaccurate without indicating this fact to the user. In our approach, for any user-specified accuracy ε>0, we approximate the kernel by a short sum of Gaussians with complex-valued exponents, and then efficiently apply the result to the input data using the unequally spaced fast Fourier transform. The resulting algorithm has computational complexity , where we evaluate the solution on an N×N grid of output points given an M×M grid of input samples. Our algorithm maintains its accuracy throughout the computational domain.

On computing distributions of products of non-negative independent random variables

Abstract We introduce a new functional representation of probability density functions (PDFs) of non-negative random variables via a product of a monomial factor and linear combinations of decaying exponentials with complex exponents. This approximate representation of PDFs is obtained for any finite, user-selected accuracy. Using a fast algorithm involving Hankel matrices, we develop a general numerical method for computing the PDF of the sums, products, or quotients of any number of non-negative independent random variables yielding the result in the same type of functional representation. We present several examples to demonstrate the accuracy of the approach.