pgm5

这部分讨论 inference 里面基本的问题，即计算 \Pr(Y \mid E = e) 这类 query，这一般可以认为等价于计算 \Pr (Y, E = e)，因为我们只需要重新 normalize 一下关于 Y 的分布就得到了需要的值，特别是像 MAP 这类 query（一般此时 Y 是 E 的补集，可以理解成为取 \max_y \Pr(Y = y \mid E = e) 只需要将这里一些 sum 换成 max 即可）也可以用类似的策略来处理。这部分我们集中考虑 variable elimination 这种策略。

不过我们应该清楚的认识到 inference 问题是“很难”的。，即便是近似解。但是一般说来我们的算法不会退化到最差的情况，因此我们觉得可以接受。

精确推断是 NPC 给定一个 BN，一个其中的随机变量 X 和可能的取值 x，决定 \Pr (X = x) 的概率是否大于零，这个问题是 NP-complete 的。

这个证明需要：说明这是一个多项式时间可验证的（给定一个 assignement 之后我们只需要按照 factor 求乘积即可）以及找到一个 NPC 说明存在一个多项式时间的归约。这个规约可以用 3-SAT，我们将 3-SAT 里面的 propositional variable 给定一个顶点 Q_k，每个从句给定一个顶点 C_k，并且将它们根据出现的关系使得 Q_k 是 C_j 的父节点，当且仅当该 variable 出现在对应从句中，且使用 deterministic CPD（导致与 3-SAT 类似的结果，即或关系）且对于 Q_k 的先验我们可以取均匀分布，我们的 X 是 C_k 的“与”，这样如果 3-SAT 成立当且仅当 \Pr(X = 1) >0。

更直接的我们可以证明计算 \Pr (X = x) 本身也是 #P-complete 的。很过分的事情不在于此，我们还有

任意的近似推断也是 NP-hard 我们称 \rho 是 \Pr (y \mid e) 的绝对误差不超过 \varepsilon 的估计，当且仅当 |P(y \mid e) - \rho| \leq \varepsilon；称 \rho 是 \Pr (y \mid e) 的相对误差不超过 \varepsilon 的估计，当且仅当 \frac{\rho}{1 + \varepsilon}\leq|P(y\mid e) - \rho| \leq \rho(1+\varepsilon)。这两个问题都是 NP-hard 的。

那么我们怎么处理这个问题呢？其实想法很简单，很直接（既然我们怎么做都不可能做得更好）。如果没有 evidence，我们相当于计算 maginal distribution \Pr(Y)，这时我们只需要去掉其他的变量就行了，比如剩余的为 X = (X_1, \ldots)，我们给定了这些变量之后一个一个加掉（积掉）他们即可。怎么加呢？我们知道如果把联合分布当成是 factor 的话，我们只需要选择与当前 X_i 相关的 factor 乘起来相加，这样就得到了这些 factor 关联的其他的 X_j 组成的 factor。这提示我们我们可以从无向图的角度来考虑：对 BN 来说我们先构造它的 moral graph，这对应于每个 CPD 对应的 factor，每个 factor 自己是一个 clique。如果给定了某个顺序，依次相加过程中消掉某些变量后可能导致几个 factor 的合并，这样得到的新 factor 仍然是一个 clique，这就会导致我们增加额外的边（fill edge）。如此消除了所有其他元素后就得到了 \Pr(Y)。

在有 evidence 的情况下问题也并没有边困难！为什么呢？我们等价于计算 \Pr(Y, E = e)，我们可以首先将 E 的部分替换成为对应的值，这导致我们在 reduced factor 上进行类似的计算即可。

这里面有两个核心的想法：

  • 由于某些 graph 的特殊结构，我们进行消元的时候不是所有的都会参与进来，往往参与的只是一小部分
  • 我们将中间的消元结果 cache，参与后续的消元就不必繁复的读取原先的 factor

另外如何选择合理的消元顺序也是一个有意思的问题。

为了将这个过程用算法实现，我们知道每个 factor 对应于一个 table，有关联的变量，我们需要有一个 table 的集合，通过变量获得相关 table 子集的方法，消去某个变量的方法。我们为每个 factor 建立一个类，其中包含数据区，若干个变量；每个变量有一个类，包含对应的维数，名称等信息。

1
2
3
4
5
6
7
8
9
10
struct drv { // for discrete random variables
  std::string name ;
  std::size_t n_vals ;
} ;
 
template <class RealType>
struct tbl {
  std::vector<drv*> var ;
  std::vector<RealType> data ;
} ;

这样我们建立一个 boost::bimap 用来检索（这是 pseudo code 了）：

1
2
3
4
5
boost::bimap<unique<tbl*>, multi<drv*> > m ;
for (int i = 0 ; i < tbls.size () ; ++ i) {
  for (int j = 0 ; j < tbls [i].var.size () ; ++ j)
    m.insert (std::make_pair (&tbls [i], tbls[i].var [j])) ;
}

