Efficient methods for Gaussian Markov random fields under sparse linear constraints
Year: 2022

Abstract:
Methods for inference and simulation of linearly constrained Gaussian Markov Random Fields (GMRF) are computationally prohibitive when the number of constraints is large. In some cases, such as for intrinsic GMRFs, they may even be unfeasible. We propose a new class of methods to overcome these challenges in the common case of sparse constraints, where one has a large number of constraints and each only involves a few elements. Our methods rely on a basis transformation into blocks of constrained versus non-constrained subspaces, and we show that the methods greatly outperform existing alternatives in terms of computational cost. By combining the proposed methods with the stochastic partial differential equation approach for Gaussian random fields, we also show how to formulate Gaussian process regression with linear constraints in a GMRF setting to reduce computational cost. This is illustrated in two applications with simulated data.

heading: 1 Introduction
text: Linearly constrained Gaussian processes (GPs) have recently gained attention, especially for Gaussian process regression where the model should obey underlying physical principle such as conservation laws or equilibrium conditions [13, 16, 17, 30, 31, 36, 37]. A well-known challenge with these models, and GPs in general, is their high computational cost for inference and prediction in the case of big data sets [2, 9, 14]. One way to reduce computational burden is to impose conditional independence assumptions. In fact, conditional independence between random variables is often explicitly or implicitly assumed in large classes of statistical models including Markov processes, hierarchical models, and graphical models. The assumption typically increases the model’s interpretability and facilitates computationally efficient methods for inference [29]. Gaussian variables with conditional independence properties are known as Gaussian Markov random fields (GMRFs), and these are widely used in areas ranging from brain imaging [24, 34, 20] to spatial statistics [7] and time series analysis [28]. GMRFs further arise as computationally efficient approximations of certain GPs [18], which is a fundamental modeling tool in both machine learning and statistics. In particular, such approximations in combination with the integrated nested Laplace approximation (INLA) methodology [29] made latent GMRFs widely used in the applied sciences [3]. GMRFs also have connections with convolutional neural networks leading to recent considerations of deep GMRFs [33].
The key feature of GMRFs that reduces computational cost is sparsity. Specifically, a GMRF
X = [X1, . . . , Xn] > ∼ N ( µ,Q−1 ) , (1)
has a sparse precision (inverse covariance) matrix Q which enables the use of sparse matrix techniques such as sparse Cholesky factorization [8] and iterative methods [35, 23, 34] for computationally efficient sampling and statistical inference [27, 28]. The sparsity is caused by conditional independence assumptions since Qij = 0 if and only if Xi and Xj are independent conditionally on all other variables in X [see 28, Chapter 2]. 35th Conference on Neural Information Processing Systems (NeurIPS 2021).
For a GMRF X, a set of k linear constraints can be formulated as
AX = b, (2)
where each row in the k × n matrix A and the vector b encodes a constraint on X. Note that observations of GMRFs can be formulated as linear constraints, where, e.g., a point observation of Xi is the simple constraint Xi = xi. These deterministic restrictions are often referred to as hard constraints [28]. If (2) is assumed to hold up to Gaussian noise, i.e., AX ∼ N (b, σ2I), then the constraints are instead referred to as soft constraints. This scenario is common when GMRFs are incorporated in hierarchical models, where the soft constraints represent noisy observations [28].
Because of the computational advantages of GMRFs, it is of interest to use GMRFs as alternatives to GPs also in the constrained case. Unfortunately, hard constraints can drastically reduce the computational benefits of GMRFs. Specifically, as we will explain in the next section, the current methods for GMRFs with hard constraints in general have a computational cost that scales cubically in the number of constraints k. This is prohibitive when there are many constraints, which is common for applications of constrained GPs. For example, physically constrained GPs such as those in [13] have k = n and any GMRF model with N observations without measurement noise will have k = N . In the first example the computational cost is cubic in n and in the second it can be cubic in N if the observations are such that A in (2) is not diagonal, which, e.g., is the case for area observations [21]. In both cases, this renders the computational benefits of the unconstrained GMRF irrelevant. Until now, there has been no methods to circumvent this problem.
Summary of contributions. The main contribution of this work is the formulation of a novel class of computationally efficient methods for linearly constrained GMRFs in tasks such as parameter estimation and simulation. These methods explore the conditional independence structure to reduce computational costs compared to traditional methods for situations with large numbers of constraints. The focus is in particular on the case, referred to as sparse hard constraints, when each constraint only involves a few number of variables so that A is sparse. This is, e.g., common for GPs under physical constraints [13] and for GMRF models with point and area observations [21]. For this case, the main idea is to perform a change of basis so that the constraints are simpler to enforce in the transformed basis. In particular, we show how the change of basis facilitates using the same efficient methods for sparse matrices that are used in the non-constrained case also for the constrained models. In order to use these models for GP regression, the methods are also generalized to GMRF models with both hard and soft constraints. An important feature of the new class of methods, that previous approaches lack, is its applicability to intrinsic GMRFs, which are improper in the sense that a set of eigenvalues of Q is zero. This makes their distributions invariant to shifts in certain directions, which is a useful property for prior distributions of Bayesian models [28, Chapter 3]. Because of this they are often used in areas such as semiparametric regression [10], medical image analysis [24, 34], and geostatistics [15, 7], and it is thus of interest to be able to enforce linear constraints also for such models. The final contribution is the derivation of GMRFs for constrained GPs, by combining the proposed methods with the stochastic partial differential equation (SPDE) approach by [18] and the nested SPDE methods by [5]. The combined approach is highly computationally efficient compared to standard covariance-based methods, as illustrated in two simulation studies.
Outline. In Section 2, the problem is introduced in more detail and the most commonly used methods for sampling and likelihood computations for GMRFs are reviewed. Section 3 introduces the new methods for GMRFs with sparse hard constraints. These methods are extended to the case with both hard and soft constraints in Section 4, followed by the GMRF methods for constrained GPs in Section 5. The methods are illustrated numerically in Section 6. A discussion closes the article, which is supported by a supplementary materials containing proofs and technical details.

