wolffd@0
|
1 function engine = frontier_inf_engine(bnet)
|
wolffd@0
|
2 % FRONTIER_INF_ENGINE Inference engine for DBNs which which uses the frontier algorithm.
|
wolffd@0
|
3 % engine = frontier_inf_engine(bnet)
|
wolffd@0
|
4 %
|
wolffd@0
|
5 % The frontier algorithm extends the forwards-backwards algorithm to DBNs in the obvious way,
|
wolffd@0
|
6 % maintaining a joint distribution (frontier) over all the nodes in a time slice.
|
wolffd@0
|
7 % When all the hidden nodes in the DBN are persistent (have children in the next time slice),
|
wolffd@0
|
8 % its theoretical running time is often similar to that of the junction tree algorithm,
|
wolffd@0
|
9 % although in practice, this algorithm seems to very slow (at least in matlab).
|
wolffd@0
|
10 % However, it is extremely simple to describe and implement.
|
wolffd@0
|
11 %
|
wolffd@0
|
12 % Suppose there are n binary nodes per slice, so the frontier takes O(2^n) space.
|
wolffd@0
|
13 % Each time step takes between O(n 2^{n+1}) and O(n 2^{2n}) operations, depending on the graph structure.
|
wolffd@0
|
14 % The lower bound is achieved by a set of n independent chains, as in a factorial HMM.
|
wolffd@0
|
15 % The upper bound is achieved by a set of n fully interconnected chains, as in an HMM.
|
wolffd@0
|
16 %
|
wolffd@0
|
17 % The factor of n arises because we need to multiply in each CPD from slice t+1.
|
wolffd@0
|
18 % The second factor depends on the size of the frontier to which we add the new node.
|
wolffd@0
|
19 % In an FHMM, once we have added X(i,t+1), we can marginalize out X(i,t) from the frontier, since
|
wolffd@0
|
20 % no other nodes depend on it; hence the frontier never contains more than n+1 nodes.
|
wolffd@0
|
21 % In a fully coupled HMM, we must leave X(i,t) in the frontier until all X(j,t+1) have been
|
wolffd@0
|
22 % added; hence the frontier will contain 2*n nodes at its peak.
|
wolffd@0
|
23 %
|
wolffd@0
|
24 % For details, see
|
wolffd@0
|
25 % "The Factored Frontier Algorithm for Approximate Inference in DBNs",
|
wolffd@0
|
26 % Kevin Murphy and Yair Weiss, UAI 01.
|
wolffd@0
|
27
|
wolffd@0
|
28 ns = bnet.node_sizes_slice;
|
wolffd@0
|
29 onodes = bnet.observed;
|
wolffd@0
|
30 ns(onodes) = 1;
|
wolffd@0
|
31 ss = length(bnet.intra);
|
wolffd@0
|
32
|
wolffd@0
|
33 [engine.ops, engine.fdom] = best_first_frontier_seq(ns, bnet.dag);
|
wolffd@0
|
34 engine.ops1 = 1:ss;
|
wolffd@0
|
35
|
wolffd@0
|
36 engine.fwdback = [];
|
wolffd@0
|
37 engine.fwd_frontier = [];
|
wolffd@0
|
38 engine.back_frontier = [];
|
wolffd@0
|
39
|
wolffd@0
|
40 engine.fdom1 = cell(1,ss);
|
wolffd@0
|
41 for s=1:ss
|
wolffd@0
|
42 engine.fdom1{s} = 1:s;
|
wolffd@0
|
43 end
|
wolffd@0
|
44
|
wolffd@0
|
45 engine = class(engine, 'frontier_inf_engine', inf_engine(bnet));
|
wolffd@0
|
46
|
wolffd@0
|
47
|
wolffd@0
|
48 %%%%%%%%%
|
wolffd@0
|
49
|
wolffd@0
|
50 function [ops, frontier_set] = best_first_frontier_seq(ns, dag)
|
wolffd@0
|
51 % BEST_FIRST_FRONTIER_SEQ Do a greedy search for the sequence of additions/removals to the frontier.
|
wolffd@0
|
52 % [ops, frontier_set] = best_first_frontier_seq(ns, dag)
|
wolffd@0
|
53 %
|
wolffd@0
|
54 % We maintain 3 sets: the frontier (F), the right set (R), and the left set (L).
|
wolffd@0
|
55 % The invariant is that the nodes in R are d-separated from L given F.
|
wolffd@0
|
56 % We start with slice 1 in F and slice 2 in R.
|
wolffd@0
|
57 % The goal is to move slice 1 from F to L, and slice 2 from R to F, so as to minimize the size
|
wolffd@0
|
58 % of the frontier at each step, where the size(F) = product of the node-sizes of nodes in F.
|
wolffd@0
|
59 % A node may be removed (from F to L) if it has no children in R.
|
wolffd@0
|
60 % A node may be added (from R to F) if its parents are in F.
|
wolffd@0
|
61 %
|
wolffd@0
|
62 % ns(i) = num. discrete values node i can take on (i=1..ss, where ss = slice size)
|
wolffd@0
|
63 % dag is the (2*ss) x (2*ss) adjacency matrix for the 2-slice DBN.
|
wolffd@0
|
64
|
wolffd@0
|
65 % Example:
|
wolffd@0
|
66 %
|
wolffd@0
|
67 % 4 9
|
wolffd@0
|
68 % ^ ^
|
wolffd@0
|
69 % | |
|
wolffd@0
|
70 % 2 -> 7
|
wolffd@0
|
71 % ^ ^
|
wolffd@0
|
72 % | |
|
wolffd@0
|
73 % 1 -> 6
|
wolffd@0
|
74 % | |
|
wolffd@0
|
75 % v v
|
wolffd@0
|
76 % 3 -> 8
|
wolffd@0
|
77 % | |
|
wolffd@0
|
78 % v V
|
wolffd@0
|
79 % 5 10
|
wolffd@0
|
80 %
|
wolffd@0
|
81 % ops = -4, -5, 6, -1, 7, -2, 8, -3, 9, 10
|
wolffd@0
|
82
|
wolffd@0
|
83 ss = length(ns);
|
wolffd@0
|
84 ns = [ns(:)' ns(:)'];
|
wolffd@0
|
85 ops = zeros(1,ss);
|
wolffd@0
|
86 L = []; F = 1:ss; R = (1:ss)+ss;
|
wolffd@0
|
87 frontier_set = cell(1,2*ss);
|
wolffd@0
|
88 for s=1:2*ss
|
wolffd@0
|
89 remcost = inf*ones(1,2*ss);
|
wolffd@0
|
90 %disp(['L: ' num2str(L) ', F: ' num2str(F) ', R: ' num2str(R)]);
|
wolffd@0
|
91 maybe_removable = myintersect(F, 1:ss);
|
wolffd@0
|
92 for n=maybe_removable(:)'
|
wolffd@0
|
93 cs = children(dag, n);
|
wolffd@0
|
94 if isempty(myintersect(cs, R))
|
wolffd@0
|
95 remcost(n) = prod(ns(mysetdiff(F, n)));
|
wolffd@0
|
96 end
|
wolffd@0
|
97 end
|
wolffd@0
|
98 %remcost
|
wolffd@0
|
99 if any(remcost < inf)
|
wolffd@0
|
100 n = argmin(remcost);
|
wolffd@0
|
101 ops(s) = -n;
|
wolffd@0
|
102 L = myunion(L, n);
|
wolffd@0
|
103 F = mysetdiff(F, n);
|
wolffd@0
|
104 else
|
wolffd@0
|
105 addcost = inf*ones(1,2*ss);
|
wolffd@0
|
106 for n=R(:)'
|
wolffd@0
|
107 ps = parents(dag, n);
|
wolffd@0
|
108 if mysubset(ps, F)
|
wolffd@0
|
109 addcost(n) = prod(ns(myunion(F, [ps n])));
|
wolffd@0
|
110 end
|
wolffd@0
|
111 end
|
wolffd@0
|
112 %addcost
|
wolffd@0
|
113 assert(any(addcost < inf));
|
wolffd@0
|
114 n = argmin(addcost);
|
wolffd@0
|
115 ops(s) = n;
|
wolffd@0
|
116 R = mysetdiff(R, n);
|
wolffd@0
|
117 F = myunion(F, n);
|
wolffd@0
|
118 end
|
wolffd@0
|
119 %fprintf('op at step %d = %d\n\n', s, ops(s));
|
wolffd@0
|
120 frontier_set{s} = F;
|
wolffd@0
|
121 end
|