This is the first 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.

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.

This first post covers exclusively regression and its connection to machine learning (ML) concepts. This is terse and more like an engineering treatment of the topic, i.e., overly practical with possibly even more hand waving than the original book and plenty of abuse of notation.

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:

[...] probability theory can be expressed in terms of two simple equations corresponding to the sum rule and the product rule. All of the probabilistic inference and learning manipulations discussed in this book, now matter how complex, amount to repeated application of these two equations.

PR & ML book, page 359.

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

A Point of Departure

Reference: Chapter 3 of PR&ML book

A good starting point is to try to exhaust a simple regression analysis by approaching it from different viewpoints and considering variations of the same general problem. The basic regression problem reads:

The word 'best' inherently implies specification of a metric which can be used to quantitatively compare the adequacy of different polynomials. The (M-1)-th degree polynomial, referred to as the model with $M$ being the model complexity, is denoted by $y$ and may be expressed as

$$y(x)=\sum _{i=0}^{M-1}{w}_i{x}^i,$$

where $w_i$ are the polynomial coefficients (or weights) that determine the behavior of the model, and $x^i$ are the basis functions. The choice of selecting the model complexity is referred to as model selection. In future posts, determintation of the model complexity will be considered, but for the time being it is assumed that this information is known a priori. As a side note, one convenient way to think of this model is to note that it corresponds to the Taylor series expansion of $y$.

The adequacy of the model can be chosen as the "closeness" of the model to the input dataset $\mathbf{x},\mathbf{t}\in {\mathbb{R}}^N$ at given observation points (training data). A particularly convenient and common measure is the (scaled) square of the $L_2$ norm

$$E(\mathbf{w};p,t)=\frac{1}{2}||\mathbf{y}(\mathbf{x};\mathbf{w})-\mathbf{t}|{|}_2^2.$$

Here, $\mathbf{y}\in {\mathbb{R}}^N$ contains all the model prediction which where convenient may be writtten as

$$\mathbf{y}=\mathbf{\Phi }\mathbf{w},$$

where

$$\Phi (\mathbf{x})=\left(\begin{array}{c}{\mathbf{\varphi }}^T(x_1)\\ {\mathbf{\varphi }}^T(x_2)\\ ⋮\\ {\mathbf{\varphi }}^T(x_N)\end{array}\right)\in {\mathbb{R}}^{N\times M}.$$

The above equation contains all of the basis functions evaluated at every input point. As alluded to, the choice of using $L_2$ norm is rather arbitrary, but a common and convenient one with interesting interpretations and properties. Other norms can be used; see e.g. the Convex Optimization book for notes on norm choice and their implications.

The "best" fit, given the error metric, is the choice of $\mathbf{w}$ which minimizes $E(\mathbf{w};\mathbf{x},\mathbf{t})$. The original probelm can thus be re-estated as:

Within the context of "machine learning" (ML), the (M-1)-th degree polynomial is referred to as the model and the input data referred to as the training set, and within the context of optimization, $E$ is known as the objective function. The above minimziation problem is a well-known unconstrained convex optimization problem seen in optimization.

Within the context of ML, the goal of the model is to make future predictions on inputs (as opposed to, for example, look for individual approximations of the coefficients themselves). For the regression problem, this then reduces to the following general steps.

1. Choose the error function (objective function) to be minimized to be $E=c||\cdot |{|}_2^2$
2. Find model coefficients $\mathbf{w}$
3. Given an input $x$, evaluate a point-wise estimate for model $y(x;\mathbf{w})$.

Sticking with the polynomial regression problem, slightly generalizing the above steps, the general framework (in surpervised learning) in ML may be expressed as

1. Choose model and objective function, $E$
2. Training: Given training set $\{\mathbf{x},\mathbf{t}\}$, find $\{\mathbf{w}\}$ such that $E$ is minimized (or maximized depending on the objective function).
3. Prediction: Evaluate the value of the model on future input data.

