Obtain Precision/Covariance Matrix for Observed Locations

[S B Lee]

Hello, 

I am currently working on deriving the covariance matrix for the spatial field evaluated at observed locations S_obs based on a fitted SPDE model defined over a set of mesh nodes S_mesh. Reproducible code is provided below. Given the projector matrix, the observed values can be expressed as:

 Z_obs = A%*% Z_mesh  

where Z_mesh denotes the latent spatial field at the mesh nodes. . Consequently, we have the distribution

  Z_obs ~ N(A%*% mu_mesh, A%*%C_mesh%*%A^T) 

where mu_mesh= E[Z_mesh] and C_mesh = Cov(Z_mesh). 

Using the Matérn covariance model implied by the SPDE formulation, the latent precision matrix Q_mesh=C_mesh^-1 can be obtained through the function inla.spde2.precision(). However, I am encountering challenges in computing the precision matrix or determinant of the projected covariance matrix A%*%C_mesh%*%A^T primarily due to rank deficiency stemming from the presence of zero-valued eigenvalues.

Would you happen to have any recommendations or best practices for dealing with this issue? 

Thank you for your time!

Code:

library(INLA)
n=500
locations<-cbind(runif(n) , runif(n))
d=2

# Covariance Function Parameters using Equations (1), (2), (4), and (5) from Lindgren and Rue (2015)
# https://www.jstatsoft.org/article/download/v063i19/838
rho<-0.2
sigma2<-1
nu<-0.5
alpha<-nu + d/2
log_kappa <-log(8*nu)/2 - log(rho)
log_tau <-  0.5* log(gamma(nu)/(gamma(alpha)*(4*pi)^(d/2)))- log(sqrt(sigma2))- nu *log_kappa
tau <- exp(log_tau)
kappa <- exp(log_kappa)

(8*nu)^(1/2)/kappa ; rho
gamma(nu)/(gamma(alpha)*(4*pi)^(d/2)*kappa^(2*nu)*tau^2) ; sigma2

# Generate Mesh
mesh<-inla.mesh.2d(locations,
                   max.edge=c(.1), cutoff=0.01)

# Define the SPDE model (Matérn)
spde <- inla.spde2.matern(mesh = mesh, alpha = nu+1)
AMat <- inla.spde.make.A(mesh, loc=locations)  

# Obtain the precision matrix (Q) and covariance matrix (inv_Q) using specified parameter values
theta <- c(log_tau, log_kappa)  # log(tau), log(kappa)
Q <- inla.spde2.precision(spde, theta = theta)
inv_Q<-solve(Q)

# Compute observation-level covariance matrix and precision matrix
covMat <- AMat %*% inv_Q %*% t(AMat)
inv_covMat<-chol2inv(chol(covMat))
eigen(covMat)$values

[Finn Lindgren]

Hi,

there is a tool for this in `fmesher`;

When Q is the precision matrix and A1 and A2 are fm_basis() or inla.spde.make.A() output matrice, 

  fm_covariance(Q, A1, A2)

computes the covariance between A1 * x and A2 * x.

If A2=A1=A, you can use the shorter form

  fm_covariance(Q, A)

The difference between that function and your code is that it doesn't evaluate the full inverse matrix solve(Q), but instead solves a linear system:

A1 %*% Matrix::solve(Q, Matrix::t(A2))

/Finn

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/b4ef2628-f82e-4837-8c3d-155b8d5f9d39n%40googlegroups.com.

--

Finn Lindgren
email: finn.l...@gmail.com

[Ben Lee]

Thank you very much for the quick response as well as pointing me towards the fm_covariance() function. 

I am sorry, but I made a blunder in the original email. I am actually trying to compute the precision matrix, not the covariance,  between  A1 * x and A2 * x as well as the corresponding determinant. So I would need the inverse of A1 %*% Matrix::solve(Q, Matrix::t(A2))A1 %*% Matrix::solve(Q, Matrix::t(A2)). 

If the columns of A are not linearly independent,  then AQ^(-1)A^t would not be invertible as this example demonstrates. Should I just add a small nugget term to the diagonal to bypass this issue?

[Finn Lindgren]

Hmmm,

I'd ask what you think you need that precision matrix for; as you note, it may be singular (e.g. if several locations lie within the same triangle) and also shouldn't be expected to have any sparsity structure; except for very specific special cases, the marginal observation precision matrix is fully dense.

The precision construction of the basis function representation of the spde models makes sense for the coefficients of the weights for the basis functions, but _not_ for the (joint marginal) distribution of the observations, whereas for those, the covariance matrix is well-defined (albeit potentially singular).

What actual question are you trying to answer? 

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/e03146b3-3566-4de0-b9fe-eb4782188783n%40googlegroups.com.