heading: 2 Standard methods for GMRFs under hard constraints
text: Hard constraints can be divided into interacting and non-interacting constraints. In the latter, (2) specifies a constraint on a subset of the variables in X so that AX = b can be written as Xc = bc, where c denotes a subset of the indices. Specifically, let Xu denote the remaining variables, then
X = [ Xc Xu ] ∼ N ([ µc µu ] , [ Qcc Qcu Quc Quu ]−1) and X|AX = b ∼ N ([ bc µu|c ] , [ 0 0 0 Q−1uu ]) ,
where µu|c = µu −Q −1 uuQuc(bc − µc). Thus, we can split the variables into two subsets and treat the unconditioned variables separately. The more difficult and interesting situation is the case of interacting hard constraints where a simple split of the variables is not possible. From now on, we will assume that we are in this scenario.
Our aim is to construct methods for sampling from the distribution of X|AX = b and for evaluating its log-likelihood function. It is straightforward to show that X|AX = b ∼ N (µ̂, Σ̂), where µ̂ = µ−Q−1A>(AQ−1A>)−1(Aµ−b) and Σ̂ = Q−1−Q−1A>(AQ−1A>)−1AQ−1. Since Σ̂ has rank n− k, likelihood evaluation and sampling is in general expensive. For example, one way is to use an eigenvalue decomposition of Σ̂ [see 28, Chapter 2.3.3]. However, this procedure is not practical since it cannot take advantage of the sparsity of Q. Also, for intrinsic GMRFs, Σ̂ and µ̂ cannot be constructed through the expressions above since Q−1 has unbounded eigenvalues.
A commonly used method for sampling under hard linear constraints, sometimes referred to as conditioning by kriging [27], is to first sample X ∼ N ( µ,Q−1 ) and then correct for the constraints by using X∗ = X−Q−1A>(AQ−1A>)−1(AX−b) as a sample from the conditional distribution. Here the cost of sampling X is CQ + SQ, where CQ denotes the computational cost of a sparse Cholesky factorization Q = R>R and SQ the cost of solving Rx = u for x given u. Adding the cost for the correction step, the total cost of the method is O(CQ + (k + 2)SQ + k3).
Let πAx(·) denote the density of Ax ∼ N (Aµ,AQ−1A>), then the likelihood of X|AX = b can be computed through the expression
π(x|Ax = b) = π∗Ax|x(b|x)π(x)
πAx(b) , (3)
where π∗Ax|x(b|x) = I(Ax = b)|AA >|−1/2 and I(Ax = b) denotes the indicator function with I(Ax = b) = 1 if Ax = b and I(Ax = b) = 0 otherwise. This result is formulated in [27] and we provide further details in the supplementary materials by showing that (3) is a density with respect to the Lebesgue measure on the level set {x : Ax = b}. The computational cost of evaluating the likelihood using this formulation is O(CQ + (k + 1)SQ + k3). Note that these methods only work efficiently for a small number of constraints, because of the term k3 in the computational costs, and for proper GMRFs. In the case of intrinsic GMRFs we cannot work with AQ−1A due to the rank deficiency of Q.