For the $L_2$ minimization, a general solution can be obtained. It is simpler for these kinds of problems to briefly reverting to the index notation (see any books on tensor algebra)

$$\begin{array}{rl}\frac{\partial E}{\partial w_m}& =\frac{\partial E}{\partial w_m}||\mathbf{y}-\mathbf{t}|{|}_2^2\\ & =\frac{\partial }{\partial w_m}\frac{1}{2}{\left({\Phi }_{ij}{w}_j-t_i\right)}^2\\ & =\left({\Phi }_{ij}{w}_j-t_i\right){\Phi }_{im}.\end{array}$$

Setting above equation to $0$ and solving for $w_j$ (or $\mathbf{w}$) and reverting back:

$$\begin{array}{r}{\mathbf{w}}_{\min }=({\mathbf{\Phi }}^T\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^T\mathbf{t}.\end{array}$$

A Statistical Viewpoint

Reference: Section 1.2.5, and Section 3.3 of PR&ML book

Assume that the observed values $\mathbf{t}$ are normally distributed whose mean is the model value $y(x,\mathbf{w})$, and which has a known fixed variance ${\beta }^{-1}>0$, namely

$$p(t|x,\mathbf{w},\beta )=\mathcal{N}(t|y(x,\mathbf{w}),{\beta }^{-1}),$$

where $\mathcal{N}(t|y(x,\mathbf{w}),{\beta }^{-1})$ is the single variable Gaussian distribution. Refer to the figure below. Note the underlying assumption: the observations $y$ are normally distributed at each known observation point. For problems where there is uncertainty in both observation points $\mathbf{x}$ and measurements $\mathbf{t}$ (e.g. both are measured using noisy sensors), this assumption may not be ideal.

To introduce the concept of the likelihood function, consider the application of Bayes' theorem to determine model weights after observing the input dataset $\mathcal{D}$:

$$\begin{array}{rl}p(\mathbf{w}|\mathcal{D})& =\frac{p(\mathcal{D}|\mathbf{w})\cdot p(\mathbf{w})}{p(\mathcal{D})}\\ & \propto p(\mathcal{D}|\mathbf{w})\cdot p(\mathbf{w}).\end{array}$$

In the above, $p(\mathcal{D}|\mathbf{w})$ is referred to as the likelihood function and expresses the probability of observing dataset $\mathcal{D}$ given model parameters $\mathbf{w}$. As a side note, the denominator can be viewed as a normalization constant such that the distribution is normalized:

$$\begin{array}{rl}p(\mathcal{D})& =\int p(\mathcal{D},\mathbf{w})d\mathbf{w}\\ & =\int p(\mathcal{D}|\mathbf{w})\cdot p(\mathbf{w})d\mathbf{w}.\end{array}$$

Returning to the regression problem, if the data is assumed to be drawn independently from the previously mentioned distribution, then the likelihood function for the full training dataset is

$$p(\mathbf{t}|\mathbf{x},\mathbf{w},\beta )=\prod _{n=0}^{N-1}\mathcal{N}(t_n|y(x_n,\mathbf{w}),{\beta }^{-1}).$$

Recall that $\log$ is concave on ${\mathbb{R}}_{++}$ (this can be verified using second-order conditions), and that that the likelihood function is a positive valued function. Equivalently therefore the minimization of the negative log or maximation of the log function can be considered, which will prove to be convenient and will be a common theme: Maximize

$$\log p(t|\mathbf{x},\mathbf{w},\beta )=-\frac{\beta }{2}\sum _{n=0}^{N-1}[y(x_n,\mathbf{w})-t_n{]}^2+\frac{N}{2}\log \beta -\frac{N}{2}\log 2\pi .$$

Maximization of the above maximum likelihood function (objective function) with respect to $\mathbf{w}$ leads to:

$${\mathbf{w}}_{ML}={\left(\sum \mathbf{\varphi }(x_n)\mathbf{\varphi }(x_n{)}^T\right)}^{-1}\sum (t_n\mathbf{\varphi }(x_n))$$

where a slight generalization over the basis function notation has been introduced namely, $\mathbf{\varphi }(x_n{)}^T:=\left(\begin{array}{ccccc}1& x& x^2& \cdots & x^n\end{array}\right)$. Maximizing with respect to $\beta$ results in:

$${\beta }_{ML}^{-1}=\frac{1}{N}\sum _{n=0}^{N-1}[y(x_n,{\mathbf{w}}_{ML})-t_n{]}^2.$$

Here the subscripts $ML$ aim to denote the Maximum Liklihood and are not to be confused with Machine Learning.

With the maximum likelihood parameters determined, a predictive distribution of new values of $x$ becomes $y(x,{\mathbf{w}}_{ML})$. In other words, the mean of the distribution is precisely as those obtained by $L_2$ minimization of the error function. While the mean is the same as before, note that now a probabilistic estimate on the model is known with both mean and variance. Alternatively, $L_2$ can be thought of as maximum likelihood of the dataset where the data is assumped to be normally distributed at each observation point.

Consider now moving towards a more Bayesian approach. The general framework will be

1. Introduce a prior over model parameters (which can include hyperparameters). The prior is generally taken to be conjugate prior of the likelihood function.
2. Evaluate a maximum likelihood function given the training set
3. Evaluate (or approximate) the (maximum) posterior distribution or the marginalized distribution.

Here, the prior over $\mathbf{w}$ is introduced assumed to be normally distributed with zero mean:

$$p(\mathbf{w}|\alpha )=\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I})={\left(\frac{\alpha }{2\pi }\right)}^{(M+1)/2}\exp \left[-\frac{\alpha }{2}{\mathbf{w}}^T\mathbf{w}\right],\alpha >0.$$

