为高斯马尔可夫随机场创建精度矩阵

我目前正在尝试为高斯马尔可夫随机场创建精度矩阵。可以说我在6x6的空间网格中有随机变量。然后我将有一个36x36的精度矩阵。为高斯马尔可夫随机场创建精度矩阵

现在假设我有3×3的邻居罩，然后我的精度矩阵将

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 
    ................................................... 
    ...................................................

等。任何人都可以建议我如何编码这个精度矩阵。我的意思是，如果我将窗口大小/邻域大小更改为5x5，那么我将拥有一个新的精度矩阵。我如何编码？其中nnbs是该元素的邻居数

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);

我试过类似的东西，但我不认为这是最有效的方法。我认为应该有更好的办法。

来源

2012-06-30 user34790

我还在等。没有答案？ – user34790

有点晚，但它可以是别人有用：你的精度矩阵是对称的Toeplitz矩阵的一种克罗内克积的线性组合：每个邻居类型对应2托普利兹矩阵的Kronecker乘积。 More info about toeplitz Matrix

实施例：要精密矩阵只与水平链路用于每个像素

写入I_n尺寸n的单位矩阵和H_{n,p}用0填充尺寸[n n]的对称的Toeplitz矩阵到处节选于第p个对角线

H_ {4,2} =

在Matlab中：

H_nonSym_n_p =托普利兹（零（N，1），[[零（1，P-1）1]的零（1，N-P）]）;

H_sym_n_p = H_nonSym + H_nonSym';

然后，如果你有一个[n m]形象，如果你想代码的每一像素的横向的邻居，你可以通过直积表达它（希望乳液状的代码将工作）：Q = - 1-N \ otimes \ H_ {n，2}。最后，为了让您的nnbs：像Q = Q - diag(sum(Q,2)) ...

现在，如果你希望其他环节，比如2水平和2个垂直链接：Q = - 1-N \ otimes \ H_ {N，2} - I_n \ otimes \ H_ {n，3} - \ H_ {n，2} \ otimes I_ {n} - \ H_ {n，3} \ otimes I_ {n}。 Q = Q - diag(sum(Q,2))

请注意，对角线邻居有点难以产生，但它仍然表示toeplitz矩阵的克罗内克积（可能这次不是对称）。

来源

2014-10-27 14:46:56

为高斯马尔可夫随机场创建精度矩阵

回答

相关问题