Laurent Demanet

Emmanuel Candes, Laurent Demanet, Lexing Ying

This paper is concerned with the fast computation of Fourier integral operators of the general form $\int_{\R^d} e^{2\piıΦ(x,k)} f(k) d k$, where $k$ is a frequency variable, $Φ(x,k)$ is a phase function obeying a standard homogeneity condition, and $f$ is a given input. This is of interest for such fundamental computations are connected with the problem of finding numerical solutions to wave equations, and also frequently arise in many applications including reflection seismology, curvilinear tomography and others. In two dimensions, when the input and output are sampled on $N \times N$ Cartesian grids, a direct evaluation requires $O(N^4)$ operations, which is often times prohibitively expensive. This paper introduces a novel algorithm running in $O(N^2 \log N)$ time, i. e. with near-optimal computational complexity, and whose overall structure follows that of the butterfly algorithm [Michielssen and Boag, IEEE Trans Antennas Propagat 44 (1996), 1086-1093]. Underlying this algorithm is a mathematical insight concerning the restriction of the kernel $e^{2\piıΦ(x,k)}$ to subsets of the time and frequency domains. Whenever these subsets obey a simple geometric condition, the restricted kernel has approximately low-rank; we propose constructing such low-rank approximations using a special interpolation scheme, which prefactors the oscillatory component, interpolates the remaining nonoscillatory part and, lastly, remodulates the outcome. A byproduct of this scheme is that the whole algorithm is highly efficient in terms of memory requirement. Numerical results demonstrate the performance and illustrate the empirical properties of this algorithm.

Laurent Demanet, Paul Hand

We address the problem of recovering an n-vector from m linear measurements lacking sign or phase information. We show that lifting and semidefinite relaxation suffice by themselves for stable recovery in the setting of m = O(n log n) random sensing vectors, with high probability. The recovery method is optimizationless in the sense that trace minimization in the PhaseLift procedure is unnecessary. That is, PhaseLift reduces to a feasibility problem. The optimizationless perspective allows for a Douglas-Rachford numerical algorithm that is unavailable for PhaseLift. This method exhibits linear convergence with a favorable convergence rate and without any parameter tuning.

Laurent Demanet, Deanna Needell, Nam Nguyen

We present a pursuit-like algorithm that we call the "superset method" for recovery of sparse vectors from consecutive Fourier measurements in the super-resolution regime. The algorithm has a subspace identification step that hinges on the translation invariance of the Fourier transform, followed by a removal step to estimate the solution's support. The superset method is always successful in the noiseless regime (unlike L1-minimization) and generalizes to higher dimensions (unlike the matrix pencil method). Relative robustness to noise is demonstrated numerically.

Laurent Demanet, Lexing Ying

This paper deals with efficient numerical representation and manipulation of differential and integral operators as symbols in phase-space, i.e., functions of space $x$ and frequency $ξ$. The symbol smoothness conditions obeyed by many operators in connection to smooth linear partial differential equations allow to write fast-converging, non-asymptotic expansions in adequate systems of rational Chebyshev functions or hierarchical splines. The classical results of closedness of such symbol classes under multiplication, inversion and taking the square root translate into practical iterative algorithms for realizing these operations directly in the proposed expansions. Because symbol-based numerical methods handle operators and not functions, their complexity depends on the desired resolution $N$ very weakly, typically only through $\log N$ factors. We present three applications to computational problems related to wave propagation: 1) preconditioning the Helmholtz equation, 2) decomposing wavefields into one-way components and 3) depth-stepping in reflection seismology.

Yunyue Elita Li, Laurent Demanet

The availability of low frequency data is an important factor in the success of full waveform inversion (FWI) in the acoustic regime. The low frequencies help determine the kinematically relevant, low-wavenumber components of the velocity model, which are in turn needed to avoid convergence of FWI to spurious local minima. However, acquiring data below 2 or 3 Hz from the field is a challenging and expensive task. In this paper we explore the possibility of synthesizing the low frequencies computationally from high-frequency data, and use the resulting prediction of the missing data to seed the frequency sweep of FWI. As a signal processing problem, bandwidth extension is a very nonlinear and delicate operation. It requires a high-level interpretation of bandlimited seismic records into individual events, each of which is extrapolable to a lower (or higher) frequency band from the non-dispersive nature of the wave propagation model. We propose to use the phase tracking method for the event separation task. The fidelity of the resulting extrapolation method is typically higher in phase than in amplitude. To demonstrate the reliability of bandwidth extension in the context of FWI, we first use the low frequencies in the extrapolated band as data substitute, in order to create the low-wavenumber background velocity model, and then switch to recorded data in the available band for the rest of the iterations. The resulting method, EFWI for short, demonstrates surprising robustness to the inaccuracies in the extrapolated low frequency data. With two synthetic examples calibrated so that regular FWI needs to be initialized at 1 Hz to avoid local minima, we demonstrate that FWI based on an extrapolated [1, 5] Hz band, itself generated from data available in the [5, 15] Hz band, can produce reasonable estimations of the low wavenumber velocity models.

