condition number of $K$ (the ratio of the largest to smallest singular value), so that preconditioning \[@chenbook\] is often essential; see @steinchenanitescufiltering for some circumstances under which one can prove that preconditioning works well.

Computing the first term in (\[score\]) requires only one solve in $K$, but the trace term requires $n$ solves (one for each column of $K_i$) for $i=1,\ldots,p$, which may be prohibitive in some circumstances. Recently, @anitescu2012mfa analyzed and demonstrated a stochastic approximation of the trace term based on the Hutchinson trace estimator \[@hutchinson\]. To define it, let $U_1,\ldots,U_N$ be i.i.d. random vectors in ${{\mathbb R}}^n$ with i.i.d. symmetric Bernoulli components, that is, taking on values 1 and $-1$ each with probability $\frac{1}{2}$. Define a set of estimating equations for $\theta$ by $$\label{ascore}
g_i(\theta,N) = \frac{1}{2}Z'K^{-1}K_iK^{-1}Z
-\frac{1}{2N}\sum_{j=1}^N
U_j'K^{-1}K_iU_j = 0$$ for $i=1,\ldots,p$. Throughout this work, $E_\theta$ means to take expectations over $Z\sim N(0,K(\theta))$ and over the $U_j$’s as well. Since $E_\theta(U_1'K^{-1}K_iU_1) = \operatorname{tr}(K^{-1}K_i)$, $E_\theta
g_i(\theta,N) = 0$ and (\[ascore\]) provides a set of unbiased estimating equations for $\theta$. Therefore, we may hope that a solution to (\[ascore\]) will provide a good approximation to the MLE. The unbiasedness of the estimating equations (\[ascore\]) requires only that the components 