Theo Wolf · Oxford

Given observations \(\mathcal D = \{x_i\}_{i=1}^{N}\) from an unknown data-generating distribution, the goal is to fit a distribution from which we can generate new samples.

\[ \mathcal D = \{x_i\}_{i=1}^{N} \overset{\mathrm{iid}}{\sim} p_{\mathrm{data}}(x), \qquad p_\theta(x) \approx p_{\mathrm{data}}(x) \quad \Longrightarrow \quad \widetilde{x} \sim p_\theta(x). \]

Parameterise the model by \(\theta\), then choose the parameters that assign the greatest likelihood to the observed dataset:

\[ \theta^\star = \underset{\theta}{\arg\max}\; \log p_\theta(\mathcal D). \]

Introduce latent variables \(z_i\) that explain the data:

\[ z_i \sim p(z_i), \qquad x_i \sim p_\theta(x_i \mid z_i), \qquad \mathcal Z = \{z_i\}_{i=1}^{N}, \] \[ p_\theta(\mathcal D,\mathcal Z) = \prod_{i=1}^{N} p_\theta(x_i \mid z_i)p(z_i). \] \[ \underbrace{p_\theta(\mathcal D)}_{\text{evidence}} = \underbrace{\int p_\theta(\mathcal D,\mathcal Z) \,d\mathcal Z}_{\text{marginalise every latent}} = \prod_{i=1}^{N} \int p_\theta(x_i \mid z_i)p(z_i)\,dz_i, \qquad d\mathcal Z = \prod_{i=1}^{N}dz_i. \] \[ \log p_\theta(\mathcal D) = \sum_{i=1}^{N} \log \underbrace{\int p_\theta(x_i \mid z_i)p(z_i)\,dz_i}_ {\text{generally intractable}}. \]

Introduce the tractable encoder \(q_\phi(z_i \mid x_i)\) to approximate the intractable posterior \(p_\theta(z_i \mid x_i)\).

\[ \begin{align*} \log p_\theta(x_i) &= \log \mathbb{E}_{q_\phi(z_i \mid x_i)} \left[\frac{p_\theta(x_i,z_i)} {q_\phi(z_i \mid x_i)}\right] \\ &\underset{\href{https://en.wikipedia.org/wiki/Jensen's_inequality} {\text{Jensen's inequality}}}{\geq} \underbrace{\mathbb{E}_{q_\phi(z_i \mid x_i)} \left[\log \frac{p_\theta(x_i,z_i)} {q_\phi(z_i \mid x_i)}\right]}_ {\displaystyle \mathcal L_{\mathrm{ELBO}}(x_i)} \\ &= \mathbb{E}_{q_\phi(z_i \mid x_i)} \left[\log p_\theta(x_i \mid z_i)\right] - D_{\mathrm{KL}}\!\left( q_\phi(z_i \mid x_i)\,\|\,p(z_i)\right). \end{align*} \] \[ \log p_\theta(\mathcal D) \geq \sum_{i=1}^{N}\mathcal L_{\mathrm{ELBO}}(x_i) \equiv \mathcal L_{\mathrm{ELBO}}(\mathcal D). \]

For a Gaussian VAE with fixed decoder variance:

\[ \begin{align*} p(z_i) &= \mathcal{N}(0,I) && \text{prior}, \\ q_\phi(z_i \mid x_i) &= \mathcal{N}\!\left(\mu_\phi(x_i), \operatorname{diag}(\sigma_\phi^2(x_i))\right) && \text{encoder}, \\ p_\theta(x_i \mid z_i) &= \mathcal{N}\!\left(f_\theta(z_i),\sigma_x^2 I\right) && \text{decoder}. \end{align*} \] \[ -\mathbb{E}_{q_\phi(z_i \mid x_i)} \!\left[\log p_\theta(x_i \mid z_i)\right] = \frac{1}{2\sigma_x^2} \mathbb{E}_{q_\phi(z_i \mid x_i)} \left[\lVert x_i-f_\theta(z_i)\rVert_2^2\right] + \frac{d_x}{2}\log(2\pi\sigma_x^2). \] \[ D_{\mathrm{KL}}\!\left( q_\phi(z_i \mid x_i)\,\|\,p(z_i)\right) = \frac{1}{2}\sum_{j=1}^{d_z} \left( \mu_{\phi,j}(x_i)^2 + \sigma_{\phi,j}(x_i)^2 - \log \sigma_{\phi,j}(x_i)^2 - 1\right). \] \[ \boxed{\begin{aligned} (\theta^\star,\phi^\star) &= \underset{\theta,\phi}{\arg\min}\; -\frac{1}{N}\mathcal L_{\mathrm{ELBO}}(\mathcal D) \\ &= \underset{\theta,\phi}{\arg\min}\; \frac{1}{N}\sum_{i=1}^{N}\Bigg[ \underbrace{\frac{1}{2\sigma_x^2} \mathbb{E}_{q_\phi(z_i \mid x_i)} \lVert x_i-f_\theta(z_i)\rVert_2^2}_ {\text{reconstruction}} + \underbrace{D_{\mathrm{KL}}\!\left( q_\phi(z_i \mid x_i)\,\|\,p(z_i)\right)}_ {\text{regularisation}}\Bigg]. \end{aligned}} \]

It is often assumed that the exponential function is the solution to the linear dynamical system.

\[ \begin{align*} \dot{x} &= Ax \\ x(0) &= x_0 \end{align*} \]

The solution is given by:

\[ x(t) = e^{At}x_0 \]

Where \(e^{At}\) is the matrix exponential, not specifically 2.718... to the power of a matrix (which is not defined), but \(e\) as the exponential function. This is defined as the infinite series:

\[ e^{A} = \sum_{n=0}^{\infty} \frac{A^n}{n!} \]

Let's start by taking small steps \(\Delta t\) and approximating the solution to the differential equation as:

\[ \begin{align*} x_{\Delta t} &= x_0 + A\Delta t x_0 \\ x_{2\Delta t} &= x_{\Delta t} + A\Delta t x_{\Delta t} \\ &= x_0 + A\Delta t x_0 + A\Delta t (x_0 + A\Delta t x_0) \\ &= x_0 + A\Delta t x_0 + A\Delta t x_0 + A^2\Delta t^2 x_0 \\ &= x_0 + 2A\Delta t x_0 + A^2\Delta t^2 x_0 \\ &= x_0 (I + 2A\Delta t + A^2\Delta t^2) \\ &= x_0 (I + A\Delta t)^2 \end{align*} \]

We can see that expanding to the \(n\)th step is:

\[ x_{n\Delta t} = x_0 (I + A\Delta t)^n \]

We set \(\Delta t = \frac{t}{n}\) and take the limit as \(n \to \infty\) (making the steps infinitely small):

\[ \lim_{n \to \infty} \Big[I+\frac{At}{n}\Big]^n \]

We can recognise the definition of the exponential function, which is a reformulation of the infinite series:

\[ e^{x} = \lim_{n \to \infty} \Big[1+\frac{x}{n}\Big]^n \]

Thus, we have shown that:

\[ x(t) = e^{At}x_0 \]

The assumptions of special relativity lead to the most famous equation in physics: \(E=mc^2\). It derives naturally from the assumption that the speed of light is the same in all frames of reference. Starting from the definition of kinetic energy:

\[ \begin{align*} \frac{dK}{dt} = v\frac{dp}{dt} = vm\frac{d}{dt}(\gamma v) \end{align*} \]

Where we have used the relativistic momentum \(p = \gamma mv\), and \(\gamma = (1-v^2/c^2)^{-1/2}\) is the Lorentz factor:

\[ \begin{align*} \frac{d(\gamma v)}{dt} &= \frac{d}{dt} \frac{v}{\sqrt{1-v^2/c^2}} \\ &= \Bigg[\Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-1/2} + \frac{v^2}{c^2} \Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-3/2}\Bigg] \frac{dv}{dt} \\ &= \Bigg[\Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-1/2}\Bigg\{1 + \frac{v^2}{c^2} \Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-1} \Bigg\} \Bigg] \frac{dv}{dt}\\ &= \Bigg[\Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-1/2}\Bigg\{\frac{1 -v^2/c^2}{1 -v^2/c^2} + \frac{v^2/c^2}{1 -v^2/c^2}\Bigg\} \Bigg] \frac{dv}{dt}\\ &= \Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-3/2} \frac{dv}{dt} \end{align*} \]

Hence:

\[ \begin{align*} \frac{dK}{dt} &= m\Bigg(1 - \frac{v^2}{c^2}\Bigg)^{-3/2} v\frac{dv}{dt} \\ &= \frac{d}{dt} \frac{m c^2}{\sqrt{1-v^2/c^2}} = \frac{d}{dt} \Big(\gamma m c^2\Big) \end{align*} \]

Where we have used the fact that:

\[ \begin{align*} \frac{d}{dt} \frac{c^2}{\sqrt{1-v^2/c^2}} = \Big(1 - \frac{v^2}{c^2}\Big)^{-3/2} v\frac{dv}{dt} \end{align*} \]

Integrating both sides, knowing that the kinetic energy is zero at rest (\(\gamma=1\) when \(v=0\)), we get:

\[ \begin{align*} K &= \gamma m c^2 + C \\ &= \gamma m c^2 - mc^2 \end{align*} \]

We can recover the classic Newtonian kinetic energy \(K = \frac{1}{2} m v^2\) by taking \(v \ll c\):

\[ \begin{align*} K &= \gamma m c^2 - mc^2 \\ &= mc^2\Bigg(\frac{1}{\sqrt{1-v^2/c^2}} - 1\Bigg) \\ &= mc^2\Bigg(1 + \frac{1}{2} \frac{v^2}{c^2} + ... - 1\Bigg) \\ &\approx \frac{1}{2} m v^2 \end{align*} \]

We can interpret \(K = \gamma m c^2 - mc^2\) as \(K = T - E\) where \(T\) is the total energy and \(E\) is the rest energy. Hence, we have:

\[ E = m c^2 \]

This one is just beautiful, links \(\pi\) and \(e\) together even though they really don't seem related.

\[ \begin{align*} \int{e^{-x^2} dx} &= \sqrt{\int{e^{-x^2} dx} \int{e^{-x^2} dx}} \end{align*} \]

Switching to a dummy variable:

\[ \begin{align*} &= \sqrt{\int{e^{-x^2} dx} \int{e^{-y^2} dy}} \\ &= \sqrt{\int\int{e^{-(x^2 + y^2)} dx\, dy}} \end{align*} \]

We can now swap to polar coordinates by using well-known identities.

\[ \begin{align*} x &= r \cos(\theta) \\ y &= r \sin(\theta) \\ dx\, dy &= r dr\, d\theta \end{align*} \]

Subbing everything in:

\[ \begin{align*} \int\int{e^{-(x^2 + y^2)} dx\, dy} &= \int\int{e^{-(r \cos(\theta))^2 - (r \sin(\theta))^2} r dr\, d\theta } \\ &= \int\int{e^{-r^2 (\cos(\theta)^2+\sin(\theta)^2 )} r dr\, d\theta} \\ &= \int\int{e^{-r^2} r dr\, d\theta} \\ &= \int_0^{2\pi} d\theta \int_0^{\infty} e^{-r^2} r dr \\ &= 2\pi \int_0^{\infty} e^{-r^2} r dr \\ &= 2\pi \int_0^{\infty} e^{-u} \frac{1}{2} du \\ &= \pi\\ \text{And so:}\\ \int{e^{-x^2} dx} &= \sqrt{\pi} \end{align*} \]

The score function trick is just so elegant, the basis of many algorithms in ML. Here, I show it off for the REINFORCE policy gradient algorithm.

Consider a distribution of states \(d(s)\) and a parameterised policy \(\pi_{\theta}(a | s)\). The objective is to maximise the expected reward under the policy.

\[ \begin{align*} J(\theta) = \mathbb{E}_{\pi_\theta, d}[R(a, s)] \end{align*} \] \[ \begin{align*} \nabla_{\theta} J(\theta) &= \nabla_{\theta} \int d(s) ds \int \pi_\theta(a | s) R(a, s) da \\ &= \int d(s) ds \int \nabla_{\theta} \pi_\theta(a|s) R(a,s) da \\ &= \int d(s) ds \int \pi_\theta(a|s) \frac{\nabla_{\theta} \pi_\theta(a|s)}{\pi_\theta(a|s)} R(a,s) da \end{align*} \]

We apply the score function trick:

\[ \nabla_{\theta} \log \pi_\theta(a|s) = \frac{\nabla_{\theta} \pi_\theta(a|s)}{\pi_\theta(a|s)} \] \[ \begin{align*} \nabla_{\theta} J(\theta) &= \int d(s) ds \int \pi_\theta(a|s) \nabla_{\theta} \log \pi_\theta(a|s) R(a,s) da \\ &= \mathbb{E}_{\pi_\theta, d} [\nabla_{\theta} \log \pi_\theta(a|s) R(a,s)] \end{align*} \]

We have moved the gradient operator inside the expectation. We can now estimate the gradient of the expected reward by sampling and therefore optimise our policy using gradient ascent.

This one is simple, elegant and a nice bridge between the frequentist and Bayesian worlds.

\[ \begin{align*} \theta_{\text{MAP}} &= \arg\max_{\theta} p(\theta | \mathcal{D})\\ &= \arg\max_{\theta} \frac{p(\mathcal{D} | \theta) p(\theta)}{p(\mathcal{D})} \\ &= \arg\max_{\theta} p(\mathcal{D} | \theta) p(\theta) \\ &= \arg\max_{\theta} \log p(\mathcal{D} | \theta) + \log p(\theta) \\ &= \arg\max_{\theta} \log \prod_{i=1}^N p(y_i | x_i, \theta) + \log p(\theta) \\ &= \arg\max_{\theta} \sum_{i=1}^N \log p(y_i | x_i, \theta) + \log p(\theta) \\ &= \arg\max_{\theta} \sum_{i=1}^N \log \mathcal{N}(y_i | \theta^T x_i, \sigma^2) + \log \mathcal{N}(\theta | 0, \alpha^2) \\ &= \arg\max_{\theta} \sum_{i=1}^N -\frac{1}{2\sigma^2} (y_i - \theta^T x_i)^2 - \frac{1}{2\alpha^2} \theta^T \theta \\ &= \arg\min_{\theta} \sum_{i=1}^N (y_i - \theta^T x_i)^2 + \frac{\sigma^2}{\alpha^2} \theta^T \theta \\ &= \arg\min_{\theta} \sum_{i=1}^N (y_i - \theta^T x_i)^2 + \lambda \theta^T \theta \end{align*} \]

Which is the well-known Ridge regression objective. \(\lambda = \frac{\sigma^2}{\alpha^2}\) is the regularisation parameter.