Leonardo Zepeda-Núñez, Adrien Scheuer, Russell J. Hewett et al.

We present a fast solver for the 3D high-frequency Helmholtz equation in heterogeneous, constant density, acoustic media. The solver is based on the method of polarized traces, coupled with distributed linear algebra libraries and pipelining to obtain an empirical online runtime $ \mathcal{O}(\max(1,R/n) N \log N)$ where $N = n^3$ is the total number of degrees of freedom and $R$ is the number of right-hand sides. Such a favorable scaling is a prerequisite for large-scale implementations of full waveform inversion (FWI) in frequency domain.

Laurent Demanet, Nam Nguyen

We consider the problem of robustly recovering a $k$-sparse coefficient vector from the Fourier series that it generates, restricted to the interval $[- Ω, Ω]$. The difficulty of this problem is linked to the superresolution factor SRF, equal to the ratio of the Rayleigh length (inverse of $Ω$) by the spacing of the grid supporting the sparse vector. In the presence of additive deterministic noise of norm $σ$, we show upper and lower bounds on the minimax error rate that both scale like $(SRF)^{2k-1} σ$, providing a partial answer to a question posed by Donoho in 1992. The scaling arises from comparing the noise level to a restricted isometry constant at sparsity $2k$, or equivalently from comparing $2k$ to the so-called $σ$-spark of the Fourier system. The proof involves new bounds on the singular values of restricted Fourier matrices, obtained in part from old techniques in complex analysis.

Laurent Demanet, Lexing Ying

This paper presents a numerical compression strategy for the boundary integral equation of acoustic scattering in two dimensions. These equations have oscillatory kernels that we represent in a basis of wave atoms, and compress by thresholding the small coefficients to zero. This phenomenon was perhaps first observed in 1993 by Bradie, Coifman, and Grossman, in the context of local Fourier bases \cite{BCG}. Their results have since then been extended in various ways. The purpose of this paper is to bridge a theoretical gap and prove that a well-chosen fixed expansion, the nonstandard wave atom form, provides a compression of the acoustic single and double layer potentials with wave number $k$ as $O(k)$-by-$O(k)$ matrices with $O(k^{1+1/\infty})$ nonnegligible entries, with a constant that depends on the relative $\ell_2$ accuracy $\eps$ in an acceptable way. The argument assumes smooth, separated, and not necessarily convex scatterers in two dimensions. The essential features of wave atoms that enable to write this result as a theorem is a sharp time-frequency localization that wavelet packets do not obey, and a parabolic scaling wavelength $\sim$ (essential diameter)${}^2$. Numerical experiments support the estimate and show that this wave atom representation may be of interest for applications where the same scattering problem needs to be solved for many boundary conditions, for example, the computation of radar cross sections.

Laurent Demanet, Xiangxiong Zhang

We provide a simple analysis of the Douglas-Rachford splitting algorithm in the context of $\ell^1$ minimization with linear constraints, and quantify the asymptotic linear convergence rate in terms of principal angles between relevant vector spaces. In the compressed sensing setting, we show how to bound this rate in terms of the restricted isometry constant. More general iterative schemes obtained by $\ell^2$-regularization and over-relaxation including the dual split Bregman method are also treated, which answers the question how to choose the relaxation and soft-thresholding parameters to accelerate the asymptotic convergence rate. We make no attempt at characterizing the transient regime preceding the onset of linear convergence.

Laurent Demanet, Pierre-David Létourneau, Nicolas Boumal et al.

This paper considers the problem of approximating the inverse of the wave-equation Hessian, also called normal operator, in seismology and other types of wave-based imaging. An expansion scheme for the pseudodifferential symbol of the inverse Hessian is set up. The coefficients in this expansion are found via least-squares fitting from a certain number of applications of the normal operator on adequate randomized trial functions built in curvelet space. It is found that the number of parameters that can be fitted increases with the amount of information present in the trial functions, with high probability. Once an approximate inverse Hessian is available, application to an image of the model can be done in very low complexity. Numerical experiments show that randomized operator fitting offers a compelling preconditioner for the linearized seismic inversion problem.

