Suppose we are given a matrix $V$ and our goal is to find non-negative matrices $W$ and $H$ such that $V \approx WH$. So we want to minimize $K(V || WH)$ (Kullback-Leibler Divergence) where $$K(V||WH) = ||V-WH||^{2} \tag{1}$$ or $$K(V||WH) = \sum_{i,j} \left[V_{i,j} \log \frac{V_{i,j}}{(WH)_{ij}}-V_{ij}+(WH)_{ij} \right] \tag{2}$$
depending on whether $\epsilon$ is normally distributed or poisson distributed respectively. Note that $\epsilon$ is noise. We want to minimize this subject to $W,H \geq 0$. For $K(V||WH)$ based on (1), we get $$H_{ak}^{t+1} = H_{ak}^{t} \left(\frac{\sum_{i}V_{ik} W_{ia}}{\sum_{i} W_{ia} \left(\sum_{b} W_{ib} H_{bk}^{t} \right)} \right)$$ and$$W_{ia}^{t+1} = H_{ia}^{t} \left(\frac{\sum_{i}H_{ak} V_{ik}}{\sum_{i} \left(\sum_{b} W_{ib} H_{bk}^{t} \right) H_{ak}} \right)$$
For $K(V||WH)$ based on (2), we get $$H_{ak}^{t+1} = H_{ak}^{t} \left(\frac{\sum_{i} \left(\frac{V_{ik}}{\sum_{b} W_{ib} H_{bk}^{t}} \right) W_{ia}}{\sum_{i} W_{ia}} \right)$$ and
$$W_{ia}^{t+1} = W_{ia}^{t} \left(\frac{\sum_{i} \left(\frac{V_{ik}}{\sum_{b} W_{ib}^{t} H_{bk}} \right) H_{ak}}{\sum_{k} H_{ak}} \right)$$
Why should these be the update rules? Are these iterations basically the best unbiased estimators for $W$ and $H$?
Added. When we factor $V$ into two non-negative matrices $W$ and $H$, we ultimately want to minimize the error between $V$ and $WH$. In other words $V = WH+\epsilon$ and we want to minimize $\epsilon$ which is the error term. We can assume that this error $\epsilon$ is normally distributed or poisson distributed. This is how we get $K(V||WH)$.