heading: 3 The basis transformation method
text: In this section, we propose the new class of methods in two steps. We first derive a change of basis in Section 3.1, and then use this to formulate the desired conditional distributions in Section 3.2. The resulting computational costs of likelihood evaluations and simulation are discussed in Section 3.3.
Before stating the results we introduce some basic notation. When working with intrinsic GMRFs the definition in (1) is inconvenient since the covariance matrix has infinite eigenvalues. Instead one can use the canonical parametrization X ∼ NC (µC ,Q) , which implies that the density of X is given by π(x) ∝ exp ( − 12x >Qx + µ>Cx )
and thus that µC = Qµ. Also, since we will be working with non-invertible matrices we will need the Moore–Penrose inverse and the pseudo determinant. We denote the Moore–Penrose inverse of a matrix B by B† and for a symmetric positive semi definite matrix M we define the pseudo determinant as |M|† = ∏ i:λi>0
λi where {λi} are the eigenvalues of M. For the remainder of this article, we will assume the following.
Assumption 1. X ∼ NC (Qµ,Q) where Q is a positive semi-definite n × n matrix with rank n− s > 0 and null-space EQ. A is a k × n matrix with rank k and we let k0 denote rank(AEQ).
Finally, we use the index notation C = (1, . . . , k) and U = (k + 1, . . . , n) so that, e.g., vC denotes the sub-vector (v1, . . . , vk)> of the first k elements in the vector v, M·,U (MU ,·) denotes the submatrix obtained by extracting the columns (rows) in M with indices in U , and MU U the submatrix obtained by extracting the columns and rows in M with indices in U .

heading: 3.1 Basis construction
text: Our main idea is to construct a basis on Rn such that the constraints are easily enforced. A key property of the basis is that A should be spanned by the first k elements of the basis. Essentially, this means that we are transforming the natural basis into one where the results for the case of non-interacting hard constraints can be used. To construct the basis, note that if USV> is the singular value decomposition (SVD) of a the matrix A, the basis V> is orthonormal and the first k rows span the image of A and the last n− k rows span the null-space of A. Now, if we let x∗ denote a vector x expressed in the basis V>, then x∗ can be transformed back to the natural basis by x = Vx∗ hence
Ax = b ⇔ USx∗ = b ⇔ [USCC 0]x∗ = b ⇔ USCC x∗C = b ⇔ x∗C = b ∗,
where b∗ = (USCC )−1b. This shows that the SVD is a natural choice of method for constructing the basis of interest, which is defined by an n× n change-of-basis matrix T = V>. In Algorithm 1 we present a simple method to build T based on the SVD. In the algorithm, id(A) is a function that returns the indices of the non-zero columns in A.
Algorithm 1 Constraint basis construction. Require: A (a k × n matrix of rank k)
1: T← In 2: D ← id(A) 3: USV> ← svd(A·,D) 4: TD,D ← V> 5: T← [T·,D T·,Dc ] 6: Return T
The computational cost of the method is dominated by that of the SVD, which is O ( k3 + k2|id(A)| ) [11, p. 493]. Clearly, the cubic scaling in the number of constraints may reduce the efficiency of any method that requires this basis construction as a first step. However, suppose that the rows of A can be split into two sub-matrices, Ã1 and Ã2, which have no common non-zero columns. Then the SVD of the two matrices can be computed separately. Suppose now that A corresponds to m such sub-matrices, then, after reordering,
A = [ Ã > 1 , . . . , Ã > m ]> where {Ãi}mi=1 represent sub-constraints such that id(Ãi) ∩ id(Ãl) = ∅ for all i and l. By replacing the SVD of Algorithm 1 by the m SVDs of the lower-dimensional matrices the computational cost is reduced to O (∑m i=1 rank(Ãi) 3 + rank(Ãi) 2|id(Ãi)| ) . This method is presented in Algorithm 2. The reordering step is easy to perform and is described in the supplementary materials.