Laurent Demanet, Alex Townsend

This paper examines the problem of extrapolation of an analytic function for $x > 1$ given perturbed samples from an equally spaced grid on $[-1,1]$. Mathematical folklore states that extrapolation is in general hopelessly ill-conditioned, but we show that a more precise statement carries an interesting nuance. For a function $f$ on $[-1,1]$ that is analytic in a Bernstein ellipse with parameter $ρ> 1$, and for a uniform perturbation level $ε$ on the function samples, we construct an asymptotically best extrapolant $e(x)$ as a least squares polynomial approximant of degree $M^*$ given explicitly. We show that the extrapolant $e(x)$ converges to $f(x)$ pointwise in the interval $I_ρ\in[1,(ρ+ρ^{-1})/2)$ as $ε\to 0$, at a rate given by a $x$-dependent fractional power of $ε$. More precisely, for each $x \in I_ρ$ we have \[ |f(x) - e(x)| = \mathcal{O}\left( ε^{-\log r(x) / \logρ} \right), \qquad\qquad r(x) = \frac{x+\sqrt{x^2-1}}ρ, \] up to log factors, provided that the oversampling conditioning is satisfied. That is, \[ M^* \leq \frac{1}{2} \sqrt{N}, \] which is known to be needed from approximation theory. In short, extrapolation enjoys a weak form of stability, up to a fraction of the characteristic smoothness length. The number of function samples, $N+1$, does not bear on the size of the extrapolation error provided that it obeys the oversampling condition. We also show that one cannot construct an asymptotically more accurate extrapolant from $N+1$ equally spaced samples than $e(x)$, using any other linear or nonlinear procedure. The proofs involve original statements on the stability of polynomial approximation in the Chebyshev basis from equally spaced samples and these are expected to be of independent interest.

Laurent Demanet, Vincent Jugnon

This note formulates a deterministic recovery result for vectors $x$ from quadratic measurements of the form $(Ax)_i \overline{(Ax)_j}$ for some left-invertible $A$. Recovery is exact, or stable in the noisy case, when the couples $(i,j)$ are chosen as edges of a well-connected graph. One possible way of obtaining the solution is as a feasible point of a simple semidefinite program. Furthermore, we show how the proportionality constant in the error estimate depends on the spectral gap of a data-weighted graph Laplacian. Such quadratic measurements have found applications in phase retrieval, angular synchronization, and more recently interferometric waveform inversion.

Laurent Demanet, Gabriel Peyré

This paper considers large-scale simulations of wave propagation phenomena. We argue that it is possible to accurately compute a wavefield by decomposing it onto a largely incomplete set of eigenfunctions of the Helmholtz operator, chosen at random, and that this provides a natural way of parallelizing wave simulations for memory-intensive applications. This paper shows that L1-Helmholtz recovery makes sense for wave computation, and identifies a regime in which it is provably effective: the one-dimensional wave equation with coefficients of small bounded variation. Under suitable assumptions we show that the number of eigenfunctions needed to evolve a sparse wavefield defined on N points, accurately with very high probability, is bounded by C log(N) log(log(N)), where C is related to the desired accuracy and can be made to grow at a much slower rate than N when the solution is sparse. The PDE estimates that underlie this result are new to the authors' knowledge and may be of independent mathematical interest; they include an L1 estimate for the wave equation, an estimate of extension of eigenfunctions, and a bound for eigenvalue gaps in Sturm-Liouville problems. Numerical examples are presented in one spatial dimension and show that as few as 10 percents of all eigenfunctions can suffice for accurate results. Finally, we argue that the compressive viewpoint suggests a competitive parallel algorithm for an adjoint-state inversion method in reflection seismology.

Leonardo Zepeda-Núñez, Laurent Demanet

We present a solver for the 2D high-frequency Helmholtz equation in heterogeneous, constant density, acoustic media, with online parallel complexity that scales empirically as $\mathcal{O}(\frac{N}{P})$, where $N$ is the number of volume unknowns, and $P$ is the number of processors, as long as $P = \mathcal{O}(N^{1/5})$. This sublinear scaling is achieved by domain decomposition, not distributed linear algebra, and improves on the $P =\mathcal{O}(N^{1/8})$ scaling reported earlier in [L. Zepeda-Núñez and L. Demanet, J. Comput. Phys., 308 (2016), pp. 347-388 ]. The solver relies on a two-level nested domain decomposition: a layered partition on the outer level, and a further decomposition of each layer in cells at the inner level. The Helmholtz equation is reduced to a surface integral equation (SIE) posed at the interfaces between layers, efficiently solved via a nested version of the polarized traces preconditioner [L. Zepeda-Núñez and L. Demanet, J. Comput. Phys., 308 (2016), pp. 347-388.]. The favorable complexity is achieved via an efficient application of the integral operators involved in the SIE.

