which are obtained by setting the gradient of the loglikelihood equal to 0. Specifically, defining $K_i =
\frac{\partial}{\partial\theta_i}K(\theta)$, the score equations for $\theta$ are given by (suppressing the dependence of $K$ on $\theta$) $$\label{score}
\tfrac{1}{2}Z'K^{-1}K_iK^{-1}Z
-\tfrac{1}{2} \operatorname{tr}\bigl(K^{-1}K_i\bigr) = 0$$ for $i=1,\ldots,p$. If these equations have a unique solution for $\theta\in\Theta$, this solution will generally be the MLE.

Iterative methods often provide an efficient (in terms of both storage and computation) way of computing solves in $K$ (expressions of the form $K^{-1}x$ for vectors $x$) and are based on being able to multiply arbitrary vectors by $K$ rapidly. In particular, assuming the elements of $K$ can be calculated as needed, iterative methods require only $O(n)$ storage, unlike matrix decompositions such as the Cholesky, which generally require $O(n^2)$ storage. In terms of computations, two factors drive the speed of iterative methods: the speed of matrix–vector multiplications and the number of iterations. Exact matrix–vector multiplication generally requires $O(n^2)$ operations, but if the data form a partial grid, then it can be done in $O(n\log n)$ operations using circulant embedding and the fast Fourier transform. For irregular observations, fast multipole approximations can be used \[@anitescu2012mfa\]. The number of iterations required is related to the 