Source author record

Laurent Demanet

Laurent Demanet appears in the imported research catalog. Authorship, coauthor and topic links are available while profile ownership is still unclaimed.

ResearcherUnclaimed source record

Catalog footprint

What is connected

20works
14topics
4close collaborators

Actions

Connect this record

Log in to claim

Research graph

See the researcher in context

Open full explorer

Inspect adjacent papers, topics, institutions and collaborators without losing the researcher page.

Building this map preview

BZPEER is loading the nearby papers, people, topics and institutions for this page.

Published work

20 published item(s)

preprint2022arXiv

Redatuming physical systems using symmetric autoencoders

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.

preprint2019arXiv

L-Sweeps: A scalable, parallel preconditioner for the high-frequency Helmholtz equation

We present the first fast solver for the high-frequency Helmholtz equation that scales optimally in parallel, for a single right-hand side. The L-sweeps approach achieves this scalability by departing from the usual propagation pattern, in which information flows in a 180 degree cone from interfaces in a layered decomposition. Instead, with L-sweeps, information propagates in 90 degree cones induced by a checkerboard domain decomposition (CDD). We extend the notion of accurate transmission conditions to CDDs and introduce a new sweeping strategy to efficiently track the wave fronts as they propagate through the CDD. The new approach decouples the subdomains at each wave front, so that they can be processed in parallel, resulting in better parallel scalability than previously demonstrated in the literature. The method has an overall O((N/p) log w) empirical run-time for N=n^d total degrees-of-freedom in a d-dimensional problem, frequency w, and p=O(n) processors. We introduce the algorithm and provide a complexity analysis for our parallel implementation of the solver. We corroborate all claims in several two- and three-dimensional numerical examples involving constant, smooth, and discontinuous wave speeds.

preprint2016arXiv

Full waveform inversion with extrapolated low frequency data

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.

preprint2016arXiv

Nested domain decomposition with polarized traces for the 2D Helmholtz equation

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.

preprint2016arXiv

Stable extrapolation of analytic functions

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.

preprint2015arXiv

A short note on the nested-sweep polarized traces method for the 2D Helmholtz equation

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.

preprint2015arXiv

The method of polarized traces for the 2D Helmholtz equation

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.

preprint2015arXiv

The recoverability limit for superresolution via sparsity

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.

preprint2014arXiv

Compressed absorbing boundary conditions via matrix probing

Absorbing layers are sometimes required to be impractically thick in order to offer an accurate approximation of an absorbing boundary condition for the Helmholtz equation in a heterogeneous medium. It is always possible to reduce an absorbing layer to an operator at the boundary by layer-stripping elimination of the exterior unknowns, but the linear algebra involved is costly. We propose to bypass the elimination procedure, and directly fit the surface-to-surface operator in compressed form from a few exterior Helmholtz solves with random Dirichlet data. The result is a concise description of the absorbing boundary condition, with a complexity that grows slowly (often, logarithmically) in the frequency parameter.

preprint2014arXiv

Scaling Law for Recovering the Sparsest Element in a Subspace

We address the problem of recovering a sparse $n$-vector within a given subspace. This problem is a subtask of some approaches to dictionary learning and sparse principal component analysis. Hence, if we can prove scaling laws for recovery of sparse vectors, it will be easier to derive and prove recovery results in these applications. In this paper, we present a scaling law for recovering the sparse vector from a subspace that is spanned by the sparse vector and $k$ random vectors. We prove that the sparse vector will be the output to one of $n$ linear programs with high probability if its support size $s$ satisfies $s \lesssim n/\sqrt{k \log n}$. The scaling law still holds when the desired vector is approximately sparse. To get a single estimate for the sparse vector from the $n$ linear programs, we must select which output is the sparsest. This selection process can be based on any proxy for sparsity, and the specific proxy has the potential to improve or worsen the scaling law. If sparsity is interpreted in an $\ell_1/\ell_\infty$ sense, then the scaling law can not be better than $s \lesssim n/\sqrt{k}$. Computer simulations show that selecting the sparsest output in the $\ell_1/\ell_2$ or thresholded-$\ell_0$ senses can lead to a larger parameter range for successful recovery than that given by the $\ell_1/\ell_\infty$ sense.

preprint2013arXiv

A parallel butterfly algorithm

