Fast Ip solution of large, sparse, linear systems: Application to seismic travel time tomography

Abstract

Least-squares inversion tends to be sensitive to small changes in the assumptions; in the statistical lexicon, it is not “robust.” Least-squares solutions of noisy, ill-conditioned systems often exhibit unphysical, high-frequency oscillations which have to be damped or smoothed in some fashion. The iteratively reweighted least squares (IRLS) algorithm provides a means of computing approximate l p solutions (1 ⩽ p). We will show that IRLS can be combined with the preconditioned conjugate gradient algorithm in order to solve large, sparse, rectangular systems of linear, algebraic equations very efficiently. Further, for p ≅ 1, the resulting algorithm is extremely stable; it seems to give reasonably smooth solutions without minute adjustment of the parameters and does not become unstable if an overly optimistic error estimate is used. It also has the ability to reject, or at least significantly diminish, the influence of outliers in the data and does not require the use of damping. In addition, we will show how to efficiently include diagonal weighting of the observations and parameters in order to incorporate a priori information into the inversion. This combined inversion algorithm has wide applicability in “real-world” inverse theory, parameter estimation, filtering, etc., where data are noisy and not uniformly reliable. Numerical examples of the l 1 inversion of noisy seismic travel time data will be shown.