Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

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.

share|improve this question
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.