We address the solution of large-scale statistical inverse problems using Bayesian inference. In this context, expensive forward models (typically partial differential equations) and high-dimensional parameter spaces (e.g., as result from discretization of continuous parameter fields) provide significant challenge for traditional black-box solution techniques. In this work, we extend the tools and insights that make feasible the solution of large-scale deterministic inverse problems (e.g., adjoints, derivatives, inexact-Newton-Krylov methods, etc.) for use in the statistical setting. Deterministic optimization is used explicitly to identify the maximum a-posteriori estimate, or MAP point, and then a low-rank approximation of an appropriately pre-conditioned Hessian operator can be used to compute covariance information for the posterior distribution. We demonstrate the approach for the problem of inferring heterogeneous wavespeeds of an elastic Earth model using synthetic seismic data.