Then, the posterior distribution of $\mathbf{w}$ is evaluated using Bayes' theorem

$$p(\mathbf{w}|\mathbf{x},\mathbf{t},\alpha ,\beta )=p(\mathbf{t}|\mathbf{x},\mathbf{w},\beta )\cdot p(\mathbf{w}|\alpha )\cdot \frac{1}{p(\mathbf{t})}.$$

Expression for the posterior distribution can be found using general posterior results (or alternatively as an exercise, a direct derivation is available here, which is one step over Exercise 3.8 in the PR&ML book which requires only directly utilizing the previously derived equations) to be

$$\begin{array}{r}p(\mathbf{w}|\mathbf{t},\mathbf{x},\beta )=\mathcal{N}\left(\mathbf{w}|\beta \mathbf{\Sigma }{\mathbf{\Phi }}^T\mathbf{t},\mathbf{\Sigma }\right),\end{array}$$

where

$${\mathbf{\Sigma }}^{-1}=\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+\alpha \mathbf{I}.$$

In the above, as will be seen in the Kernel Methods section, $\mathbf{\Phi }\in {\mathbb{R}}^{N\times M}$ is referred to as the design matrix and $\mathbf{\Phi }{\mathbf{\Phi }}^T\in {\mathbb{R}}^{M\times M}$ as the Gram matrix.

Alternatively, the maximum posterior (MAP) can be found by directly expanding the negative log of the above posterior distribution:

$$-\log p(\mathbf{w}|\mathbf{x},\mathbf{t},\alpha ,\beta )=\frac{\beta }{2}\sum _n(t_n-y_n{)}^2+\frac{\alpha }{2}{\mathbf{w}}^T\mathbf{w}+C$$

where $C$ represents combination of the normalization factors. The maximization of above equation with respect to $\mathbf{w}$ leads to

$$\left(\frac{\alpha }{\beta }\mathbf{I}+\sum _n\mathbf{\varphi }(x_n)\mathbf{\varphi }(x_n{)}^T\right)\mathbf{w}=\sum _n{t}_n\mathbf{\varphi }(x_n)$$

