Bayesian Learning_ELBO, Variational Inference, and EM Algorithm

Reference:

https://mbernste.github.io/posts/elbo/

https://mbernste.github.io/posts/variational_inference/

https://mbernste.github.io/posts/em/

https://mbernste.github.io/posts/gmm_em/

Theodoridis S. Machine learning: a Bayesian and optimization perspective[M]. Academic press, 2015.

Latent variable model

• We posit that our observed data $x$ is a realization from some random variable $X$.
• We posit the existence of another random variable $Z$ where $X$ and $Z$ are distributed according to a joint distribution $p(X,Z;\theta )$ where $\theta$ parameterizes the distribution.
• Our data is only a realization of $X$, not $Z$, and therefore $Z$ remains unobserved (i.e. latent).

There are two predominant tasks that we may be interested in accomplishing:

1. Given some fixed value for $\theta$, compute the posterior distribution $p(Z|X;\theta )$ [Can be solved by variational inference]
2. Given that $\theta$ is unknown, find the maximum likelihood estimate of $\theta$ [Can be solved by EM]

Both variational inference and EM rely on the ELBO.

The Evidence Lower Bound (ELBO)

What is the ELBO?

To understand the evidence lower bound, we must first understand what we mean by “evidence”: it is just a name given to the likelihood function evaluated at a fixed $\theta$
$$\text{evidence}:=\log p(x;\theta )\text{\hspace{1em}(ELBO.1)}$$

Why is this quantity called the “evidence”?

Intuitively, if we have chosen the right model $p$ and $\theta$, then we would expect that the marginal probability of our observed data $x$, would be high. Thus, a higher value of $\log p(x;\theta )$ indicates, in some sense, that we may be on the right track with the model $p$ and parameters $\theta$ that we have chosen. That is, this quantity is “evidence” that we have chosen the right model for the data.

If we happen to know (or posit) that $Z$ follows some distribution denoted by $q$, s.t.
$$p(x,z;\theta ):=p(x|z;\theta )q(z)\text{\hspace{1em}(ELBO.2)}$$
Then the evidence lower bound is just a lower bound on the evidence that makes use of the known $q$. Specifically,
$$\log p(x;\theta )\ge E_{Z\sim q}\left[\log \frac{p(x,Z;\theta )}{q(Z)}\right]\text{\hspace{1em}(ELBO.3)}$$
where the ELBO is simply the right-hand side of the above equation:
$$ELBO:=E_{Z\sim q}\left[\log \frac{p(x,Z;\theta )}{q(Z)}\right]\text{\hspace{1em}(ELBO.4)}$$

Derivation

Jensen’s Inequality: if $X$ is a random variable and $\phi$ is a convex (concave) function, then
$$\phi (E[X])\le (\ge )E[\phi (X)]$$
Since $\log (\cdot )$ is a concave function,
$$\begin{array}{rl}{E}_{Z\sim q}\left[\log \frac{p(x,Z;\theta )}{q(Z)}\right]& \le \log \left[E_{Z\sim q}\left(\frac{p(x,Z;\theta )}{q(Z)}\right)\right]\\ & =\log \left[\int q(Z)\frac{p(x,Z;\theta )}{q(Z)}dz\right]\\ & =\log p(x;\theta )\end{array}$$

The gap between the evidence and the ELBO

It turns out that the gap between the evidence and the ELBO is precisely the Kullback Leibler divergence between $q(Z)$ and $p(Z|x;\theta )$.
$$\text{evidence}-\text{ELBO}:=\log p(x;\theta )-E_{Z\sim q}\left[\log \frac{p(x,Z;\theta )}{q(Z)}\right]=KL(q(Z)\Vert p(Z|x;\theta ))\text{\hspace{1em}(ELBO.5)}$$
This fact forms the basis of the [variational inference algorithm] for approximate Bayesian inference.

Derivation
$$\begin{array}{rl}\log p(x;\theta )-E_{Z\sim q}\left[\log \frac{p(x,Z;\theta )}{q(Z)}\right]& =\int q(Z)\log p(x;\theta )dz-\int q(Z)\log \frac{p(x,Z;\theta )}{q(Z)}dz\\ & =\int q(Z)\log \frac{p(x;\theta )q(Z)}{p(x,Z;\theta )}dz\\ & =\int q(Z)\log \frac{q(Z)}{p(Z|x;\theta )}dz\\ & =KL(q(Z)\Vert p(Z|x;\theta ))\end{array}$$

