We consider the computation of u(t) = exp(−tA)ϕ using rational Krylov subspace reduction for 0 ≤ t < ∞, where u(t), ϕ ∈ R N and 0 < A = A * ∈ R N×N . The objective of this work is the optimization of the shifts for the rational Krylov subspace (RKS). We consider this problem in the frequency domain and reduce it to a classical Zolotaryov problem. The latter yields an asymtotically optimal solution with real shifts. We also construct an infinite sequence of shifts yielding a nested sequence of the RKSs with the same (optimal) Cauchy-Hadamard convergence rate. The effectiveness of the developed approach is demonstrated on an example of the three-dimensional diffusion problem for Maxwell's equation arising in geophysical exploration.

We compute u(t) = exp(−tA)ϕ using rational Krylov subspace reduction for 0 ≤ t < ∞, where u(t), ϕ ∈ R N and 0 < A = A * ∈ R N×N. A priori optimization of the rational Krylov subspaces for this problem was considered in [V. Druskin, L. Knizhnerman, and M. Zaslavsky, SIAM J. Sci. Comput., 31 (2009), pp. 3760-3780]. There was suggested an algorithm generating sequences of equidistributed shifts, which are asymptotically optimal for the cases with uniform spectral distributions. Here we develop a recursive greedy algorithm for choice of shifts taking into account nonuniformity of the spectrum. The algorithm is based on an explicit formula for the residual in the frequency domain allowing adaptive shift optimization at negligible cost. The effectiveness of the developed approach is demonstrated in an example of the three-dimensional diffusion problem for Maxwell's equation arising in geophysical exploration. We compare our approach with the one using the above-mentioned equidistributed sequences of shifts. Numerical examples show that our algorithm is able to adapt to the spectral density of operator A. For examples with near-uniform spectral distributions, both algorithms show the same convergence rates, but the new algorithm produces superior convergence for cases with nonuniform spectra.

We estimate the wave speed in the acoustic wave equation from boundary measurements by constructing a reduced-order model (ROM) matching discrete time-domain data. The state-variable representation of the ROM can be equivalently viewed as a Galerkin projection onto the Krylov subspace spanned by the snapshots of the time-domain solution. The success of our algorithm hinges on the data-driven Gram-Schmidt orthogonalization of the snapshots that suppresses multiple reflections and can be viewed as a discrete form of the Marchenko-Gel'fand-Levitan-Krein algorithm. In particular, the orthogonalized snapshots are localized functions, the (squared) norms of which are essentially weighted averages of the wave speed. The centers of mass of the squared orthogonalized snapshots provide us with the grid on which we reconstruct the velocity. This grid is weakly dependent on the wave speed in traveltime coordinates, so the grid points may be approximated by the centers of mass of the analogous set of squared orthogonalized snapshots generated by a known reference velocity. We present results of inversion experiments for one-and two-dimensional synthetic models.

The motivation of this work is an inverse problem for the acoustic wave equation, where an array of sensors probes an unknown medium with pulses and measures the scattered waves. The goal of the inversion is to determine from these measurements the structure of the scattering medium, modeled by a spatially varying acoustic impedance function. Many inversion algorithms assume that the mapping from the unknown impedance to the scattered waves is approximately linear. The linearization, known as the Born approximation, is not accurate in strongly scattering media, where the waves undergo multiple reflections before they reach the sensors in the array. Thus, the reconstructions of the impedance have numerous artifacts. The main result of the paper is a novel, linear-algebraic algorithm that uses a reduced order model (ROM) to map the data to those corresponding to the single scattering (Born) model. The ROM construction is based only on the measurements at the sensors in the array. The ROM is a proxy for the wave propagator operator, that propagates the wave in the unknown medium over the duration of the time sampling interval. The output of the algorithm can be input into any off-the-shelf inversion software that incorporates state of the art linear inversion algorithms to reconstruct the unknown acoustic impedance.

We introduce a novel nonlinear imaging method for the acoustic wave equation based on datadriven model order reduction. The objective is to image the discontinuities of the acoustic velocity, a coefficient of the scalar wave equation from the discretely sampled time domain data measured at an array of transducers that can act as both sources and receivers. We treat the wave equation along with transducer functionals as a dynamical system. A reduced order model (ROM) for the propagator of such system can be computed so that it interpolates exactly the measured time domain data. The resulting ROM is an orthogonal projection of the propagator on the subspace of the snapshots of solutions of the acoustic wave equation. While the wavefield snapshots are unknown, the projection ROM can be computed entirely from the measured data, thus we refer to such ROM as data-driven. The image is obtained by backprojecting the ROM. Since the basis functions for the projection subspace are not known, we replace them with the ones computed for a known smooth kinematic velocity model. A crucial step of ROM construction is an implicit orthogonalization of solution snapshots. It is a nonlinear procedure that differentiates our approach from the conventional linear imaging methods (Kirchhoff migration and reverse time migration -RTM). It resolves all dynamical behavior captured by the data, so the error from the imperfect knowledge of the velocity model is purely kinematic. This allows for almost complete removal of multiple reflection artifacts, while simultaneously improving the resolution in the range direction compared to conventional RTM.