这样给定了 elimination 顺序后，如 std::vector<drv*>，有

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
boost::shared_ptr<tbl>
eliminate (const std::vector<drv*>& order) {
  // assume we have m
  std::vector<boost::shared_ptr<tbl> > aux_tbls ; // for interim factors
  for (int i = 0 ; i < order.size () ; ++ i) {
    // find in bimap all related factors
    std::pair<..., ...> ranges = m.right.find (order[i]) ;
    // call sum-product
    boost::shared_ptr<tbl> new_factor = sum_product (ranges.first, ranges.second, order[i]) ;
    // remove factors from bimap
    // ...
    // add new to bimap
    // ...
    // save it
    aux_tbls.push_back (new_factor) ;
  }
  return aux_tbls[aux_tbls.size () - 1] ;
}

这里的 sum_product 的实现大致如下

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
boost::shared_ptr<tbl>
sum_product (Iterator begin, Iterator end, drv* a) {
  boost::shared_ptr<tbl> t = new tbl () ;
  std::size_t s = 1 ;
  // add vars to t and compute size
  for (Iterator iter = begin ; iter != end ; ++ iter) {
    for (int i = 0 ; i < (*iter).var.size () ; ++ i) {
      drv* cur = (*iter).var [i] ;
      if (a != cur) {
        t.var.push_back [cur] ;
        s *= cur->n_vals ;
      }
    }
  }
  t -> data.resize (s) ;
  // building the table
  // foreach value
  {
    // compute the product from each of the input tables
    // sum up and place the sum into output table
  }
  return t ;
}

嗯大致如此了。不知道使用 metaprogramming 会不会有点好处（还是纯粹折腾？）那么很重要的一点就是我们需要理解清楚这个算法的“瓶颈”在哪里？其实很容易看出来就是我们中间产生的那个 table 里面尺寸最大的，为了填充它，我们将耗费 O(d^k) 的代价（假定每个 r.v.s 取 d 个值，这个 table 是 k 维的），这反应出来图的性质其实是对应的 clique size – 1。如果算上 fill edge，我们称这个 graph 为某个消元序下的 induced graph，这里 k + 1 其实就是最大 clique 的 size 了（我们称 k 为此图的 width，或者所谓 induced width）。这刻画了我们的消元算法的时间复杂度。

下面的定理告诉我们，最优消元序也不好弄 -,-b

定理 给定图，界 K，确定是否存在不超过 K 的 induced width 是 NPC。

为此一般选择的策略都是 heuristic 的，通常有

  • min-neighbors 与 weighted 版本：weight 一般是 neighbor factor 含有变量个数的乘积
  • min-fill 与 weighted 版本：每条边的 weight 是两个顶点的 weight 的乘积

注意所有的 induced graph 都是 chordal，我们可以证明每个 chordal graph 也对应着一个消元序（即不会额外引入 fill edge）。这个消元序可以使用 maximal cardinality 算法获得（每次选择 marked neighbors 最多的点，并把当前点 mark，按此顺序的逆序即可）。注意因为我们并不知道消元序，所以我们无法直接使用这个算法，而如果对于非 chordal graph 这个算法获得的结果往往不如前面的。

除了以上的策略以外，还有一种称为 conditioning 的策略（video 中并没有 cover），其解决的核心问题是消元法可能产生的 table 尺寸过大而无法放入内存。那 conditioning 的 idea 是什么呢？其实也很简单就是直接拿某些东西进行累加：

\displaystyle\Pr (Y) = \sum_{u \in \mathrm{Val}(U)} \Pr(Y, U = u)

这样做的意义在于，我们只需要遍历所有可能的取值，而不必要存放其中的 factor，这样我们等价于每次用那部分 factor计算出来的一项完成结果中的一部分，然后遍历时累加最后得到结果，这明显没减少总共的计算次数。（与消元法的差异在于消元时每次都先加好了，最后再乘，所以只是顺序上的调整）这里复杂度分析也是借助对应 induced graph（没大看明白有什么用）。

在某些情况下 conditioning 可以得到加速：

  • conditioning 是从外面向里走，而消元是从里向外走（直观上来看处理连加的策略），因此很多时候我们可以借助消元简化再来做 conditioning；
  • 存在独立的部分，这时候网络就 decompose 了，conditioning 在独立部分的时候其实什么都可以不做

有了这两种基本的求解方法后，我们来看两个具体带有结构的 CPD 上如何简化。一个例子是 noisy OR，对于 OR 这类操作，简化的想法是将原因分开进行 OR 而不要一次全部 OR 起来，这样做会增大 conditioning 的代价。另一个例子是 context-specific independency，这时我们使用了 tree 或者 rule 类型的表示，这时我们就应该利用这些结构来改进 conditioning。（感觉这类型见的不多，后面有空详细看例子吧！）

——————-
Therefore Abimelech rose early in the morning, and called all his servants, and told all these things in their ears: and the men were sore afraid.

posted @ 2014-06-07 22:12  prml  阅读(281)  评论(0)  编辑  收藏  举报