where $\frac{\alpha }{\beta }\mathbf{I}$ is positive definite which naturally appears and plays the role of a regularization parameter to avoid overfitting the data. Note that this is precisely the same result obtained for the mean of the weights earlier. This type of regularization could and is sometimes included in an adhoc manner in the previously encountered regression problem.

The coefficients maximizing this objective function, denoted ${\mathbf{w}}_{ML}$ as before, can be used to evaluate the mean of the probability distribution.

Bayesian Regression

Reference: Section 1.2.6, 3.3.2 of PR&ML book

As pointed out in previous example, the posterior distribution is not directly used. The dependence on the model weights $\mathbf{w}$ can be removed if it is integrated out by evaluating the marginal distribution.

$$\begin{array}{rl}p(t_{n+1}|\mathbf{t},\alpha ,\beta )& =\int \mathcal{N}(t_{n+1}|y(x,\mathbf{w}),{\beta }^{-1})\mathcal{N}(\mathbf{w}|\mathbf{t},\alpha ,\beta )d\mathbf{w}\\ & =\int \mathcal{N}(t_{n+1}|\mathbf{\varphi }(x{)}^T\mathbf{w},{\beta }^{-1})\mathcal{N}(\mathbf{w}|{\mathbf{\mu }}_w,{\mathbf{S}}_N)d\mathbf{w}\\ & =\mathcal{N}(t_{n+1}|{\mathbf{\mu }}_w^T\mathbf{\varphi },{\sigma }_N^2),\end{array}$$

where,

$${\sigma }_N^2={\beta }^{-1}+{\mathbf{\varphi }}^T\mathbf{S}\mathbf{\varphi }.$$

The last equation for the marginal distribution follows by either expanding the integrand and separating and integrating out the $\mathbf{w}$ components which turn out to be Gaussian, or directly using the general results which is accomplished using the general approach.

One of the caveats considered thus far has been introduction of fixed hyperparameters $\alpha ,\beta$. A full Bayesian treatment would consider introduction of prior distributions over these parameters as well as the model weights. Doing so however will no longer result in anlytical solutions for the posterior and predictive distributions, and instead approximate solutions to the corresonding distributions must be saught. More on this later.

Kernel Methods

Reference: 3.3.3, Chapter 6 of PR&ML book

In preceding section the primary focus was on restricting the basis functions to polynomial basis functions.

Consider the MAP solution considered before, i.e. the minimization of the objective function

$$f(\mathbf{w}):={\mathbf{w}}^T\mathbf{w}+\frac{\beta }{\alpha }{\left(\mathbf{\Phi }\mathbf{w}-\mathbf{t}\right)}^T\left(\mathbf{\Phi }\mathbf{w}-\mathbf{t}\right),$$

whose minimizer was found to be

$${\mathbf{w}}_{ML}=-\frac{\beta }{\alpha }{\mathbf{\Phi }}^T\left(\mathbf{\Phi }\mathbf{w}-\mathbf{t}\right).$$

Now, define $\mathbf{a}$ as

$$\mathbf{a}:=-\frac{\beta }{\alpha }(\mathbf{\Phi }\mathbf{w}-\mathbf{t})$$

and note that, although not explicitely specified, $\mathbf{a}$ itself is a function of $\mathbf{w}$. The minimzier can be thus be implicitely expressed as

$$\mathbf{w}={\mathbf{\Phi }}^T\mathbf{a}$$

and the original objective function can be re-written as

