combined to approximate likelihoods.

We consider methods based on composite likelihoods to be the main competitor to solving (\[ascore\]). The approximate loglikelihoods described by @vecchia [@steinchiwelty; @caragea] can all be written in the following form: for some sequence of pairs of matrices $(A_j,B_j)$, $j=1,\ldots,q$, all with $n$ columns, at most $n$ rows and full rank, $$\label{composite}
\sum_{j=1}^q \log f_{j,\theta}(A_j
Z\mid B_j Z),$$ where $f_{j,\theta}$ is the conditional Gaussian density of $A_j Z$ given $B_j Z$. As proposed by @vecchia and @steinchiwelty, the rank of $B_j$ will generally be larger than that of $A_j$, in which case the main computation in obtaining (\[composite\]) is finding Cholesky decompositions of the covariance matrices of $B_1Z,\ldots,
B_q Z$. For example, @vecchia just lets $A_j Z$ be the $j$th component of $Z$ and $B_j Z$ some subset of $Z_1,\ldots,Z_{j-1}$. If $m$ is the largest of these subsets, then the storage requirements for this computation are $O(m^2)$ rather than $O(n^2)$. Comparable to increasing the number of $U_j$’s in the randomized algorithm used here, this approach can be updated to obtain a better approximation of the likelihood by increasing the size of the subset of $Z_1,\ldots,Z_{j-1}$ to condition on when computing the conditional density of $Z_j$. However, 