matrix $A$, $${\operatorname{tr}}(A)=E_U\bigl(U'AU\bigr) \quad\mbox{and}\quad
{\operatorname{diag}}(A)=E_U(U\circ AU),$$ where the expectation $E_U$ is taken over $U$, a random vector with i.i.d. symmetric Bernoulli components. One can unbiasedly estimate $\mathcal{I}(\theta)$ and $\mathcal{J}(\theta)$ by $$\label{Iij}
{\widehat}{\mathcal{I}}_{ij}(\theta)=\frac{1}{2N_2}\sum
_{k=1}^{N_2}U_k'W^iW^jU_k$$ and $$\begin{aligned}
\label{Jij2} {\widehat}{\mathcal{J}}_{ij}(\theta) & = & \frac{1}{N_2}\sum
_{k=1}^{N_2}U_k'W^iW^jU_k
+\frac{1}{N_2}\sum_{k=1}^{N_2}U_k'W^i
\bigl(W^j\bigr)'U_k
\nonumber\\[-8pt]\\[-8pt]
&&{} -2\sum_{\ell=1}^n \Biggl[
\frac{1}{N_2}\sum_{k=1}^{N_2}
\bigl(U_k\circ W^iU_k\bigr)
\Biggr]_{\ell} \Biggl[\frac{1}{N_2}\sum_{k=1}^{N_2}
\bigl(U_k\circ W^jU_k\bigr)
\Biggr]_{\ell}.
\nonumber\end{aligned}$$ Note that here the set of vectors $U_k$ need not be the same as that in (\[ascore\]) and that $N_2$ may not be the same as $N$, the number of $U_j$’s used to compute the estimate of $\theta$. Evaluating ${\widehat}{\mathcal{I}}(\theta)$ and ${\widehat}{\mathcal
{J}}(\theta)$ requires linear solves since $W^iU_k=K^{-1}(K_iU_k)$ and $(W^i)'U_k=K_i(K^{-1}U_k)$. Note that one can also unbiasedly estimate $\mathcal{J}_{ij}(\theta)$ as the sample covariance of $U'_kW^iU_k$ and $U'_kW^jW_k$ for $k=1,\ldots,N$, but (\[Jij2\]) directly exploits properties of symmetric Bernoulli variables (e.g., $U^2_{1j}=1$). Further study would be needed to see when each approach is preferred.

Linear solver {#sec4.1}
-------------

We consider an iterative solver for solving a set of linear equations $Ax=b$ for a symmetric positive definite matrix $A\in{{\mathbb R}}^{n\times n}$, given a right-hand vector $b$. Since the matrix $A$ (in our case the covariance matrix) is symmetric positive definite, the conjugate gradient algorithm is naturally used. 