Efficient methods for Gaussian Markov random fields under sparse linear constraints
Authors: David Bolin
Jonas Wallin
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian processes
GMRF
Markov random fields
Gaussian processes regression
linear constraints

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.

Efficient methods for Gaussian Markov random fields under sparse linear constraints
David Bolin
Jonas Wallin Department of Statistics, Lund University jonas wwwlwalllll.. meallindstat.lng lu.se
King Abdullah University of Science and Technology david. bolinkkaust....mmmmmmm edu.
Abstract
Methods for inference and simulation of linearly constrained Gaussian Markov Random Fields (GMRF) are computationally prohibitive when the number 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 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 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 GMRF setting to reduce computational cost. This is illustrated in two applications with simulated data.
Introduction
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 cquilibrium conditions [33 16 17 30 31 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 [4] One way to reduce computational burden is to impose conditional independence assumptions. In fact, conditional independence between random variables is often explicitly implicitly assumed 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
widely used in areas ranging from brain imaging [24] 34 20 to spatial statistics ( and time series analysis (28). GMRFs further arise as computationally efficient approximations of certain GPs (18), which is fundamental modeling tool in both machine learning and statistics. In particular, such approximations
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 [X,...,, X,] ~ (m,, Q-1),
has sparse precision (inverse covariance) matrix Q which enables the 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 indepen- dence assumptions since Qij = if and only if X, and X, are independent conditionally on all other
variables in X [see 28, Chapter 2].
35th Conference on Ncural Information Processing Systems (NcurIPS 2021).

---

For GMRF X, a set of & linear constraints can be formulated as
where each row in the x n matrix A and the vector encodes a constraint on X. Note that observations of GMRFs can be formulated as linear constraints, where, e.g.. a point observation of X; is the simple constraint X; = Tj. These deterministic restrictions are often referred to as hard constraints [28) If (2) is assumed to hold up to Gaussian noise, i..., AX A/b,, o I²I), 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, 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 computational cost that scales cubically in the number of constraints k. This is prohibitive when there are many constraints, which 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 and in the second it can be cubic in N if the observations are such that in 22) is not diagonal, which, e.g., is the case for area observations (1 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 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 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 of eigenvalues of Q is zero. This makes their distributions invariant to shifts in certain directions, which is useful property for prior distributions of Bayesian models (28 Chapter 3]. Because of this
intrinsic GMRFs, which are improper in the sense that set
they are often used in areas such as semiparametric regression [10). medical image analysis (24, [34].
and geostatistics [15] 7, and 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 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.
2 Standard methods for GMRFs under hard constraints
Hard constraints can be divided into interacting and non-interacting constraints. In the latter, specifies a constraint on a subset of the variables in X so that AX b can be written as X, = be. where denotes subset of the indices. Specifically, let X., denote the remaining variables, then
X. X
[e mu
Qoc
Qcu Quu
and X|AX = -bb-b~a b~ (
be [Mlcc]
x
~
N
Que
Q
2

---

where Pule = M. - Quc -cp). 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 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 b and for evaluating its log-likelihood function. It straightforward show that X b N Σ), where ς = μ - Q A (AQ A AT)- - b) and Σ = Q - QA'--'A' A (AQ A 'AT)- AQ Since Σ has rank
k. likelihood evaluation and sampling is in general expensive. For example, one way is to use 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, E and j cannot be constructed through the expressions above since Q has unbounded eigenvalues.
commonly used method for sampling under hard linear constraints, sometimes referred to as conditioning by kriging (27 is to first sample X N (M,, Q) and then correct for the constraints by using X = X Q Q A (AQ A (AX -0-) b) as a sample from the conditional distribution. Here the cost of sampling X is Cq Sq. where Co denotes the computational cost of sparse Cholesky factorization Q = R R and Sq the cost of solving Rx = for x given u. Adding the cost for the correction step, the total cost of the method is O(Cq + (k: + 2)SQ + k*).
Let Axx denote the density of Ax (( AQ A ). then the likelihood of X AX- b can be computed through the expression
T Axxx (b(x)()
x(x[Ax=b) b)
Tx (b)
where T xxxxx((xx = - 1- and I(Ax = b) denotes the indicator function with [[Ax= = b) = 1 Ax =: b and I(Ax b) = O 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 OCq + (A: + 1)8Q + ).)).
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- -1 A due to the rank deficiency of Q.
3 The basis transformation method
In this section, we propose the new class of methods in two steps. We first derive a change of basis in [3.I 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 GMRFss the definition in (1) is inconvenient since the covariance matrix has infinite eigenvalues. Instead one can use the canonical parametrization X ~ NC (Mc, Q), which implies that the density of X is given by (xx x exp (- Qx + xxx) and thus that fcc = Qu.. 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 matrix B by B' and for symmetric positive semi definite matrix M we define the pseudo determinant as M = II 1= so A where { } are the eigenvalues of M. For the remainder of this article, will assume the following.
Assumption 1. X No (Q/,, Q) where Q positive semi-definite n x n matrix with rank 8 > 0 and null-space EQ. A is k matrix with rank k and we let ko denote rank(A (AEq).
Finally, use the index notation 8 -1,...,,k k) and 2 (k + .....,,) n) that, e.g., Ve denotes the sub-vector (vi,..., vx,) T of the first k elements in the vector v, M. w (M denotes the submatrix obtained by extracting the columns (rows) in M with indices in 2, and M M the submatrix obtained by extracting the columns and rows in M with indices in 2.
3

---

3.1 Basis construction
Our main idea to construct a basis on R" such that the constraints are easily enforced. 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 the matrix A, the basis V is orthonormal and the first & rows span the image of A and the last k rows span the null-space of A. Now, if we let x" denote vector expressed in the basis V then x" can be transformed back to the natural basis by Vx" hence
Ax + USx* = b
( 0jx² =b
US = b + X b"
where b" = (US This shows that the SVD a natural choice of method for constructing the basis of interest, which is defined by an n x n change-of matrix T V
In Algorithm we present a simple method to build T based on the SVD. In the al- gorithm, id(A)
function that returns the indices of the non-zero columns in A. The computational cost of the method is dominated by that of the SVD, which is O (k* (A)) [11 p. 493]. Clearly, the cubic scaling the number of constraints may reduce the efficiency of any method that re- quires this basis construction as a first step. However, suppose that the rows of A can be split into two sub-matrices, A and Ag. which have no common non-zero columns. Then the SVD of the two matrices can be computed sep- arately. Suppose now that A corresponds to m such sub-matrices then, after reordering, A [A, ] where ! repre-
Algorithm 1 Constraint basis construction.
Require: A (a k x matrix of rank k)
I: T < In
2: D +- id(A)
3: USV
sud( (..))
4: TD. V
5: T [T.D T.p=
6: Return T
sent consttratiss such that id(A,) n id(A,) = for all and I. By replacing the SVD of Algorithm by the SVDs of the lower-dimensional matrices the computational cost is reduced to O Σ i= rank(A, A + rank( (ank(A,) [id(( A.)). This method is presented in Algorithm 2 The
reordering step is easy to perform and is described the supplementary materials.
3.2 Conditional distributions
Using the change of basis from the previous subsection, we can now derive alternative formulations of the distributions AXX 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 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" No (Q'p² p² Q"), where μ Tp 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 l matrix of rank k)
10: Th:(k+u),D (V
1: T < In
11:
h h u
2: Djull < id(A)
12:
if D > u then
3: h
13:
TR(4+DI-u),D + (V 1 1+ DD- -ut1
(w+1):)D)D
4: 1+
14:
5: Reorder that A A
A m
15: end if
6: for == m. do
16: end for
7:
D id,
17: T --
T T.pp (L)).
8: USV + sed As)
18: Return T
9:
+ (U)

---

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 vector and by multiplying with T from the left and T from the right for matrix.
We use the matrix therefore has inverse H
(AT))
Co which is equal to US from the SVD in Algorithm and U Finally, recall that b* -.
Theorem 1. Under Assumption 7 it
follows that
Tbx(b)= (2T) AA 11/2 exp (-2(0-- m&e) Q (b" -
(b-
(b--
m
where Q/ew QEe - Qeo (Q w/ Qwe and 1 1= Qu 1-
Q is positive definite we have Q = Qee - Vee (Qwr) Qie and we can then replace 1e with Q the expression of (b). Theorem 2. Under Assumption it follows that
x|AX No (Axiatl Qx1s) I(AX =b);
(4)
where Qx Ty Qe Tw is positive semi-definite with rank (k: ko) and T T with μi" Mis Q 2se (b' m.)
b'
Note that Qx A O, which implies that the right side of 4 is (possibly intrinsic) density with
respect to Lebesgue measure on the level set {x Ax = b). Further, note that Q T Eg = 0, which implies that X is improper the span of Eo
If Q is positive definite we get the following corollary.
Corollary 1. Under Assumption / with O. we have
XAA= (i. I(AX (IE))
(5)
where E = T x Tw, is positive semi-definite matrix of rank k and jı T i
with μ pb (Quw) Que (b' pb)
3.3 Sampling and likelihood evaluations
The standard method for sampling a GMRF X ~ N (14, Q-1) is to first compute the Cholesky factor R of Q, then sample Z . (0, I), and finally set
x= '.
(6)
To sample X = b we use this method in combination with Theorem as shown Algorithm
3
The cost of using the algorithm for sampling, and for computing the expectation of X in Theorem dominated by Cq given that T has been pre-computed. Similarly, the cost for evaluating the likelihood in Theorem 1 dominated by the costs of the Cholesky factors Co w + Cq + CAAT T.
These costs are not directly comparable to costs the methods from Section 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 (n²). Thus, to understand the computational cost must understand the sparsity structure of the transformed matrix. To that end, first note that only the rows id(A) (A) in Q" will have a sparsity structure that different from that in Q. In general, the variables involved for the ith constraint, id( AX)), will in the constrained distribution share all their neighbors. This implies that if E id(A), then [ > 0 QQ,,1 >
and we might have |Q, > if Of ERcid(A) Qigg > O. This provides a worst-case scenario for the amount of non-zero elements in Q" where see that the sparsity of the constraints is important.
5

---

GMRFs under hard and soft constraints
As previously mentioned. one can view observations of GMRF hard constraints. In many cases, these observations are assumed to be taken under Gaussian measurement noise, which can be seen soft constraints on the GMRF. It is therefore common have models with both soft and hard constraints (e.g., a model with noisy observations of a field with local sum-to-zer constraints). Here, we extend the methods the previous section this case. Specifically, we consider the following hierarchical model
X ~ No (Qu, Q), subject to AX =
Y ~~~ W of
where Y I represent noisy observations of the linear combinations BX of E R", with m S n. satisfying Assumption| D and B is an m X matrix with rank m. To deal with this type of models we present two results in this section. First, Theorem 3 how 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 - to denote quantities for distributions conditionally on the observa- tions Y y. We also use the notation from Theorem 2 and additionally introduce B BT and y = y BT We start by deriving the likelihood, (yjax(y(b) which is needed for inference. Theorem 3. For the model in 7 one has
(y(b) -x 2
y
T + uu
T.T Quu D2 a ))
exp
(2T) 0 1
oy
where co > 0, and
Q Qu = Q +
(B ) B
oy
(8)
Pi = Q ai
(Quw + og
(Bx) 1 )r)
The computational cost of evaluating the likelihood is Co. theorem introduces the distribution of X given the event (AX- = b, y), which for example is needed when the model is used for prediction.
+ +o " The following
Theorem 4. For X in (7) one has X|{AX = b, b,Y (R, Q) where Q T T Tw,
and ja T [i Here Q' and aad are given in 8 Further, let EQuuw be the rull space of
Qie v then rank Q) = - kk- + rank B ) )
Since the distribution in the theorem is normal distribution,
sample from X given the event (AX b, Y T=yy yy} using sparse Cholesky factorization as shown Algorithm
Constrained Gaussian processes and the SPDE approach
Gaussian processes and random fields typically specified in terms of their mean and covariance functions. However, a problem with any novariance-based Gaussian model the computational cost for inference and simulation. Several authors have proposed solutions 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,
o²
(h)
(xh) (sh),
1e 2 0,
which is the most popular covariance model is spatial statisties, 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 solution the SPDE
(k² x oW,
(9)

---

Algorithm 3 Sampling X N (H,
Algorithm 4 Sampling X N N Q Q-1 subject to
subject to AX- b.
AX b and = where Y (BX,
Require: A, b, Q, pa, T
Require: A, b, Q, p. T. B,
I: 8 + nrooo A)
I: E 1: nroww ()
2 t (nrow(A)) AA +1) (A)
2: V + (nrow A) + 1) ncol(A)
Q" f- TQT
3: Q" + TQT
4: R + chol( QQ
4: B + BT
5: b' t solve (AT)) )(6,)))
5: R < chol chol(Quv 1 (B" B*)
(B")
m² ---
6: b" solve (AT)) leee
solve(R. T, Qere (b" -T Txıl))
7: y y y - BT
7: Sample Z N (0, Iww)
8: X -- -lbb, solve m" +Z)]
8: m² + solve
R Q w T wh+
9: + '
(B²) y Svr (b" -Teu))
10: Return X
9: Sample Z N /,,00, I0,
10: X + bb, solve R, 1mm + 2)]
11: X +- T
12: Return
where the exponent a' is related to v via the relation a = = 1 d/2, A 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 (2) 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 XA (s) = Σx=1 (s), where {x)) are piecewise linear basis functions induced by a triangulation of the domain, and the vector X with all weights X; is 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 GMRE In particular, when a = 2 and á = 1, the precision matrix of X is Q (R&C + G)C (x²C + G), where C is a diagonal matrix with diagonal (s)(sss. and G is a sparse matrix with elements Giz = (s))s)
elements Ca
Clearly, linear constraint (s) can be written
constraint X. For example, if X (s) is observed at a location in 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 $1 $1...., 5%k in the domain, we can write Y (Xx(s)).... with (Ay), = (s,).. ) 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- baased approach in terms of computational cost.
(sx)) = Ay X where Ay is k x n matrix
Through the nested SPDE approach in (5) one can also construct computationally efficient repre- sentations of differentiated Gaussian Matérn fields like (ss) (v V (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 (s Σm= (s) where now U (0, A Q A). Here Ay is sparse matrix representing the directional derivative and Q is the precision matrix of the GMRF representation of (s) (5). If we introduce X ~O,, Q-1). we may write U = A X, and we can thus enforce a restriction on the directional derivative of X as a linear restriction Ar = b. As an example, v = (1, and
results in the restriction Bal (s) + 0 X(s) = 0, or in other words that the field 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.
7

---

Likelihood, = 2
Sampling
100
10
Evataatoon tonnttrell
Samping time
0.1
0.01
GMRF new GMRE cid
Covarance Basis build
GMRF olt, a-2 GMRF 0k6,
GMRF new, a-2
0.001
GMRF
0.01
1000
2000
3000
4000
5000
1000
2000
3000
4000
5000
k
Likelihood, a =4
Largest connected set
100
Eveluston thee 65
10
spee Serga conneccert
0.1
0.01
GMRF new GMRF cld
Covarance
0.001
Basis build
1000
2000
3000
4000
5000
1000 2000
3000
4000
5000
k
Figure I: Average computation times (based on 10 replications) for one likelihood evaluation (left) and one sample from X|AyX y (top right) of the Matérn model with functions of the number of observations includes the time needed construct the basis matrix T, which also is shown separately and which depends on the largest set of connected observations (bottom right).
2 and 4 as
The computation times for the new GMRF method
Numerical illustrations
In this section we present two applications. In both cases, timings are obtained through R 25 implementations, available in the R package (6). run on an iMac Pro computer with a 3.2 GHz Intel Xeon processor.
6.1 Observations as hard constraints
Suppose that we have an SPDE approximation (x(s) of Gaussian Matérn field X (s)), as described above with D [0, 1]2 and triangulation for the GMRF approximation that is based uniform mesh with 100 x 100 nodes in D. We consider the costs of sampling (s(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 x² 0.5, & 1 and a 2 or a 4.
In the top left panel of Figure Twe we show the computation time for one log ligeliiolooodd evaluation using the standard method from Section 2 on the GMRF of weights for the basis expansion of (s)s) The panel also shows the corresponding computation time for the computation times are evaluated for different values of k. where for each & the observation locations sampled uniformly over triangles, and uniformly within triangles, under the restriction that there can be only one observation per triangle, which guarantees that ArQ A has full rank. In each iteration, the values of x² and á that are evaluated in the likelihood are sampled from a uniform distribution [1, 2). The curves shown in the figure are computed 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 (s) is a Gaussian Matérn field, which means that we evaluate the log likelihood !(Y) of k k-dimensional (O, Σ) distribution without using any sparsity properties. This is done by calculating the Cholesky factor R of Σ and then computing {(Y) E log R. Y R -R-TRT- Y..
method from Section 3 The
An alternative method for this example would be to adapt the mesh to 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 do not explore it here.

---

The covariance- prrreance-baeed 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 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 [ 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 A >> 5000. However, it should be noted that the construction of the basis needed for the new method only has to be done once, 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
needed.
The specific values of x and & have no effect on the computational costs for any of the methods. The only parameter that affects the results is a 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 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 and 3 are shown with similar results.
In the top right panel of Figure ! we show the time needed to sample (s)ss conditionally on the observations X y. i... to simulate from Ay = 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 x² and & that sampled from a uniform distribution on [L. 2]. The results are again shown for 2 and a = 4, and the results for = 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.
6.2 Gaussian process regression with linear constraints
We now consider an application to GP regression from (13), in which one is given noisy observations
Y, = f(s.) + Eq, with E, (0, I²I) of a bivariate function f = (f1, fe) R² + IR² where (i(s) = e aas1" sin($182) 82) - $1 c0s(8182) 52)) and fe(s) = e as1 * $1 (1 - as2 sin $1 $2)). The goal is to use Gaussian process regression reconstruct f, under the assumption that we know that it is divergence-free, i.e., J., + 8 = O. 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 σ₂ = 10 and generate 50 observations at randomly selected locations in (0, 4) [0, 4] and predict the function at 2 20 regularly spaced locations in the square. Independent Matérn priors with a = 4 are assumed for fı and fz and the covariance- based approach by (13) is taken as 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] x [-2, 6] (the extension 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 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 for the SPDE model as a function of n. based on these 50 sinulations, 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 comparable RMSE as long as n is large enough.

---

SPDE method
8
SPDE methed
Meeni
Covoriance -basdd method
RMMGE
Covariance -based method
Preede soom
8
2000
3000
4000
2000
5000
7003
8000
300 1000 1500
2000
Figure 2: Left: Average RMSE with corresponding 95% pointwise confidence band for 50 reconstruc- tions of f based SPDE method with different number of basis functions (red), the corresponding RMSE for the GP model (black) and the RMSE for estimating f by a constant equal to the mean of the observations (blue). Right: Average computation times for the prediction of as functions of the number of observations (m) with envelopes showing the range of times for the 50 replications.
We next fix n = 3600 and consider the time it takes to compute a prediction given the estimated parameters. In the right panel of Figure we show this computation time as function of the number of observations, m, for the baseline method and for the SPDE-based method. Also here we see that the covariance-baset 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 > 600.
7 Discussion
We have proposed new methods for GMRFs under linear constraints, which can greadl reduce computational costs for models with large number of constraints. In addition, we showed how combine these methods with the SPDE approach to allow for computationally efficient linearly constrained Gaussian process regression. Some recent papers 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 that the methods are only efficient the constraints are sparse. For instance, our proposed method does not work for sum-to-zero constraint Σx X; 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- processoog step as described Section
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. combining, the proposed methods with those in (35) thus another interesting topic for future work. further question is whether Algorithm 2 can be improved in terms of computational cost. We have used simple SVD implementation and for cases when the construction of the basis dominates the total cost could worth exploring whether more modern iterative methods such as [ [2 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
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 Thittle-MMténn models (4). Our proposed methods are directly applicable to these models 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 is solcly concerned with improving the performance of existing methods.
10

---

References
[1] Zeyuan Allen-Zhr and Yuanzhi Li. "LazySVD: Even Faster SVD Decomposition Yet Without Agonizing Pain". In: Advances in Neural Information Processing Systems. Ed. by D. Lee, M. Sugiyama, U. Luxburg, I. Guyon, and R. Gamett. Vol. 29. Curran Asso- ciates, Inc., 2016. URL: https prw.eedingee.s...cooo neurips. cc peeppaper// 2016 file/ c661968308591
hBt.://f2ccoocccccccccccccaaa www.ccccaaaa2aacaapppeeeo Paper.pdff
[2] Yoann Altmann, Stephen McLaughlin, and Nicolas Dobigeon. "Sampling from multivariate Gaussian distribution truncated on a simplex: A review". In: 2014 IEEE Workshop on Statistical Signal Processing (SSP). 2014, pp. 1 113-116. DOI: 10. 1109/SSP 2014. 6884588
[3] Haakon Bakka, Håvard Rue, Geir-Arne Fuglstad, Andrea Riebler, David Bolin, Janine Illian, Elias Krainski, Daniel Simpson, and Finn Lindgren. "Spatial modeling with R-INLA: a review". In: Wiley Interdiscip. Rev. Comput, Stat. 10.6 (2018), e c1443, 24. ISSN: 1939-5108.
DOT: 10. 1002/0icciis.. 1443
[4] David Bolin and Kristin Kirchner. "The rational SPDE approach for Gaussian random fields with general smoothness" In: J. Comput. Graph. Starist, 29.2 (2020), pp. 274-285, ISSN: 1061-8600. DOI: 10. 1000/10018800000
2019. 1665537
[5] David Bolin and Finn Lindgren. "Spatial models generated by nested stochastic partial differ- ential equations, with an application to global ozone mapping" In: Ann. Appl. Stat. 5.1 (2011),
pp. 523-550. ISSN: 1932-6157. DOI: 10. 1214/10-AOASA
[6] David Bolin and Jonas Wallin. CB: Contraint basis construction for GMRFs. Version 1.0.
2021. URL: https: /estps:/ggithnn wht.b.tttt......o/gg com/ hoam//Jnswaaann////
[7] David Bolin, Johan Lindström, Lars Eklundh, and Finn Lindgren. "Fast estimation of spatially dependent temporal vegetation trends using Gaussian Markov random fields". In: Comput. Statist. Data Anal. 53.8 (2009), pp. 2885-2896. ISSN: 0167-9473. DOT: 10. 1016/ csda
2008.09.017
[8] Yanqing Chen, Timothy A. Davis, William W. Hager, and Sivasankaran Rajamanickam. "Algorithm 887: CHOLMOD, superodal sparse Cholesky factorization and update/downdate" In: ACM Truns. Math. Software 35.3 (2008), Art. 22, 14. ISSN: 0098-3500. DOI: 10. 1 1145/ 20rg/10.11 1011/1111199900 1391995
1391989. 1391995 URL: /ttht:ttt//aatt
[9] Yulai Cong, Bo Chen, and Mingyuan Zhou. "Fast simulation of hyperplane hyperplane-truncatee trunccttd multi- variate normal distributions". In: Bayesian Anal. 12.4 (2017), pp. 1017-1037.
ISSN: 1936-0975.
DO1: 1et.1111117BA000
[10] Ludwig Fahrmeir and Stefan Lang. "Bayesian semiparametric regression analysis of mul- ticategorical time-space data". In: vol. 53. 1. Nonlinear non-Gaussian models and related URI.:
filtering methods (Tokyo, 2000). 2001, pp. 1 11-30. DOI: 10.
https://doi.
1023 1017904188111
[11] Gene H. Golub and Charles F. Van Loan. Matrix computations. Fourth. Johns Hopkins Stud- ies in the Mathematical Sciences. Johns Hopkins University Press, Baltimore, MD, 2013, 1-4214-0794-9; 978-1-4214-0859-1
pp. xiv+756. ISBN:
[12] Peter Guttorp and Tilmann Gneiting. "Studies in the history of probability and statistics. XLIX. On the Matérn correlation family". In: Biometrika 93.4 (2006), pp. 989-995 ISSN: 0006-3444
101:3/ www10/0099me..... http:/ttmmeett//////ee 989
[13] Carl Jidling. Niklas Wahlström, Adrian Wills, and Thomas B Schön. "Linearly constrained Gaussian processes". In: Advances in Neural Information Processing Systems 30. Ed. by I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett. Curran Associates, Inc., 2017, pp. 1215-1224.
[14] Carl Jidling, Johannes Hendriks, Niklas Wahlström, Alexander Gregg, Thomas B. Schön, Christopher Wensrich, and Adrian Wills. "Probabilistic modelling and reconstruction of strain". In: Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 436 (2018), pp. 141-155. ISSN: 0168-583X. DOI: htt:httss////doorrg
10.1011 1016/j nimb. 2018. 08.051
[15]
Hans R. Küinsch. "Intrinsic autoregressions and related models on the two-dimensional lattice". In: Biometrika 74.3 (1987), pp. 517-524. ISSN: 0006-3444, DOI: 10. w093///imet/74/////
517
URL: https: / /tt. wttp:////........ooommeeeee/////
1093/

---

[16] Markus
"Algorithmic Linearly Constrained Gaussian Processes". In: Ad- vances Neural Information Processing Systems 31. Ed. by Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Desa-Bianchi, and R. Garnett. Curran Associates, Inc., 2018, pp. 2137-2148.
[17] Markus Henge-Hegnmmnn.. "Linearly Constrained Gaussian Processes with Boundary Con- ditions". In: Proceedings of The 24th International Conference on Artificial Intelligence and Statistics. Ed. by Arindam Banerjee and Kenji Fukumizu. Vol. 130. Proceedings of Machine Learning Research. PMLR, 2021, pp. 1090-1098. URL.: htR.:/htp::// ha00langg--hegrmaannii htmm
.mlr press/
[18] Finn Lindgren, Hävard Rue, and Johan Lindström. "An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach". In: J. R. Stat. Soc. Ser. B Stat. Methodol. 73.4 (201 12) With discussion and a reply by the authors, pp. 423-498. DOI: [0. 1 1 11444777866 2011 0077000777 .x
[19] Per-Gunnar Martinsson and Joel A. Tropp. "Randomized numerical linear algebra: Foun- dations and algorithms". In: Acta Numerica 29 (2020), 403-572 DOI: 10 1017 /
50962492920000021
[20] Amanda F. Mejia, Yu Yue, David Bolin, Finn Lindgren, and Martin A. Lindquist. "A Bayesian general linear modeling approach to cortical surface IMRI data analysis". In: J. Amer Statist, Assoc. 115.530 (2020), pp. 501-520 ISSN: 0162 1459. DOr: 10. 1080 0 0002145921499.. 2019 1611582 URL: httt:::tt://// / org/ 10. 1080/01621459 2019 1611582
[21] Paula Moraga, Susamna M. Cramb, Kerrie L. Mengersen, and Marcello Pagano. "A geostatisti- cal model for combined analysis of point-level and area-level data using INLA and SPDE". In: Spar. Star. 21 part A (2017), pp. 27-41. DOE: 10. 1016/ spasta. 2017 04. 006 URI.:
https://dol: orgoorg/o//oo 10. 1016/ www.sspttttaaa...aaa...pa 20. 006
[22]
Cameron Musco and Christopher Musco. "Randomized Block Krylov Methods for Stronger and Faster Approximate Singular Value Decomposition"" In: Advances in Neural Information Processing Systems. Ed. by C. Cortes, N. Lawrence. D. Lee, M. Sugiyama, and R. Garnett.
Vol. 28. Curran Associates, Inc., 2015. URL: htwphhttts://rcceeddiiee nwwrips. neurips.cc/pppper/ Paf
[23] George Papandreou and Alan Yuille. "Gaussian sampling by local perturbations". In: Advances in Neural Information Processing Systems. Ed. by J. Lafferty, C. Williams, Shawe-Taylor, R. Zemel, and A. Culotta. Vol. 23. Curran Associates. Inc., 2010. URL.: https / proceedings neurips cc paper / 2010 / file / pdf. Trujilo-Barreto, and Karl J Friston. "Bayesian fMRI time series analysis with spatial priors". In: Neurolmage 24.2 (2005), pp. 350-362, D01: 10. 1016/j
[24] William D Penny, Nelson J Trujillo Barreto,
neuroimage. 2004. 08. 034
[25] R Core Team. R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing. Vienna, Austria, 2020. URL: URL: httpss://tw. R- Rrojerect. org/ [26] Carl Edward Rasmussen and Christopher K. Williams. Gaussian Processes for Machine Learning. Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA, 2006, pp. xviii+248. ISBN: 978-0-262-18253-4
[27] Håvard Rue. "Fast sampling of Gaussian Markov random fields". In: J. R. Stat. Soc. Ser B Stat. Methodol. 63.2 (2001), pp. 325-338. ISSN: 1369-7412. DOL: 10 10..111111667-9998 91711988880088888888888888888888
[28] Hävard Rue and Leonhard Held. Gaussian Markov random fields. Vol. 104, Monographs on Statistics and Applied Probability. Theory and applications. Chapman & Hall/CRC, Boca
Raton, FL, 2005, pp. xii+263. ISBN: 978-1-58488-432-3; 1-58488-432-0, DO1: 10. 1201/ 9780203492024
[29] Håvard Rue, Sara Martino, and Nicolas Chopin. "Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations" In: R. Stat. Soc. Ser.
B Stat. Methodol. 71.2 (2009), pp. 319-392. ISSN: 1369-7412. DOI: 10 1111 / j. 1467
9868. 2008. 0008800000.
[30] Mathieu Salzmann and Raquel Urtasun. "Implicitly Constrained Gaussian Process Regression for Monocular Non- Rig-Rigid Pose Estimation".. In: Advances in Newral Information Processing Systems 23. Ed. by J. D. Lafferty, C. K. I. Williams, J. Shawe-Taylor R. S. Zemel, and A. Culotta, Curran Associates, Inc., 2010, pp. 2065-2073.
12

---

[31] Simo Särkkä. "Linear Operators and Stochastic Partial Differential Equations in Gaussian Process Regression". In: Artificial Neural Networks and Machine Learning ICANN 2011. Ed. by Timo Honkela, Wlodzislaw Duch, Mark Girolami, and Samuel Kaski. Berlin, Heidelberg: Springer Berlin Heidelberg, 2011, 1, pp. 151-158. ISBN: 978-3-642-21738--2
[32] Michael Scheuerer and Martin Schlather. "Covariance models for divergence- and curl-free random vector fields". In: Stoch. Models 28.3 (2012), pp. 433 431-451. ISSN: 1532-6349. DOI: wtw.ddo.....g.gooo org/ 10. 1000000332349..20022 2012
10. N0.1080/15326349 20122233322022.... 2t..2.0222.2.... 699756
URL:
699756
[33] Per Sidén and Fredrik Lindsten. "Deep Gaussian Markov Random Fields". In: Proceedings of the 37th International Conference on Machine Learning. Ed. by Hal Daumé IIII and Aarti Singh. Vol. 119. Proceedings Machine Learning Research. PMLR, 2020, pp. 8916 8926.
[34] Per Sidén, Anders Eklund, David Bolin, and Mattias Villani. FFast Bayesian whole-brain IMRI analysis with spatial 3D priors". In: Neurolmage 146 (2017), pp. 211-225. ISSN: 1053-
8119. DOI: https:tpp... htt.//do..ooo/ org 10. 1016/ neuroimage 2016 11 040 URL: https
//ww.sciencedirect. http://eeccccccccceeeeertrrttttt http:///wcaaaaaaaaaa...///////aa
[35] Daniel Simpson, Ian W. Turner, and Anthony Pettitt. "Sampling from Gaussian Markov random fields conditioned on linear constraints". In: ANZIAM J. 48.(C) (2008), C1041-C1053 (2009). ISSN: 1446-1811.
[36] Arno Solin, Manon Kok, Niklas Wahlström, Thomas Schön, and Simo Särkkä. "Mod- eling and Interpolation of the Ambient Magnetic Field by Gaussian Processes". In: IEEE Transactions on Robotics 34.4 (2018), pp. 1112-1127. DOI: 10. 1109/TRO 2018. 2830326
[37] Niklas Wahlström, Manon Kok, Thomas B. Schön, and Fredrik Gustafsson. "Modeling mag- netic fields using Gaussian processes". In: 2013 IEEE International Conference on Acous- tics, Speech and Signal Processing. 2013, pp. 3522-3526.
DOI: 10.
2013
6638813
13


Reviews:
[{"review_id": "_5SZgUKnAX", "reviewer": "Reviewer_rVma", "summary": "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. ", "questions": "", "limitations": "", "rating": 7, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "This is a neat contribution touching upon computational spatial statistics and Gaussian Process modelling. The main approach is quite simple yet potentially very useful in GP applications, as illustrated at the end of the paper. The paper is written in a quite formal way, which I personally apreciate, but could make the reading not so easy to part of ist potential readership. To take a simple example, while it is of interest to tackle the intrinsic case it might be valuable to motivate more in which circumstances it could be of practical use. Also, there is room for better motivating the relevance of the sparse constraint assumption: when it practice is this assumption met / reasonable to assume vs when is it not; also in examples in could be worth stressing more to what extent the constraints are sparse and how the matrix A could actually be refactored in simpler form. What if the constraints are not exactly sparse (but block diagonal up to very small terms outside of blocks)? \n\nA few additional comments along the pages follow: \n•\tOne speaks of fields but from Equation (1) one works with a random vector (n coordinates). \n•\tI am not fond oft he hat notation as it suggests estimation. Better use tilda for instance? \n•\tIn what sense is Q^-1 unbounded? You mean Q^-1 has 0 eigenvalues, right? \n•\tWhat is k0 in Assumption 1? Is it related to something else or named on the spot? \n•\tShould curly C and U be vectors rather than sets as the order matters? \n•\tUnless the number of authors is very large I would avoid „et al.“ in the references \n\nPost-discussion: thank you to the authors for their responses. I maintain my score. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Ok", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2.5"}, "extra_text": {}}, {"review_id": "L3p94_ogKf0", "reviewer": "Reviewer_VGtS", "summary": "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.\n", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Originality: The efficient GMRF methods introduced in this work perform transformation of basis that separates constrained and unconstrained subspaces using SVD, thereby bypassing sparsity challenges and allowing for effective enforcement of the constraints in the transformed space. Nevertheless, the idea of precision matrix decomposition through SVD, or many other decomposition methods for that matter has been extensively studied. In that sense, this work appears to be a combination of well-known practices with sparse gaussian random field models. In other words, the main originality lies in using matrix decomposition in a linearly constrained GMRF setting, however, the idea of precision matrix decomposition in general is well-studied. Also, related studies involving efficient GMRFs based on SVD/eigenvalue/Cholesky decomposition do not seem to be given credit.\n\nQuality: The proposed GMRF methods are (1) introduced primarily for GMRFs with sparse hard (linear) constraints, (2) extended to cases of both hard and soft constraints, and (3) derived and adapted for constrained Gaussian processes. Each variant is technically sound and well supported by theoretical claims involving reasonable assumptions on the positive semi-definiteness of the precision matrix Q (typically assumed in similar settings) and likelihood derivations. Also, using the change of basis, alternative formulations of the distributions of AX and X|AX = b are derived, suitable for efficient sampling and likelihood calculation.\n\nClarity: The paper is generally well written and organized. The notation is clear, not difficult to follow, and consistent throughout the paper.\n\nSignificance: The problem considered in this paper is of considerable importance as efficient sampling and inference with GMRFs as well as constrained Gaussian processes find application in a broad range of areas. On the other hand, I believe that the narrow (i.e., quite specific) settings of the experimental design (particularly in the data generation process) limit the contribution’s significance.\n\nUPDATE AFTER THE REBUTTAL: In their response, the authors have addressed the majority of the comments that I have raised. The authors have also conducted the additional assessment of the GMRF variants under different levels of sparsity that I suggested, summarized the results in the response and will reflect them in the appendix of the paper. Certain points around (1) the contribution's significance and (2) the relation between the SVD-based constraint basis construction and other decomposition methods remain unresolved; however, the authors claim that further clarifications on these points will be included in the revised version of the paper. Therefore, I decided to increase my initially assigned score from 5 to 6 (Marginally above the acceptance threshold).", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors claim that the rows of A can be split into m sub-matrices (blocks) such that each pair of blocks have no common non-zero columns; and then the SVD sub-matrices can be computed separately. I am wondering if the authors have investigated other approaches for the same purpose (e.g., block decomposition methods such as the block SVD power method)? If so, I would encourage the authors to briefly discuss this in their response.\n\nThe specifics of the covariance-based method are not explained nor is it disclosed which particular algorithm was used. Moreover, no comparison to sparse inverse covariance estimation methods has been performed. In that regard, the proposed GMRF methods are compared solely to other exact methods. In some cases, an approximation might exhibit sufficient performance, while being more efficient.\n\nA more thorough experimental evaluation is needed. The main limitations in this regard are summarized as follows:\nThere is an analysis on how the number of observations or basis functions affect the RMSE and the sampling/prediction time. However, there is no analysis on how any of the proposed GMRF variants performs under different levels of sparsity. Has such an experiment been conducted? If so, I would encourage the authors to elaborate on their observations.\n\nToo many assumptions are made in the synthetic data generation process. I am wondering if the experiments with GMRF with hard constraints have been repeated for other values of $k^2$, $\\phi$, $\\alpha$  (other than the [0.5, 1, 2] setting)? Moreover, in the experiments with the constrained Gaussian processes, $\\alpha$ and $\\sigma_2$ are simply fixed to $0.01$ and $10^{-4}$, respectively, with no explanation.\n\nNo experiments on real-world regression datasets have been conducted, particularly in the  case of the GMRF variant for Gaussian process regression.\n", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "10"}, "extra_text": {}}, {"review_id": "1VHnF46_Y4", "reviewer": "Reviewer_2z31", "summary": "This paper proposes a basis transformation method to speed up inference in Gauss--Markov random fields under certain sparse constraints. \n\n", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "\n- Typo? (Line 74) \n\nI believe the first minus sign should be a plus in the conditional mean expression. \n\n- sum-to-zero constraints a violation of assumptions? (Line 191)\n\nIt would seem like a sum-to-zero constraint would violate the assumptions on the constraint matrix from Sec 3.1?\n\n- y-axis in Fig 1. \n\nIt might be preferable to use a logarithmic scaling of the y-axis for readability.\n\n- Computational time decreasing with number of observations (Line 269)\n\nThis would require further explanation.\n\n   ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "yes", "review_text": "", "extra_scores": {"ethical_concerns": "no concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "1.5"}, "extra_text": {}}, {"review_id": "1QC6VcCiTfQ", "reviewer": "Reviewer_Ws93", "summary": "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. \n\nThe 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.\n\n-----------------------\nThe 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.", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "**Writing**\n\nI found the paper to be globally well written and the exposition easy to follow. One shortcoming in my view is that the method could be better motivated (I had to look at the references to convince myself that the problem tackled by the paper had a practical significance). In its current form the paper only specifies that \"hard constraints have a cost that scales cubically with the number of constraints\" but since this is the core issue that is addressed in this work a little hand holding would be appreciated. In my opinion providing more explanations about why the issue arises and a better positioning of the paper with respect to the prior work (ie how this is currently dealt with) would significantly improve the paper quality.\n\n**Methodology**\n\nThe proposed methodology is sound but there is a chance that the assumptions that are added on the sparsity patterns of $A$ when introducing the method could be used in the first place to cut down the complexity of the sampling and of the likelihood evaluation. I believe this could for example be done with the sparse inverse subset method as described in [1].\n\n**Significance**\n\nThe feeling that this contribution is of interest to the community would be much stronger if the authors were to give a better description of how the problem they address is currently tackled (or unsolved) and what are the shortcomings. At the moment it seems that they mainly compare a naive baseline which results in impressive results but does not seem to be a fair comparison. For example, my understanding is that using observations as nodes for the triangulation would be a perfectly valid way to tackle the problem in Experiment 6.1.\n\n**Experiments**\n\nThree remarks:\n - Since the proposed method loses its interest when $A$ cannot be split into \"independent rows\", it would be interesting to increase the number of constraints up to the point where this becomes the case.\n - Could an explanation of why the evaluation time decreases with the number of constraints be included in the result analysis?\n - Is the computation of $T$ included in the timings that are reported?\n\n\n[1] Durrande, Nicolas, et al. \"Banded matrix operators for Gaussian Markov models in the automatic differentiation era.\" The 22nd International Conference on Artificial Intelligence and Statistics. PMLR, 2019.\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "No specific comments", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}]

Author rebuttal:

Meta-review: Following reviewers discussion, I recommend accepting the paper under the condition that the authors will explain that its main contributions are applications for the projective clustering problem  in probability/statistics, and cite the relevant existing solutions from different fields (coresets in computational geometry, subspace clustering in DB, dictionary learning in signal processing). It should also be clarified that the suggested solution is given as a simple and possibly inefficient solver of projective clustering for the special case where no noise exists (the points lie on the subspace).
In this case, the paper can serve as an interesting bridge between the CS, ML and statistics community.
    

Decision: Accept (Poster)