The butterfly algorithm is a fast algorithm which approximately evaluates a discrete analogue of the integral transform \int K(x,y) g(y) dy at large numbers of target points when the kernel, K(x,y), is approximately low-rank when restricted to subdomains satisfying a certain simple geometric condition. In d dimensions with O(N^d) quasi-uniformly distributed source and target points, when each appropriate submatrix of K is approximately rank-r, the running time of the algorithm is at most O(r^2 N^d log N). A parallelization of the butterfly algorithm is introduced which, assuming a message latency of αand per-process inverse bandwidth of β, executes in at most O(r^2 N^d/p log N + βr N^d/p + α)log p) time using p processes. This parallel algorithm was then instantiated in the form of the open-source DistButterfly library for the special case where K(x,y)=exp(i Φ(x,y)), where Φ(x,y) is a black-box, sufficiently smooth, real-valued phase function. Experiments on Blue Gene/Q demonstrate impressive strong-scaling results for important classes of phase functions. Using quasi-uniform sources, hyperbolic Radon transforms and an analogue of a 3D generalized Radon transform were respectively observed to strong-scale from 1-node/16-cores up to 1024-nodes/16,384-cores with greater than 90% and 82% efficiency, respectively.

preprint2013arXiv

Eventual linear convergence of the Douglas Rachford iteration for basis pursuit

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.

preprint2013arXiv

Stable optimizationless recovery from phaseless linear measurements

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.

preprint2013arXiv

Super-resolution via superset selection and pruning

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.

preprint2013arXiv

Velocity estimation via registration-guided least-squares inversion

This paper introduces an iterative scheme for acoustic model inversion where the notion of proximity of two traces is not the usual least-squares distance, but instead involves registration as in image processing. Observed data are matched to predicted waveforms via piecewise-polynomial warpings, obtained by solving a nonconvex optimization problem in a multiscale fashion from low to high frequencies. This multiscale process requires defining low-frequency augmented signals in order to seed the frequency sweep at zero frequency. Custom adjoint sources are then defined from the warped waveforms. The proposed velocity updates are obtained as the migration of these adjoint sources, and cannot be interpreted as the negative gradient of any given objective function. The new method, referred to as RGLS, is successfully applied to a few scenarios of model velocity estimation in the transmission setting. We show that the new method can converge to the correct model in situations where conventional least-squares inversion suffers from cycle-skipping and converges to a spurious model.

preprint2012arXiv

Sublinear randomized algorithms for skeleton decompositions

Let $A$ be a $n$ by $n$ matrix. A skeleton decomposition is any factorization of the form $CUR$ where $C$ comprises columns of $A$, and $R$ comprises rows of $A$. In this paper, we consider uniformly sampling $ł\simeq k \log n$ rows and columns to produce a skeleton decomposition. The algorithm runs in $O(ł^3)$ time, and has the following error guarantee. Let $\norm{\cdot}$ denote the 2-norm. Suppose $A\simeq X B Y^T$ where $X,Y$ each have $k$ orthonormal columns. Assuming that $X,Y$ are incoherent, we show that with high probability, the approximation error $\norm{A-CUR}$ will scale with $(n/ł)\norm{A-X B Y^T}$ or better. A key step in this algorithm involves regularization. This step is crucial for a nonsymmetric $A$ as empirical results suggest. Finally, we use our proof framework to analyze two existing algorithms in an intuitive way.

preprint2011arXiv

Conditioning bounds for traveltime tomography in layered media

This paper revisits the problem of recovering a smooth, isotropic, layered wave speed profile from surface traveltime information. While it is classic knowledge that the diving (refracted) rays classically determine the wave speed in a weakly well-posed fashion via the Abel transform, we show in this paper that traveltimes of reflected rays do not contain enough information to recover the medium in a well-posed manner, regardless of the discretization. The counterpart of the Abel transform in the case of reflected rays is a Fredholm kernel of the first kind which is shown to have singular values that decay at least root-exponentially. Kinematically equivalent media are characterized in terms of a sequence of matching moments. This severe conditioning issue comes on top of the well-known rearrangement ambiguity due to low velocity zones. Numerical experiments in an ideal scenario show that a waveform-based model inversion code fits data accurately while converging to the wrong wave speed profile.

preprint2011arXiv

Matrix probing and its conditioning

When a matrix A with n columns is known to be well approximated by a linear combination of basis matrices B_1,..., B_p, we can apply A to a random vector and solve a linear system to recover this linear combination. The same technique can be used to recover an approximation to A^-1. A basic question is whether this linear system is invertible and well-conditioned. In this paper, we show that if the Gram matrix of the B_j's is sufficiently well-conditioned and each B_j has a high numerical rank, then n {proportional} p log^2 n will ensure that the linear system is well-conditioned with high probability. Our main application is probing linear operators with smooth pseudodifferential symbols such as the wave equation Hessian in seismic imaging. We demonstrate numerically that matrix probing can also produce good preconditioners for inverting elliptic operators in variable media.

preprint2011arXiv

Matrix probing: a randomized preconditioner for the wave-equation Hessian

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.

preprint2008arXiv

Compressive Wave Computation

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.