assume the mean is known to be 0 throughout this work. For the application to ozone data in Section \[sec6\], we avoid modeling the mean by removing the monthly mean for each pixel. The simulations in Section \[sec5\] all first preprocess the data by taking a discrete Laplacian, which filters out any mean function that is linear in the coordinates, so that the results in those sections would be unchanged for such mean functions. The loglikelihood is then, up to an additive constant, given by $${{\mathcal}L}(\theta) = -\tfrac{1}{2}Z'K(\theta)^{-1}Z -
\tfrac{1}{2}\log\det\bigl\{K(\theta)\bigr\}.$$ If $K$ has no exploitable structure, the standard direct way of calculating ${{\mathcal}L}(\theta)$ is to compute the Cholesky decompositon of $K(\theta)$, which then allows $Z'K(\theta)^{-1}Z$ and $\log\det\{K(\theta)\}$ to be computed quickly. However, the Cholesky decomposition generally requires $O(n^2)$ storage and $O(n^3)$ computations, either of which can be prohibitive for sufficiently large $n$.

Therefore, it is worthwhile to develop methods that do not require the calculation of the Cholesky decomposition or other matrix decompositions of $K$. If our goal is just to find the maximum likelihood estimate (MLE) and the corresponding Fisher information matrix, we may be able to avoid the computation of the log determinants by considering the score equations, 