heading: 3.2 Conditional distributions
text: Using the change of basis from the previous subsection, we can now derive alternative formulations of the distributions of AX and X|AX = b which are suitable for sampling and likelihood-evaluation. There are two main results in this section. The first provides an expression of the density of AX that allows for computationally efficient likelihood evaluations for observations AX = b. The second formulates the conditional distribution for X|AX = b in a way that allows for efficient sampling of X given observations AX = b. To formulate the results, let T = CB(A) be the output of Algorithm 1 or Algorithm 2 and X∗ = TX which, under Assumption 1, has distribution X∗ ∼ NC (Q∗µ∗,Q∗) , where µ∗ = Tµ and Q∗ = TQT>. Henceforth we use stars to denote
Algorithm 2 CB(A). Constraint basis construction for non-overlapping subsets of constraints.
Require: A (a k × n matrix of rank k) 1: T← In 2: Dfull ← id(A) 3: h← 1 4: l← k + 1 5: Reorder so that A = [ Ã > 1 . . . Ã > m
]> 6: for i = 1 : m do 7: D ← id(Ãi) 8: USV> ← svd(Ãi) 9: u← ncol(U)
10: Th:(h+u),D ← ( V> ) 1:u,D 11: h← h+ u+ 1 12: if |D| > u then 13: Tl:(l+|D|−u),D ← ( V> ) (u+1):|D|,D 14: l← l + |D| − u+ 1 15: end if 16: end for 17: T← [ T·,1:|Dfull| (In)·,Dcfull
] 18: Return T
quantities such as means and precisions in the transformed space. Note that we can move from the transformed space to the original space by multiplying with T> for a vector and by multiplying with T> from the left and T from the right for a matrix.
We use the matrix H = ( AT> ) CC
, which is equal to USCC from the SVD in Algorithm 1, and therefore has inverse H−1 = S−1CC U
>. Finally, recall that b∗ = H−1b. Theorem 1. Under Assumption 1 it follows that
πAX(b) = |Q∗C |U | † 2
(2π)k/2|AA>|1/2 · exp
( −1
2 (b∗ − µ∗C ) > Q∗C |U (b
∗ − µ∗C ) ) ,
where Q∗C |U = Q ∗ CC −Q ∗ CU (Q ∗ U U ) † Q∗U C and |Q ∗ C |U | † 2 = |Q| 12 |Q∗U U |− 1 2 .
If Q is positive definite we have Q∗C |U = Q ∗ CC −Q ∗ CU (Q ∗ U U ) −1 Q∗U C and we can then replace |Q∗C |U | † 2 with |Q∗C |U | 1 2 in the expression of πAX(b).
Theorem 2. Under Assumption 1 it follows that X|AX = b ∼ NC ( QX|bµ̃,QX|b ) I (AX = b) , (4)
where QX|b = T > U ,Q ∗ U U TU , is positive semi-definite with rank n− s− (k− k0) and µ̃ = T >µ̃∗ with µ̃∗ = [ b∗
µ∗U −Q ∗† U U Q ∗ U C (b ∗ − µ∗C )
] .
Note that QX|bA = 0, which implies that the right side of (4) is a (possibly intrinsic) density with respect to Lebesgue measure on the level set {x : Ax = b}. Further, note that Q∗U U TU E0 = 0, which implies that X is improper on the span of TU E0.
If Q is positive definite we get the following corollary. Corollary 1. Under Assumption 1 with s = 0, we have
X|AX = b ∼ N ( µ̃, Σ̃ ) I (AX = b) , (5)
where Σ̃ = T>U , (Q ∗ U U ) −1 TU , is a positive semi-definite matrix of rank n − k and µ̃ = T>µ̃∗ with µ̃∗ = [ b∗
µ∗U − (Q ∗ U U ) −1 Q∗U C (b ∗ − µ∗C )
] .

heading: 3.3 Sampling and likelihood evaluations
text: The standard method for sampling a GMRF X ∼ N ( µ,Q−1 ) is to first compute the Cholesky factor R of Q, then sample Z ∼ N (0, I), and finally set
X = µ+ R−1Z. (6)
To sample X|AX = b we use this method in combination with Theorem 2 as shown in Algorithm 3. The cost of using the algorithm for sampling, and for computing the expectation of X in Theorem 2, is dominated by CQ∗U U given that T has been pre-computed. Similarly, the cost for evaluating the likelihood in Theorem 1 is dominated by the costs of the Cholesky factors CQ∗U U + CQ + CAA> . These costs are not directly comparable to costs of the methods from Section 2 since they involve operations with the transformed precision matrix Q∗ = TQT> which may have a different, and often denser, sparsity structure than Q. In fact if T is dense the method will not be practically useful since even the construction of Q∗ would be O ( n2 ) . Thus, to understand the computational cost we must understand the sparsity structure of the transformed matrix. To that end, first note that only the rows id(A) in Q∗ will have a sparsity structure that is different from that in Q. In general, the variables involved for the ith constraint, id(AiX), will in the constrained distribution share all their neighbors. This implies that if i ∈ id(A), then |Q∗i,j | > 0 if |Qi,j | > 0 and we might have |Q∗i,j | > 0 if ∑ k∈id(A) |Q ∗ k,j | > 0. This provides a worst-case scenario for the amount of non-zero elements in Q∗, where we see that the sparsity of the constraints is important.

