This is the second part of series of notes taken to summarize notes taken while reading the Pattern Recognition and Machine Learning (Springer, 2006) book and includes simple examples with implementations for demonstration including. As in the first post, the goal of these posts is not to be rigorous, comprehensive, or provide any new or advanced information. The target audience was first and foremost myself, but these notes are made available for anyone else who happens to find it useful, including

• Individuals curious enough to want to explore pattern recognition and machine learning, but have no previous exposure to it.
• Individuals who had a course on it and wish to refresh their memory.
• People interested in the ideas and equations, but not interested in only learning recipes and/or only learning how to use a library or toolkit.

Some contents are hidden away from default view inside modal dialogs where no page redirection takes place. These links are underlined (see e.g. here). Either clicking outside the window, clicking the close button, or clicking the back button in the browser window will close window.

Basic background in mathematics including: basic calculus, basic linear algebra, exposure to optimization, and basic statistics is assumed; see suggested readings at the end of this page. If your statistics is rusty, fear not, as the author of the PR & ML book writes:

All snippets are done using python 3 using only numpy and matplotlib add-on packages.

Evidence Approximation

Reference: Section 3.5 of PR&ML book

In the treatment of the regression problem in the first post model complexity and hyperparameters were assumed to be be known ahead of time. Extension of the approach described can be used to learn these paramters during the learning phase.

In a full Bayesian treatment of the problem, prior probabilities would be introduced over the hyperparameters. This approach will be revisited later, but for now consider the maximization of the marginalized likelihood function. Recall that the likelihood function of the data set was defined to be

$$\begin{array}{rl}p(\mathbf{t}|\mathbf{w},\beta )& =\prod _n^N\mathcal{N}(t_n|{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_n),{\beta }^{-1})\\ & =C_1\exp \{-\frac{\beta }{2}\sum _n(t_n-{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_n){)}^2\}\\ & =C_1\exp \{-\frac{\beta }{2}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}^2\},\end{array}$$

where

$$C_1={\left(\frac{\beta }{2\pi }\right)}^{N/2}$$

is the Gaussian distribution normalization constant. Also, the prior distribution introduced over model weights was given by

$$\begin{array}{rl}p(\mathbf{w},\alpha )& :=\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I})\\ & =C_2\exp \{-\frac{\alpha }{2}{\mathbf{w}}^T\mathbf{w}\},\end{array}$$

where

$$C_2={\left(\frac{\alpha }{2\pi }\right)}^{M/2}$$

is the Gaussian normalization constant for the prior distribution. Then, the maginalized likelihood function, given the training set alone, is

$$\begin{array}{rl}p(\mathbf{t}|\alpha ,\beta )& =\int p(\mathbf{t}|\mathbf{w},\beta )p(\mathbf{w},\alpha )d\mathbf{w}\\ & =C_1{C}_2\int \exp \left\{-\frac{1}{2}\left(\beta ||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}^2+\alpha {\mathbf{w}}^T\mathbf{w}\right)\right\}d\mathbf{w}.\end{array}$$

The integral can be evaluated the same way as in Part 1 of the discussion, namely expanding and completing the square, or using the general results of evaluating the joint distribution (see notes in Part 1). Consider expanding it out once again as an exercise. Expanding the exponential term excluding the preceding factor gives

$$\begin{array}{rl}\beta ||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}^2+\alpha {\mathbf{w}}^T\mathbf{w}& =\beta (\mathbf{t}-\mathbf{\Phi }\mathbf{w}{)}^T(\mathbf{t}-\mathbf{\Phi }\mathbf{w})+\alpha {\mathbf{w}}^T\mathbf{w}\\ & =\beta {\mathbf{t}}^T\mathbf{t}-\beta {\mathbf{t}}^T\mathbf{\Phi }\mathbf{w}-\beta {\mathbf{w}}^T{\mathbf{\Phi }}^T\mathbf{t}+{\mathbf{w}}^T(\underset{=:\mathbf{A}\in {\mathbb{R}}^{M\times M}}{\underbrace{\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+\alpha \mathbf{I}}})\mathbf{w}\\ & ={\mathbf{w}}^T\mathbf{A}\mathbf{w}-2\beta {\mathbf{w}}^T{\mathbf{\Phi }}^T\mathbf{t}+\beta {\mathbf{t}}^T\mathbf{t}\\ & =(\mathbf{w}-\mathbf{m}{)}^T\mathbf{A}(\mathbf{w}-\mathbf{m})-{\mathbf{m}}^T\mathbf{A}\mathbf{m}+\beta {\mathbf{t}}^T\mathbf{t}.\end{array}$$

The last line follows by completing the square, and $\mathbf{m}$ is defined as

$$\mathbf{m}:=\beta {\mathbf{A}}^{-T}{\mathbf{\Phi }}^T\mathbf{t}.$$

As a side note, with a little advance insight, the last two terms can be rewritten as (as are in the PR&ML book)

$$-{\mathbf{m}}^T\mathbf{A}\mathbf{m}+\beta {\mathbf{t}}^T\mathbf{t}=\beta ||\mathbf{t}-\mathbf{\Phi }\mathbf{m}|{|}^2+\alpha {\mathbf{m}}^T\mathbf{m}=:2E(\mathbf{m}).$$

Isolating the $\mathbf{w}$ terms for integration, it follows that the original integral can be written as

$$\begin{array}{rl}p(\mathbf{t}|\alpha ,\beta )& =C_1{C}_2\int \exp \left\{-\frac{1}{2}\left(\beta ||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}^2+\alpha {\mathbf{w}}^T\mathbf{w}\right)\right\}d\mathbf{w}\\ & =C_1{C}_2\exp \left\{-E(\mathbf{m})\right\}\int \exp \left\{-\frac{1}{2}(\mathbf{w}-\mathbf{m}{)}^T{\mathbf{A}}^{-1}(\mathbf{w}-\mathbf{m})\right\}d\mathbf{w}\\ & =C_1{C}_2\exp \left\{-E(\mathbf{m})\right\}(2\pi {)}^{M/2}|\mathbf{A}{|}^{-1/2}\\ & ={\left(\frac{\beta }{2\pi }\right)}^{N/2}{\left(\frac{\alpha }{2\pi }\right)}^{M/2}(2\pi {)}^{M/2}|\mathbf{A}{|}^{-1/2}\exp \left\{-E(\mathbf{m})\right\}\\ & ={\left(\frac{\beta }{2\pi }\right)}^{N/2}{\alpha }^{M/2}|\mathbf{A}{|}^{-1/2}\exp \left\{-E(\mathbf{m})\right\}\end{array}$$

In above, the property of the unit integral property of the multivariate Gaussian distribution is used. As before, it may be convenient to consider the log of the marginal likelihood function:

$$\log p(\mathbf{t}|\alpha ,\beta )=\frac{N}{2}\left(\log (\beta )-\log (2\pi )\right)+\frac{M}{2}\log (\alpha )-\frac{1}{2}\log (|\mathbf{A}|)-E(\mathbf{m}).$$

Note that in the above expression $\mathbf{A}=\mathbf{f}(\mathbf{x},\beta ,\alpha )$ and that varying $M$ can have an effect on the $|\mathbf{A}|$ value.

Consider now the maximization of the log-likelihood with respect to the hyperparameters $\alpha$ and $\beta$. Recall that if ${\lambda }_i$ are the eigenvalues of the symmetric matrix ${\mathbf{A}}^s$, then ${\lambda }_i+\alpha$ are the eigenvalues of ${\mathbf{A}}^s+\alpha \mathbf{I}$. The log-likelihood can thus be expressed as

$$\begin{array}{rl}\log p(\mathbf{t}|\alpha ,\beta )& =\frac{N}{2}\left(\log (\beta )-\log (2\pi )\right)+\frac{M}{2}\log (\alpha )-\frac{1}{2}\log (|\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+\alpha \mathbf{I}|)-E(\mathbf{m})\\ & =\frac{N}{2}\left(\log (\beta )-\log (2\pi )\right)+\frac{M}{2}\log (\alpha )-\frac{1}{2}\log \left\{\prod _i\{{\lambda }_i+\alpha \}\right\}-E(\mathbf{m})\\ & =\frac{N}{2}\left(\log (\beta )-\log (2\pi )\right)+\frac{M}{2}\log (\alpha )-\frac{1}{2}\sum _i\log \left\{{\lambda }_i+\alpha \right\}-E(\mathbf{m}),\end{array}$$

where ${\lambda }_i$ are the eigenvalues of $\beta {\mathbf{\Phi }}^T\mathbf{\Phi }$. Optimizing with respect to $\alpha$ requires

$$\frac{M}{2\alpha }-\frac{1}{2}\sum _i\frac{1}{{\lambda }_i+\alpha }-\frac{1}{2}{\mathbf{m}}^T\mathbf{m}=0$$

or,

$$\begin{array}{rl}\alpha {\mathbf{m}}^T\mathbf{m}& =M-\alpha \sum _i^M\frac{1}{{\lambda }_i+\alpha }\\ & =\sum _i\frac{{\lambda }_i}{{\lambda }_i+\alpha }=:\gamma ,\end{array}$$

leading to an implicit expression for $\alpha$

$$\begin{array}{r}\alpha =\frac{\gamma (\alpha )}{{\mathbf{m}}^T\mathbf{m}}.\end{array}$$ $$\begin{array}{r}\alpha =\frac{\gamma (\alpha )}{{\mathbf{m}}^T\mathbf{m}}.\end{array}$$

Consider next maximization of the log-likelihood with respect to $\beta$. First, note that if the eigenvalues of $\beta {\mathbf{A}}^s$ are ${\lambda }_i$, then the eigenvalues of ${\mathbf{A}}^s$ are ${\hat{\lambda }}_i={\beta }^{-1}{\lambda }_i$. Then,

$$\frac{\partial {\lambda }_i}{\partial \beta }={\hat{\lambda }}_i={\beta }^{-1}{\lambda }_i.$$

and the minimizer of the log-likelihood requires

$$\begin{array}{rl}\frac{\partial }{\partial \beta }\log p(\mathbf{t}|\alpha ,\beta )& =\frac{N}{2\beta }-\frac{1}{2\beta }\sum \frac{{\lambda }_i}{{\lambda }_i+\alpha }-\frac{1}{2}||\mathbf{t}-\mathbf{\Phi }\mathbf{m}|{|}^2=0.\end{array}$$

Rearranging,

$${\beta }^{-1}=\frac{1}{N-\gamma }||\mathbf{t}-\mathbf{\Phi }\mathbf{m}|{|}^2,$$

which is again an implicit expression for $\beta$.

Expectation Maximization

Reference: Section 9.3 of PR&ML book

In Bayesian treatment of the problem in the previous post, the hyperparameters in the prior distribution were assumed to be known a priori. The previous section showed an example of these hyperparameters can be learned from the training data, e.g. by maximizing the (log) likelihood function. This section looks at an alterantive approach of learning the hyperparameters including by maximizing the expected joint distribution through what is known as Expectation Maximization (EM).

Expectation maximization, is a more general concept as will be seen later, and as the name suggests, is an algorithm for finding the maximum likelihood solutions for models having latent variables given model parameters. Within the context of Bayesian treatment of the regression problem, the latent variables are the weights, and the model parameters are the hyperparameters $\alpha$ and $\beta$.

Specifically, consider the full observation dataset denoted by $\mathbf{X}$:

$$\mathbf{X}=\left(\begin{array}{c}{\mathbf{x}}_1^T\\ {\mathbf{x}}_2^T\\ ⋮\\ {\mathbf{x}}_N^T\end{array}\right)$$

where $\mathbf{x}$ denotes the n-th observation point. Similarly, latent (unobserved) values are denoted by $\mathbf{Z}$ (e.g. weights) and the model parameters denoted by $\mathbf{\theta }$ (e.g. model mean and covariance).

The log likelihood function of the observed data is expressed as marginization of the joint distribution of the observed and latent variable distribution.

$$\log p(\mathbf{X}|\mathbf{\theta })=\log \sum _{\mathbf{Z}}p(\mathbf{X},\mathbf{Z}|\mathbf{\theta })$$

The data set $\{\mathbf{X},\mathbf{Z}\}$ is referred to as the complete data set, and data set $\{\mathbf{X}\}$ as the incomplete data set.

Once again, within the context of regression, the latent variables $\mathbf{Z}$ are the weighing functions and are not known in advance. Under EM, since the complete data set is (generally) not known, the expected value of the log likelihood under the posterior distribution of the latent variable is considered. Namely, the problem statement may be expressed as

The Expectation Maximization (EM) algorithm involves staggering individual components of the log likelihood function using previous values of $\mathbf{\theta }$ which is guaranteed to increase the incomplete data log likelihood (unless it is already at a local maximum).

EM Bayesian Regression

Reference: Section 9.3.4 of PR&ML book

For the linear regression problem, following the EM algorithm above, implicit expressions can be obtained for the hyperparamters assuming previous, posterior and prior distributions.

For the linear regression problem, the latent (unobserved) variables, $\mathbf{Z}$ are the model weights $\mathbf{w}$, and model parameters $\mathbf{\theta }$ are the hyperparameters $\alpha ,\beta$. The complete (joint) data log-likelihood is given by

$$\log p(\mathbf{t},\mathbf{w}|\alpha ,\beta )=\log p(\mathbf{t}|\mathbf{w},\beta )+\log p(\mathbf{w},\alpha )$$

where, as before,

$$\begin{array}{rl}p(\mathbf{t}|\mathbf{w},\beta )& =\prod _n\mathcal{N}(t_n|{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_n),{\beta }^{-1})\\ \log p(\mathbf{t}|\mathbf{w},\beta )& =\sum _n\log \mathcal{N}(t_n|{\mathbf{w}}^T\mathbf{\varphi }({\mathbf{x}}_n),{\beta }^{-1})\\ & =\frac{N}{2}\log (\beta )-\frac{N}{2}\log (2\pi )-\beta \frac{1}{2}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2\\ p(\mathbf{w}|\alpha )& =\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I}).\end{array}$$

Recall that the posterior distribution within the Bayesian treatment was taken to be

$$p(\mathbf{w}|\mathbf{t})=\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N),$$

where

$$\begin{array}{rl}{\mathbf{S}}_N^{-1}& =\alpha \mathbf{I}+\beta {\mathbf{\Phi }}^T\mathbf{\Phi }\\ {\mathbf{m}}_N& =\beta {\mathbf{S}}_N{\mathbf{\Phi }}^T\mathbf{t}.\end{array}$$

The model parameters are evaluated by maximizing the log-likelihood of the joint distribution:

$$\begin{array}{rl}\mathbb{E}[\log p(\mathbf{t},\mathbf{w}|\alpha ,\beta )]& =\sum _{\mathbf{Z}}\log p(\mathbf{t},\mathbf{w}|\alpha ,\beta )p(\mathbf{w}|\mathbf{t},\alpha ,\beta )\\ & =\sum _{\mathbf{Z}}\log p(\mathbf{t},\mathbf{w}|\alpha ,\beta )\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N)\\ & =\frac{N}{2}\log \left(\frac{\beta }{2\pi }\right)-\frac{\beta }{2}\sum _n\mathbb{E}[(t_n-{\mathbf{w}}^T{\mathbf{\varphi }}_n{)}^2]-\frac{\alpha }{2}\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]+\frac{M}{2}\log \left(\frac{\alpha }{2\pi }\right).\end{array}$$

The expectation $\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]$ can be evaluated as

$$\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]={\mathbf{m}}_N^T{\mathbf{m}}_N+\text{tr}({\mathbf{S}}_N).$$

Optimizing with respect to $\alpha$ leads to

$$\begin{array}{rl}\alpha & =\frac{M}{\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]}\\ & =\frac{M}{{\mathbf{m}}_N^T{\mathbf{m}}_N+\text{tr}({\mathbf{S}}_N)}.\end{array}$$

Similarly, optimizing with respect to $\beta$ leads to

$$\frac{N}{\beta }-\sum _n\mathbb{E}[(t_n-{\mathbf{w}}^T{\mathbf{\varphi }}_n{)}^2]=0$$

or,

$$\beta =\frac{N}{{\sum }_n\mathbb{E}[(t_n-{\mathbf{w}}^T{\mathbf{\varphi }}_n{)}^2]}.$$

The expectation can be evaluated by using the result for the expectation $\mathbb{E}[\mathbf{w}{\mathbf{w}}^T]$ with respect to the distribution $\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N)$.

$$\begin{array}{rl}\mathbb{E}[(t_n-{\mathbf{w}}^T{\mathbf{\varphi }}_n{)}^2]& =t_n^2+({\mathbf{m}}_N{\mathbf{m}}_N^T+{\mathbf{S}}_N):{\mathbf{\varphi }}_n{\mathbf{\varphi }}_n^T-2t_n{\mathbf{m}}_N^T{\mathbf{\varphi }}_n\end{array}$$

leading to

$${\beta }^{-1}=\frac{1}{N}\sum _n\left(t_n^2+({\mathbf{m}}_N{\mathbf{m}}_N^T+{\mathbf{S}}_N):{\mathbf{\varphi }}_n{\mathbf{\varphi }}_n^T-2t_n{\mathbf{m}}_N^T{\mathbf{\varphi }}_n\right).$$

Note that the above expressions for $\alpha$ and $\beta$ are implicit as before and their solution requires some iterations.

Variational Inference

Reference: Section 10.1 of PR&ML book

Consider now a full Bayesian models where a prior distribution is introduced over all variables including the hyperparameters. Observation dataset will be denoted by $\mathbf{X}$ and latent variables (including hyperparameters which will be assumed to be stochastic) will be denoted by $\mathbf{Z}$. The probablistic model will specify the joint distribution of the observed and the latent variables, i.e. $p(\mathbf{X},\mathbf{Z})$, and the objective is to find (an approximation to) the posterior distribution $p(\mathbf{Z},\mathbf{X})$ as well as the model evidence $p(\mathbf{X})$.

It will be instructive to note that log marginal distribution $\log p(\mathbf{X})$ can be expressed as

$$\begin{array}{rl}\log p(\mathbf{X})& =\mathcal{L}(q)+\text{KL}(q||p)\\ \mathcal{L}(q)& :=\int q(\mathbf{Z})\log \left(\frac{p(\mathbf{X},\mathbf{Z})}{q(\mathbf{Z})}\right)d\mathbf{Z}\\ \text{KL}(q||p)& :=-\int q(\mathbf{Z})\log \left(\frac{p(\mathbf{Z}|\mathbf{X})}{q(\mathbf{Z})}\right)d\mathbf{Z},\end{array}$$

where $\text{KL}(q||p)\ge 0$ is the Kullback-Leibler divergence, and $\mathcal{L}(q)$ is a lower bound on $\log p(\mathbf{X})$.

A key result in variational inference results from considering the $q$ distribution to be separable as

$$q(\mathbf{Z})=\prod _i^M{q}_i({\mathbf{Z}}_i),$$

where ${\mathbf{Z}}_i$ are disjoint. The lower bound $\mathcal{L}$ under this distribution becomes

$$\begin{array}{rl}\mathcal{L}(q)& =\int \left(\prod _i{q}_i\right)\left[\log (p(\mathbf{X},\mathbf{Z}))-\sum _i\log q_i\right]d\mathbf{Z}\\ & =\int q_j\underset{=:{\mathbb{E}}_{i\ne j}[\log p(\mathbf{X},\mathbf{Z})]}{\underbrace{\left(\int \log p(\mathbf{X},\mathbf{Z})\prod _{i\ne j}{q}_{i}d{\mathbf{Z}}_i\right)}}d{\mathbf{Z}}_j-\int q_j\log q_{j}d{\mathbf{Z}}_j-\int \sum _{i\ne j}{q}_i\log q_{i}d{\mathbf{Z}}_i\\ & =\int q_j{\mathbb{E}}_{i\ne j}[\log p(\mathbf{X},\mathbf{Z})]d{\mathbf{Z}}_j-\int q_j\log q_{j}d{\mathbf{Z}}_j-\int \sum _{i\ne j}{q}_i\log q_{i}d{\mathbf{Z}}_i,\end{array}$$

where $q_i:=q({\mathbf{Z}}_i)$. From the above, it can be found that maximization of the lower bound with respect to a single factor $q_j$ requires that

The above expression is an implicit expression which shows that the log of the minimizer of the lower bound, $\mathcal{L}$ of the marginalized log distribution $\log p(\mathbf{X})$ is obtained by considering the expectation of the log of the joint distribution over all variables with respect to all factors $\{q_i\}$ for $i\ne j$.

Variational Linear Regression

Reference: Section 10.3 of PR&ML book

Reconsider the linear regression problem where the likelihood and weight priors are repeated for convenience:

$$\begin{array}{rl}p(\mathbf{t}|\mathbf{w})& =\prod _n\mathcal{N}(t_n|{\mathbf{w}}^T{\mathbf{\varphi }}_n,{\beta }^{-1})\\ p(\mathbf{w}|\alpha )& =\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I}).\end{array}$$

This time, priors are introduced over $\alpha$ and $\beta$. Since the conjugate prior for a Gaussian distribution with known mean is a Gamma distribution, $p(\alpha )$ and $p(\beta )$ are taken to be

$$\begin{array}{rl}p(\alpha )& =\text{Gam}(\alpha |a_0,b_0)\\ p(\beta )& =\text{Gam}(\beta |c_0,d_0).\end{array}$$

The marginalized distribution over the latent variables $\mathbf{Z}=(\mathbf{w},\alpha ,\beta )$ decomposed into Kullback-Leibler divergence and lower bound on

$$p(\mathbf{X})=\mathcal{L}(q)+\text{KL}(q||p)$$

with the maximum optimum value determined using

$$\begin{array}{rl}\log q_j^{\ast }& ={\mathbb{E}}_{i\ne j}[\log p(\mathbf{X},\mathbf{Z})]+\text{const}\\ & =\int \log p(\mathbf{X},\mathbf{Z})\prod _{i\ne j}{q}_{i}d{\mathbf{Z}}_i+\text{const}.\end{array}$$

More specifically, in this context, the joint distribution is

$$\begin{array}{rl}p(\mathbf{X},\mathbf{Z})& =p(\underset{\mathbf{X}}{\underbrace{\mathbf{t}}},\underset{\mathbf{Z}}{\underbrace{\mathbf{w},\alpha ,\beta }})\\ & =p(\mathbf{t}|\mathbf{w})p(\mathbf{w}|\alpha ,\beta )p(\alpha )p(\beta ).\end{array}$$

Following the discussion above, the distribution $q(\mathbf{Z})$, with abuse of notation, is assumed to be factorable as

$$\begin{array}{rl}q(\mathbf{Z})& =q(\mathbf{w},\alpha ,\beta )\\ & =q(\mathbf{w})q(\alpha )q(\beta ).\end{array}$$

The optimal individual components of the above expression are determined using the above equation:

$$\begin{array}{rl}\log q^{\ast }(\alpha )& ={\mathbb{E}}_{q_i,i\ne j}[\log p(\mathbf{X},\mathbf{Z})]+C\\ & =\int \log p(\mathbf{t},\mathbf{w},\alpha )\prod _{i\ne j}{q}_{i}d{\mathbf{Z}}_i\\ & =\int \log \left[p(\mathbf{t}|\mathbf{w},\beta )p(\mathbf{w}|\alpha ,\beta )p(\alpha )p(\beta )\right]q(\mathbf{w})q(\beta )d\mathbf{w}d\beta \\ & =\int \log p(\alpha )q(\mathbf{w})q(\beta )d\mathbf{w}d\beta +\int \log p(\mathbf{w}|\alpha )q(\mathbf{w})q(\beta )d\mathbf{w}d\beta +\int \log p(\mathbf{t}|\mathbf{w})q(\mathbf{w})q(\beta )d\mathbf{w}d\beta \\ & =\log p(\alpha )+{\mathbb{E}}_{\mathbf{w}}[\log p(\mathbf{w}|\alpha ,\beta )]+C\\ & =-\log \Gamma (a_0)+a_0\log b_0+(a_0-1)\log \alpha -b_0\alpha -\frac{M}{2}\log 2\pi +\frac{M}{2}\log \alpha -\frac{\alpha }{2}{\mathbb{E}}_{\mathbf{w}}[{\mathbf{w}}^T\mathbf{w}]+C\\ & =(a_0-1)\log \alpha -b_0\alpha +\frac{M}{2}\log \alpha -\frac{\alpha }{2}\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]+C_2\\ & =\left(a_0+\frac{M}{2}-1\right)\log \alpha -\left(b_0+\frac{1}{2}\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]\right)\alpha +C_2.\end{array}$$

In above, $C,C_2$ denote any terms that are not a function of $\alpha$. Comparing the above expression with the log of the Gamma distribution indicates that the resulting expression (as expected due to the choice of the prior) is a Gamma distribution given by

$$\begin{array}{rl}{q}^{\ast }(\alpha )& =\text{Gam}(\alpha |a_N,b_N)\\ a_N& =a_0+\frac{M}{2}\\ b_N& =b_0+\frac{1}{2}{\mathbb{E}}_{\mathbf{w}}[{\mathbf{w}}^T\mathbf{w}]\\ & =b_0+\frac{1}{2}\int {\mathbf{w}}^T\mathbf{w}q(\mathbf{w})d\mathbf{w}.\end{array}$$

Repeating the same process for $q(\mathbf{w})$:

$$\begin{array}{rl}\log q^{\ast }(\mathbf{w})& ={\mathbb{E}}_{q_i,i\ne j}[\log p(\mathbf{X},\mathbf{Z})]+C\\ & =\int \log p(\mathbf{t},\mathbf{w},\alpha ,\beta )\prod _{i\ne j}{q}_{i}d{\mathbf{Z}}_i\\ & =\int \log \left[p(\mathbf{t}|\mathbf{w})p(\mathbf{w}|\alpha ,\beta )p(\alpha )p(\beta )\right]q(\alpha )q(\beta )d\alpha d\beta \\ & =\int p(\mathbf{t}|\mathbf{w})q(\alpha )q(\beta )d\alpha d\beta +\int \log p(\mathbf{w}|\alpha )q(\alpha )q(\beta )d\alpha d\beta +\int \log p(\alpha )q(\alpha )q(\beta )d\alpha d\beta \\ & =-\frac{1}{2}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2\mathbb{E}[\beta ]-\frac{1}{2}{\mathbf{w}}^T\mathbf{w}\int \alpha q(\alpha )q(\beta )d\alpha d\beta +\int \log p(\alpha )q(\alpha )q(\beta )d\alpha d\beta \\ & =-\frac{1}{2}\mathbb{E}[\beta ]||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2-\frac{1}{2}{\mathbf{w}}^T\mathbf{w}{\mathbb{E}}_{\alpha }[\alpha ]+C\\ & =-\frac{1}{2}\left(\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+{\mathbb{E}}_{\alpha }[\alpha ]\mathbf{I}\right)+\beta {\mathbf{w}}^T{\mathbf{\Phi }}^T\mathbf{t}-\frac{\beta }{2}{\mathbf{t}}^T\mathbf{t}+C\\ & =-\frac{1}{2}{\mathbf{w}}^T\left(\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+{\mathbb{E}}_{\alpha }[\alpha ]\mathbf{I}\right)\mathbf{w}+\beta {\mathbf{w}}^T{\mathbf{\Phi }}^T\mathbf{t}+C_2\\ & =-\frac{1}{2}(\mathbf{w}-{\mathbf{m}}_N{)}^T{\mathbf{S}}_N^{-1}(\mathbf{w}-{\mathbf{m}}_N)+C_3\end{array}$$

where the $C,C_2,C_3$ are terms which are not a function of $\mathbf{w}$, and the last line follows by completing the square with

$$\begin{array}{rl}{\mathbf{S}}_N^{-1}& =\left({\mathbb{E}}_{\alpha }[\alpha ]\mathbf{I}+\mathbb{E}[\beta ]{\mathbf{\Phi }}^T\mathbf{\Phi }\right)\\ {\mathbf{m}}_N& =\mathbb{E}[\beta ]{\mathbf{S}}_N{\mathbf{\Phi }}^T\mathbf{t}.\end{array}$$

The corresponding distribution therefore is

$$q^{\ast }(\mathbf{w})=\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N).$$

Finally,

$$\begin{array}{rl}\log q^{\ast }(\beta )& =\int \log p(t|\mathbf{w},\alpha ,\beta )q(\alpha )q(\mathbf{w})d\alpha d\mathbf{w}+\int \log p(\mathbf{w}|\alpha )q(\alpha )q(\mathbf{w})d\alpha d\mathbf{w}\\ & +\int \log p(\alpha )q(\alpha )q(\mathbf{w})d\alpha d\mathbf{w}+\int \log p(\beta )q(\alpha )q(\mathbf{w})d\alpha d\mathbf{w}.\end{array}$$

Evaluating each integral is straightforward:

$$\begin{array}{rl}\int \log p(\mathbf{t}|\mathbf{w},\alpha ,\beta )q(\alpha )q(\mathbf{w})d\alpha d\mathbf{w}& =\int \sum _n^N\log \left\{\frac{\beta }{(2\pi {)}^{1/2}}\exp \left[-\frac{\beta }{2}(t_n-{\mathbf{w}}^T{\mathbf{\varphi }}_n{)}^2\right]\right\}q(\alpha )q(\mathbf{w})d\alpha d\mathbf{w}\\ & =N\log \beta -\frac{\beta }{2}{\mathbb{E}}_{\mathbf{w}}||\mathbf{t}-\mathbf{\Phi }\mathbf{w}|{|}_2^2+C\\ & =N\log \beta -\frac{\beta }{2}\left[{\mathbf{t}}^T\mathbf{t}-2{\mathbf{t}}^T\mathbf{\Phi }\mathbb{E}[\mathbf{w}]+{\mathbf{\Phi }}^T\mathbf{\Phi }:\mathbb{E}[\mathbf{w}{\mathbf{w}}^T]\right]+C.\end{array}$$

Since $p(\alpha )$ and $p(\beta )$ are both Gamma distributions, evaluation of ${\mathbb{E}}_{\mathbf{w},\alpha }[p(\beta )]$ is identical to ${\mathbb{E}}_{\mathbf{w},\beta }[p(\alpha )]$

$${\mathbb{E}}_{\mathbf{w},\alpha }[p(\beta )]=(c_0-1)\log \beta -d_0\beta +C_2.$$

Combining the terms, and ignoring the integrals resulting in terms independent of $\beta$ gives

$$\begin{array}{rl}\log q^{\ast }(\beta )& =(c_0+N-1)\log \beta -\left(d_0+{\mathbf{t}}^T\mathbf{t}-2{\mathbb{E}}_{\mathbf{w}}[{\mathbf{w}}^T]{\mathbf{\Phi }}^T\mathbf{t}+{\mathbf{\Phi }}^T\mathbf{\Phi }:{\mathbb{E}}_{\mathbf{w}}[\mathbf{w}{\mathbf{w}}^T]\right)\beta +C_3.\end{array}$$

Similar to $q^{\ast }(\alpha )$, this indicates that the optimum distribution is once again a Gamma distribution:

$$\begin{array}{rl}{q}^{\ast }(\beta )& =\text{Gam}(\alpha |c_N,d_N)\\ c_N& =c_0+N\\ d_N& =d_0+{\mathbf{t}}^T\mathbf{t}-2{\mathbb{E}}_{\mathbf{w}}[\mathbf{w}]{\mathbf{\Phi }}^T\mathbf{t}+{\mathbf{\Phi }}^T\mathbf{\Phi }:{\mathbb{E}}_{\mathbf{w}}[\mathbf{w}{\mathbf{w}}^T].\end{array}$$

With $q(\mathbf{w})$, $q(\alpha )$, and $q(\beta )$ distributions determined (Gaussian and Gamma distributions), the expected values shown in above expressions can be evaluated using standard results of Gamma distribution first moment, and previous expressions for Gaussian distributions:

$$\begin{array}{rl}{\mathbb{E}}_{\alpha }[\alpha ]& =\frac{a_N}{b_N}\\ {\mathbb{E}}_{\beta }[\beta ]& =\frac{c_N}{d_N}\\ {\mathbb{E}}_{\mathbf{w}}[{\mathbf{w}}^T\mathbf{w}]& ={\mathbf{m}}_N^T{\mathbf{m}}_N+\text{tr}({\mathbf{S}}_N).\\ {\mathbb{E}}_{\mathbf{w}}[\mathbf{w}{\mathbf{w}}^T]& ={\mathbf{m}}_N{\mathbf{m}}_N^T+{\mathbf{S}}_N.\\ {\mathbb{E}}_{\mathbf{w}}[\mathbf{w}]& ={\mathbf{m}}_N.\end{array}$$

Predictive Distribution

For completeness, the predictive distribution is considered as in Part 1. Recall that the goal is, given an observation point $x$, to make future predictions on the value of $t$. The model weights are marginalized as before

$$\begin{array}{rl}p(t|\mathbf{x},\mathbf{t},\alpha ,\beta )& =\int p(t|\mathbf{x},\mathbf{w},\beta )p(\mathbf{w}|\mathbf{t},\alpha ,\beta )d\mathbf{w}\\ & \approx \int p(t|\mathbf{x},\mathbf{w},\beta )q(\mathbf{w}|\mathbf{t},\alpha ,\beta )d\mathbf{w}\\ & =\int \mathcal{N}(t|{\mathbf{w}}^T\mathbf{\varphi }(\mathbf{x}),{\beta }^{-1})\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N)d\mathbf{w}\\ & =\mathcal{N}\left(t|{\mathbf{m}}_N^T\mathbf{\varphi }(\mathbf{x}),{\beta }^{-1}+\mathbf{\varphi }(x{)}^T{\mathbf{S}}_N\mathbf{\varphi }(\mathbf{x})\right)\end{array}$$

where the last equality follows using the same approach(es) discussed in Part 1. Note that the integration of the distributions with the priors $p(\alpha )$ and $p(\beta )$ can quickly become intractable, and thus have been left out, and are assumed to have been learned from the training set.

Summary

Exercises

• Show that $\mathbb{E}[{\mathbf{w}}^T\mathbf{w}]={\mathbf{m}}_N^T{\mathbf{m}}_N+\text{tr}({\mathbf{S}}_N)$ with respect to the distribution $\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N)$.

  Solution

• Show that $\mathbb{E}[\mathbf{w}{\mathbf{w}}^T]={\mathbf{m}}_N{\mathbf{m}}_N^T+{\mathbf{S}}_N$ with respect to the distribution $\mathcal{N}(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N)$.

  Solution

• Show that if ${\lambda }_i$ are the eigenvalues of the symmetric matrix ${\mathbf{A}}^s$, then ${\lambda }_i+\alpha$ are the eigenvalues of ${\mathbf{A}}^s+\alpha \mathbf{I}$.

  Solution

• The covariance matrix $\Sigma$ of the Gaussian distribution may be taken to be symmetric.

  Solution

• Show that

  $${\mathbf{w}}^T\mathbf{A}\mathbf{w}-2\beta {\mathbf{w}}^T{\mathbf{\Phi }}^T\mathbf{t}+\beta {\mathbf{t}}^T\mathbf{t}=(\mathbf{w}-\mathbf{m}{)}^T\mathbf{A}(\mathbf{w}-\mathbf{m})-{\mathbf{m}}^T\mathbf{A}\mathbf{m}+\beta {\mathbf{t}}^T\mathbf{t}.$$ Solution

• Show that $p(\mathbf{X})$ can be expressed as

  $$\begin{array}{rl}\log p(\mathbf{X})& =\mathcal{L}(q)+\text{KL}(q||p)\\ \mathcal{L}(q)& :=\int q(\mathbf{Z})\log \left(\frac{p(\mathbf{X},\mathbf{Z})}{q(\mathbf{Z})}\right)d\mathbf{Z}\\ \text{KL}(q||p)& :=-\int q(\mathbf{Z})\log \left(\frac{p(\mathbf{Z}|\mathbf{X})}{q(\mathbf{Z})}\right)d\mathbf{Z},\end{array}$$ Solution

• Show that the conjugate prior of the Gaussian likelihood function with known variance ${\sigma }^2$ is Gaussian.

  Solution

References

Suggested general reference books:

• Bishop, Christopher M. Pattern Recognition and Machine Learning. Springer, 2006.

  Base reference for content of this page, though certainly not the best book that I have read.

• Boyd, Stephen, and Lieven Vandenberghe. Convex Optimization. Cambridge university press, 2004.

  A classic and a must-read introductory book to (convex) optimization for everyone regardless of their background. PDF copy is also available on Standford webpage as of this writing.

• Greenberg, Michael D. Foundations of Applied Mathematics. Courier Corporation, 2013.

  Good reference book or refresher on applied math.

• Bulmer, Michael George. Principles of Statistics. Courier Corporation, 1979.

  Very easy read and short refresher on basics of statistics.