Variational Inference

Why variational inference?

Variational inference is a paradigm for estimating a posterior distribution when computing it explicitly is intractable.

Assume that we have a model that involves hidden random variables $Z$, observed random variables $X$, and some posited probabilistic model over the hidden and the observed random variables $P(X,Z)$. The goal is to compute the posterior distribution $P(Z|X)$.

Ideally, we would do so by using Bayes theorem:
$$p(z|x)=\frac{p(x|z)p(z)}{p(x)}$$
In practice, it is often difficult to compute $p(z|x)$ via Bayes theorem because the denominator $p(x)$ does not have a closed form. Usually, the denominator $p(x)$ can be only be expressed as an integral that marginalizes over $z$: $p(x)=\int p(x,z)dz$. In such scenarios, we’re often forced to approximate $p(z|x)$ rather than compute it directly. Variational inference is one such approximation technique.

Details

Instead of computing $p(z|x)$ exactly via Bayes theorem, variational inference attempts to find another distribution $q(z)$ that is ‘close’ to $p(z|x)$, where the ‘closeness’ is measured by the KL-divergence.
$$KL(q(Z)\Vert p(Z|x))=\int q(Z)\log \frac{q(Z)}{p(Z|x)}dz=E_{Z\sim q}\left[\log \frac{q(Z)}{p(Z|x)}\right]\text{\hspace{1em}(VI.1)}$$
We restrict the search for $q(z)$ to a family of surrogate distributions over $Z$, called the variational distribution family, denoted by the set of distributions $\mathcal{Q}$. Each member of $\mathcal{Q}$ is characterized by the values of a set of variational parameters $\varphi$.

Then we can formulate the optimization problem as:
$$\hat{q}:=q(z|\hat{\varphi })=\arg \underset{q}{\min }KL(q(Z)\Vert p(Z|x))\text{\hspace{1em}(VI.2)}$$
From $(ELBO.5)$, we know that minimizing the KL-divergence is equivalent to maximizing the ELBO:
$$\hat{q}:=q(z|\hat{\varphi })=\arg \underset{q}{\max }ELBO(q)\text{\hspace{1em}(VI.3)}$$
where $ELBO(q):=E_{Z\sim q}[\log p(x,Z)]-E_{Z\sim q}[\log p(Z)]$.

Why is the method called ‘variational’ inference?

The term “variational” in “variational inference” comes from the mathematical area of the calculus of variations. The calculus of variations is all about optimization problems that optimize functions of functions (called functionals).

More specifically, let’s say we have some set of functions $\mathcal{F}$ where each $f\in \mathcal{F}$ maps items from some set $A$ to some set $B$. That is,
$$f:A\to B$$
Let’s say we have some function $g$ that maps functions in $\mathcal{F}$ to real numbers $\mathbb{R}$. That is,
$$g:\mathcal{F}\to \mathbb{R}$$
Then, we may wish to solve an optimization problem of the form:
$$\arg \underset{f\in \mathcal{F}}{\max }g(f)$$
This is precisely the problem addressed in the calculus of variations. In the case of variational inference, the functional, $g$, that we are optimizing is the ELBO. The set of functions, $\mathcal{F}$, that we are searching over is the set of measurable functions in the variational family, $\mathcal{Q}$.

Expectation-Maximization (EM)

