This is the third 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 and second posts, 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. Also, ffmpeg package needs to be installed to create the videos.

Linear Discriminants

Reference: Section 4.1.5

Part 1 and Part 2 of ML&PR notes considered supervised learning applied to regression problems. This post takes a step in a different direction and considers the application of supervised learning to classification problems. Specifically, consider an input training dataset $\mathbf{x}$ where each data point is assigned to known target classes. Then the classification problem can be stated as follows.

As in Part 1 this problem can be approached by optimizing an objective function. Consider first the restriction to a simpler problem, namely the linear classification model for a two-class problem ($K=2$. The hyperplane (decision surface) defined by

$$y(\mathbf{x})={\mathbf{w}}^T\mathbf{x}+w_0$$

may be used to implicitely subdivide the the two halfspaces, where each halfspace identifies a separate class. Specifically, referring to the figure below, note that $||\mathbf{w}|{|}^{-1}y(x)$ provides a signed measure of the distance from the decision surface. The two classes can thus be implicitely defined through the relationship

$$\begin{array}{rl}{C}_1& =\{\mathbf{x}|y(\mathbf{x})\ge 0\}\\ C_2& =\{\mathbf{x}|y(\mathbf{x})<0\},\end{array}$$

while the decision surface corresponds to $y=0$.

The orientation of the supporting halfplane may be determined by considering the minimization of an appropriate objective function as was considered in Part 1 and Part 2 of the posts covering treatment of regression problems. Specifically, consider the objective function defined by

$$\begin{array}{rl}E& :=\frac{1}{2}||\mathbf{y}(\mathbf{x})-\mathbf{t}|{|}_2^2\\ y_n({\mathbf{x}}_n)& :={\mathbf{w}}^T{\mathbf{x}}_n+w_0\\ t_n& :=\{\begin{array}{ll}+\frac{N}{N_1}& {\mathbf{x}}_n\in {\mathcal{C}}_1\\ -\frac{N}{N_2}& {\mathbf{x}}_n\in {\mathcal{C}}_2\end{array}\end{array}$$

In above $t_n$ denotes the target class of point $n$. Note that under this choice, the sum of the target value will be zero when all points are classified correctly:

$$\begin{array}{rl}\sum _n{t}_n& =\sum _{n\in {\mathcal{C}}_1}{t}_n+\sum _{n\in {\mathcal{C}}_2}{t}_n\\ & =\frac{N}{N_1}{N}_1-\frac{N}{N_2}{N}_2\\ & =0.\end{array}$$

Consider now the minimization of the objective function: the minimization of the objective function due to $w_0$ requires that

$$\sum _n\left({\mathbf{w}}^T\mathbf{x}+w_0\right)=\sum _n{t}_n=0$$

or

$$\begin{array}{rl}{w}_0& ={\mathbf{w}}^T\underset{:=\mathbf{m}}{\underbrace{\frac{1}{N}\sum _n{x}_n}}\\ & ={\mathbf{w}}^T\mathbf{m},\end{array}$$

where $\mathbf{m}$ is understood to be the mean of the total dataset. If the individual class means are defined to be

$${\mathbf{m}}_k=\frac{1}{N_k}\sum _{n\in {\mathcal{C}}_k}{\mathbf{x}}_n,k\in \{1,2\}$$

then

$$\begin{array}{rl}\mathbf{m}& =\frac{1}{N}\sum _n{\mathbf{x}}_n\\ & =\sum _k\sum _{n\in {\mathcal{C}}_k}{\mathbf{x}}_n\\ & =\frac{1}{N}\left(N_1{\mathbf{m}}_1+N_2{\mathbf{m}}_2\right).\end{array}$$

Similarly, minimization with respect to $\mathbf{w}$ requires that

$$\sum _n\left({\mathbf{w}}^T{\mathbf{x}}_n-\underset{w_0}{\underbrace{{\mathbf{w}}^T\mathbf{m}}}-t_n\right){\mathbf{x}}_n=0.$$

or equivalently

$$\sum _n\left({\mathbf{x}}_n{\mathbf{x}}_n^T-{\mathbf{x}}_n\mathbf{m}\right)\mathbf{w}=\sum _n{t}_n{\mathbf{x}}_n.$$

After some work, the above equation can be re-expressed as

$$\left({\mathbf{S}}_W+\frac{N_1{N}_2}{N}{\mathbf{S}}_B\right)\mathbf{w}=N\cdot ({\mathbf{m}}_1-{\mathbf{m}}_2)$$

where,

$$\begin{array}{rl}{\mathbf{S}}_W& :=\sum _{n\in {\mathcal{C}}_k,k\in \{1,2\}}({\mathbf{x}}_n-{\mathbf{m}}_k)({\mathbf{x}}_n-{\mathbf{m}}_k{)}^T\\ {\mathbf{S}}_B& :=({\mathbf{m}}_2-{\mathbf{m}}_1)({\mathbf{m}}_2-{\mathbf{m}}_1{)}^T.\end{array}$$

It is worth noting that ${\mathbf{S}}_b\mathbf{w}$ is in the direction of ${\mathbf{m}}_2-{\mathbf{m}}_1$ and thus we re-express previous equations as

$$\mathbf{w}\propto {\mathbf{S}}_W^{-1}({\mathbf{m}}_1-{\mathbf{m}}_2).$$

The orientation of the supporting hyperplane can therefore be explicitely learned from the data. Specifically, note that the previously described Training/Learning phase followed by Prediction phase is still applicable:

1. Training: Given input training set $\mathbf{X}\in {\mathbb{R}}^{NtimesD}$ and corresponding target classes ${\mathcal{C}}_k$, determine the direction of the hyperplane (model weights) through

   $$\begin{array}{rl}\mathbf{w}& \propto {\mathbf{S}}_W^{-1}({\mathbf{m}}_1-{\mathbf{m}}_2).\\ w_0& ={\mathbf{w}}^T\mathbf{m}\end{array}$$

   Note that there are no iterations involved in the training phase. Also, note that the in-class means can be evaluated sequentially, and that ${\mathbf{S}}_W\in {\mathbb{R}}^{D\times D}$.

2. Prediction: For a new observation $x$, compute $y$

   $$y(\mathbf{x})={\mathbf{w}}^T\mathbf{x}+w_0$$

   and assign a class based its value.

Though may not be immediately obvious, it will be shown below that this result corresponds to the the maximization of the separation of the projected class means ${\mathbf{m}}_k$ onto $\mathbf{w}$ while the in-class variance $s_k^2$ is minimized.

Fisher's Discriminant

Consider now an alternative viewpoint of the same problem, namely the Fisher's linear discriminant:

The numerator of the Fisher's criteria provides a measure of the (projected) class mean separation, while the denominator provides a measure of the (summed) in-class variance (see figure below). Maximization of Fisher's criteria thus maximizes the separation of the projected class means ${\mathbf{m}}_k$ onto $\mathbf{w}$ while the in-class variance $s_k^2$ is minimized.

Note that the objective function is to be optimized in terms of the model weights $\mathbf{w}$. Using the definitions above, the objective function can be explicitely written in terms of the model weights as

$$J(\mathbf{w})=\frac{\mathbf{w}{\mathbf{S}}_B\mathbf{w}}{\mathbf{w}{\mathbf{S}}_w\mathbf{w}}$$

where

$$\begin{array}{rl}{\mathbf{S}}_B& :=({\mathbf{m}}_2-{\mathbf{m}}_1)({\mathbf{m}}_2-{\mathbf{m}}_1{)}^T\\ {\mathbf{S}}_W& :=\sum _{n\in {\mathcal{C}}_1}({\mathbf{x}}_n-{\mathbf{m}}_1)({\mathbf{x}}_n-{\mathbf{m}}_1{)}^T+\sum _{n\in {\mathcal{C}}_2}({\mathbf{x}}_n-{\mathbf{m}}_2)({\mathbf{x}}_n-{\mathbf{m}}_2{)}^T.\end{array}$$

The maximization of the above objective function is straightforward. Taking the derivative with respect to $\mathbf{w}$ and rearranging results in

$$D_{\mathbf{w}}J(\mathbf{w})=\frac{{\mathbf{S}}_B\mathbf{w}}{{\mathbf{w}}^{\mathbf{T}}{\mathbf{S}}_W\mathbf{w}}-\frac{{\mathbf{w}}^T{\mathbf{S}}_B\mathbf{w}}{({\mathbf{w}}^{\mathbf{T}}{\mathbf{S}}_W\mathbf{w}{)}^2}{\mathbf{S}}_W\mathbf{w}=0$$

or

$${\mathbf{S}}_W\mathbf{w}=\frac{({\mathbf{w}}^T{\mathbf{S}}_W\mathbf{w})}{({\mathbf{w}}^T{\mathbf{S}}_B\mathbf{w})}{\mathbf{S}}_B\mathbf{w}.$$

Using the definiton of ${\mathbf{S}}_B$, note that ${\mathbf{S}}_B\mathbf{w}$ scales in the direction of $\alpha \cdot ({\mathbf{m}}_w-{\mathbf{m}}_1)$. Also, the leading coefficient term in the above equation dependent of $\mathbf{w}$ can be neglected since the objective is to find the orientation of $\mathbf{w}$. Hence, we can find an expression for the model weights as

$$\begin{array}{rl}\mathbf{w}& =\alpha {\mathbf{S}}_W^{-1}{\mathbf{S}}_B\mathbf{w}\\ & =\hat{\alpha }{\mathbf{S}}_W^{-1}({\mathbf{m}}_2-{\mathbf{m}}_1)\end{array}$$

where $\alpha$ can be treated as a scaling constant. Note that, as alluded to earlier, this result is the same as the one obtained for least squares solution considered in the preceeding section.

Probabilistic Generative Models

Reference: Section 4.2

Consider now taking a probabilistic approach to classification such that it is generalizable to more than two classes. Assume a joint distribution of classification and datapoint given by $p(\mathbf{x},{\mathcal{C}}_k)$. Assuming prior class distribution $p({\mathcal{C}}_k)$, the posterior distribution can be expressed as

$$p({\mathcal{C}}_1|\mathbf{x})=\frac{p(\mathbf{x}|{\mathcal{C}}_k)p({\mathcal{C}}_k)}{{\sum }_{k}p(\mathbf{x}|{\mathcal{C}}_k)p({\mathcal{C}}_k)}.$$

By using

$$a_k=\log \left(p(\mathbf{x}|{\mathcal{C}}_k)p({\mathcal{C}}_k)\right),$$

the above expression can be expressed as

$$p({\mathcal{C}}_k|\mathbf{x})=\frac{\exp \{a_k\}}{{\sum }_m\exp \{a_m\}}.$$

Now consider assuming that the class conditional distributions are Gaussian with shared covariance given by

$$p(\mathbf{x}|{\mathcal{C}}_k)=\mathcal{N}(\mathbf{x}|{\mathbf{\mu }}_k,\mathbf{\Sigma }),$$

and consider finding the paramters by maximum likelihood solution. Let $t_{nk}$ denote the target class for the $n$-th point where ${\mathbf{t}}_n=(t_1,t_2,\cdots ,t_K)$ is zero except for the correct class. Also, let the prior class probablities $p({\mathcal{C}}_k)={\pi }_k$. The likelihood function is given by

$$p(\mathbf{t},\mathbf{x}|{\pi }_k,{\mathbf{\mu }}_k,\mathbf{\Sigma })=\prod _n^N\prod _k^K{\pi }_k\mathcal{N}({\mathbf{x}}_n|{\mathbf{\mu }}_k,\mathbf{\Sigma }{)}^{t_{nk}}$$

As before, consider the maximization of the log likelihood. Before doing that, recall that prior probabilities must sum to one which can be included into the likelihood function as a equality constraint through introduction of a Lagrange multiplier $\lambda$

$$\log p(\mathbf{t},\mathbf{x}|{\pi }_k,{\mathbf{\mu }}_k,\mathbf{\Sigma })=\sum _n^N\sum _k^K{t}_{nk}\left(\log ({\pi }_k)+\log \mathcal{N}({\mathbf{x}}_n|{\mathbf{\mu }}_k,\mathbf{\Sigma })\right)+\lambda \left(\sum _k\left({\pi }_k\right)-1\right).$$

The priors can then be found to be ${\pi }_k=\frac{N_k}{N}$ by maximizating the constrained log likelihood with respect to ${\pi }_k$.

Consider now the maximzation of the log likelihood with respect to the model means:

$$\begin{array}{rl}\frac{\partial E}{\partial {\mathbf{\mu }}_k}& =\frac{\partial }{\partial {\mathbf{\mu }}_k}\sum _n\sum _k-\frac{1}{2}{t}_{nk}\left(-{\mathbf{x}}_n^T{\mathbf{\Sigma }}^{-1}{\mathbf{\mu }}_k-{\mathbf{\mu }}_k^T{\mathbf{\Sigma }}^{-1}{\mathbf{x}}_n+{\mathbf{\mu }}_k^T{\mathbf{\Sigma }}^{-1}{\mathbf{x}}_n+{\mathbf{\mu }}_k^T{\mathbf{\Sigma }}^{-1}{\mathbf{\mu }}_k\right)\\ & =\sum _n\sum _k-\frac{1}{2}{t}_{nk}\left(-{\mathbf{x}}_n^T{\mathbf{\Sigma }}^{-1}+{\mathbf{\mu }}_k^T{\mathbf{\Sigma }}^{-1}\right)=0.\end{array}$$

By rearranging the above expression, an explicit expression for the model mean is found to be

$${\mathbf{\mu }}_k^T=\frac{1}{N_k}\sum _k{t}_{nk}{\mathbf{x}}_n^T.$$

Finally, consider the maximization of the objective function with respect to the covariance matrix. Before proceeding, it's worth mentioning two identities:

1. $$\begin{array}{rl}\frac{\partial }{\partial \mathbf{A}}\left(\mathbf{A}{\mathbf{A}}^{-1}\right)& =0\\ & =\mathbf{A}\frac{\partial {\mathbf{A}}^{-1}}{\partial \mathbf{A}}+{\mathbf{A}}^{-1}\end{array}$$

   or

   $$\frac{\partial {\mathbf{A}}^{-1}}{\partial \mathbf{A}}=-{\mathbf{A}}^{-2}$$
2. $$\frac{\partial }{\partial \mathbf{A}}\log |\mathbf{A}|=\frac{1}{|\mathbf{A}|}|\mathbf{A}|{\mathbf{A}}^{-1}$$

   See Convex Optimization section A4.1 (or a book on tensor algebra) for proof.

The full log likelihood including all terms including the covariance matrix is

$$p:=\sum _n\sum _k-\frac{1}{2}\log |\mathbf{\Sigma }|-\frac{1}{2}(\mathbf{x}-{\mathbf{\mu }}_k{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{x}-{\mathbf{\mu }}_k)+C$$

Then the gradient with respect to the covariance matrix, utilizing the identities listed above, becomes

$$\sum _n\sum _k\left[{\mathbf{\Sigma }}^{-1}-(\mathbf{x}-{\mathbf{\mu }}_k{)}^T{\mathbf{\Sigma }}^{-2}(\mathbf{x}-{\mathbf{\mu }}_k)\right]t_{nk}=0$$

Rearranging, multiplying through by ${\mathbf{\Sigma }}^2$ and rearranging again gives

$$\mathbf{\Sigma }=\frac{1}{N}\sum _n\sum _k{t}_{nk}({\mathbf{x}}_n-{\mathbf{\mu }}_k)({\mathbf{x}}_n-{\mathbf{\mu }}_k{)}^T.$$

Summarizing, given Gaussian class-conditional probability densities $p(\mathbf{x}|{\mathcal{C}}_k)=\mathcal{N}(\mathbf{x}|{\mathbf{\mu }}_k,\mathbf{\Sigma })$ with shared covariance and class priors $p({\mathcal{C}}_k)=\pi$, the posterior distributions are given by

$$\begin{array}{rl}p({\mathcal{C}}_k|\mathbf{x})& =\frac{\exp a_k}{{\sum }_j{a}_j}\\ a_k& =\log p(\mathbf{x}|{\mathcal{C}}_k)p({\mathcal{C}}_k)\\ {\pi }_k& =\frac{N_k}{N}\\ {\mathbf{\mu }}_k& =\frac{1}{N_k}\sum _n{t}_{nk}{\mathbf{x}}_n\\ \mathbf{\Sigma }& =\frac{1}{N}\sum _n\sum _k{t}_{nk}({\mathbf{x}}_n-{\mathbf{\mu }}_k)({\mathbf{x}}_n-{\mathbf{\mu }}_k{)}^T\end{array}$$

The previous training followed by prediction paradigm then holds:

1. Training: Given input data set, and class count $K$ compute the $\pi ,{\mathbf{\mu }}_k,\mathbf{\Sigma }$.
2. Prediction: Given an observation $\mathbf{x}$, compute the class posterior distribution $p({\mathcal{C}}_k|\mathbf{x})$ using above expression for $k\in [1,K]$. The point is assigned to the value which has the highest posterior class probability.

Logistic Regression

Reference: Section 4.3

In the preceding section, the class posterior distributions were found by assuming class-conditional probablities and priors over the class conditions, and using maximum likelihood to determine paramters defining these probabilities. Alternatively, one may consider working directly with the functional form of the posterior distribution (and defining decision boundaries).

As indicated to above, the ${\mathcal{C}}_k$ class posterior distribution is a normalized exponential

$$\begin{array}{rl}{y}_k(\mathbf{x})& :=p({\mathcal{C}}_k|\mathbf{x})\\ & =\sigma ({\mathbf{w}}_k^T\mathbf{x})\\ & =\frac{\exp a_k}{{\sum }_j\exp a_j}.\end{array}$$

Note that $y_k(\mathbf{x})\in (0,1)$.

Consider now the special case where the activation coefficients $a_k$ are given by

$$a_k={\mathbf{w}}^T\mathbf{x}$$

which resembles the linear decision surface seen in the preceeding sections and alluded to in the preceeding example.

Without introducing priors over the distributions, consider now the maximization of the likelihood function

$$p(\mathbf{T}|{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_K)=\prod _n^N\prod _k^{K}p(\mathcal{C}|{\mathbf{x}}_n{)}^{t_{nk}}$$

where $\mathbf{T}\in {\mathbb{Z}}_0^{N\times K}$ contains target variables with elements $t_{nk}$. As before, consider maximization of the negative log likelihood, namely

$$\begin{array}{rl}E({\mathbf{w}}_1,\cdots ,{\mathbf{w}}_K)& :=-\log p(\mathbf{T}|{\mathbf{w}}_1,{\mathbf{w}}_2,\cdots ,{\mathbf{w}}_K)\\ & =-\sum _n\sum _k{t}_{nk}\log y_{nk}\\ & =\boxed{-\sum _n\sum _k{t}_{nk}\left(a_k-\log \sum _j{a}_j\right)}\end{array}$$

Note that this is a nonlinear function of the model weights, and the maximization of the log likelihood requires at least its first derivative. Note that

$$\begin{array}{rl}\frac{\partial y_k}{\partial a_m}& =\frac{\partial }{\partial a_m}\frac{\exp \{a_k\}}{{\sum }_j\exp \{a_j\}}\\ & =\frac{{\delta }_{km}\exp \{a_k\}}{{\sum }_j\exp \{a_j\}}-\frac{\exp \{a_k\}}{({\sum }_j\exp \{a_j\}{)}^2}\exp \{a_j\}\\ & =y_k\cdot ({\delta }_{km}-y_j),\end{array}$$

and thus

$${\nabla }_{{\mathbf{w}}_j}E({\mathbf{w}}_1,\cdots ,{\mathbf{w}}_K)=\sum _n(y_{nj}-t_{nj}){\mathbf{x}}_n$$

Support Vector Machines (SVM)

Reference: Chapter 7

This is a good place to take a look at SVMs. Unlike the logistic regression approach considered above which relied on a probabilistic approach for arriving at the objective function, we'll need to start out with a geometric viewpoint.

The support vector machine is inherently for two-class classfication problems (in practice it is used for multi-class problems as well). Consider once again the classification of dataset linearly separable dataset. As before, we consider finding a hyperplane that "best" separates the two datasets. In SVM, we choose "best" to be a hyperplane that maximizing the orthogonal distance between the hyperplanes and the closest dataset to the hyperplane from each dataset.

Recall that $y$ provides a signed distance measure of the a data-point to the hyperplane. Without going into detail (see Section 7.1 of PRML book), we can choose to construct the weights such that $t_n{y}_n=1$ for points that are closest to the hyperplane, which then requires that under for all points

$$t_n{y}_n\ge 1,$$

where $1-t_n{y}_n\ge 0$ implies a point incorrectly classified.

With some handwaving, we can then set the objective function as the classfication error given by

$$E(\mathbf{W})=\sum _n\max (0,1-t_n{y}_n).$$

Despite SVM inherently being a two-class classifier, we can heuristicly extend it for multi-class classification. One approach, commonly known as "one-vs-all" does this by constructing the error function such that a non-zero error is introduced when a wrong class is associated to a given data point.

$$\boxed{E(\mathbf{W})=\sum _n\sum _{j\ne y_n}\max (s_{nj}-s_{y_n}+1,0)},$$

where $y_n$ is the known correct class index for given datapoint under consideration.

Consider the pseudo-example below:

Neural Networks

Reference: Chapter 5

In the preceeding section, direct use of the functional form of the posterior distribution (and defining decision boundaries) was considered. Specifically, discriminants and the decision boundaries were chosen to be a linear function of the "activations". Consider now revisiting Neural Networks which were briefly introduced in Part 1 with application to regression.

Consider a more generalize form of the two-layer network:

$$\begin{array}{r}{y}_k=\sigma \left(\sum _j^M{w}_{kj}^{(2)}h(\sum _j^D{w}_{ji}^{(1)}{x}_i+w_{j0}^{(1)})+w_{k0}^{(2)}\right),k\in \{1,\cdots ,K\}\end{array}$$

In above,

$M$ : Model complexity

$D$: Dimension of dataset (e.g. 2 dataset)

$K>1$: Number of training class datasets

$w_{ji}^{(1)}\in R^{M\times D}$ : First layer model weights

$w_{j0}^{(1)}\in R^M$ First layer model biases

$w_{kj}^{(2)}\in R^{K\times M}$ Second layer model weights

$w_{k0}^{(2)}\in R^K$ Second layer model biases

$h(\cdot )$ (First layer) activation function commonly taken to be sigmoid or $\tanh$ functions

$\sigma (\cdot )$ Second layer activation function, e.g. a logistic sigmoid function

Also, within the context of classification, $y_k$ represents the posterior probability of class ${\mathcal{C}}_k$, while ${\mathbf{x}}_{nk}$ denotes the $n$-th data point belonging to class ${\mathcal{C}}_k$.

Recall that the log likelihood function from the previous was

$$\begin{array}{r}E(\mathbf{w})=-\sum _n^N\sum _k^K{t}_{nk}{\log }_k{y}_k({\mathbf{x}}_n,\mathbf{w})\end{array}$$

Let

$$\begin{array}{rl}h(a)& :=\tanh (a)\\ \sigma (a)& :=\frac{1}{1+\exp \{-x\}}\end{array}$$

where component-wise operation is understood where necessary. Solution to the resulting nonlinear optimization problem requires at minimum a evaluation of the first derivative of the model outputs with respect to the model weights. Using expressions for the derivative of the model output with respect to the model weights, the first derivatives of the objective function with respect to model weights can be computed as

$$\frac{\partial E}{\partial \mathbf{w}}=-\sum _n^N\sum _k^K{t}_{nk}\frac{1}{y_{nk}}\frac{\partial y_{nk}}{\partial {\hat{w}}_{mp}}$$

where $\hat{w}$ represents the model weights including biases.

Summary

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.