Toward the Future of Global Waveform Tomography: Scalable Gauss-Newton and Alternative Data-Partitioning Schemes

Scott French, Vedran Lekic, and Barbara Romanowicz


With increasingly fine-scale parameterization of global tomographic models comes growth of the linear systems that underlie most iterative inversion schemes (e.g. SEMum: Lekic and Romanowicz, 2010). As assembly and solution of these systems becomes intractable on modern workstations, distributed-memory parallel solvers provide an attractive alternative. We summarize the development of a parallel solver tailored to the regularized Gauss-Newton scheme often used in global tomography.

In a separate effort, we explore alternative data-partitioning schemes for waveform inversion using a synthetic dataset. The SPICE tomographic benchmark was developed to evaluate the effects of inversion methodology on the resolving power of tomographic images (Qin et al., 2008). We discuss our use of the simpler SPICE ``preliminary'' waveform dataset to examine data-partitioning schemes different from those employed in Berkeley global tomographic efforts to date.

A parallel Gauss-Newton solver for waveform tomography

In waveform tomography, we seek to invert a dataset of long-period waveforms $\mathbf{d}$ for a model of earth structure $\mathbf{m}$. The non-linear dependence of the wavefield upon $\mathbf{m}$ requires structural sensitivity to be linearized and $\mathbf{m}$ iteratively estimated. The generalized least-squares formalism is popular choice for such problems (e.g. Tarantola, 2005), and gives rise to the following regularized Gauss-Newton scheme:

\begin{displaymath}\left( \mathbf{G}_k^T \mathbf{C}_d^{-1} \mathbf{G}_k + \epsil...
...mathbf{d}_k + \epsilon^2 \mathbf{C}_m^{-1} \delta \mathbf{m}_p \end{displaymath}

where $\mathbf{m}_k$ and $\mathbf{m}_{k+1}$ are current and updated model estimates, $\mathbf{m}_p$ and $\mathbf{C}_m$ are a priori model and covariance matrix estimates, $\mathbf{C}_d$ the data covariance matrix, and $\mathbf{G}_k$ the current waveform Jacobian matrix. Also, for clarity: $\delta\mathbf{m}_k = \mathbf{m}_{k+1} - \mathbf{m}_k$, $\delta\mathbf{d}_k = \mathbf{d}- g(\mathbf{m}_k)$, and $\delta\mathbf{m}_p = \mathbf{m}_p - \mathbf{m}_k$, where $g(\cdot)$ is the forward modeling operator.

While $\mathbf{G}_k$ possesses a sizable zero component, it is generally non-sparse. Thus, the $m$-by-$m$ matrix $\mathbf{G}_k^T \mathbf{C}_d^{-1} \mathbf{G}_k + \epsilon^2 \mathbf{C}_m^{-1}$, where $m$ is the dimension of $\mathbf{m}$, is likely fully dense. In the case of the SEMum model of Lekic and Romanowicz (2010), where $m \simeq 38500$, wall-clock times of approximately one day were observed during assembly and sequential solution of the above system using MATLAB on a modern workstation. Further, repeated solution is required at each linearized iteration ( $k\rightarrow k+1$) in order to evaluate model sensitivity to regularization. Thus, parallel solution for the regularized Gauss-Newton update described above is a desirable capability.

A parallel solver mimicking the functionality of the sequential code was developed and validated, utilizing the Scalable Linear Algebra PACKage (ScaLAPACK) - named for its attractive weak-scaling behavior on large dense problems such as these (Choi et al., 1996). A robust parallel implementation of LU decomposition with partial pivoting was used. Benchmark tests were performed for realizations of the above system arising in the development of SEMum, and wall-clock solve times for a range of processor grids are shown in Table 2.1.

Table 2.1: Mean wall-clock times for a series of validation solves.
Processor Cores 4 16 64
Mean Wall-clock Time (s) 4388 1612 464.7

Perhaps most importantly, this savings in compute-time allows for faster and more detailed enumeration of model sensitivity to regularization - hopefully leading to greater flexibility in scheduling of compute resources for far heavier workloads, such as wavefield modeling with the Spectral Element Method (SEM: e.g. Komatitsch and Vilotte, 1998).

Synthetic tests of alternative data-partitioning schemes

The simple isotropic model underlying the SPICE ``preliminary'' benchmark contains three large ($r>2500$ km) and five small ($r\simeq1000$ km) columnar $V_s$ perturbations of $\pm 2.5$% and $\pm 6$% respectively, from which $V_p$ and $\rho$ are scaled (Qin et al., 2008). These structures have roughly Gaussian profile and extend from the surface to the CMB. The associated SEM synthetic dataset is valid to 50 s period, which we further bandpass filtered with corners at 80 and 250 s and cutoff at 60 and 400 s. Following the waveform partitioning and weighting scheme of Li and Romanowicz (1996), first and second orbit Rayleigh wave fundamental mode and overtone wavepackets were selected from vertical component waveforms for all 27 synthetic events recorded at the 256 stations in the dataset.

While maintaining a 1D crust, we inverted for mantle structure using the path-average approximation (PAVA: Woodhouse and Dziewonski, 1984) and non-linear asymptotic coupling theory (NACT: Li and Romanowicz, 1995) for fundamental mode and overtone data, respectively. We partitioned the inversion at $\sim300$ km, below which the resolving power of fundamental mode data rapidly decreases, and adjusted the relative weight assigned to the fundamental mode and overtone datasets in each partition. We found that we could more effectively extend model sensitivity into the transition zone over simply merging the datasets with a uniform weighting. The cubic spline model parametrization was kept continuous across the partition so that the model would remain smooth. This represents a methodology broadly consistent with Berkeley tomographic studies to date.

To explore alternative approaches to optimizing depth sensitivity, we proposed a simple multi-frequency scheme. Overlapping filters with cutoff periods of 60/140 s, 120/200 s, and 180/260 s were applied to SPICE vertical component waveforms, which were further weighted according to path redundancy as in Li and Romanowicz (1996), as well as maximum waveform amplitude. As the filtered waveforms are dominated by fundamental mode surface waves, we chose to invert full waveforms using PAVA only. During the inversion, sensitivities of each of the filtered datasets were independently weighted to maximize smoothness of total sensitivity with depth. Unsurprisingly, resolution of structure through the transition zone, best constrained by overtone surface waves, was notably degraded relative to the wavepacket-partitioned inversion. A far more complex multi-frequency scheme (e.g. AMI: Lebedev et al., 2005) would likely be required in order to compete with the wavepacket-based scheme.

As approximate forward modeling schemes are supplanted by SEM in waveform tomography (e.g. SEMum: Lekic and Romanowicz, 2010), mitigation of computational cost takes on new importance. One proposed solution is the source-stacking technique (Capdeville et al., 2005). Linearity of the elastic wave equation with respect to source implies that SEM waveforms stacked over individual sources are equivalent to those resulting from simultaneous triggering of sources in a single simulation. Thus, inversion of source-stacked traces may provide a novel avenue for waveform tomography at vastly reduced cost. As a naïve test, we sought to invert the SPICE preliminary benchmark using PAVA for both the sensitivity kernels and as a stand-in for SEM.

Vertical component SPICE synthetic waveforms were bandpassed with corners at 120 and 250 s (cutoff at 100 and 400 s), aligned on event origin time, weighted by scalar moment magnitude, and stacked over all 27 events. PAVA sensitivity kernels and synthetics were similarly weighted and stacked. Using these approximate methods alone, the structure of the SPICE model became recognizable as early as the second inversion iteration (see Figure 2.30). Of course, while these results are highly preliminary and do not reflect realistic data reliability and noise conditions, they reiterate the promise of the source stacking method. Overcoming the shortcomings present in realistic datasets, as well as developing methods for optimized event-station clustering, are of primary concern moving forward.

Figure 2.30: Depth snapshot at 50 km after stacked-waveform inversion iteration 2, PAVA only.
\begin{figure}\epsfig{file=sfrench10_1_1.eps, width=8cm}\end{figure}


This research was supported by National Science Foundation grant NSF/EAR-0738284.


Capdeville, Y., Y. Gung, and B. Romanowicz, Towards global earth tomography using the spectral element method: a technique based on source stacking, Geophys. J. Int, 162, 541-554, 2005.

Choi, J., J. Demmel, I. Dhillon, J. Dongarra, S. Ostrouchov, A. Petitet, K. Stanley, D. Walker, and R.C. Whaley, ScaLAPACK: a portable linear algebra library for distributed memory computers - design issues and performance, Comp. Phys. Comm., 97, 1-15, 1996.

Komatitsch, D. and J. Vilotte, The spectral element method: an efficient tool to simulate the seismic responds of 2D and 3D geological structures, Bull. Seism. Soc. Am., 88, 368-392, 1998.

Lebedev, S., G. Nolet, T. Meier, and R.D. van der Hilst, Automated multimode inversion of surface and S waveforms, Geophys. J. Int., 162, 951-964, 2005.

Lekic, V. and B. Romanowicz, Inferring upper mantle structure by full waveform tomography with the Spectral Element Method, Submitted to Geophys. J. Int., 2010.

Li, X.D. and B. Romanowicz, Comparison of global waveform inversions with and without considering cross-branch modal coupling, Geophys. J. Int., 121, 695-709, 1995.

Li, X.D. and B. Romanowicz, Global mantle shear-velocity model developed using nonlinear asymptotic coupling theory, Geophys. J. R. Astr. Soc., 101, 22,245-22,272, 1996.

Qin, Y., Y. Capdeville, V. Maupin, J.-P. Montagner, S. Lebedev, and E. Beucler, SPICE Benchmark for global tomographic methods, Geophys. J. Int., 175, 598-616, 2008.

Tarantola, A., Inverse problem theory and methods for model parameter estimation, SIAM, 2005.

Woodhouse, J.H. and A.M. Dziewonski, Mapping the upper mantle: Three dimensional modeling of earth structure by inversion of seismic waveforms, J. Geophys. Res, 89, 5953-5986, 1984.

Berkeley Seismological Laboratory
215 McCone Hall, UC Berkeley, Berkeley, CA 94720-4760
Questions or comments? Send e-mail:
© 2007, The Regents of the University of California