[Ben Lee]

Thank you for the clarification and the quick response! 

I’m working with a spatial dataset and wanted to evaluate the multivariate normal distribution using a covariance or precision matrix derived from a Matérn covariance function. I was curious whether I could use the SPDE approach to efficiently obtain the precision matrix (at the observed locations) and compute its determinant or eigenvalues. The main goal is to see if this would be more computationally efficient than fitting the full model without any approximations.

Thank you for your time. 

You received this message because you are subscribed to a topic in the Google Groups "R-inla discussion group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/r-inla-discussion-group/4goZIg8ots4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAHC4gza09EWYTggDE%2BLoVmqdNbc4XarvJ%3Dnu_NJgqxb5SRgd_g%40mail.gmail.com.

[Ben Lee]

My apologies for the double message, but I had one other thought.

If I we were to build a mesh such that the mesh consists of only the observed locations, we could use the  Q precision matrix directly in evaluating the multivariate normal distribution. Then we would not have to worry about the singularity issue as previously mentioned. That being said, is it possible to generate a mesh where the mesh nodes are the locations themselves? If there are additional mesh vertices, in addition to the observation locations, I could just remove the corresponding rows in the resulting Q matrix. 

[Finn Lindgren]

Hi Ben,

I see. The problem is that the common intuition for how to compute things that is based on _starting_ with the covariance matrix isn't valid when what we have more direct access to is a sparse precision matrix. Nearly all expressions can be reformulated (usually using the "matrix inversion lemma") into expressions that do not actually involve the marginal observation covariance matrix. In short, when we work with sparse precision matrices, we rarely need to compute covariance matrices, and most quantities can be evaluated using a different set of standard computational expressions than those based on covariances.

For example, to evaluate the density for the observations, conditionally on the covariance parameters, one can rewrite the expression into functions of the field precision matrix Q_x, the conditional observation precision Q_e (typically a diagonal matrix) and the data-conditional precision for the field, Q_{x|y}, which is equal to Q_x +A' Q_e A, when A is the observation matrix (from fm_basis() or inla.spde.make.A()). The details can be seen in Section 3.1 of my notes on how to compute with Gaussian Markov random fields (a condensation of the relevant methods from the Rue&Held 2005 book on Gaussian Markov Random Fields) available at https://webhomes.maths.ed.ac.uk/~flindgre/cuso2019/gmrf.pdf

One of the few cases where we might explicitly want to have access to elements of the covariance matrix for y, S_y = Q_y^{-1}  is to access the pointwise marginal variances at the mesh nodes.

These can be computed by selected/partial inversion of Q_x, plus the noise variances;

  y_var = diag(inla.qinv(Q_x)) + diag(inla.qinv(Q_e))

Approximate variance evaluation can then be done with A %*% y_var, or (A %*% y_var^0.5)^2.

As long as A' A has non-sparsity that's contained in the sparsity of Q_x, the variances of A %*% x can also be obtained using tricks with the partial inverses; some of that is also written in the gmrf.pdf document.

To specifically answer your last question: It is possible to construct the mesh so it includes specific locations of interest (but this isn't always a numerically beneficial approach), but

you can _not_ obtain the joint precision matrix for _just_ those points by removing the other rows and columns from the matrix; doing that would only give you the _conditional_ precision matrix, i.e. the inverse of the _conditional_ covariance matrix, given all the other points. This can be seen either from the matrix inversion lemma or from block matrix inversion calculations;

If the covariance is

  S = [S_{AA}, S_{AB}; S_{BA}, S_{BB}]

and the precision is

  Q = [Q_{AA}, Q_{AB}; Q_{BA}, Q_{BB}],

then

  Q_{BB}^{-1} = S_{BB} - S_{BA} S_{AA}^{-1} S_{AB},

which is the conditional covariance of block B, given block A.

Conversely, the _marginal_ precision matrix for block B is

  S_{BB}^{-1} = Q_{BB} - Q_{BA} Q_{AA}^{-1} Q_{AB},

but there are few cases where calculating this is useful.

Finn

To view this discussion, visit https://groups.google.com/d/msgid/r-inla-discussion-group/f1ca392e-0b19-4bfc-844c-547f1fdd4f79n%40googlegroups.com.

[Ben Lee]

Hi Finn, 

My apologies for the late response here. Thank you very much for the detailed and thoughtful comments and suggestions!

This is very promising for my project, and I will look into the details provided in your notes. Since we are able to compute Q_x, Q_{x|y} and Q_e easily, this would help with getting through the computational bottleneck. I also appreciate your feedback on my suggestion of using the locations as mesh nodes. 

Again, thank you very much!

best, 

Ben