Maximizing the likelihood has been widely used for estimating the unknown covariance parameters of spatial Gaussian processes. However, evaluating and optimizing the likelihood function can be computationally intractable, particularly for large number of (possibly) irregularly spaced observations, due to the need to handle the inverse of ill-conditioned and large covariance matrices. Extending the "inversion-free" method of Anitescu, Chen and Stein [1], we investigate a broad class of covariance parameter estimation based on inversion-free surrogate losses and block diagonal approximation schemes of the covariance structure. This class of estimators yields a spectrum for negotiating the trade-off between statistical accuracy and computational cost. We present fixed-domain asymptotic properties of our proposed method, establishing √ n-consistency and asymptotic normality results for isotropic Matern Gaussian processes observed on a multi-dimensional and irregular lattice. Simulation studies are also presented for assessing the scalability and statistical efficiency of the proposed algorithm for large data sets.MSC 2010 subject classifications: Primary 62M30, 62M40; secondary 60G15. We are grateful to Mihai Anitescu and Michael Stein for valuable discussions and suggestions. 1 arXiv:1811.12602v3 [math.ST] 7 Jul 2019 2. Preconditioning and inversion-free surrogate loss 2.1. Gaussian processes observed on irregular lattices Consider a zero mean, real valued, and stationary Gaussian process G on domain D, where D is a bounded subset of R d such as [0, 1]d . The dependence structure of G is typically parametrized by a variance parameter φ 0 > 0 and a (correlation) range parameter ρ 0 . Specifically, if G is a geometric anisotropic process on D, then there are a fully known covariance function K and a matrix ρ 0 ∈ R d×d such thatThe objective is to estimate the microergodic parameters of the covariance function, given n measurements from one realization of G at locations D n = {s 1 , . . . , s n } ⊂ D. Throughout the paper, we assume that ρ 0 2 2 .However throughout the paper and for simplifying the theoretical analysis, we only consider the case of non-overlapping bins. It will also be assumed that w i = 1 for all i ∈ {1, . . . , b n }.Remark 3.2. It is informative to take an alternative viewpoint of the LIF objective function in (3.2) as corresponding to a block diagonal approximation of the covariance matrix. Interestingly, as a consequence of the asymptotic theory developed in the next section, this approximation does not affect the asymptotic estimation rate, but it can substantially help to speed up the computation. The block diagonal approximation of K n,m (ρ) corresponding to partitioning scheme B, to be denoted by K B n,m (ρ), can be described as follows. Choose any s, s ∈ D n , and let t, t denote the index of the elements in B containing s and s , i.e., s ∈ B t and s ∈ B t . The entries of K B n,m (ρ) can be equivalently represented by K B n,m (ρ) s,s = [K n,m (ρ)] s,s 1 {t=t } .Remark 4.3. It has been discus...