  Gaussian Markov random fields are used in a large number of disciplines in
machine vision and spatial statistics. The models take advantage of sparsity in
matrices introduced through the Markov assumptions, and all operations in
inference and prediction use sparse linear algebra operations that scale well
with dimensionality. Yet, for very high-dimensional models, exact computation
of predictive variances of linear combinations of variables is generally
computationally prohibitive, and approximate methods (generally interpolation
or conditional simulation) are typically used instead. A set of conditions are
established under which the variances of linear combinations of random
variables can be computed exactly using the Takahashi recursions. The ensuing
computational simplification has wide applicability and may be used to enhance
several software packages where model fitting is seated in a maximum-likelihood
framework. The resulting algorithm is ideal for use in a variety of spatial
statistical applications, including \emph{LatticeKrig} modelling, statistical
downscaling, and fixed rank kriging. It can compute hundreds of thousands exact
predictive variances of linear combinations on a standard desktop with ease,
even when large spatial GMRF models are used.