Dmitry Batenkov, Laurent Demanet, Hrushikesh N. Mhaskar

Soft extrapolation refers to the problem of recovering a function from its samples, multiplied by a fast-decaying window and perturbed by an additive noise, over an interval which is potentially larger than the essential support of the window. A core theoretical question is to provide bounds on the possible amount of extrapolation, depending on the sample perturbation level and the function prior. In this paper we consider soft extrapolation of entire functions of finite order and type (containing the class of bandlimited functions as a special case), multiplied by a super-exponentially decaying window (such as a Gaussian). We consider a weighted least-squares polynomial approximation with judiciously chosen number of terms and a number of samples which scales linearly with the degree of approximation. It is shown that this simple procedure provides stable recovery with an extrapolation factor which scales logarithmically with the perturbation level and is inversely proportional to the characteristic lengthscale of the function. The pointwise extrapolation error exhibits a Hölder-type continuity with an exponent derived from weighted potential theory, which changes from 1 near the available samples, to 0 when the extrapolation distance reaches the characteristic smoothness length scale of the function. The algorithm is asymptotically minimax, in the sense that there is essentially no better algorithm yielding meaningfully lower error over the same smoothness class. When viewed in the dual domain, the above problem corresponds to (stable) simultaneous de-convolution and super-resolution for objects of small space/time extent. Our results then show that the amount of achievable super-resolution is inversely proportional to the object size, and therefore can be significant for small objects.

Augustin Cosse, Laurent Demanet

This paper studies the problem of deterministic rank-one matrix completion. It is known that the simplest semidefinite programming relaxation, involving minimization of the nuclear norm, does not in general return the solution for this problem. In this paper, we show that in every instance where the problem has a unique solution, one can provably recover the original matrix through two rounds of semidefinite programming relaxation with minimization of the trace norm. We further show that the solution of the proposed semidefinite program is Lipschitz-stable with respect to perturbations of the observed entries, unlike more basic algorithms such as nonlinear propagation or ridge regression. Our proof is based on recursively building a certificate of optimality corresponding to a dual sum-of-squares (SOS) polynomial. This SOS polynomial is built from the polynomial ideal generated by the completion constraints and the monomials provided by the minimization of the trace. The proposed relaxation fits in the framework of the Lasserre hierarchy, albeit with the key addition of the trace objective function. We further show how to represent and manipulate the moment tensor in favorable complexity by means of a hierarchical low-rank decomposition.

Pawan Bharadwaj, Matthew Li, Laurent Demanet

This paper considers physical systems described by hidden states and indirectly observed through repeated measurements corrupted by unmodeled nuisance parameters. A network-based representation learns to disentangle the coherent information (relative to the state) from the incoherent nuisance information (relative to the sensing). Instead of physical models, the representation uses symmetry and stochastic regularization to inform an autoencoder architecture called SymAE. It enables redatuming, i.e., creating virtual data instances where the nuisances are uniformized across measurements.

Matthew Li, Laurent Demanet, Leonardo Zepeda-Núñez

We propose an end-to-end deep learning framework that comprehensively solves the inverse wave scattering problem across all length scales. Our framework consists of the newly introduced wide-band butterfly network coupled with a simple training procedure that dynamically injects noise during training. While our trained network provides competitive results in classical imaging regimes, most notably it also succeeds in the super-resolution regime where other comparable methods fail. This encompasses both (i) reconstruction of scatterers with sub-wavelength geometric features, and (ii) accurate imaging when two or more scatterers are separated by less than the classical diffraction limit. We demonstrate these properties are retained even in the presence of strong noise and extend to scatterers not previously seen in the training set. In addition, our network is straightforward to train requiring no restarts and has an online runtime that is an order of magnitude faster than optimization-based algorithms. We perform experiments with a variety of wave scattering mediums and we demonstrate that our proposed framework outperforms both classical inversion and competing network architectures that specialize in oscillatory wave scattering data.

Matthew Li, Laurent Demanet, Leonardo Zepeda-Núñez

