each block. For spatial data, one would typically group spatially contiguous observations within blocks. How to block for space–time data is less clear. The results here focus on the variance of the randomized trace approximation. @avron obtain bounds on the probability that the approximation error is less than some quantity and note that these results sometimes give rankings for various randomized trace approximations different from those obtained by comparing variances.

Computational aspects {#sec4}
=====================

Finding $\theta$ that solves the estimating equations (\[ascore\]) requires a nonlinear equation solver in addition to computing linear solves in $K$. The nonlinear solver starts at an initial guess $\theta^0$ and iteratively updates it to approach a (hopefully unique) zero of (\[ascore\]). In each iteration, at $\theta^i$, the nonlinear solver typically requires an evaluation of $g(\theta^i,N)$ in order to find the next iterate $\theta^{i+1}$. In turn, the evaluation of $g$ requires employing a linear solver to compute the set of vectors $K^{-1}Z$ and $K^{-1}U_j$, $j=1,\ldots,N$.

The Fisher information matrix $\mathcal{I}(\theta)$ and the matrix $\mathcal{J}(\theta)$ contain terms involving matrix traces and diagonals. Write ${\operatorname{diag}}(\cdot)$ for a column vector containing the diagonal elements of a matrix and $\circ$ for the Hadamard (elementwise) product of matrices. For any real 