pgm5

这部分讨论 inference 里面基本的问题,即计算 Pr(Y mid E = e) 这类 query,这一般可以认为等价于计算 Pr (Y, E = e),因为我们只需要重新 normalize 一下关于 Y 的分布就得到了需要的值,特别是像 MAP 这类 query(一般此时 YE 的补集,可以理解成为取 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_kC_j 的父节点,当且仅当该 variable 出现在对应从句中,且使用 deterministic CPD(导致与 3-SAT 类似的结果,即或关系)且对于 Q_k 的先验我们可以取均匀分布,我们的 XC_k 的“与”,这样如果 3-SAT 成立当且仅当 Pr(X = 1) >0

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

任意的近似推断也是 NP-hard 我们称 
hoPr (y mid e) 的绝对误差不超过 varepsilon 的估计,当且仅当 |P(y mid e) - 
ho| leq varepsilon;称 
hoPr (y mid e) 的相对误差不超过 varepsilon 的估计,当且仅当 frac{
ho}{1 + varepsilon}leq|P(ymid e) - 
ho| leq 
ho(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 是什么呢?其实也很简单就是直接拿某些东西进行累加:

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

原文地址:https://www.cnblogs.com/focus-ml/p/3775438.html