I am trying to understand an experiment from this paper, specifically Section 5.2.
In the paper, they propose a new algorithm for computing the log-determinant of sparse matrices, and in section 5 they test it on a dataset that they generate.
They want to test it on a synthetic dataset, so they create a sparse matrix of size 5000x5000 whose precision matrix (the inverse of the covariance matrix) is parameterized by $\rho = -0.22$. According to the paper, each node has 4 neighbors with partial correlation of $\rho$. Then, they use a Gibbs sampler to take one sample from the multivariate gaussian distribution that is described by the matrix J. On this sample, I compute the log-likelihood as: $$\log(x|\rho) = \log\det J(\rho) - x^TJ(\rho)x + G$$. and I plot the values as in Figure 2.
If my understanding correct, they do evaluate the log-likelihood given only one sample? I understand the plot in figure 2 is the following formula above, which is calculated only for one sample? I usually compute the log-likelihood on a dataset, not just on a single sample.
My question is the following: what exactly is the meaning $\rho$, and how do I create $J(\rho)$ and sample from it? (i.e. with a python package? otherwise, what is the algorithm?)?
I think the underlying assumption is that the $\log\det J(\rho)$ for two different samples of $J(\rho)$ is the same, why?
I actually went to look to the much-cited referenced book , which is a very good book about GMRF, but I have not found any clear link between a single parameter $\rho$ and the matrix they generate. The parameterization of GMRFs is described in Section 2.7, page 87. There, a single parameter is never used, and the parameter space is actually described by a 2 dimensional vector $\Theta$:
$$ \pi(x|\Theta) \propto exp(\frac{-\theta_1}{2}\sum_{i\approx j} (x_i - x_j)^2 - \frac{\theta_2}{2}\sum_i x_i^2 )$$ But probably they are referring to a different matrix.
Update Actually, I think the precision matrix $J(\rho)$ which describe the interaction between 4 neighbors is just a band matrix i.e. a matrix with multiple diagonals. In this case (I imagine) with 2 upper and 2 lower diagonals, all filled with $-0.22$, and just $1$ on the main diagonal.
Still, how can I sample from the distribution described by the precision matrix? Should I invert it and obtain the covariance matrix of the data and then sample from it? If yes, below is the code to do it. It might be useful for someone to see the code we can use to sample from this GMRF, assuming $\vec(0)$ mean and a matrix dimension of 30.
import numpy as np
def precision(k, numero):
return np.diag(np.repeat(1, k)) + np.diag(np.repeat(numero, k-1), -1) + np.diag(np.repeat(numero, k-2), -2) + np.diag(np.repeat(numero, k-1), 1) + np.diag(np.repeat(numero, k-2), 2)
J = precision(30, -0.22)
Sigma = np.linalg.inv(J)
np.random.multivariate_normal(np.repeat(0, 30), Sigma)
