# Thinkwee's Blog

Too Stupid to Give Up Learning

0%

Deep-Bayes 2018 Summer Camp的习题 填不动了，就到这吧

# Deep|Bayes summer school. Practical session on EM algorithm

• 第一题就是应用EM算法还原图像，人像和背景叠加在一起，灰度值的概率分布形式已知，设计人像在背景中的位置为隐变量，进行EM迭代推断。
• 具体说明在官网和下面的notebook注释中有，实际上公式已经给出，想要完成作业就是把公式打上去，可以自己推一下公式。

One of the school organisers decided to prank us and hid all games for our Thursday Game Night somewhere.

Let's find the prankster!

When you recognize him or her, send: * name * reconstructed photo * this notebook with your code (doesn't matter how awful it is :)

privately to Nadia Chirkova at Facebook or to info@deepbayes.ru. The first three participants will receive a present. Do not make spoilers to other participants!

Please, note that you have only one attempt to send a message!

### Data

We are given a set of $$K$$ images with shape $$H \times W$$.

It is represented by a numpy-array with shape $$H \times W \times K$$:

(100, 200, 1000)

Example of noisy image:

[255. 255.  41. 255.   0.  51. 255.  15. 255.  59.   0.   0.   0. 255.
0. 255. 255.   0. 175.  74. 184.   0.   0. 150. 255. 255.   0.   0.
148.   0. 255. 181. 255. 255. 255.   0. 255. 255.  30.   0.   0. 255.
0. 255. 255. 206. 234. 255.   0. 255. 255. 255.   0. 255.   0. 255.
0. 255. 255. 175.  30. 255.   0.   0. 255.   0. 255.  48.   0.   0.
0. 232. 162. 255.  26.   0.   0. 255.   0. 255. 173. 255. 255.   0.
0. 255. 255. 119.   0.   0.   0.   0.   0.   0. 255. 255. 255. 255.
0.   0. 248.   5. 149. 255.   0. 255. 255. 255.   0. 108.   0.   0.
255.   0. 255. 255. 255.   0.   0. 193.  79.   0. 255.   0.   0.   0.
233. 255.   0.  65. 255. 255. 255.   0. 255.   0.   0.   0. 255.  58.
226. 255.   0. 242. 255. 255.   0. 255.   4. 255. 255.  97. 255.   0.
0. 255.   0. 255.   0.   0.   0. 255.   0.  43. 219.   0. 255. 255.
255. 166. 255.   0. 255.  42. 255.  44. 255. 255. 255. 255. 255. 255.
255. 255.  28.   0.  52. 255.  81. 104. 255. 255.  48. 255. 255. 255.
102.  25.  30.  73.]

### Goal and plan

Our goal is to find face $$F$$ ($$H \times w$$).

Also, we will find: * $$B$$: background ($$H \times W$$) * $$s$$: noise standard deviation (float) * $$a$$: discrete prior over face positions ($$W-w+1$$) * $$q(d)$$: discrete posterior over face positions for each image (($$W-w+1$$) x $$K$$)

Implementation plan: 1. calculating $$log\, p(X \mid d,\,F,\,B,\,s)$$ 1. calculating objective 1. E-step: finding $$q(d)$$ 1. M-step: estimating $$F,\, B, \,s, \,a$$ 1. composing EM-algorithm from E- and M-step

### Implementation

#### 1. Implement calculate_log_probability

For $$k$$-th image $$X_k$$ and some face position $$d_k$$: $p(X_k \mid d_k,\,F,\,B,\,s) = \prod_{ij} \begin{cases} \mathcal{N}(X_k[i,j]\mid F[i,\,j-d_k],\,s^2), & \text{if}\, (i,j)\in faceArea(d_k)\\ \mathcal{N}(X_k[i,j]\mid B[i,j],\,s^2), & \text{else} \end{cases}$

Important notes: * Do not forget about logarithm! * This implementation should use no more than 1 cycle!

OK

#### 2. Implement calculate_lower_bound

$\mathcal{L}(q, \,F, \,B,\, s,\, a) = \sum_k \biggl (\mathbb{E} _ {q( d_k)}\bigl ( \log p( X_{k} \mid {d}_{k} , \,F,\,B,\,s) + \log p( d_k \mid a)\bigr) - \mathbb{E} _ {q( d_k)} \log q( d_k)\biggr)$

Important notes: * Use already implemented calculate_log_probability! * Note that distributions $$q( d_k)$$ and $$p( d_k \mid a)$$ are discrete. For example, $$P(d_k=i \mid a) = a[i]$$. * This implementation should not use cycles!

OK

#### 3. Implement E-step

$q(d_k) = p(d_k \mid X_k, \,F, \,B, \,s,\, a) = \frac {p( X_{k} \mid {d}_{k} , \,F,\,B,\,s)\, p(d_k \mid a)} {\sum_{d'_k} p( X_{k} \mid d'_k , \,F,\,B,\,s) \,p(d'_k \mid a)}$

Important notes: * Use already implemented calculate_log_probability! * For computational stability, perform all computations with logarithmic values and apply exp only before return. Also, we recommend using this trick: $\beta_i = \log{p_i(\dots)} \quad\rightarrow \quad \frac{e^{\beta_i}}{\sum_k e^{\beta_k}} = \frac{e^{(\beta_i - \max_j \beta_j)}}{\sum_k e^{(\beta_k- \max_j \beta_j)}}$ * This implementation should not use cycles!

OK

#### 4. Implement M-step

$a[j] = \frac{\sum_k q( d_k = j )}{\sum_{j'} \sum_{k'} q( d_{k'} = j')}$ $F[i, m] = \frac 1 K \sum_k \sum_{d_k} q(d_k)\, X^k[i,\, m+d_k]$ $B[i, j] = \frac {\sum_k \sum_{ d_k:\, (i, \,j) \,\not\in faceArea(d_k)} q(d_k)\, X^k[i, j]} {\sum_k \sum_{d_k: \,(i, \,j)\, \not\in faceArea(d_k)} q(d_k)}$ $s^2 = \frac 1 {HWK} \sum_k \sum_{d_k} q(d_k) \sum_{i,\, j} (X^k[i, \,j] - Model^{d_k}[i, \,j])^2$

where $$Model^{d_k}[i, j]$$ is an image composed from background and face located at $$d_k$$.

Important notes: * Update parameters in the following order: $$a$$, $$F$$, $$B$$, $$s$$. * When the parameter is updated, its new value is used to update other parameters. * This implementation should use no more than 3 cycles and no embedded cycles!

OK

#### 5. Implement EM_algorithm

Initialize parameters, if they are not passed, and then repeat E- and M-steps till convergence.

Please note that $$\mathcal{L}(q, \,F, \,B, \,s, \,a)$$ must increase after each iteration.

OK

### Who is the prankster?

To speed up the computation, we will perform 5 iterations over small subset of images and then gradually increase the subset.

If everything is implemented correctly, you will recognize the prankster (remember he is the one from DeepBayes team).

Run EM-algorithm:

And this is the background:

### Optional part: hard-EM

If you have some time left, you can implement simplified version of EM-algorithm called hard-EM. In hard-EM, instead of finding posterior distribution $$p(d_k|X_k, F, B, s, A)$$ at E-step, we just remember its argmax $$\tilde d_k$$ for each image $$k$$. Thus, the distribution q is replaced with a singular distribution: $q(d_k) = \begin{cases} 1, \, if d_k = \tilde d_k \\ 0, \, otherwise\end{cases}$ This modification simplifies formulas for $$\mathcal{L}$$ and M-step and speeds their computation up. However, the convergence of hard-EM is usually slow.

If you implement hard-EM, add binary flag hard_EM to the parameters of the following functions: * calculate_lower_bound * run_e_step * run_m_step * run_EM

After implementation, compare overall computation time for EM and hard-EM till recognizable F.