$$\begin{array}{rl}f(\mathbf{w})& ={\left(\mathbf{\Phi }\mathbf{w}-\mathbf{t}\right)}^T\left(\mathbf{\Phi }\mathbf{w}-\mathbf{t}\right)+\frac{\alpha }{\beta }{\mathbf{w}}^T\mathbf{w}\\ & ={\left(\mathbf{\Phi }{\mathbf{\Phi }}^T\mathbf{a}-\mathbf{t}\right)}^T\left(\mathbf{\Phi }{\mathbf{\Phi }}^T\mathbf{a}-\mathbf{t}\right)+\frac{\alpha }{\beta }{\mathbf{a}}^T\mathbf{\Phi }{\mathbf{\Phi }}^T\mathbf{a}\\ & ={\mathbf{a}}^T\mathbf{K}\mathbf{K}\mathbf{a}-{\mathbf{a}}^T\mathbf{K}\mathbf{t}-{\mathbf{t}}^T\mathbf{K}\mathbf{a}+{\mathbf{t}}^T\mathbf{t}+\frac{\alpha }{\beta }\mathbf{\alpha }\mathbf{K}{\mathbf{a}}^T,\end{array}$$

where $\mathbf{\Phi }$ is known as the design matrix, and $\mathbf{K}:=\mathbf{\Phi }{\mathbf{\Phi }}^T\in {\mathbb{R}}^{N\times N}$ is symmetric and is known as the Gram matrix. Note that the individual components of the Gram matrix can be written as

$$K_{nm}=\mathbf{\varphi }(x_n{)}^T\mathbf{\varphi }(x_m)=:k(x_n,x_m),$$

where $x_n$ is the $n$th observation point, and $k(x,x^{\prime })$ is referred to as the kernel function.

Minimization of the objective function in terms of $\mathbf{a}$ leads to

$$\underset{\mathbf{a}}{\mathrm{argmin}}f(\mathbf{a})=:{\mathbf{a}}_{min}={\left(\mathbf{K}+\frac{\alpha }{\beta }\mathbf{I}\right)}^{-1}\mathbf{t},$$

and the target values can therefore be written as

$$\begin{array}{rl}y(x)& ={\mathbf{w}}^T\mathbf{\varphi }\\ & ={\mathbf{a}}_{min}^T\mathbf{\Phi }\mathbf{\varphi }\\ & ={\mathbf{t}}^T{\left(\mathbf{K}+\frac{\alpha }{\beta }\mathbf{I}\right)}^{-1}\mathbf{\Phi }\mathbf{\varphi }\\ & ={\mathbf{t}}^T{\left(\mathbf{K}+\frac{\alpha }{\beta }\mathbf{I}\right)}^{-1}\mathbf{k}(x)\\ & ={\mathbf{k}}^T(x){\left(\mathbf{K}+\frac{\alpha }{\beta }\mathbf{I}\right)}^{-1}\mathbf{t}\\ & =\sum _{n}k(x,x_n){\hat{t}}_n,\end{array}$$

where $\mathbf{k}\in {\mathbb{R}}^N$ with components

$$\begin{array}{rl}\mathbf{k}(x):=\mathbf{\Phi }\mathbf{\varphi }& =\left(k(x_1,x),k(x_2,x),\cdots ,k(x_N,x)\right).\end{array}$$

Side note: this form of the solution is not new and was already encountered in a previous example where the model weight prior was assumed to be Gaussian and zero mean and isotropic covariance. Repeating the results from the example:

$$\begin{array}{rl}y(x)& ={\mathbf{m}}_N^T\mathbf{\varphi }(x)\\ & ={\mathbf{\varphi }}^T(x)\left(\beta {\mathbf{S}}_N{\mathbf{\Phi }}^T\mathbf{t}\right)\\ & =\beta {\mathbf{\varphi }}^T{\mathbf{S}}_N{\mathbf{\Phi }}^T\mathbf{t}\\ & =\sum _n\beta {\mathbf{\varphi }}^T{\mathbf{S}}_N\mathbf{\varphi }(x_n)t_n\\ & =\sum _n\hat{k}(x,x_n)t_n.\end{array}$$

The second to last line follows since from linear algebra it is know that a matrix-vector multiplication can is linear combination of the columns of the matrix scaled by each row of the multiplying vector.