Expectation-maximization (EM) is a popular algorithm for performing maximum-likelihood estimation of the parameters in a [latent variable model](# Latent variable model).

Despite not observing $z$, we wish to find the value for $\theta$ that makes $x$ most likely under our model. Since we have not observed $z$, our likelihood function must marginalize over $z$:
$$l(\theta ):=\log p(x;\theta )=\log \int p(x,z;\theta )dz\text{\hspace{1em}(EM.1)}$$
In practice, maximizing this function over $\theta$ may be difficult to do analytically due to this integral. For such situations, the EM algorithm may provide a method for computing a local maximum of this function with respect to $\theta$.

Description of EM

The EM algorithm alternates between two steps: an expectation-step (E-step) and a maximization-step (M-step). Throughout execution, the algorithm maintains an estimate of the parameters $\theta$ that is updated on each iteration. As the algorithm progresses, this estimate of $\theta$ converges to a local maximum of $l(\theta )$.

Let $t$ denote a given step of the algorithm. On the $t$-th step, we’ll denote ${\theta }_t$ as the $t$-th estimate of $\theta$.

• E-Step

  • Compute the conditional probability $p(z|x;{\theta }_t)$

  • Then formulate the Q-function:
    $$Q_t(\theta ):=E_{Z|x,{\theta }_t}[\log p(x,z;\theta )]=\int p(z|x;{\theta }_t)\log p(x,z;\theta )dz\text{\hspace{1em}(EM.2)}$$

• M-Step

  Assign ${\theta }_{t+1}$ to be the value that maximizes the Q-function that was formulated on the previous E-step:
  $${\theta }_{t+1}:=\arg \underset{\theta }{\max }{Q}_t(\theta )\text{\hspace{1em}(EM.3)}$$

Note:

The core of formulating an EM algorithm for a particular model is calculating $p(z|x;{\theta }_t)$. If the model admits an easy calculation of this function, then the EM algorithm may be a good a choice.

We repeat the E-Step and M-Step until the difference between ${\theta }_t$ and ${\theta }_{t+1}$ is small, indicating convergence. We also note that ${\theta }_0$ can be set to an arbitrary value – usually set randomly. Since EM finds a local maximum of $l(\theta )$, it is often wise to try many different ${\theta }_0$’s.

EM is coordinate ascent on the evidence lower bound

What is coordinate ascent?

Given a function $f(a,b)$ that takes two arguments $a$ and $b$, the idea of coordinate ascent is that we will iteratively fix $a$ to some value $\hat{a}$ and then choose for $b$ the value $\hat{b}$ that maximizes $f(\hat{a},b)$. Given $\hat{b}$, we then re-assign to $a$ the value that maximizes $f(a,\hat{b})$. This procedure is repeated until convergence.

In our setting, we do not yet know which value to use for $\theta$, nor do we know the distribution $q$ for $Z$. Thus, the ELBO becomes a function of both $q$ and $\theta$:

$$F(q,\theta ):=E_{Z\sim q}\left[\log \frac{p(x,Z;\theta )}{q(Z)}\right]\text{\hspace{1em}(EM.4)}$$
EM is just coordinate ascent on the ELBO. In the E-step, we fix $\theta$ and solve for $q$. In the M-step, we fix $q$ and solve for $\theta$.

• Fix $\theta$ to ${\theta }_t$ (the current value for $\theta$) and solve for $q$:
  $$\begin{array}{rl}{\mathrm{argmax}}_{q}F\left(q,{\theta }_t\right)& ={\mathrm{argmax}}_q{E}_{z\sim q}\left[\log \frac{p(Z,x;\theta )}{q(Z)}\right]\\ & ={\mathrm{argmin}}_q\mathrm{K}\mathrm{L}\left(q(z)\Vert p\left(z|x;{\theta }_t\right)\right)\\ & =p\left(z|x;{\theta }_t\right)\end{array}$$
  [[Make use of Eq. $(ELBO.5)$]](# The gap between the evidence and the ELBO)

  Notably, computing $p\left(z|x;{\theta }_t\right)$ is exactly what is needed to formulate the Q-function in the E-Step.

• Next, if we hold $q$ fixed to $q_t:=p\left(z|x;{\theta }_t\right)$, we see that maximizing $F\left(q_t,\theta \right)$ with respect to $\theta$ is equivalent to maximizing the Q-function:
  $$\begin{array}{rl}{\mathrm{argmax}}_{\theta }F\left(q_t,\theta \right)& ={\mathrm{argmax}}_{\theta }{E}_{z\sim q}\left[\log \frac{p(Z,x;\theta )}{p\left(z|x;{\theta }_t\right)}\right]\\ & ={\mathrm{argmax}}_{\theta }{E}_{Z|x,{\theta }_t}[\log p(x,z;\theta )]\end{array}$$
  The last line in the above derivation is simply maximizing the Q-function - exactly the M-Step.

Visualization

• Coordinate ascent on the evidence lower bound
  

• EM algorithm - overview

We can observe that EM is in line with a more general class of optimization algorithms, known as minorize-maximization (MM) methods. A surrogate function that minorizes the cost function (lower bound), which is easier to be maximized, is employed. Then maximizing this lower bound iteratively pushes the cost function to a local maximum.

Intuition behind the Q-function

• A generalization of the likelihood function

  In the case in which $Z$ is a discrete random variable, the Q-function can be interpreted as a “generalized” version of the complete data likelihood function.

  Let us say that both $X$ and $Z$ are observed where $X=x$ and $Z=z$. Then, we could write the complete data likelihood as
  $$\log p(x,z;\theta ):=\sum _{z^{\prime }\in \mathcal{Z}}\mathbb{I}\left(z=z^{\prime }\right)\log p\left(x,z^{\prime };\theta \right)$$
  where $\mathcal{Z}$ is the support of $Z$'s probability mass function and $\mathbb{I}\left(z=z^{\prime }\right)$ is an indicator that equals one when $z^{\prime }=z$ and equals zero otherwise.

  Now recall the Q-function:
  $$Q_t(\theta ):=\sum _{z^{\prime }}p(z^{\prime }|x;{\theta }_t)\log p(x,z^{\prime };\theta )$$
  We note that the difference between these two functions is that the indicator variable in the complete data likelihood $\mathbb{I}\left(z=z^{\prime }\right)$ is replaced by the probability $p(z^{\prime }|x;{\theta }_t)$. This leads us to viewing the Q-function as a sort of generalization of the complete data likelihood. That is, $p(z^{\prime }|x;{\theta }_t)$ acts as a weight for its corresponding term in the summation, and measures our current certainty that the hidden data is equal to $z^{\prime }$.

• A likelihood function over a hypothetical dataset

  Another way to view the Q-function is as a complete data likelihood function over a hypothetical dataset where this hypothetical dataset is generated from a distribution that depends on our current best guess of $\theta$ (i.e. ${\theta }_t$ at the $t$-th iteration).

  That is, let us assume that we generate a very large set of values for $Z$ drawn i.i.d. from the distribution specified by $p(z|x;{\theta }_t)$:
  $$z_1^{\prime },z_2^{\prime },\cdots ,z_n^{\prime }\sim p(z|x;{\theta }_t)$$
  For each $z_i^{\prime }$, we create a new ‘sub’-dataset $(x,z_i^{\prime })$. Note that the observed data $x$ is duplicated in each of these sub-datasets. We then merge all of the datasets $(x,z_1^{\prime }),\cdots ,(x,z_n^{\prime })$ to form our new hypothetical dataset. Then, the likelihood over these data is
  $$l^{\prime }(\theta ):=\prod _{i=1}^{n}p(x,z_i^{\prime };\theta )$$
  Now we will show that maximizing the Q-function is equivalent to maximizing this likelihood function over the hypothetical data.

  Let Z ′ : = { z 1 ′ , z 2 ′ , ⋯   , z n ′ } Z':=\{z_1',z_2',\cdots,z_n'\} Z′:={z1′​,z2′​,⋯,zn′​}. Now let $c(z)$ be the cardinality of { z ′ ∈ Z ′ ∣ z ′ = z } \{z'\in Z'|z'=z\} {z′∈Z′∣z′=z}. That is, $c(z)$ is the number of items in $Z^{\prime }$ that equal $z$. As $n$ approaches infinity then
  $$\frac{c(z)}{{\sum }_{z^{\ast }}c(z^{\ast })}=p(z|x;{\theta }_t)$$
  Finally we can formulate our maximum likelihood problem on this hypothetical data:
  $$\begin{array}{rl}{\mathrm{argmax}}_{\theta }{l}^{\prime }(\theta )& :={\mathrm{argmax}}_{\theta }\prod _{i=1}^{n}p\left(x,z_i^{\prime };\theta \right)\\ & ={\mathrm{argmax}}_{\theta }\prod _{z\in \mathcal{Z}}p(x,z;\theta {)}^{c(z)}\\ & ={\mathrm{argmax}}_{\theta }\sum _{z\in Z}\frac{c(z)}{{\sum }_{z^{\ast }}c\left(z^{\ast }\right)}\log p(x,z;\theta )\\ & ={\mathrm{argmax}}_{\theta }\sum _{z\in Z}p\left(z|x;{\theta }_t\right)\log p(x,z;\theta )\end{array}$$
  The final line in the above calculation is the Q-function!

Gaussian Mixture Model (GMM)

Definition

Mixture modelling: Model the unknown PDF, $p(\mathbf{x})$, as a linear combination of different distributions, i.e.,
$$p(\mathbf{x})=\sum _{k=1}^K{\alpha }_{k}p(\mathbf{x}|k)\text{\hspace{1em}(GMM.1)}$$
where ${\alpha }_k$ is the parameter weighting the specific contributing PDF, $p(\mathbf{x}|k)$. To guarantee that $p(\mathbf{x})$ is a PDF, the weighting parameters must be nonnegative and add to one (${\sum }_{k=1}^K{\alpha }_k=1$).

The physical interpretation of Eq. $(GMM.1)$ is that we are given a set of $K$ distributions, $p(\mathbf{x}|k),k=1,2,\cdots ,K$. Each observation ${\mathbf{x}}_n,n=1,2,\cdots ,N$, is drawn from one of these $K$ distributions, but we are not told from which one. All we know is a set of parameters, ${\alpha }_k,k=1,2,\cdots ,K$, each one providing the probability that a sample has been drawn from the corresponding PDF, $p(\mathbf{x}|k)$.

It can be shown that for a large enough number of mixtures, $K$, and appropriate choice of the involved parameters, one can approximate arbitrarily close any continuous PDF.

Mixture modelling is a typical task involving latent variables; that is, the labels $k$ of the PDF from which an obtained observation has originated. In practice, each $p(\mathbf{x}|k)$ is chosen from a known PDF family, parameterized via a set of parameters, and $(GMM.1)$ can be rewritten as
$$p(\mathbf{x};\mathbf{\Theta })=\sum _{k=1}^K{\alpha }_{k}p(\mathbf{x}|k;{\mathbf{\Theta }}_k)\text{\hspace{1em}(GMM.2)}$$
and the task is to estimate $({\alpha }_k,{\mathbf{\Theta }}_k),k=1,2,\cdots ,K$, based on a set of observations ${\mathbf{x}}_n,n=1,2,\cdots ,N$.

If we choose the PDFs $p(\mathbf{x}|k;{\mathbf{\Theta }}_k)$ to be Gaussian, then we obtain the GMM:
$$p(\mathbf{x};\mathbf{\Theta })=\sum _{k=1}^K{\alpha }_k\varphi (\mathbf{x};{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)\text{\hspace{1em}(GMM.3)}$$
where $\varphi$ is the PDF of the multivariate Gaussian distribution:
$$\varphi (\mathbf{x};\mathbf{\mu },\mathbf{\Sigma }):=\frac{1}{(2\pi {)}^{\frac{N}{2}}\det (\mathbf{\Sigma }{)}^{\frac{1}{2}}}\exp \left[-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })\right]$$
and
Θ : = { μ 1 , ⋯   , μ K , Σ 1 , ⋯   , Σ K , α 1 , ⋯   , α K } \mathbf \Theta:=\{\boldsymbol \mu_1,\cdots,\boldsymbol \mu_K,\boldsymbol \Sigma_1,\cdots,\boldsymbol \Sigma_K,\alpha_1,\cdots,\alpha_K \} Θ:={μ1​,⋯,μK​,Σ1​,⋯,ΣK​,α1​,⋯,αK​}
Here’s an example density function for a two-dimensional GMM with three Gaussians (i.e. $K=3$):

A model for data clustering

GMM’s provide a framework for finding clusters within data points ${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N$. To cluster the data, we make a very strong assumption that our data points were samples from a GMM with $K$ Gaussians. Unfortunately, we do not know the GMM’s parameters $\mathbf{\Theta }$, nor do we know which Gaussian generated each data point.

If we denote $z_1,\cdots ,z_N$ to be the label of ${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N$, then ${\mathbf{x}}_1,\cdots ,{\mathbf{x}}_N$ constitutes the observed data and $z_1,\cdots ,z_N$ constitutes the latent data. The situation we are facing can be depicted in the figure below:

Parameter estimation for such a problem naturally lends itself to be treated via the [EM algorithm](# Expectation-Maximization (EM)).

Maximum-likelihood estimation for GMM’s via the EM algorithm

• E-Step

  • Compute the conditional probability $p(Z_n=k|{\mathbf{x}}_n;{\mathbf{\Theta }}_{k,t})$
    $$\begin{array}{rl}p(Z_n=k|{\mathbf{x}}_n;{\mathbf{\Theta }}_{k,t})& =\frac{p({\mathbf{x}}_n|Z_n=k;{\mathbf{\Theta }}_{k,t})p(Z_n=k)}{{\sum }_{k=1}^{K}p({\mathbf{x}}_n|Z_n=k;{\mathbf{\Theta }}_{k,t})p(Z_n=k)}\\ & =\frac{\varphi ({\mathbf{x}}_n|{\mathbf{\mu }}_{k,t},{\mathbf{\Sigma }}_{k,t}){\alpha }_{k,t}}{{\sum }_{k=1}^K\varphi ({\mathbf{x}}_n|{\mathbf{\mu }}_{k,t},{\mathbf{\Sigma }}_{k,t}){\alpha }_{k,t}}\end{array}\text{\hspace{1em}(GMM.4)}$$
    For simplicity, denote ${\gamma }_{n,k,t}\triangleq p(Z_n=k|{\mathbf{x}}_n;{\mathbf{\Theta }}_{k,t})$.

  • Then formulate the Q-function from $(GMM.3)$ and $(GMM.4)$:
    $$Q_t(\mathbf{\Theta }):=\sum _{n=1}^N\sum _{k=1}^K{\gamma }_{n,k,t}\log [{\alpha }_k\varphi ({\mathbf{x}}_n;{\mathbf{\mu }}_k,{\mathbf{\Sigma }}_k)]\text{\hspace{1em}(GMM.5)}$$

• M-Step

  Assign ${\mathbf{\Theta }}_{t+1}$ to be the value that maximizes the Q-function that was formulated on the previous E-step:
  $${\mathbf{\Theta }}_{t+1}:=\arg \underset{\mathbf{\Theta }}{\max }{Q}_t(\mathbf{\Theta })\text{s.t. }\sum _{k=1}^K{\alpha }_{k,t}=1\text{\hspace{1em}(GMM.6)}$$

  The solution to this optimization problem is given by
  $$\begin{array}{rl}{\alpha }_{k,t+1}& :=\frac{1}{n}\sum _{i=1}^n{\gamma }_{n,k,t}\\ {\mathbf{\mu }}_{k,t+1}& :=\frac{1}{{\sum }_{n=1}^N{\gamma }_{n,k,t}}\sum _{n=1}^N{\gamma }_{n,k,t}{\mathbf{x}}_n\\ {\mathbf{\Sigma }}_{k,t+1}& :=\frac{1}{{\sum }_{n=1}^N{\gamma }_{n,k,t}}\sum _{i=1}^n{\gamma }_{n,k,t}\left({\mathbf{x}}_n-{\mathbf{\mu }}_{k,t}\right){\left({\mathbf{x}}_n-{\mathbf{\mu }}_{k,t}\right)}^T\end{array}\text{\hspace{1em}(GMM.7)}$$
  Note:

  • ${\alpha }_{k,t+1}$ can be solved using Lagrange multipliers. The value of $\lambda$ can be determined by making use of ${\sum }_{k=1}^K{\gamma }_{n,k,t}=1$.
  • $\frac{\partial }{\partial s}(\mathbf{x}-\mathbf{s}{)}^T\mathbf{W}(\mathbf{x}-\mathbf{s})=-2{\mathbf{W}}^{-1}(\mathbf{x}-\mathbf{s})$
  • $\frac{\partial }{\partial \mathbf{X}}\log |\det (\mathbf{X})|={\left({\mathbf{X}}^T\right)}^{-1}$
  • $\frac{\partial }{\partial \mathbf{X}}{\mathbf{a}}^T{\mathbf{X}}^{-1}\mathbf{b}=-{\left({\mathbf{X}}^{-1}\right)}^T\mathbf{a}{\mathbf{b}}^T{\left({\mathbf{X}}^{-1}\right)}^T$

  The complete derivation can be found here.

Visualization: