# Creating a precision matrix for Gaussian markov random field

I am currently trying to create a precision matrix for a Gaussian markov random field. Lets say I have random variables in a spatial grid of 6x6. Then I will have a precision matrix of 36x36.

Now suppose that I have a neighbor hood of 3x3, then my precision matrix will be

Q= nnbs[1] -1      0        0    0     0    -1.......0
-1      nnbs[2] -1       0    0     0     0 ......0
0       -1      nnbs[3]  -1   0     0     0 ......0
...................................................
...................................................


and so on. Can anyone suggest me how can I code this precision matrix. I mean if I change the window size/neighborhood size to 5x5, then I will have a new precision matrix. How can I code this?

rows=20;
columns=20;

%Random initialization
data=zeros(1000,3);
index=1;
value=-1;

%3x3 neighborhood
%For each element the neighbors are accessible within 1 hop so neighbors=1
neighbors=1;

for i=1:rows
for j=1:columns

for k=1:neighbors
%same row right
if j+k <= columns
data(index,1) = (i-1)*columns+j;
data(index,2) = ((i-1)*columns) + (j+k);
data(index,3) = value;
index=index+1;
end

%same row left
if j-k >= 1;
data(index,1) = (i-1)*columns+j;
data(index,2) = ((i-1)*columns) + (j-k);
data(index,3) = value;
index=index+1;
end
end

%row below -> bottom left right
for k=i+1:i+neighbors
if k <= rows
%bottom
data(index,1) = (i-1)*columns+j;
data(index,2) = (k-1)*columns + j;
data(index,3) = value;
index=index+1;

for l=1:neighbors
%right
if j+l <= columns
data(index,1) = (i-1)*columns+j;
data(index,2) = ((k-1)*columns) + (j+1);
data(index,3) = value;
index=index+1;
end

%left
if j-l >= 1;
data(index,1) = (i-1)*columns+j;
data(index,2) = ((k-1)*columns)+(j-1);
data(index,3) = value;
index=index+1;
end
end

end

end

%row above top left right
for k=i-1:i-neighbors
if k >= 1
%top
data(index,1) = (i-1)*columns+j;
data(index,2) = ((k-1)*columns) +j;
data(index,3) = value;
index=index+1;

for l=1:neighbors
%right
if j+l <= columns
data(index,1) = (i-1)*columns+j;
data(index,2) = ((k-1)*columns) + (j+1);
data(index,3) = value;
index=index+1;
end

%left
if j-k >= 1;
data(index,1) = (i-1)*columns+j;
data(index,2) = ((k-1)*columns) + (j-1);
data(index,3) = value;
index=index+1;
end
end
end
end
end
end

%Get the values for the diagonal elements(which is equal to the number of
%neighbors or absolute sum of the nondiagonal elements of the corresponding
%row)

diagonal_values = zeros(rows*columns,3);
for i=1:rows*columns
pointer=find(data(:,1) == i);
diag_value=abs(sum(data(pointer,3)));
diagonal_values(i,1) = i;
diagonal_values(i,2) = i;
diagonal_values(i,3) = diag_value;
end

data(index:index+rows*columns-1,:)=diagonal_values(:,:);

Q = sparse(data(:,1), data(:,2), data(:,3), rows*columns, rows*columns);


Currently I have written this code. But I don't think it is that efficient.

-
Hm, let me try to understand. What is nnbs[.]. Why do you have the '-1?'s? –  Hauke Strasdat Jul 3 '12 at 10:45
In general, the following is true. If you have $n$ noded/variables in your MRF, the precision matrix is $n\times n$. So far, so good. Now, we know that the entry $(i,j)$ is $0$ if node $i$ is NOT connected to node $j$. This will give you the sparseness pattern of your precision matrix. –  Hauke Strasdat Jul 3 '12 at 10:48
Are you trying to generate a random precision matrix given the window, neighborhood values? If you are assuming certain graphical structures, you can code easily.. –  tatterdemalion Jul 9 '12 at 15:23
No it's not random precision matrix. It depends on the graph itself. If I am considering a window of 3x3, then for each center node I will have a window of 3x3 or 8 neighbors, which are connected. So these neighbors will have a non zero component in the precision matrix and so on. The above code works for variable window sizes like 3x3,5x5 and so on. But I don't think the code is efficient. There should have been a simpler way to do it –  user34790 Jul 9 '12 at 23:14