Gaussian Processes

A random process $y(x)$ is a Gaussian process if it is defined as a probability distribution over $y(x)$ such that the its values evaluated at arbitrary points jointly have a Gaussian distribution. In such processes, the joint distribution is over $N$ variables is completely specified using the mean and the covariance.

It may be convenient to assume that the random processes have zero mean, and hence being able to completely specify them through the covariance, $k(x,x^{\prime })$. Now let $x,x^{\prime }\in [a,b]$. Then a zero-mean Gaussian process $y$ can be expressed as

$$p(y)=N(x|0,k(x,x^{\prime })).$$

We can consider kernel functions directly, e.g.

$$\begin{array}{rl}{k}_{\text{gauss}}(x,x^{\prime })& =\exp \{-\frac{1}{2}{\sigma }^{-2}||x-x^{\prime }|{|}^2\}\\ k_{\text{exp}}(x,x^{\prime })& =\exp \{-\theta |x-x^{\prime }|\}\end{array}$$

The two figures below show sampling 5 Gaussian zero-mean processes with Gaussian and Exponential kernels (covariancces) over a prespecified range.

In order to relate this back to the regression model previously encountered, consider revisiting the linear regression model where the target values are expressed as linear combination of weights and basis functions and re-consider the zero mean isotropic Gaussian prior introduced over the weights:

$$\begin{array}{rl}\mathbf{y}(x)& =\mathbf{\Phi }\mathbf{w}\\ p(\mathbf{w}|\alpha )& =\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I}).\end{array}$$

Considering the mean and covariance of the target values, we have

$$\begin{array}{rl}\mathbb{E}[\mathbf{y}]& =\mathbf{\Phi }\mathbb{E}[\mathbf{w}]=0\\ \text{cov}[\mathbf{y}]& =\mathbb{E}[\mathbf{y}{\mathbf{y}}^{\mathbb{T}}]\\ & =\mathbf{\Phi }\mathbb{E}[\mathbf{w}{\mathbf{w}}^T]{\mathbf{\Phi }}^T\\ & =\mathbf{\Phi }\int \mathbf{w}{\mathbf{w}}^T\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I})d\mathbf{w}{\mathbf{\Phi }}^T\\ & ={\alpha }^{-1}\mathbf{\Phi }{\mathbf{\Phi }}^T=:\hat{\mathbf{K}}\end{array}$$

Last line of equality follows a simple integration by parts if the unit integral property of the Gaussian is used. Also this line defines the kernel function for a linear regression model basis with zero mean isotropic Gaussian prior defined over the weights.

If the target value distribution conditioned on the model values is assumed to be isotropic Gaussian, i.e.

$$p(\mathbf{t}|\mathbf{y})=\mathcal{N}(\mathbf{t}|\mathbf{y},{\beta }^{-1}\mathbf{I})$$

While restricting the model being a zero-mean Gaussian process namely

$$p(\mathbf{y})=\mathcal{N}(\mathbf{y}|0,\mathbf{K}),$$

then the marginalized distribution of the model trained on the data set can be written as

$$\begin{array}{rl}p(\mathbf{t})& =\int p(\mathbf{t}|\mathbf{y})p(\mathbf{y})d\mathbf{y}\\ & =\mathcal{N}(\mathbf{t}|0,\mathbf{C}),\end{array}$$

where

$$C(x_n,x_m)=k(x_n,x_m)+{\beta }^{-1}{\delta }_{nm},$$

and ${\delta }_{nm}$ is the Kronecker delta, $k$ is a kernel function, and the last line follows as before by either using the generalized results for predictive distributions or following process similar to the marginal distribution evaluation previously mentioned.

Using the trained results to evalute predictive values for new observations is straightforward. Consider the joint distribution

