canonmargs.m

Gaussian belief propagation code in matlab.

% [B, BC] = CANONMARGS(G,LXX,LXY,Y) computes the marginal probabilities on
% individual nodes and pairs of adjacent nodes given a factored Gaussian
% distribution with canonical parameterization. G is structure with fields
% G.A and G.E. G.A is an N x N adjacency matrix, where N in the number of
% nodes in the graphical model, such that G.A(i,j) > 0 if and only if nodes
% i and j are adjacent. The value A(i,j) = u specifies the index of the
% factor (see the parameter LXX below for more information). Since this is
% an undirected graphical model, A(i,j) = A(j,i).

% LXX is a cell array of length M, where M is the number of edges in the
% graph. LXX{u} is the precision matrix (also called the "inverse covariance
% matrix") between nodes i and j, where 
%
%        i = G.E(u,1)
%        j = G.E(u,2)
%
% LXX{u} is an F x F matrix, where F is the dimension of the random
% variables in the model. We assume that the precision matrices are
% symmetric, so there is no change if we switch i and j. The precision
% matrix of the joint of Xi and Xj is
%
%                    | I    Lxx |
%                    | Lxx'  I  |
%
% Y is an F x N matrix, where entry Y(:,i) is the ith observation. LXY{i} is
% the precision matrix between hidden state Xi and observation Yi. The
% joint on X and Y is given by
%
%                    | I    Lxy |
%                    | Lxy'  I  |
%
% The output B is an N x 1 cell array, where each entry B{i} is the marginal
% distribution of Xi conditioned on all the observations. BC is an M x 1
% cell array, and each entry B{u} is the marginal distribution on (Xi,Xj)
% conditioned on the observations, where i = G.E(u,1) j = G.E(u,2), as
% before.

function [b, bc] = canonmargs (g, Lxx, Lxy, y)
  
  % Get the number of nodes.
  n = length(Lxy);
  
  % Get the number of edges.
  m = length(Lxx);
  
  % Get the number of dimensions for each hidden state and observation.
  F = size(y,1);
  
  % Build the joint precision (inverse covariance) matrix over all the
  % observations and hidden states. The matrix over the joint p(x,y) is 
  %    | Jxx  Jxy |
  %    | Jxy'  I  |
  nf  = n*F;
  Jxx = speye(nf);
  Jxy = sparse(nf,nf);
  
  % This is correct only because the covariance matrices are symmetric.
  for u = 1:m
    i          = g.E(u,1);
    j          = g.E(u,2);
    xi         = extract(i,F);
    xj         = extract(j,F);
    Jxx(xi,xj) = Lxx{u};
    Jxx(xj,xi) = Lxx{u};
  end
  
  % Add the node potentials to the covariance matrix.
  for i = 1:n
    xi         = extract(i,F);
    Jxy(xi,xi) = Lxy{i};
  end
  
  % Condition on the observations.
  [nu L] = canoncond(zeros(nf,1),Jxx,Jxy,y(:));
  
  % Now that we have the conditional, let's compute the node marginals.
  b = momentpar(F,n);
  for i = 1:n
    xi               = extract(i,F);
    [b(i).mu b(i).S] = marginal(xi,nu,L,n,F);
  end
  
  % Compute the edge marginals.
  bc = momentpar(2*F,m);
  for u = 1:m
    i  = g.E(u,1);
    j  = g.E(u,2);
    xi = [extract(i,F) extract(j,F)];
    [bc(u).mu bc(u).S] = marginal(xi,nu,L,n,F);
  end

function x = extract (i, F)
  x = F*(i-1) + [1:F];
  
function xj = remainder (xi, F, n)
  A     = ones(1,n*F);
  A(xi) = 0;
  xj    = find(A);
  
function [mi, Si] = marginal (xi, nu, L, n, F)
  
  xj = remainder(xi,n,F);
  if length(xj)
    
    Lii = L(xi,xi);
    Ljj = L(xj,xj);
    Lij = L(xi,xj);
    
    % Compute the marginal.
    [ni Li] = canonmarg(nu(xi),nu(xj),Lii,Ljj,Lij);
  else
    
    % We don't need to do any marginalization.
    ni = nu;
    Li = L;
  end
  
  % Convert the canonical parameters to the moment parameters.
  [mi Si] = canontomoment(ni,full(Li));