heading: 4 GMRFs under hard and soft constraints
text: As previously mentioned, one can view observations of a GMRF as hard constraints. In many cases, these observations are assumed to be taken under Gaussian measurement noise, which can be seen as soft constraints on the GMRF. It is therefore common to have models with both soft and hard constraints (e.g., a model with noisy observations of a field with local sum-to-zero constraints). Here, we extend the methods of the previous section to this case. Specifically, we consider the following hierarchical model

heading: X ∼ NC (Qµ,Q) , subject to AX = b,
text: Y ∼ N ( BX, σ2Y I ) ,
(7)
where Y ∈ Rm represent noisy observations of the linear combinations BX of X ∈ Rn, withm ≤ n, satisfying Assumption 1, and B is an m× n matrix with rank m. To deal with this type of models we present two results in this section. First, Theorem 3 shows how to compute the likelihood of the model. Second, the result in Theorem 4 can be used to efficiently compute the mean of X given the constraints and to sample from it.
We use the hat notation – like Q̂ – to denote quantities for distributions conditionally on the observations Y = y. We also use the notation from Theorem 2 and additionally introduce B∗ = BT> and y∗ = y −BT>C b
∗. We start by deriving the likelihood, πY|AX(y|b), which is needed for inference. Theorem 3. For the model in (7) one has
πY|AX(y|b) = σ−mY |Q ∗ U U | † 2
(2π) c0 |Q̂ ∗ U U | † 2
exp ( −1
2
[ y∗Ty∗
σ2Y + µ̃∗>U Q ∗ U U µ̃ ∗ U − µ̂ ∗> U Q̂ ∗ U U µ̂ ∗ U
]) ,
where c0 > 0, and
Q̂ ∗ U U = Q ∗ U U + 1
σ2Y (B∗U ) > B∗U ,
µ̂∗U = Q̂ ∗† U U ( Q∗U U µ̃ ∗ U + 1
σ2Y (B∗U )
> y∗ ) .
(8)
The computational cost of evaluating the likelihood is CQ̂∗U U + SQ̂∗U U + CQ∗U U . The following theorem introduces the distribution of X given the event {AX = b,Y = y}, which for example is needed when the model is used for prediction.
Theorem 4. For X in (7) one has X|{AX = b,Y = y} ∼ N ( µ̂, Q̂ ) where Q̂ = T>U ,Q̂ ∗ U U TU ,
and µ̂ = T> [ b∗
µ̂∗U
] . Here Q̂ ∗ U U and µ̂ ∗ U are given in (8). Further, let EQ∗U U be the null space of
Q∗U U then rank(Q̂) = n− s− (k − k0) + rank(B ∗ U EQ∗U U ).
Since the distribution in the theorem is a normal distribution, we can sample from X given the event {AX = b,Y = y} using sparse Cholesky factorization as shown in Algorithm 4.