$$p(t_{N+1},{\mathbf{t}}_N)=\mathcal{N}((t_{N+1},{\mathbf{t}}_N)|0,\left(\begin{array}{cc}{\mathbf{C}}_N& \mathbf{k}\\ {\mathbf{k}}^T& c\end{array}\right)),$$

where $\mathbf{k}=(k(x_1,x_{N+1}),k(x_2,x_{N+1}),\cdots )$ and $c=k({\mathbf{x}}_{N+1},{\mathbf{x}}_{N+1})+{\beta }^{-1}$. Then the conditional distribution distribution $p(t_{n+1}|{\mathbf{t}}_N)$ can be found using similar approaches as before to be

$$p(t_{n+1}|{\mathbf{t}}_N)=\mathcal{N}(t_{n+1}|{\mathbf{k}}^T{\mathbf{C}}_N^{-1}\mathbf{t},c-{\mathbf{k}}^T{\mathbf{C}}_N^{-1}\mathbf{k}).$$

Note that as in the previous example, generally there is no distinguishable training phase which can be separated from the prediction phase. Specifically, the $\mathbf{k}$ vector must be constructed using the training data at each evaluation point. If the basis functions are limited to a finite set, then the original equations encountered for Bayesian treatment of the regression problem is recovered. Additionally, note that ${\mathbf{C}}_N^{-1}\in R^{N\times N}$ involves inversion of an $N\times N$ matrix which computationally grows in a cubic fasion with training set size. The advantage of the Gaussian process viewpoint however is that covariance functions which can only be expressed in terms of infinite dimensional basis functions.

Neural Networks

Reference: Chapter 5 of PR&ML book

The basic neural network can be described as a series of functional transformations. The word network refers to the existence of a (nonlinear) coupling of the weights between basis functions.

The approximation at an input point $x$ is expressed in terms of the model's network coefficients at different "layers". For a two-layer network, following the regression theme, the model is given by

$$\begin{array}{rl}{a}_j^{(1)}& :=w_{j1}^{(1)}x+w_{j0}^{(1)}j=1,2,\cdots ,M\\ z_j& :=h(a_j^{(1)})\\ a^{(2)}& :=\sum _{j=1}^M{w}_j^{(2)}{z}_j+w_0^{(2)}\\ y(x)& =\sigma (a^{(2)})\\ & =\sigma \left(\sum _{j=0}^{M-1}{w}_j^{(2)}h(w_{j1}^{(1)}x+w_{j0}^{(1)})+w_{j0}^{(2)}\right)\end{array}$$

where ${\mathbf{w}}_1^{(1)},{\mathbf{w}}_0^{(1)}\in {\mathbb{R}}^M$, and $\mathbf{a}\in {\mathbb{R}}^M$. The superscripts ${}^{(1)}$ and ${}^{(2)}$ denotes the first and second layers number of the neural network, respectively, and $M\ge 1$ is the model complexity. Following the earlier notation, $w_j^{(1)}$ are referred to as the weights and $w_{j0}^{(1)}$ referred to as the model biases (of the first layer). Additionally $a_j$ are referred to as the activations and $h$ is called activation function.

The above expressions can be made less intimating by rewriting them as

$$\begin{array}{rl}{\mathbf{a}}^{(1)}& ={\mathbf{W}}^{(1)}\mathbf{\varphi }(\mathbf{x})\\ {\mathbf{W}}^{(1)}& =\left(\begin{array}{ccc}& w_{10}& w_{11}\\ & w_{20}& w_{21}\\ & ⋮& ⋮\\ & w_{M0}& w_{M1}\end{array}\right)\in R^{M\times 2}\\ {\mathbf{\varphi }}^T(x)& =\left(\begin{array}{cc}1& x\end{array}\right)\\ {\hat{\mathbf{z}}}^T& =\left(\begin{array}{ccccc}1& h(a_1)& h(a_2)& \cdots & h(a_M)\end{array}\right)\in {\mathbb{R}}^{M+1}\\ a^{(2)}& ={\hat{\mathbf{z}}}^T{\mathbf{w}}^{(2)}\\ y(x)& =\sigma (a^{(2)}).\end{array}$$

