@article {daskalakis2019SIAMJISasr,
	title = {Accelerating Sparse Recovery by Reducing Chatter},
	journal = {SIAM Journal on Imaging Sciences},
	volume = {13},
	number = {3},
	year = {2020},
	note = {(SIAM Journal on Imaging Sciences)},
	month = {07},
	pages = {1211{\textendash}1239},
	abstract = {Compressive Sensing has driven a resurgence of sparse recovery algorithms with l_1-norm minimization. While these minimizations are relatively well understood for small underdetermined, possibly inconsistent systems, their behavior for large over-determined and inconsistent systems has received much less attention. Specifically, we focus on large systems where computational restrictions call for algorithms that use randomized subsets of rows that are touched a limited number of times. In that regime, l_1-norm minimization algorithms exhibit unwanted fluctuations near the desired solution, and the Linear Bregman iterations are no exception. We explain this observed lack of performance in terms of chatter, a well-known phenomena observed in non-smooth dynamical systems, where intermediate solutions wander between different states stifling convergence. By identifying chatter as the culprit, we modify the Bregman iterations with chatter reducing adaptive element-wise step lengths in combination with potential support detection via threshold crossing. We demonstrate the performance of our algorithm on carefully selected stylized examples and a realistic seismic imaging problem involving millions of unknowns and matrix-free matrix-vector products that involve expensive wave-equation solves.},
	keywords = {chatter, inconsistent linear systems, Kacmarz, linearized Bregman dynamical systems, non-smooth dynamics, sparsity promotion},
	doi = {10.1137/19M129111X},
	url = {https://slim.gatech.edu/Publications/Public/Journals/SIAMJournalOnImagingSciences/2020/daskalakis2019SIAMJISasr/daskalakis2019SIAMJISasr.pdf},
	author = {Emmanouil Daskalakis and Felix J. Herrmann and Rachel Kuske}
}
@article {silva2018alr,
	title = {Applications of low-rank compressed seismic data to full waveform inversion and extended image volumes},
	journal = {Geophysics},
	year = {2019},
	note = {Accepted on January 22, 2019.},
	abstract = {Conventional oil and gas fields are increasingly difficult to explore and image, resulting in the call for more complex wave-equation based inversion algorithms that require dense long-offset samplings. Consequently, there is an exponential growth in the size of data volumes and prohibitive demands on computational resources. In this work, we propose a method to compress and process seismic data directly in a low-rank tensor format, which drastically reduces the amount of storage required to represent the data. We first outline how seismic data exhibits low-rank structure in a particular transform-domain, which can be exploited to compress the dense data in one extremely storage-efficient tensor format when the data is fully sampled. In the more realistic case of missing data, we can use efficient interpolation techniques to approximate the fully sampled volume in compressed form. In either case, once we have our data represented in its compressed tensor form, we design an algorithm to extract source or receiver gathers directly from the compressed parameters. This extraction process can be done on-the-fly directly on the compressed data, in a full waveform inversion context, and does not require scanning through the entire dataset in order to form shot gathers. To the best of our knowledge, this work is one of the first major contributions to working with seismic data applications directly in the compressed domain without reconstructing the entire data volume. We use a stochastic inversion approach, which works with small subsets of source experiments at each iteration, further reducing the computational and memory costs of full waveform inversion. We also demonstrate how this data compression and extraction technique can be applied to forming full subsurface image gathers through probing techniques.},
	keywords = {common shot/receiver gather, extended image volume, full waveform inversion, low rank, tensor},
	doi = {10.1190/geo2018-0116.1},
	url = {https://slim.gatech.edu/Publications/Public/Journals/Geophysics/2019/silva2018alr/silva2018alr.pdf},
	author = {Curt Da Silva and Yiming Zhang and Rajiv Kumar and Felix J. Herrmann}
}
@article {dasilva2017uls,
	title = {A unified 2D/3D large scale software environment for nonlinear inverse problems},
	journal = {ACM Transactions on Mathematical Software},
	year = {2019},
	note = {Accepted on January 27, 2019.},
	abstract = {Large scale parameter estimation problems are some of the most computationally demanding problems. An academic researcher{\textquoteright}s domain-specific knowledge often precludes that of software design, which results in software frameworks for inversion that are technically correct, but not scalable to realistically-sized problems. On the other hand, the computational demands of the problem for realistic problems result in industrial codebases that are geared solely for performance, rather than comprehensibility or flexibility. We propose a new software design that bridges the gap between these two seemingly disparate worlds. A hierarchical and modular design allows a user to delve into as much detail as she desires, while using high performance primitives at the lower levels. Our code has the added benefit of actually reflecting the underlying mathematics of the problem, which lowers the cognitive load on user using it and reduces the initial startup period before a researcher can be fully productive. We also introduce a new preconditioner for the Helmholtz equation that is suitable for fault-tolerant distributed systems. Numerical experiments on a variety of 2D and 3D test problems demonstrate the effectiveness of this approach on scaling algorithms from small to large scale problems with minimal code changes.},
	keywords = {large scale, Matlab, Optimization, PDE-constrained inversion},
	url = {https://slim.gatech.edu/Publications/Public/Journals/ACMTOMS/2019/dasilva2017uls/dasilva2017uls.html},
	author = {Curt Da Silva and Felix J. Herrmann}
}
@article {fang2017uqfip,
	title = {Uncertainty quantification for inverse problems with weak partial-differential-equation constraints},
	journal = {Geophysics},
	volume = {83},
	number = {6},
	year = {2018},
	note = {(Geophysics)},
	pages = {R629-R647},
	abstract = {In a statistical inverse problem, the objective is a complete statistical description of unknown parameters from noisy observations in order to quantify uncertainties of the parameters of interest. We consider inverse problems with partial-differential-equation-constraints, which are applicable to a variety of seismic problems. Bayesian inference is one of the most widely-used approaches to precisely quantify statistics through a posterior distribution, incorporating uncertainties in observed data, modeling kernel, and prior knowledge of the parameters. Typically when formulating the posterior distribution, the partial-differential-equation-constraints are required to be exactly satisfied, resulting in a highly nonlinear forward map and a posterior distribution with many local maxima. These drawbacks make it difficult to find an appropriate approximation for the posterior distribution. Another complicating factor is that traditional Markov chain Monte Carlo methods are known to converge slowly for realistically sized problems. In this work, we relax the partial-differential-equation-constraints by introducing an auxiliary variable, which allows for Gaussian deviations in the partial-differential-equations. Thus, we obtain a new bilinear posterior distribution consisting of both data and partial-differential-equation misfit terms. We illustrate that for a particular range of variance choices for the partial-differential-equation misfit term, the new posterior distribution has fewer modes and can be well-approximated by a Gaussian distribution, which can then be sampled in a straightforward manner. Since it is prohibitively expensive to explicitly construct the dense covariance matrix of the Gaussian approximation for intermediate to large-scale problems, we present a method to implicitly construct it, which enables efficient sampling. We apply this framework to two-dimensional seismic inverse problems with 1,800 and 92,455 unknown parameters. The results illustrate that our framework can produce comparable statistical quantities to those produced by conventional Markov chain Monte Carlo type methods while requiring far fewer partial-differential-equation solves, which are the main computational bottlenecks in these problems.},
	keywords = {acoustic, FWI, UQ, weak-constraint},
	doi = {10.1190/geo2017-0824.1},
	url = {https://slim.gatech.edu/Publications/Public/Journals/Geophysics/2018/fang2017uqfip/fang2017uqfip.html},
	author = {Zhilong Fang and Curt Da Silva and Rachel Kuske and Felix J. Herrmann}
}
@article {kumar2016bls,
	title = {Beating level-set methods for 5D seismic data interpolation: a primal-dual alternating approach},
	journal = {IEEE Transactions on Computational Imaging},
	year = {2017},
	note = {(published online in IEEE Transactions on Computational Imaging)},
	month = {04},
	abstract = {Acquisition cost is a crucial bottleneck for seismic workflows, and low-rank formulations for data interpolation allow practitioners to {\textquoteright}fill in{\textquoteright} data volumes from critically subsampled data acquired in the field. Tremendous size of seismic data volumes required for seismic processing remains a major challenge for these techniques. Residual-constrained formulations require less parameter tuning when the target noise floor is known. We propose a new approach to solve residual constrained formulations for interpolation. We represent the data volume in a compressed manner using low-rank matrix factors, and build a block-coordinate algorithm with constrained convex subproblems that are solved with a primal-dual splitting scheme. The develop optimization framework works on the whole seismic temporal frequency slices and does not require windowing or non-trivial sorting of seismic data. The new approach is competitive with state of the art level-set algorithms that interchange the role of objectives with constraints. We use the new algorithm to successfully interpolate a large scale 5D seismic data volume (upto 1010 data-points), generated from the geologically complex synthetic 3D Compass velocity model, where 80\% of the data has been removed. We also develop a robust extension of the primal-dual approach to deal with the outliers (or noise) in the data.},
	keywords = {alternating minimization, matrix completion, primal-dual splitting, seismic data, seismic trace interpolation},
	doi = {10.1109/TCI.2017.2693966},
	url = {https://slim.gatech.edu/Publications/Public/Journals/IEEETransComputationalImaging/2017/kumar2016bls/kumar2016bls.pdf},
	author = {Rajiv Kumar and Oscar Lopez and Damek Davis and Aleksandr Y. Aravkin and Felix J. Herrmann}
}
@article {kumar2014GEOPemc,
	title = {Efficient matrix completion for seismic data reconstruction},
	journal = {Geophysics},
	volume = {80},
	number = {5},
	year = {2015},
	note = {(Geophysics)},
	month = {09},
	pages = {V97-V114},
	abstract = {Despite recent developments in improved acquisition, seismic data often remain undersampled along source and receiver coordinates, resulting in incomplete data for key applications such as migration and multiple prediction. We have interpreted the missing-trace interpolation problem in the context of matrix completion (MC), and developed three practical principles for using low-rank optimization techniques to recover seismic data. Specifically, we strive for recovery scenarios wherein the original signal is low rank and the subsampling scheme increases the singular values of the matrix. We use an optimization program that restores this low-rank structure to recover the full volume. Omitting one or more of these principles can lead to poor interpolation results, as we found experimentally. In light of this theory, we compensate for the high-rank behavior of data in the source-receiver domain by using the midpoint-offset transformation for 2D data and a source-receiver permutation for 3D data to reduce the overall singular values. Simultaneously, to work with computationally feasible algorithms for large-scale data, we use a factorization-based approach to MC, which significantly speeds up the computations compared with repeated singular value decompositions without reducing the recovery quality. In the context of our theory and experiments, we also find that windowing the data too aggressively can have adverse effects on the recovery quality. To overcome this problem, we carried out our interpolations for each frequency independently while working with the entire frequency slice. The result is a computationally efficient, theoretically motivated framework for interpolating missing-trace data. Our tests on realistic 2D and 3D seismic data sets show that our method compares favorably in terms of computational speed and recovery quality with existing curvelet- and tensor-based techniques.},
	keywords = {2D, 3D, algorithm, interpolation, low-rank, signal processing},
	doi = {10.1190/geo2014-0369.1},
	url = {https://slim.gatech.edu/Publications/Public/Journals/Geophysics/2015/kumar2014GEOPemc/kumar2014GEOPemc.pdf},
	url2 = {http://library.seg.org/doi/abs/10.1190/geo2014-0369.1},
	author = {Rajiv Kumar and Curt Da Silva and Okan Akalin and Aleksandr Y. Aravkin and Hassan Mansour and Ben Recht and Felix J. Herrmann}
}
@article {dasilva2014htuck,
	title = {Optimization on the Hierarchical Tucker manifold - applications to tensor completion},
	journal = {Linear Algebra and its Applications},
	volume = {481},
	year = {2015},
	note = {(Linear Algebra and its Applications)},
	month = {09},
	pages = {131-173},
	abstract = {In this work, we develop an optimization framework for problems whose solutions are well-approximated by Hierarchical Tucker (HT) tensors, an efficient structured tensor format based on recursive subspace factorizations. By exploiting the smooth manifold structure of these tensors, we construct standard optimization algorithms such as Steepest Descent and Conjugate Gradient for completing tensors from missing entries. Our algorithmic framework is fast and scalable to large problem sizes as we do not require SVDs on the ambient tensor space, as required by other methods. Moreover, we exploit the structure of the Gramian matrices associated with the HT format to regularize our problem, reducing overfitting for high subsampling ratios. We also find that the organization of the tensor can have a major impact on completion from realistic seismic acquisition geometries. These samplings are far from idealized randomized samplings that are usually considered in the literature but are realizable in practical scenarios. Using these algorithms, we successfully interpolate large-scale seismic data sets and demonstrate the competitive computational scaling of our algorithms as the problem sizes grow.},
	keywords = {differential geometry, Gauss{\textendash}Newton, hierarchical tucker tensors, low-rank tensor, Riemannian manifold optimization, tensor completion},
	doi = {10.1016/j.laa.2015.04.015},
	url = {https://slim.gatech.edu/Publications/Public/Journals/LinearAlgebraAndItsApplications/2015/dasilva2014htuck/dasilva2014htuck.pdf},
	url2 = {http://www.sciencedirect.com/science/article/pii/S0024379515002530},
	author = {Curt Da Silva and Felix J. Herrmann}
}
