canongbp.m

Gaussian belief propagation code in matlab.

% [B, BC] = CANONGBP(G,LXX,LXY,Y,VERBOSE,SCHED) runs belief propagation
% on an undirected graphical model with Gaussian factors following a
% canonical parameterization. See CANONMARGS for explanation of the first
% five inputs and the two outputs. 
%
% VERBOSE and SCHED are optional parameters. SCHED is a user-defined message
% passing schedule. It is only useful to specify a schedule when the graph
% is not a tree. The default message passing schedule will converge to the
% correct marginals as long as the graph does not contain any cycles. See
% MPP for more information on how to specify a message passing schedule.

function [b, bc] = canongbp (g, Lxx, Lxy, y, verbose, sched)
  
  % Get the number of nodes in the graph.
  n = length(Lxy);

  % Get the number of edges in the graph.
  e = length(Lxx);
  
  % Get the number of features.
  F = size(y,1);
  
  % Check to see if the user set verbose. If not, set it to the default.
  if nargin < 5
    verbose = 0;
  end
  
  % If a message passing schedule was not specified, create one.
  if nargin < 6
    sched = mpp(g.A > 0);
  end
  
  % Map from a pair of edges to a message index.
  T = g.A + (tril(g.A) > 0)*e;

  % Initialize the messages. Note that if the schedule is correct, we
  % will never use the initialized messages.
  M = canonpar(F,e*2);
  
  % Compute the node potentials (more precisely, the mean and
  % precision). The precision of the node potentials is identity, so we
  % don't need to compute it.
  nu = zeros(F,n);
  for i = 1:n
    nu(:,i) = -Lxy{i}*y(:,i);
  end

  I = eye(F);

  % Repeat for each message.
  for m = 1:length(sched)
    
    % The message is from i to j.
    i = sched(m).i;
    j = sched(m).j;

    if verbose
      disp(sprintf('%d -> %d', i, j));
    end
    
    % Get i's neighbours that are not j. (We don't want to integrate over
    % the message from j to i.)
    g.A(j,i) = 0;
    ks       = find(g.A(:,i)' > 0);
    g.A(j,i) = 1;
    
    % Multiply the node potential and the messages (k,i).
    nu0 = nu(:,i);
    L0  = I;
    for k = ks
      t   = T(k,i);
      nu0 = nu0 + M(t).nu;
      L0  = L0  + M(t).L;
    end

    % Get the index of the message (i,j).
    t = T(i,j);

    % Get the index of the edge (i,j).
    u = g.A(i,j);
    
    % This is the message update.
    M(t).L  = -Lxx{u}*inv(L0)*Lxx{u};
    M(t).nu = -Lxx{u}*inv(L0)*nu0;
  end
  clear m
  
  % Now that the messages have converged (assuming the graph is a tree),
  % we compute the marginal beliefs.
  b = momentpar(F,n);
  for i = 1:n
    
    % Start with the node potential.
    ni = nu(:,i);
    Li = I;

    % Get the neighbours of i.
    ks = find(g.A(:,i)' > 0);

    % Multiply by the messages coming into node i.
    for k = ks
      t  = T(k,i);
      ni = ni + M(t).nu;
      Li = Li + M(t).L;
    end
    
    % Convert the canonical parameters to the moment parameters.
    [b(i).mu b(i).S] = canontomoment(ni,Li);
  end

  % Compute the edge marginals.
  bc = momentpar(2*F,e);
  for u = 1:e
    i = g.E(u,1);
    j = g.E(u,2);
    
    % Get all the neighbours of i except j.
    g.A(j,i) = 0;
    ks       = find(g.A(:,i)' > 0);
    g.A(j,i) = 1;
    
    % Multiply the node potential and the messages coming into node i.
    ni = nu(:,i);
    Li = I;
    for k = ks
      t  = T(k,i);
      ni = ni + M(t).nu;
      Li = Li + M(t).L;
    end
  
    % Get all the neighbours of j except i.
    g.A(i,j) = 0;
    ks       = find(g.A(:,j)' > 0);
    g.A(i,j) = 1;
    
    % Multiply the node potential and the messages coming into node i.
    nj = nu(:,j);
    Lj = I;
    for k = ks
      t  = T(k,j);
      nj = nj + M(t).nu;
      Lj = Lj + M(t).L;
    end

    Lu  = [ Li     Lxx{u};
            Lxx{u}   Lj   ];
    [bc(u).mu bc(u).S] = canontomoment([ni; nj],Lu);
  end
end