heading: 5 Constrained Gaussian processes and the SPDE approach
text: Gaussian processes and random fields are typically specified in terms of their mean and covariance functions. However, a problem with any covariance-based Gaussian model is the computational cost for inference and simulation. Several authors have proposed solutions to this problem, and one particularly important solution is the GMRF approximation by [18]. This method is applicable to GPs and random fields with Matérn covariance functions,
r(h) = σ2
Γ(ν)2ν−1 (κh)νKν(κh), h ≥ 0,
which is the most popular covariance model is spatial statistics, inverse problems and machine learning [12, 26]. The method relies on the fact that a Gaussian random field X(s) on Rd with a Matérn covariance function can be represented as a solution to the SPDE
(κ2 −∆)α2 X = φW, (9)
Algorithm 3 Sampling X ∼ N ( µ,Q−1 ) subject to AX = b. Require: A,b,Q,µ,T
1: C ← 1 : nrow(A) 2: U ← (nrow(A) + 1) : ncol(A) 3: Q∗ ← TQT> 4: R← chol(Q∗U U ) 5: b∗ ← solve (( AT> ) CC ,b ) 6: m∗ ← µU− solve(R>,Q∗U C (b
∗ −TCµ)) 7: Sample Z ∼ N (0, IU U ) 8: X∗ ← [b∗, solve(R,m∗ + Z)]>
9: X← T>X∗ 10: Return X
Algorithm 4 Sampling X ∼ N ( µ,Q−1 ) subject to
AX = b and Y = y where Y ∼ N ( BX, σ2Y I ) .
Require: A,b,Q,µ,T,y,B, σ2Y 1: C ← 1 : nrow(A) 2: U ← (nrow(A) + 1) : ncol(A) 3: Q∗ ← TQT> 4: B∗ ← BT> 5: R← chol(Q∗U U + 1σ2Y (B ∗) > B∗)
6: b∗ ← solve (( AT> ) CC ,b ) 7: y∗ ← y −BT>C b ∗
8: m∗ ← solve ( R>,Q∗U U TU µ+
1 σ2Y
(B∗) > y∗ −Q∗U C (b ∗ −TCµ) ) 9: Sample Z ∼ N (0, IU U )
10: X∗ ← [b∗, solve(R,m∗ + Z)]> 11: X← T>X∗ 12: Return X
where the exponent α is related to ν via the relation α = ν + d/2, ∆ is the Laplacian,W is Gaussian white noise on Rd, and φ is a constant that controls the variance of X . The GMRF approximation by [18] is based on restricting (9) to a bounded domain D, imposing homogeneous Neumann boundary conditions on the operator, and approximating the solution via a finite element method (FEM). The resulting approximation is Xh(s) = ∑n i=1Xiϕi(s), where {ϕi(s)} are piecewise linear basis functions induced by a triangulation of the domain, and the vector X with all weights Xi is a centered multivariate Gaussian distribution. This can be done for any α > d/2 [4], but the case α ∈ N is of particular importance since X then is a GMRF. In particular, when α = 2 and φ = 1, the precision matrix of X is Q = (κ2C + G)C−1(κ2C + G), where C is a diagonal matrix with diagonal elements Cii = ∫ ϕi(s)ds, and G is a sparse matrix with elements Gij = ∫ ϕi(s)ϕj(s)ds.
Clearly, a linear constraint on Xh(s) can be written as a constraint on X. For example, if Xh(s) is observed at a location in a given triangle, it creates a linear constraint on the three variables in X corresponding to the corners of the triangle. Thus, if we draw some observation locations s1, . . . , sk in the domain, we can write Y = (Xh(s1), . . . , Xh(sk))> = AY X where AY is a k × n matrix with (AY )ij = ϕj(si). A model where a Gaussian Matérn fields is observed without measurement noise can therefore be handled efficiently by combining the SPDE approach with the methods from Section 3. The next section contains a simulation study that compares this combined approach with a standard covariance-based approach in terms of computational cost.
Through the nested SPDE approach in [5] one can also construct computationally efficient representations of differentiated Gaussian Matérn fields like U(s) = (v>∇)X(s), where v>∇ is the directional derivative in the direction given by the vector v and X(s) is a sufficiently differentiable Matérn field. A FEM approximation of this model can be written as Uh(s) = ∑n i=1 Uiϕi(s) where now U ∼ N (0,AUQ−1A>U ). Here AU is a sparse matrix representing the directional derivative and Q is the precision matrix of the GMRF representation of X(s) [5]. If we introduce X ∼ N (0,Q−1), we may write U = AUX, and we can thus enforce a restriction on the directional derivative of X as a linear restriction AUX = b. As an example, v = (1, 1)> and b = 0 results in the restriction ∂ ∂ s1
X(s) + ∂∂ s2X(s) = 0, or in other words that the field is divergence-free. In the next section we use this in combination with the methods in Section 4 to construct a computationally efficient Gaussian process regression under linear constraints.

heading: 6 Numerical illustrations
text: In this section we present two applications. In both cases, timings are obtained through R [25] implementations, available in the CB R package [6], run on an iMac Pro computer with a 3.2 GHz Intel Xeon processor.