Neural network models therefore describe nonlinear models and will generally require iterative solutions.

Recall that model training in the case of least squares was an unconstrained minimization of the $L_2$ norm of the error function. This process is repeated for the model training of neural networks, i.e. minimizer of

$$\begin{array}{r}E(\mathbf{x}):=\frac{1}{2}||\mathbf{y}(\mathbf{x})-\mathbf{t}|{|}_2^2\end{array}$$

is saught, where $\mathbf{y}$ denotes the model values at each observation point given by above expressions, and $\mathbf{t}$ are the observation values. The key difference here is that the model is nonlinear and the minimization of the objective function would need to iterated.

Summary

Key findings for each approach is tersely summarized in the table below meant to be used only as a refresher of the overall approach.

Approach	Model | Likelihood	Objective | Prior | Conditional	Post. Mean	Post. Var
Regression	$\mathbf{y}=\mathbf{\Phi }\mathbf{w}$	$1/2||\mathbf{y}-\mathbf{t}|{|}_2^2$	$({\mathbf{\Phi }}^T\mathbf{\Phi }{)}^{-1}{\mathbf{\Phi }}^T\mathbf{t}$	-
MAP	$$\prod _{n=0}^{N-1}\mathcal{N}(t_n|y(x_n,\mathbf{w}),{\beta }^{-1}).$$	$\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I})$	${\left(\frac{\alpha }{\beta }\mathbf{I}+{\mathbf{\Phi }}^T\mathbf{\Phi }\right)}^{-1}{\mathbf{\Phi }}^T\mathbf{t}$	-
Bayesian Regression	$\mathcal{N}(\mathbf{t}|\mathbf{y}={\mathbf{\varphi }}^T\mathbf{w},{\beta }^{-1}\mathbf{I})$	$\mathcal{N}(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I})$	$\mathbb{E}[p(t|\cdots )]=\beta \mathbf{\Sigma }{\mathbf{\Phi }}^T\mathbf{t}$	${\mathbf{\Sigma }}^{-1}=\beta {\mathbf{\Phi }}^T\mathbf{\Phi }+\alpha \mathbf{I}$
Nadaray-Watson (Using Gaussian kernel)	$y(x)=\frac{{\sum }_{n}g(x-x_n)t_n}{{\sum }_{n}g(x-x_n)}$	-	$\frac{{\sum }_n\exp \left\{-\frac{1}{2{\sigma }^2}(x-x_n{)}^2\right\}{t}_n}{{\sum }_m\exp \left\{-\frac{1}{2{\sigma }^2}(x-x_m{)}^2\right\}}$	${\sigma }^2-y^2$
Gaussian Process	$p(y)=\mathcal{N}(\mathbf{y}|0,\mathbf{K})$	$p(\mathbf{t}|\mathbf{y})=\mathcal{N}(\mathbf{t}|\mathbf{y},{\beta }^{-1}\mathbf{I})$	$m(x_i)=\mathbf{k}(\mathbf{x},x_i{)}^T\mathbf{a}$	${\sigma }^2(x_i)=c-\mathbf{k}(\mathbf{x},x_i{)}^T{\mathbf{C}}_N^{-1}\mathbf{k}(\mathbf{x},x_i)$ $c=k(x_{N+1},x_{N+1})+{\beta }^{-1}$
Neural Network	$\sigma \left({\sum }_{j=0}^{M-1}{w}_j^{(2)}h(w_{j1}^{(1)}x+w_{j0}^{(1)})+w_{j0}^{(2)}\right)$	$1/2||\mathbf{y}-\mathbf{t}|{|}_2^2$	Requires iterations	-

There are various other topics like sparse kernel machines, approximate inference, etc. that will be addressed in future posts.

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.