We introduce an end-to-end deep learning architecture called the wide-band butterfly network (WideBNet) for approximating the inverse scattering map from wide-band scattering data. This architecture incorporates tools from computational harmonic analysis, such as the butterfly factorization, and traditional multi-scale methods, such as the Cooley-Tukey FFT algorithm, to drastically reduce the number of trainable parameters to match the inherent complexity of the problem. As a result WideBNet is efficient: it requires fewer training points than off-the-shelf architectures, and has stable training dynamics, thus it can rely on standard weight initialization strategies. The architecture automatically adapts to the dimensions of the data with only a few hyper-parameters that the user must specify. WideBNet is able to produce images that are competitive with optimization-based approaches, but at a fraction of the cost, and we also demonstrate numerically that it learns to super-resolve scatterers in the full aperture scattering setup.

Leonardo Zepeda-Núñez, Laurent Demanet

We present a solver for the 2D high-frequency Helmholtz equation in heterogeneous acoustic media, with online parallel complexity that scales optimally as $\mathcal{O}(\frac{N}{L})$, where $N$ is the number of volume unknowns, and $L$ is the number of processors, as long as $L$ grows at most like a small fractional power of $N$. The solver decomposes the domain into layers, and uses transmission conditions in boundary integral form to explicitly define "polarized traces", i.e., up- and down-going waves sampled at interfaces. Local direct solvers are used in each layer to precompute traces of local Green's functions in an embarrassingly parallel way (the offline part), and incomplete Green's formulas are used to propagate interface data in a sweeping fashion, as a preconditioner inside a GMRES loop (the online part). Adaptive low-rank partitioning of the integral kernels is used to speed up their application to interface data. The method uses second-order finite differences. The complexity scalings are empirical but motivated by an analysis of ranks of off-diagonal blocks of oscillatory integrals. They continue to hold in the context of standard geophysical community models such as BP and Marmousi 2, where convergence occurs in 5 to 10 GMRES iterations.

Leonardo Zepeda-Núñez, Laurent Demanet

We present a variant of the solver in Zepeda-Núñez and Demanet (2014), for the 2D high-frequency Helmholtz equation in heterogeneous acoustic media. By changing the domain decomposition from a layered to a grid-like partition, this variant yields improved asymptotic online and offline runtimes and a lower memory footprint. The solver has online parallel complexity that scales \emph{sub linearly} as $\mathcal{O} \left( \frac{N}{P} \right)$, where $N$ is the number of volume unknowns, and $P$ is the number of processors, provided that $P = \mathcal{O}(N^{1/5})$. The variant in Zepeda-Núñez and Demanet (2014) only afforded $P = \mathcal{O}(N^{1/8})$. Algorithmic scalability is a prime requirement for wave simulation in regimes of interest for geophysical imaging.

Emmanuel Candes, Laurent Demanet, Lexing Ying

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πi Φ(x,ξ)} a(x,ξ) \hat{f}(ξ) \mathrm{d}ξ$ at points given on a Cartesian grid. Here, $ξ$ is a frequency variable, $\hat f(ξ)$ is the Fourier transform of the input $f$, $a(x,ξ)$ is an amplitude and $Φ(x,ξ)$ is a phase function, which is typically as large as $|ξ|$; 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 πi Φ(x,ξ)} a(x,ξ)$ 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.

Laurent Demanet, Wilhelm Schlag

We make a detailed numerical study of the spectrum of two Schroedinger operators L_- and L_+ arising in the linearization of the supercritical nonlinear Schroedinger equation (NLS) about the standing wave, in three dimensions. This study was motivated by a recent result of the second author on conditional asymptotic stability of solitary waves in the case of a cubic nonlinearity. Underlying the validity of this result is a spectral condition on the operators L_- and L_+, namely that they have no eigenvalues nor resonances in the gap (a region of the positive real axis between zero and the continuous spectrum,) which we call the gap property. The present numerical study verifies this spectral condition, and further shows that the gap property holds for NLS exponents of the form 2*beta + 1, as long as beta_* < beta <= 1, where beta_* = 0.913958905 +- 1e-8. Our strategy consists of rewriting the original eigenvalue problem via the Birman-Schwinger method. From a numerical analysis viewpoint, our main contribution is an efficient quadrature rule for the kernel 1/|x-y| in R^3, i.e., provably spectrally accurate. As a result, we are able to give similar accuracy estimates for all our eigenvalue computations. We also propose an improvement of the Petviashvili's iteration for the computation of standing wave profiles which automatically chooses the radial solution.