heading: 6.1 Observations as hard constraints
text: Suppose that we have an SPDE approximation Xh(s) of a Gaussian Matérn field X(s), as described above with D = [0, 1]2 and a triangulation for the GMRF approximation that is based on a uniform mesh with 100× 100 nodes in D. We consider the costs of sampling Xh(s) conditionally on k point observations without measurement noise, and of log-likelihood evaluations for these observations. In both cases, the observations are simulated using the parameters κ2 = 0.5, φ = 1 and α = 2 or α = 4.
In the top left panel of Figure 1 we show the computation time for one log-likelihood evaluation using the standard method from Section 2 on the GMRF X of weights for the basis expansion of Xh(s). The panel also shows the corresponding computation time for the new method from Section 3. The computation times are evaluated for different values of k, where for each k the observation locations are sampled uniformly over triangles, and uniformly within triangles, under the restriction that there can be only one observation per triangle, which guarantees that AY Q−1A>Y has full rank. In each iteration, the values of κ2 and φ that are evaluated in the likelihood are sampled from a uniform distribution on [1, 2]. The curves shown in the figure are computed as averages of 10 repetitions for each value of k. As a benchmark, we show the time it takes to evaluate the log-likelihood assuming that X(s) is a Gaussian Matérn field, which means that we evaluate the log-likelihood `(Y) of a k-dimensionalN (0,Σ) distribution without using any sparsity properties. This is done by calculating the Cholesky factor R of Σ and then computing `(Y) = − ∑k i=1 logRii − 1 2Y >R−1R−>Y.
An alternative method for this example would be to adapt the mesh to the observation locations, so that one only has observations at mesh nodes. This is, however, in general not a good solution as it might require using poor meshes (from a numerical stability point of view) or using meshes with a large number of nodes, increasing the computational cost, and hence we do not explore it here.
The covariance-based method is the fastest up to approximately k = 1000 observations, since the problem then is too small for sparsity to be beneficial, whereas the new method is fastest for k > 1000. If k is large relative to the total size of the GMRF, constraints may no longer be local since many constraints will be interacting. As illustrated in the bottom right panel of Figure 1, this causes the largest set of connected observations to increase rapidly for k > 4000. The computation time for creating T will roughly scale cubically in the size of the largest set of connected observations, which is the main reason for the large increase in the computation time for the new method when k > 5000. However, it should be noted that the construction of the basis needed for the new method only has to be done once, so even though the computation time for the construction of T is high when the largest set of connected observations becomes large, the method will still be efficient for maximum likelihood estimation using numerical optimization, where several likelihood evaluations are needed.
The specific values of κ and φ have no effect on the computational costs for any of the methods. The only parameter that affects the results is α since it changes the sparsity structure of the GMRF where a larger value results in precision matrices with more non-zero elements. The bottom left panel in Figure 1 shows the results for α = 4, where we can note that the standard GMRF method is slower than the covariance-based method for all values of k, whereas our proposed method still performs the best for large values of k. In the supplementary materials, corresponding results for α = 1 and α = 3 are shown with similar results.
In the top right panel of Figure 1 we show the time needed to sample Xh(s) conditionally on the observations AY X = y, i.e., to simulate from X|AY X = y. Both the old method (conditioning by kriging) and the new method (using (6) with mean and precision from Theorem 2) are shown. We do not show the covariance-based method since it is much slower than both GMRF methods. Also here the displayed values are averages of 10 repetitions for each value of k, and for each repetition the simulation is performed using values of κ2 and φ that sampled from a uniform distribution on [1, 2]. The results are again shown for α = 2 and α = 4, and the results for α = 1, 3 are presented in the supplementary materials. The results are similar to those for likelihood evaluations, where the new method is much faster for large numbers of observations.

heading: 6.2 Gaussian process regression with linear constraints
text: We now consider an application to GP regression from [13], in which one is given noisy observations Yi = f(si) + εi, with εi ∼ N (0, σ2eI) of a bivariate function f = (f1, f2) : R2 → R2 where f1(s) = e
−as1s2(as1 sin(s1s2)− s1 cos(s1s2)) and f2(s) = e−as1s2(s2 sin(s1s2)−as2 sin(s1s2)). The goal is to use Gaussian process regression to reconstruct f , under the assumption that we know that it is divergence-free, i.e., ∂∂ s1 f + ∂ ∂ s2
f = 0. We thus want to improve the regression estimate by incorporating this information in the GP prior for f . This can be done as in [32, 37, 13] by encoding the information directly in the covariance function, or by imposing the restriction through a hard constraint at each spatial location through the nested SPDE approach. This is done by setting B = AY and A = AU in (7) where the matrices AY and AU are defined in Section 5.
Following [13], we choose a = 0.01 and σ2 = 10−4 and generate 50 observations at randomly selected locations in [0, 4]× [0, 4] and predict the function f at N2 = 202 regularly spaced locations in the square. Independent Matérn priors with α = 4 are assumed for f1 and f2 and the covariancebased approach by [13] is taken as a baseline method. As an alternative, we consider the SPDE approximation of the Matérn priors, with n basis functions obtained from a regular triangulation of an extended domain [−2, 6] × [−2, 6] (the extension is added to reduce boundary effects). To be able to use Algorithm 2, we only enforce the divergence constraint at every third node for the SPDE model. This procedure can be seen as an approximation of the divergence operator that will converge to the true operator when the number of basis functions increases.
The parameters of the baseline model and of the SPDE model are estimated using maximum likelihood, where the likelihood for the SPDE model is computed using Theorem 3. The function f is then reconstructed using the posterior mean of the GP given the data, which is calculated using Theorem 4. This experiment is repeated for 50 randomly generated datasets. In the left panel of Figure 2 we show the average root mean squared error (RMSE) for the reconstruction of f for the SPDE model as a function of n, based on these 50 simulations, together with the corresponding RMSE of the baseline method. The shaded region for the SPDE model is a pointwise 95% confidence band. One can see that the SPDE model gives a comparable RMSE as long as n is large enough.
We next fix n = 3600 and consider the time it takes to compute a prediction of f given the estimated parameters. In the right panel of Figure 2 we show this computation time as a function of the number of observations, m, for the baseline method and for the SPDE-based method. Also here we see that the covariance-based method is the fastest for small numbers of observations, whereas the GMRF method (that has a computational cost that scales with the number of basis functions of the SPDE approximation rather than with the number of observations) is fastest whenever m > 600.

heading: 7 Discussion
text: We have proposed new methods for GMRFs under linear constraints, which can greatly reduce computational costs for models with a large number of constraints. In addition, we showed how to combine these methods with the SPDE approach to allow for computationally efficient linearly constrained Gaussian process regression. Some recent papers on constrained GPs, such as [9], consider methods that are similar in spirit to those we have developed here. However, to the best of our knowledge, our methods are the first to account for sparsity, which is crucial for GMRFs.
Clearly the proposed methods will not be beneficial if the number of constraints is small. Another limitation is that the methods are only efficient if the constraints are sparse. For instance, our proposed method does not work for a sum-to-zero constraint ∑n i=1Xi = 0, which is commonly used in hierarchical models to ensure identifiability. An interesting topic for future research is to handle problems where both sparse and dense constraints are included. In that case one could combine the proposed method with a conditioning by kriging approach where the dense constraints are handled in a post-processing step as described in Section 2.
We have only considered exact methods in this work, but if one is willing to relax this requirement, an interesting alternative is the iterative Krylov subspace methods by [35]. Comparing, or combining, the proposed methods with those in [35] is thus another interesting topic for future work. A further question is whether Algorithm 2 can be improved in terms of computational cost. We have used a simple SVD implementation and for cases when the construction of the basis dominates the total cost it could worth exploring whether more modern iterative methods such as [1, 22] could be applicable. How that could be done is not clear since most such methods provide a low rank solution whereas we need the full matrix V>. A related idea, when the exact methods presented here are not suitable, is to use subspace embeddings or sketch-based methods [19] to approximate the linear constraints. It is, however, an open question if one can get good results through such an approach, and how one could describe the accuracy of the corresponding approximated distribution.
The SPDE approach also allows for more flexible non-stationary covariance structures like the generalized Whittle–Matérn models [4]. Our proposed methods are directly applicable to these models in the Markov case (with integer smoothness), and can also be extend to the case with general smoothness by combining the constraint basis with the rational approximations in [4].
Finally, we can think of no potential negative societal impacts that this work may have, given that it is solely concerned with improving the performance of existing methods.

Review summaries:
This paper introduces efficient methods for GP inference under the Markov assumption when conditioning on sparse linear combinations. The approach relies on astute changes of basis. It assumes that the matrix A representing the linear mapping involved in the conditioning can be represented in block diagonal form, enabling numerical savings in the resulting SVD. Decomposing the finite dimensional space of interest appropriately also enables enjoying elegant and practical results, as illustrated by their application to GP/GMRF modelling under hard and soft (sparse) constraints. A noticeable fact about the presented approach ist hat it also applies to intrinsic GMRFs. 

The paper proposes a method for evaluating the likelihood and for sampling from the posterior of Gaussian Markov random field $X$ given some linear constraints $X | AX = b$ where $A$ is a sparse matrix. The method is based on the SVD of $A$, which may be computed for a small cost if one can identify rows of $A$ that don't share columns with non-zero elements.   The method is illustrated on two toy examples, where it is shown that it provides massive computational gains. It is however unclear if the competitors considered in the benchmark are the actual alternatives a practitioner would consider.  ----------------------- The author's rebuttal makes me confident that the shortcomings I've identified have been understood by the authors and that they are in a good position to address them at the camera ready stage. I'm thus upgrading my score from 5 to 6.

This paper formulates a class of computationally efficient methods for Gaussian Markov Random field (GMRF) models with linear constraints. The main idea behind the methods relies on a transformation of basis such that the constraints are easier to enforce upon transformation. The methods outperform existing alternative formulations of GMRFs under scenarios involving large numbers of sparse constraints. Finally, a GMRF for linearly constrained Gaussian process regression is derived by utilizing the proposed methods in conjunction with the stochastic partial differential equation (SPDE) approach.