text:
\section{General Data}

\begin{tabularx}{\textwidth}{ | L | c | c | }
  \hline
  \textbf{Parameter}                    & \textbf{Value}   & \textbf{Reference} \\ \hline
  \endfirsthead
  \hline
  \textbf{Parameter}                    & \textbf{Value}   & \textbf{Reference} \\ \hline
  \endhead
  Length                                & 8.28 m           & \cite{Cessna172_InformationManual,Janes20042005} \\ \hline
  Wingspan                              & 11.00 m          & \cite{Cessna172_InformationManual,Janes20042005} \\ \hline
  Height                                & 2.72 m           & \cite{Cessna172_InformationManual,Janes20042005} \\ \hline
  Wheelbase                             & 1.63 m           & \cite{Janes20042005} \\ \hline
  Wheel track                           & 2.53 m           & \cite{Janes20042005} \\ \hline
  Wing area                             & 16.17 m\textsuperscript{2} & \cite{Cessna172_InformationManual,Janes20042005} \\ \hline
  Mean aerodynamic chord                & 1.49 m           & \cite{Cessna172_InformationManual} \\ \hline
  Wing airfoil                          & NACA 2412        & \cite{Janes20042005} \\ \hline
  Horizontal tail area                  & 2.00 m\textsuperscript{2} & \cite{Janes20042005} \\ \hline
  Horizontal tail airfoil at root (Cessna 177) & NACA 0012 & \cite{Torenbeek1982} \\ \hline
  Horizontal tail airfoil at tip (Cessna 177)  & NACA 0009 & \cite{Torenbeek1982} \\ \hline
  Vertical tail area                    & 1.04 m\textsuperscript{2} & \cite{Janes20042005} \\ \hline
  Vertical tail airfoil at root         & NACA 0009        & \cite{NASA-CR-158974} \\ \hline
  Vertical tail airfoil at tip          & NACA 0006        & \cite{NASA-CR-158974} \\ \hline
  Ailerons deflection limit             & up 20\degree, down 15\degree & \cite{EASA-IM-A-051} \\ \hline
  Ailerons area (total)                 & 1.70 m\textsuperscript{2} & \cite{Janes20042005} \\ \hline
  Elevator deflection limit             & up 28\degree, down 23\degree & \cite{EASA-IM-A-051} \\ \hline
  Elevator area (including trim tab)    & 1.35 m\textsuperscript{2} & \cite{Janes20042005} \\ \hline
  Elevator trim tab deflection limit    & up 22\degree, down 19\degree & \cite{EASA-IM-A-051} \\ \hline
  Rudder deflection limit               & $\pm$17.7\degree & \cite{EASA-IM-A-051} \\ \hline
  Flaps area                            & 1.98 m\textsuperscript{2} & \cite{Janes20042005} \\ \hline
  Flaps deflection limit                & 30\degree        & \cite{EASA-IM-A-051} \\ \hline
  Standard empty weight                 & 754 kg           & \cite{Cessna172_InformationManual} \\ \hline
  Total fuel tanks capacity             & 212 l            & \cite{Cessna172_InformationManual,EASA-IM-A-051} \\ \hline
  Downwash angle derivative with respect to the aircraft angle of attack & 0.25 & \cite{NASA-CR-2605} \\ \hline
  Engine manufacturer                   & Textron Lycoming & \cite{Cessna172_InformationManual} \\ \hline
  Engine model                          & IO-360-L2A       & \cite{Cessna172_InformationManual} \\ \hline
  Engine rated power (at 2~700 RPM)     & 134.2 kW         & \cite{Cessna172_InformationManual,Cessna172_TextronLycomingOperatorsManual} \\ \hline
  Engine displacement                   & 5.9 l            & \cite{Cessna172_InformationManual} \\ \hline
  Engine height                         & 0.631 m          & \cite{Cessna172_TextronLycomingOperatorsManual} \\ \hline
  Engine width                          & 0.848 m          & \cite{Cessna172_TextronLycomingOperatorsManual} \\ \hline
  Engine length                         & 0.757 m          & \cite{Cessna172_TextronLycomingOperatorsManual} \\ \hline
  Engine standard dry weight            & 126.1 kg         & \cite{Cessna172_TextronLycomingOperatorsManual} \\ \hline
  Specific fuel consumption (at 2~200 RPM) & 253.4 g/(kW$\cdot$h) & \cite{Cessna172_TextronLycomingOperatorsManual} \\ \hline
  Propeller manufacturer                & McCauley         & \cite{Cessna172_InformationManual} \\ \hline
  Propeller model                       & 1A170E/JHA7660   & \cite{Cessna172_InformationManual} \\ \hline
  Number of blades                      & 2                & \cite{Cessna172_InformationManual} \\ \hline
  Propeller diameter                    & 1.93 m           & \cite{Cessna172_InformationManual} \\ \hline
  \caption{General Data}
\end{tabularx}

ext: tex

max_stars_repo_name: marek-cel/mscsim-docs

max_issues_repo_name: marek-cel/mscsim-docs

max_forks_repo_name: marek-cel/mscsim-docs


text:
% Information Theory
%
% Brief introduction to the bits of information theory essential for understanding BLAST

\subsection{Information Theory}
\begin{frame}
  \frametitle{Bits of Information}
  \begin{itemize}
    \item<1-> The amount of information in a thing (e.g. a message, or picture) can be measured in \emph{bits}, and represented as $H$ (entropy)
    \item<1-> One bit = one "yes/no" or "1/0" answer
    \item<2-> The information associated with a probability $p$ can be represented in bits: \\
              $H(p) = -\text{log}_2(p)$
    \item<2-> For a series of $n$ probabilities $\{p_i,$\ldots$, p_n\}$ we have: \\
              $H = -\sum^{n}_{i} p_i \text{log}_2 (p_i)$
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Information in Random Alphabets}
  \begin{itemize}
    \item DNA and protein alphabets have multiple symbols, so have different entropies (are not equally informative)
    \item Assuming random sequences (no symbol bias):
    \begin{itemize}
      \item For a coin, $p(H) = p(T) = 0.5 \implies H(coin) = 1\text{bit}$
      \item For DNA, $p(A) = p(C) = p(G) = p(T) = 0.25 \implies H(\text{DNA}) = 2\text{bits}$
      \item For protein, $p(A) = p(C) = $\ldots$ = 0.05 \implies H(\text{protein}) = 4.32\text{bits}$
    \end{itemize}
    \item Twice as much information in protein than DNA alphabet
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Information in Biased Alphabets}
  \begin{itemize}
    \item Deviation from random sequence reduces entropy
    \begin{itemize}
      \item For DNA, $p(A) = p(C) = p(G) = p(T) = 0.25 \implies H(\text{DNA}) = 2\text{bits}$
      \item For DNA: $p(A) = 0.4; p(C) = 0.3; p(G) = 0.2; p(T) = 0.1 \implies H(\text{DNA}) = 1.85\text{bits}$
    \end{itemize}
    \item Biological sequences/databases deviate from randomness
  \end{itemize}
\end{frame}

\begin{frame}
  \frametitle{Log-odds Scores}
  \begin{itemize}
    \item<1-> How often does a symbol substitution occur, compared to random expectation?
    \item<1-> If $p(\text{Met}) = 0.01$, $p(\text{Leu}) = 0.1$, we would expect a Met-Leu substitution with frequency $p(\text{Met})p(\text{Leu}) = 0.001$
    \item<2-> Observe Met-Leu substitution with frequency $q(\text{Met-Leu}) = 0.002$
    \item<2-> Odds ratio is $\frac{q(\text{Met-Leu})}{p(\text{Met})p(\text{Leu})} = 2$
    \item<3-> Log odds ratio is $\text{log}_2(\frac{q(\text{Met-Leu})}{p(\text{Met})p(\text{Leu})}) = 1\text{bit}$
  \end{itemize}
\end{frame}

max_stars_repo_name: peterjc/Teaching-Intro-to-Bioinf

max_issues_repo_name: peterjc/Teaching-Intro-to-Bioinf

max_forks_repo_name: peterjc/Teaching-Intro-to-Bioinf


text:
\chapter{Linear Models for Classification}
\label{Linear Models for Classification}
The goal in classification is to take an input vector $\vec{x}$ and to assign it to one of $K$ discrete classes $\mathcal{C}_k$ where $K = 1,...,K$.The input space is thereby divided into \textbf{decision regions} whose boundaries are called \textbf{decision boundaries} or \textbf{decision surfaces}.Data sets whose classes can be separated exactly by linear decision surfaces are said to be \textbf{linearly separable}.

In Chapter 1,we identified three distinct approaches to the classification problem.The simplest involves constructing a \textbf{discriminant function} that directly assigns
each vector x to a specific class. A more powerful approach, however, models the \textbf{conditional probability distribution} $p(C_k|x)$ in an \textbf{inference} stage, and then subsequently uses this distribution to make optimal \textbf{decisions}. By separating inference and decision, we gain numerous benefits, as discussed in \ref{chapter:Introduction}. There are two different approaches to determining the conditional probabilities $p(C_k|x)$. One technique is to model them \textbf{directly}, for example by representing them as parametric models and then optimizing the parameters using a training set. Alternatively, we can adopt a generative approach in which we model the \textbf{class-conditional densities} given by p(x|Ck), together with the \textbf{prior probabilities} p(Ck) for the classes, and then we compute the required posterior probabilities using Bayes’ theorem
\begin{equation}
p(\mathcal{C}_k|\vec{x}) = \dfrac{p(\vec{x}\mathcal{C}_k)p(\mathcal{C}_k)}{p(\vec{x})}
\end{equation}

In the linear regression models considered in Chapter 3, the model prediction
y(x,w) was given by a linear function of the parameters w. In the simplest case,
the model is also linear in the input variables and therefore takes the form y(x) =
$\vec{w}^Tx+w_0$, so that y is a real number. For classification problems, however, we wish
to predict discrete class labels, or more generally posterior probabilities that lie in
the range (0, 1). To achieve this, we consider a generalization of this model in which
we transform the linear function of w using a nonlinear function $f( · )$ so that
\begin{equation}
y(\vec{x}) = f(\vec{w}^T\vec{x} + w_0)
\end{equation}
In the machine learning literature $f(\cdot)$ is known as an \textbf{activation function},whereas its inverse is called a \textbf{link function} in the statistics literature.The decision surfaces correspond to y(x) = constant, so that $\vec{w}^T\vec{x} + w_0 = constant$ and hence the decision surfaces are linear functions of $\vec{x}$, even if the function $f(·)$ is nonlinear. For this reason, the class of models described by (4.3) are called \textbf{generalized linear models} (McCullagh and Nelder,1996).However,in contrast to the models used for regression,they are no longer linear in the parameters due to the presence of the nonlinear  function $f(\cdot)$.
\section{Discriminant Functions}
A discriminant is a function that takes an input vector $\vec{x}$ and assigns it to be one of $K$ classes,denoted  $\mathcal{C}_k$.
\subsection{Two classes}
Linear discriminant function
\begin{align}
y(\vec{x}) = \vec{x}^T\vec{x}+w_0
\end{align}
where $\vec{w}$ is called a \textbf{weight vector},and $w_0$ is $bias$.The negative of the bias is sometimes called \textbf{threshold}.

The value of $y(\vec{x})$ gives a signed measure of the perpendicular distance $r$ of point $\vec{x}$ from the decision surface.Consider an arbitrary point $\vec{x}$ and let $\vec{x}_\perp$ be its orthogonal projection onto the decision surface,then
\begin{align}
\vec{x} &= \vec{x}_\perp + r\dfrac{\vec{w}}{\parallel\vec{w}\parallel} \\
\end{align}
Multiplying both sides by $\vec{w}^T$ and adding $\vec{w}_0$,
\begin{align}\label{eqn:distance from point to hyperplane}
y(\vec{x}) &= y(\vec{x}_{\perp}) + r\dfrac{\vec{w}\cdot\vec{w}}{\parallel\vec{w}\parallel} \\
&= r\parallel\vec{w}\parallel\\
\therefore
r &= \dfrac{y(\vec{x})}{\parallel\vec{w}\parallel}
\end{align}

\subsection{Multiple classes}
Consider the extension of linear discriminants to $ K >2$ classes,building a $K$-class discriminants by combining a number of two-class discriminant functions on which leads to some serious difficulties(Duda and Hart,1973).

Consider the use of $K-1$ classifier each of which solves a two-class problem of separating points int particular class $\mathcal{C}_k$ from points not in that class.This is known as a \textbf{one-versus-the-rest} classifier.An alternative is to introduce $K(K-1)/2$ binary discriminant functions,one for every possible pair classes.This is known as a \textbf{one-versus-one} classifier.Both of these two approaches run into the problem of ambiguous regions.We can avoid these difficulties by considering a single $K$-class discriminant comprising $k$ linear functions of the form
\begin{equation}
y_k(\vec{x}) = \vec{w_k}^T\vec{x}+w_{k0}
\end{equation}
and then assigning a point $\vec{x}$ to class $\mathcal{C}_k$ if $y_k(\vec{x}) > y_j(\vec{x})$ for all $j \neq k$.The decision boundary between class $\mathcal{C}_k$ and $\mathcal{C}_j$ is therefore given by $y_k(\vec{x}) = y_j(\vec{x})$ and hence corresponds to a $(D-1)$-dimensional hyperplane defined by
\begin{equation}
(\vec{w_k}-\vec{w_j})^T + (w_{k0}-w_{j0}) = 0
\end{equation}
The decision regions of such a discriminant are always singly connected and convex.

\subsection{Least squares for classification}
Even as a discriminant function(where we use it to make decisions directly and dispense with any probabilistic interpretation)it suffers from some severe problems.Least-squares solutions lack robustness to outliers,and this applies equally to the classification application.Additional points(outliers) produce a significant change in the location of the decision boundary.The sum-of-squares error function penalizes predictions that are 'too correct' in that they lie a long way on the correct side of the decision boundary.More than lack of robustness,least-squares solutions may gives poor results on classification problems.

The failure of least squares lies in the fact that it corresponds to maximum likelihood under the assumption of a Gaussian conditional distribution,whereas binary target vectors clearly have a distribution that is far from Gaussian(a Bernoulli).

\subsection{Fisher's linear discriminant}
\subsubsection{representation}
One way to view a linear classification model is in terms of dimensionality reduction.Consider first the case of two classes,and suppose we take the $D-dimensional$ input vector $\vec{x}$ and project it down to one dimension using
\begin{equation}\label{eqn:Fisher LDA projection}
y = \vec{w}^T\vec{x}.
\end{equation}
If we place a threshold on $y$ and classify $y \geq -w_0$ as class $\mathcal{C_1}$,and otherwise class $\mathcal{C_2}$,then we obtain our standard linear classifier discussed in the previous section.
\subsubsection{evaluation}
By adjusting the components of the weight vector $\vec{w}$,we can select a projection that maximizes the class separation.To begin with,consider a two-class problem in which there are $N_1$ points of class $\mathcal{C_1}$ and $N_2$ points of class $\mathcal{C_2}$,so that the mean vectors of the two classes are given by
\begin{equation}
\vec{m_1} = \dfrac{1}{N_1} \sum_{n\in \mathcal{C_1}}{\vec{x_n}},\vec{m_1} = \dfrac{1}{N_2}\sum_{n\in \mathcal{C_2}}\vec{x_n}
\end{equation}
The simplest measure of the separation of the classes,when projected onto $\vec{w}$,is the separation of the projected  class means.This suggests that we might choose $\vec{w}$ so as to maximize
\begin{equation}
m_2 - m_1 = \vec{x}^T(\vec{m_2}-\vec{m_1})
\end{equation}
where\begin{equation}
m_k = \vec{w}^T\vec{m_k}
\end{equation}
is the mean of the projected data from class $\mathcal{C_k}$.However this expression can be makde arbitrarily large simply by increasing the magnitude of $\vec{w}$.To solve this,we could constrain $\vec{2}$ to have unit length,so that $\sum_{i}{w_i^2}=1$.Using a Lagrange multiplier to perform the constrained maximization,we then find that $\vec{w} \propto (\vec{m_2} - \vec{m_1})$.There is still a problem with this approach that the projected data have considerable overlap for strongly nondiagonal convariances of the class distributions.

\subsubsection{optimization}
The idea proposed by Fisher is to \textbf{maximize a function that will give a large separation between the projected class means while also giving a small variance within each class,thereby minimizing the class overlap}.The projection formula \ref{eqn:Fisher LDA projection} transforms the set of labelled data points in $\vec{x}$ into a labelled set in the one-dimensional space $y$.The within-class variance of the transformed data from class $\mathcal{C_k}$ is therefore given by
\begin{equation}
s_k^2 = \sum_{n\in \mathcal{C_k}}(y_n-m_k)^2
\end{equation}
where\begin{equation}
y_n = \vec{w}^T\vec{x_n}
\end{equation}
We can define the total within-class variance for the whole data set to be simply $s_1^2+s_2^2$.The Fisher criterion is defined to be the ratio of the between-class variance to the within-class variance and is given by
\begin{equation}
J(\vec{w}) = \dfrac{(m_2-m_1)^2}{s_1^2+s_2^2}
\end{equation}

We can make the dependence on \vec{w} explicit by rewrite the Fisher criterion in the form
\begin{equation}\label{eqn:Fisher criterion}
J(\vec{w}) = \dfrac{\vec{w}^T\vec{S_B}\vec{w}}{\vec{w}^T\vec{S_W}\vec{w}}
\end{equation}
where $\vec{S_B}$ is the $between-class$ covariance matrix,given by
\begin{equation}\label{eqn:Fisher between-class covariance matrix}
\vec{S_B} = (\vec{m_2}-\vec{m_1})(\vec{m_2}-\vec{m_1})^T
\end{equation}
so
\begin{equation}
\vec{w}^T\vec{S_B}\vec{w} = \vec{w} ^T(\vec{m_2} - \vec{m_1})(\vec{m_2}-\vec{m_1})^T\vec{w} = (m_2-m_1)(m_2-m_1)
\end{equation}
and $\vec{S_W}$ is the total $within-class$ covariance matrix,given by
\begin{equation}
\vec{S_W} = \sum_{\mathcal{C}}\sum_{n\in \mathcal{C}_1}(\vec{x_n} - \vec{m_1})(\vec{x_n}-\vec{m_1})^T
\end{equation}
so
\begin{equation}
\vec{w}^T\vec{S_W}\vec{w} = \sum_{\mathcal{C}_k}\sum_{n\in \mathcal{C}_k}{\vec{w}^T(\vec{x_n}-\vec{m_k})(\vec{x_n}-\vec{m_k})^T} = \sum_{\mathcal{C}_i}\sum_{n \in \mathcal{C}_k}(y_n-m_k)^2
\end{equation}
Differentiating \ref{eqn:Fisher criterion} with respect to $\vec{w}$,we find that $J(\vec{w})$ is maximized when
\begin{equation}
{(\vec{w}^T\vec{S_B}\vec{w})\vec{S_W}\vec{w}} = (\vec{w}^T\vec{S_W}\vec{w})\vec{S_B}\vec{w}
\end{equation}
Rewriting this as
\begin{equation}\label{eqn:Fisher generalized eigenvalue}
\vec{S_B}\vec{w} = \lambda\vec{S_W}\vec{w}
\end{equation}
where
\begin{equation}
\lambda = \dfrac{(\vec{w}^T\vec{S_B}\vec{w})}{(\vec{w}^T\vec{S_W}\vec{w})}
\end{equation}
which is called a \textbf{generalized eigenvalue} problem solution.

\begin{proof}
\begin{eqnarray}\because\begin{cases}
	\dfrac{\partial}{\partial x}\dfrac{f(x)}{g(x)}
	= \dfrac{f'g-fg'}{g^2}, \text{where } f^\prime
	= \dfrac{\partial f(x)}{\partial x} \text{ and } g'
	= \dfrac{\partial g(x)}{\partial x}	\\
	\dfrac{\partial}{\partial \vec{x}} \vec{x}^T\vec{A}\vec{x} = (\vec{A}+\vec{A}^T)\vec{x}
\end{cases}	\\
\therefore\end{eqnarray}
\begin{eqnarray}
\dfrac{\partial}{\partial \vec{w}}J(\vec{w}) & =  \dfrac{2\vec{S_B}\vec{w}(\vec{w}^T\vec{S_W}\vec{w})-(\vec{w}^T\vec{S_B}\vec{w})\vec{S_W}\vec{w}}{(\vec{w}^T\vec{S_W}\vec{w})^2}
\end{eqnarray}
Setting the derivative to zero,we can get the result.
\end{proof}
In particular,since
\begin{equation}
\lambda\vec{w} = \vec{S_W}(\vec{m_2}-\vec{m_1})(m_2-m_1)
\end{equation}
Multiply both sides of \ref{eqn:Fisher generalized eigenvalue} by $\vec{S_W}^{-1}$ we then obtain
\begin{equation}\label{eqn:Fisher's linear discriminant}
\vec{w} \propto \vec{S_W}^{-1}(\vec{m_2}-\vec{m_1})
\end{equation}
Note that if the within-class covariance is isotropic,so that $\vec{S_W}$ is proportional to the unit matrix($\vec{S_W} \propto \vec{I}$),we find that $\vec{w}$ is proportional to the difference vector of the class means.

The result \ref{eqn:Fisher's linear discriminant} is known as $Fisher's linear discriminant$,although strictly it is not a discriminant but rather a specific choice of direction for projection of the data down to one dimension.However,the projected data can subsequently be used to construct a discriminant,by choosing a threshold $y_0$ so that we classify a new point as belonging to $\mathcal{C}_1$ if $y(\vec{x}) \geq y_0$ and can classify it as belonging to $\mathcal{C}_2$ otherwise.

\subsection{Extensions to higher dimensions and multiple classes}
We can extend the above idea to multiple classes,and to higher dimensional subspaces,by finding a projection $matrix \vec{W}$ which maps from $D$ to $L$ so as to maximize
\begin{equation}
J(\vec{W}) = \dfrac{|\vec{W}\vec{\Sigma_B}\vec{W}^T|}{|\vec{W}\Sigma_W\vec{W}^T|}
\end{equation}
where
\begin{eqnarray}
\Sigma_B \triangleq \sum_{c}{\dfrac{N_c}{N}(\vec{m_c}-\vec{m})(\vec{m_c}-\vec{m})^T} \\
\Sigma_W \triangleq \sum_{c}\dfrac{N_c}{N}\Sigma_c \\
\Sigma_c \triangleq \sum_{i:y_i=c}(\vec{x_i}-\vec{m_c})(\vec{x_i}-\vec{m_c})^T
\end{eqnarray}
The solution can be shown to be
\begin{equation}
\vec{W} = \Sigma_W^{-\frac{1}{2}}\vec{U}
\end{equation}
where $\vec{U}$ are the $\vec{L}$ leading eigenvectors of $\Sigma_W^{-\frac{1}{2}}\Sigma_B\Sigma_W^{-1\frac{1}{2}}$,assuming $\Sigma_W$ is non-singular.(If it is singular,we can first perform PCA on all the data).
\subsection{Probabilistic interpretation of FLDA *}
\subsection{Relation to least squares}
\subsection{Fisher's discriminant for multiple classes}
\subsection{The perceptron algorithm}
Another example of a linear discriminant model is the perceptron of Rosenblatt(1962),which corresponds to a two-class transformation to give a feature vector $\phi(\vec{x})$,and this is then used to construct a generalized linear model of the form
\begin{equation}
y(\vec{x}) = f(\vec{w}^T\vec{\phi}(\vec{x}))
\end{equation}
where the nonlinear activation function $f(\cdot)$ is given by a step function of the form
\begin{equation}
f(a) = \begin{cases}
+1,a \geq 0 \\
-1,a \leq 0
\end{cases}
\end{equation}

\subsubsection{error function}
A natural choice of error function would be the total number of misclassified patterns,which does not lead to a simple learning algorithm because the error is a piecewise constant function of $\vec{w}$,with discontinuities wherever a change in $\vec{w}$ causes the decision boundary to move across one of the data points.Methods based on gradient can't be applied.An alternative error function is known as the \textbf{perceptron criterion},given by
\begin{equation}
E_P(\vec{w}) = -\sum_{n\in \mathcal{M}}{\vec{w}^T\vec{\phi_n}(\vec{x})t_n}
\end{equation}
where we use $t\in \{-1,+1\}$ coding scheme, $\mathcal{M}$ denotes the set of all misclassified patterns.The perceptron criterion associates zero error with any pattern that is correctly classified,whereas for a misclassified pattern $\vec{x}$ it tries to minimize the quantity $-\vec{w}^T\vec{\phi_n}t_n$.The perceptron learning rule is not guaranteed to reduce the total error function at each stage.

However,the \textbf{perceptron convergence theorem} states that if there exists an exact solution(linearly separable),then the perceptron learning algorithm is guaranteed to find an exact solution if a finite number of steps.

Aside from difficulties with the learning algorithm,the perceptron does not provide probabilistic outputs,nor does it generalize readily to $K>2$ classes.The most important limitation,however,arises from the fact that it is based on linear combinations of fixed basis functions.

\subsubsection{optimization}
We now apply the stochastic gradient descent algorithm to this error function.The change in the weight vector $\vec{w}$ is then given by
\begin{equation}
\vec{w}^{\tau+1} = \vec{w}^{\tau} - \eta\nabla E_P(\vec{w}) = \vec{w}^{\tau} + \eta\phi_n t_n
\end{equation}
where $\eta$ is the learning rate parameter and $\tau$ is an integer that indexes the steps of the algorithm.Because the perceptron function $y(\vec{x},\vec{w})$ is unchanged if we multiply $\vec{w}$ by a constant,so we can set the learning rate $\eta$ equal to 1 without of generality.

\section{Probabilistic Generative Models}
We turn next to a probabilistic view of classification and show how models with linear decision boundaries arise from simple assumptions about the distribution of the data,adopting a generative approach in which we model the class-conditional densities $p(\vec{x}|\mathcal{C}_k)$,as well as the class priors $p(\mathcal{c}_k)$,and then use these to compute posterior probabilities $p(\mathcal{C}_k|\vec{x})$ through Bayes' theorem.

\subsection{Two classes}
Consider first of all the case two classes.The posterior probability for class $\mathcal{C}_1$ can be written as
\begin{align}
p(\mathcal{C}_1|\vec{x}) &= \dfrac{p(\vec{x}|\mathcal{C}_1)p(\mathcal{C}_1)}{p(\vec{x}|\mathcal{C}_1)p(\mathcal{C}_1)+p(\vec{x}|\mathcal{C}_2)p(\mathcal{C}_2)} \\
&= \dfrac{1}{1+\exp(-a)} = \sigma(a)
\end{align}
where we have defined
\begin{equation}
a = \ln \dfrac{p(\vec{x}|\mathcal{C}_1)p(\mathcal{C}_1)}{p(\vec{x}|\mathcal{C}_2)p(\mathcal{C}_2)}
\end{equation}
and $\sigma(a)$ is the \textbf{logistic sigmoid} function defined by
\begin{equation}
\sigma(a) = \dfrac{1}{1+\exp(-a)}
\end{equation}
The term 'sigmoid' means S-shaped,sometimes called a 'squashing function' because it maps the whole real axis into a finite interval.It satisfies symmetry property
\begin{equation}
\sigma(-a) = 1 - \sigma(a)
\end{equation}
The inverse of the logistic sigmoid is given by
\begin{equation}
a = \ln(\dfrac{\sigma}{1-\sigma})
\end{equation}
and is known as the \textbf{logit} function.It represents the log of the ratio of probabilities $ln[p(\mathcal{C}_1|\vec{x})/p(\mathcal{C}_2|\vec{x})]$ for the two classes,also known as the \textbf{log odds}.
\subsection{Multiple classes}
For the case of $K>2$ classes,we have posterior
\begin{align}
p(\mathcal{C}_k|\vec{x})
&=\dfrac{p(\vec{x}|\mathcal{C}_k)p(\mathcal{C}_k)}{\sum_{j}{p(\vec{x}|\mathcal{C}_j)p(\mathcal{C}_j)}} \\
&=\dfrac{\exp(a_k)}{\sum_j\exp(a_j)}
\end{align}
which is known as the \textbf{normalized exponential} and can be regarded as a multiclass generalization of the logistic sigmoid.Here the quantities $a_k$ are defined by
\begin{equation}
a_k=\ln p(\vec{x}|\mathcal{C}_k)p(\mathcal{C}_k)
\end{equation}
The normalized function is also known as the \textbf{softmax function},as it represents a smoothed version of the 'max' function because,if $a_k \gg a_j$ for all $j\neq k$,then $p(\mathcal{C}_k|\vec{x}) \simeq 1$ and,$p(\mathcal{C}_j|\vec{x}) \simeq 0$.

The following investigate the consequences of choosing specific forms for the class-conditional densities,continuous and discrete inputs.
\subsection{Continuous inputs}
Assume that the class-conditional densities are \textbf{Gaussian} sharing the \textbf{same covariance matrix} and then explore the resulting form for the posterior probabilities.Thus the density(pdf) for class $\mathcal{C}_k$ is given by
\begin{equation}
p(\vec{x}|\mathcal{C}_k) = \dfrac{1}{(2\pi)^{D/2}}\dfrac{1}{|\vec{\Sigma}|^{1/2}}
\exp\{-\dfrac{1}{2}(\vec{x}-\vec{\mu}_k)^T\vec{\Sigma}^{-1}(\vec{x}-\vec{\mu}_k)  \}
\end{equation}
In the case of two classes,we have
\begin{align}
a(\vec{x}) &= \log\dfrac{p(\vec{x}|\mathcal{C}_1)p(\mathcal{C}_1)}
{p(\vec{x}|\mathcal{C}_2)p(\mathcal{C}_2)} \\
&=-\dfrac{1}{2}(\vec{x}-\vec{\mu}_1)^T\vec{\Sigma}^{-1}(\vec{x}-\vec{\mu}_1) +
\dfrac{1}{2}(\vec{x}-\vec{\mu}_2)^T\vec{\Sigma}^{-1}(\vec{x}-\vec{\mu}_2) + \log\dfrac{p(\mathcal{C}_1)}{p(\mathcal{C}_2)} \\
&=(\vec{\mu_1}-\vec{\mu_2})^T\vec{\Sigma}^{-1}\vec{x}-
\dfrac{1}{2}\vec{\mu_1}^T\vec{\Sigma}^{-1}\vec{\mu_1} +
\dfrac{1}{2}\vec{\mu_2}^T\vec{\Sigma}^{-1}\vec{\mu_2} + \log\dfrac{p(\mathcal{C}_1)}{\mathcal{C}_2} \\
&=\vec{w}^T\vec{x}+w_0
\end{align}
where
\begin{align}
\vec{w} &= \vec{\Sigma}^{-1}(\vec{\mu_1}-\vec{\mu_2}) \\
w_0 &=  -\dfrac{1}{2}\vec{\mu_1}^T\vec{\Sigma}^{-1}\vec{\mu_1} +
\dfrac{1}{2}\vec{\mu_2}^T\vec{\Sigma}^{-1}\vec{\mu_2} + \log\dfrac{p(\mathcal{C}_1)}{\mathcal{C}_2}
\end{align}
We see that the quadratic terms in $\vec{x}$ from the exponents of the Gaussian densities have cancelled (due to the common covariance matrices) leading to a linear function of $\vec{x}$ in the argument of logistic sigmoid.The prior $p(\mathcal{C}_k)$ enter only through the bias parameter $w_0$ so that it have effect of the parallel contours of constant posterior probability.

For the general case of $K$ cases we have
\begin{align}
a_k(\vec{x}) &= \log p(\vec{x}|\mathcal{C}_k)p(\mathcal{C}_k)  \\
&=\vec{\mu_k}^T\vec{\Sigma}^{-1}\vec{x} - \dfrac{1}{2}{\vec{\mu_k}}^T\vec{\Sigma}^{-1}\vec{\mu_k} + \vec{x}^T\vec{\Sigma}^{-1}\vec{x}+
 \log p(\mathcal{C}_k) + const
\end{align}
Notice that $\vec{x}\vec{\Sigma}^{-1}\vec{x}$ is cancelled in the exponent of the softmax function.Therefore we can write
\begin{equation}
a_k(\vec{x}) = \vec{w_k}^T\vec{x}+w_{k0}
\end{equation}
where
\begin{align}
\vec{w_k} &= \vec{\Sigma}^{-T}\vec{\mu_k} \\
w_{k0} &= -\dfrac{1}{2}\vec{\mu_k}^T\vec{\Sigma}^{-1}\vec{\mu_k} + \ln p(\mathcal{C}_k)
\end{align}
We see that the $a_k(\vec{x})$ are again linear functions of $\vec{x}$ due to the shared covariances,otherwise we obtain quadratic functions of $\vec{x}$,giving rise to a \textbf{quadratic discriminant}.

\subsection{Maximum likelihood solution}
Having specified a parametric functional form for the class-conditional densities $p(\vec{x}|\mathcal{C}_k)$,we can then determine the values of the parameters,together with the prior class probabilities $p(\mathcal{C}_k)$,using \textbf{maximum likelihood}.

Consider first the case of two classes,each having a Gaussian class-conditional density with a shared covariance matrix,and suppose we have a data set $\{\vec{x_n},\vec{t_n}\}$ where $n=1,...N$.Here $t_n=1$ denotes class $\mathcal{C}_1$ and $0$ denotes $\mathcal{C}_2$.Denote the prior class probability $p(\mathcal{C}_1) = \pi$ so that $p(\mathcal{C}_2) =  1- \pi$.For a data point $\vec{x_n}$ from class $\mathcal{C}_1$,we have $t_n=1$ and hence
\begin{equation}
p(\vec{x_n},\mathcal{C}_1) = p(\mathcal{C}_1)p(\vec{x}|\mathcal{C}_1) =
\pi \mathcal{N}(\vec{x_n}|\vec{\mu_1},\vec{\Sigma})
\end{equation}
For class $\mathcal{C}_2$ ,and hence
\begin{equation}
p(\vec{x_n},\mathcal{C}_2) = p(\mathcal{C}_1)p(\vec{x}|\mathcal{C}_2) =
\pi \mathcal{N}(\vec{x_n}|\vec{\mu_2},\vec{\Sigma})
\end{equation}
Thus the likelihood is given by
\begin{equation}
p(\vec{t}|\pi,\vec{\mu_1},\vec{\mu_2}) = \prod_{n=1}^{N}
[\pi \mathcal{N}(\vec{x_n}|\vec{\mu_1},\vec{\Sigma})]^{t_n}
[(1-\pi)\mathcal{N}(\vec{x_n}|\vec{\mu_2},\vec{\Sigma})]^{1-t_n}
\end{equation}
where $\vec{t} = (t_1,...,t_N)^T$.
Maximize the log likelihood
\begin{equation}
\log p(\vec{t}|\pi,\vec{\mu_1},\vec{\mu_2}) =
\sum_{n=1}^{N}\{t_n\log \pi + t_n\log \mathcal{N(\vec{x_n}|\vec{\mu_1},\vec{\Sigma})}+
(1-t_n)\log(1-\pi) + (1-t_n)\log \mathcal{N}(\vec{x_n}|\vec{\mu_2},\vec{\Sigma})  \}
\end{equation}
The terms that depend on $\pi$ are
\begin{equation}
\sum_{n=1}^{N}\{t_n\log\pi + (1-t_n)\ln(1-\pi)  \}
\end{equation}
Setting the derivative with respect to $\pi$ equal to zero and rearranging,we obtain
\begin{align}
&\sum_{n=1}^{N}\{t_n\dfrac{1}{\pi} +(t_n-1)\dfrac{1}{1-\pi}  \} = 0 \\
\Rightarrow & \pi\sum_{n=1}^{N}(t_n-1) -\pi\sum_{n=1}^{N}t_n + \sum_{n=1}^{N}t_n = 0 \\
\Rightarrow & \pi = \dfrac{1}{N}\sum_{n=1}^{N}t_n = \dfrac{N_1}{N}=\dfrac{N_1}{N_1+N_2}
\end{align}
where $N_k$ denotes the total number of data points in class $\mathcal{C}_k$.Thus the maximum likelihood estimate for $pi$ is simply the fraction of points in class as expected,which can be easily generalized to the multiclass case where ... .

Now consider the maximization with respect to $\vec{\mu_1}$.The log likelihood function those terms that depend on $\vec{\mu_1}$
\begin{equation}
\sum_{n=1}^{N}t_n\ln\mathcal{N}(\vec{x_n}|\vec{\mu_1},\Sigma) =
-\dfrac{1}{2}\sum_{n=1}^{N}t_n(\vec{x_n} - \vec{\mu_1})^T\vec{\Sigma}^{-1}(\vec{x_n}-\vec{\mu_1})+ const
\end{equation}
Setting the derivative with respect to $\vec{\mu_1}$ to zero and arranging ,we obtain
\begin{equation}
\vec{\mu_1} = \dfrac{1}{N_1}t_n\vec{x_n}
\end{equation}
which is simply the mean of all the input vectors $\vec{x_n}$ assigned to corresponding class.
Similarly,
\begin{equation}
\vec{\mu_2} = \dfrac{1}{N_2}t_n\vec{x_n}
\end{equation}
Finally,consider the maximum likelihood solution for the shared covariance matrix $\Sigma$.We have
\begin{align}
&-\dfrac{1}{2}\sum_{n=1}^{N}t_n\ln|\Sigma| -
\dfrac{1}{2}\sum_{n=1}^{N}t_n(\vec{x_n} - \vec{\mu_1})^T\Sigma^{-1}(\vec{x_n}-\vec{\mu_1}) \\
&-\dfrac{1}{2}\sum_{n=1}^{N}(1-t_n)\ln|\Sigma| - \dfrac{1}{2}\sum_{n=1}^{N}(1-t_n)(\vec{x_n}-\vec{\mu_2})^T\Sigma^{-1}(\vec{x_n}-\vec{\mu_2})\\
&= -\dfrac{N}{2}\ln|\Sigma| -\dfrac{N}{2}Tr\{\Sigma^{-1}\vec{S}\}
\end{align}
where we have defined
\begin{align}
\vec{S} &= \dfrac{N_1}{N}\vec{S_1} + \dfrac{N_2}{N}\vec{S_2} \\
\vec{S_1}&= \dfrac{1}{N_1}\sum_{n\in \mathcal{C}_1}{(\vec{x_n}-\vec{\mu_1})(\vec{x_n}-\vec{\mu_1})} \\
\vec{S_2}&= \dfrac{1}{N_2}\sum_{n\in \mathcal{C}_2}{(\vec{x_n}-\vec{\mu_2})(\vec{x_n}-\vec{\mu_2})}
\end{align}
Using the standard result for the maximum likelihood solution for a Gaussian distribution,we see that
\begin{equation}
\vec{\Sigma} = \vec{S}
\end{equation}
which represents a weighted average of the covariance matrices associated with each of the two classes separately.\footnote{?}

This result is easily extended to the $K$ class problem.Note that the approach of fitting Gaussian distributions to the classes is not robust to outliers,because the maximum likelihood estimation of a Gaussian is not robust.

\subsection{Discrete features}
Assume binary feature value $x_i\in \{0,1\}$ is Bernoulli distributed,thus we have class-conditional distributions of the form
\begin{equation}
p(\vec{x}|\mathcal{C}_k) = \prod_{i=1}^{D}\mu_{ki}^{x_i}(1-\mu_{ki})^{1-x_i}
\end{equation}
which contains $D$ independent parameters for each class,giving
\begin{equation}
a_k(\vec{x}) = \sum_{i=1}^{D}\{x_i\ln\mu_{ki}+(1-x_i)\ln(1-\mu_ki)\} + \ln p(\mathcal{C}_k)
\end{equation}
which again are linear functions of the input values $x_i$.Analogous results are obtained for discrete variables each of which can take $M>2$ states.

\subsection{Exponential family}
As we have see so far,for both Gaussian distributed and discrete inputs,the posterior class probabilities are given by generalized linear models with logistic sigmoid($K=2$ classes) or softmax($K\gg 2$ classes) activation functions.These are particular cases
of a more general result obtained by assuming that the class-conditional densities
p(x|Ck) are members of the exponential family of distributions.
Using the form  \ref{sec:exponential-family} for members of the exponential family, we see that the distribution of x can be written in the form
\begin{equation}
p(\vec{x}|\vec{\lambda_k}) = h(\vec{x})g(\vec{\lambda}_k)
\exp\{\vec{\lambda}_k^T\vec{u}(\vec{x}) \}
\end{equation}
We restrict attention to the subclass for which $\vec{u}(\vec{x})  = \vec{x}$.Making use of \ref{eqn:density scale parameter} to introduce a scaling parameter $s$,then
\begin{equation}\label{eqn:restricted exponential family}
p(\vec{x}|\vec{\lambda_k},s) = \dfrac{1}{s}h(\vec{x}/s)g(\vec{\lambda}_k)
\exp\{\dfrac{1}{s}\vec{\lambda}_k^T\vec{u}(\vec{x}) \}
\end{equation}
Note that each class have its own parameter vector $\vec{\lambda}_k$ but share the same scale parameter $s$.

For the two-class problem,the posterior is given by a logistic sigmoid acting on a linear function
\begin{equation}
a(\vec{x}) = (\vec{\lambda}_1-\vec{\lambda}_2)^T\vec{x} + \ln g(\vec{\lambda}_1) - \ln g(\vec{\lambda}_2) +\ln p(\mathcal{C}_1)-\ln p(\mathcal{C}_2)
\end{equation}
Similarly,for the $K$-class problem,
\begin{equation}
a_k(\vec{x}) = \vec{\lambda}_k^T\vec{x}+\ln g(\vec{\lambda}_k)+\ln p(\mathcal{C}_k)
\end{equation}
and so again is a linear function of $\vec{x}$.

\section{Probabilistic Discriminative Models}
We have seen that the posterior probability for the two-class classification and multiclass case can be written as a logistic and softmax function respectively for a wide choice of class-conditional distributions $p(\vec{x}|\mathcal{C}_k)$.For a specific class-conditional densities $p(\vec{x}|\mathcal{C}_k)$,we use maximum likelihood to determine the parameters of the densities as well as the
class priors $p(\mathcal{C}_k)$ and use Bayes' theorem to find the posterior class probabilities.However,an alternative approach is to use the functional form of the generalized linear model explicitly and to determine its parameters directly by using maximum likelihood,such as \textbf{iterative reweighted least squares} or IRLS.

The indirect approach to finding the parameters of a generalized linear model by fitting class-conditional densities  and class priors separately and then applying Bayes' theorem,represents an example of \textbf{gernerative modelling}，because we could take such a model and generate synthetic data by drawing values of $\vec{x}$ from the marginal distribution $p(\vec{x})$.In the direct approach,we maximize a likelihood function defined through the conditional distribution $p(\mathcal{C}_k|\vec{x})$,which represents a form of \textbf{discriminative} training.

\subsection{Fixed basis functions}
So far,we have considered classification models that work directly with the original input vector $\vec{x}$.However,all of the algorithms can equally applicable if we first make a fixed  nonlinear transformation of the inputs using a vector of basis functions $\phi(\vec{x})$.The resulting decision boundaries will be linear in the feature space $\phi$,corresponding to nonlinear decision boundaries in the original $\vec{x}$ space.One of the basis function is typically set to a constant,say $\phi_0(\vec{x})=1$,so that the corresponding parameter $w_0$ plays the role of a bias.

Note that nonlinear transformations cannot remove class overlap between the class-conditional densities $p(\vec{x}|\mathcal{C}_k)$.

\subsection{Logistic regression}
\subsubsection{representation}
We begin with two-class classification.The posterior probability of class $\mathcal{C}_1$ can be written as a logistic sigmoid acting on a linear function of the feature vector $\phi$ so that
\begin{equation}
p(\mathcal{C}_1|\phi) = y(\phi) = \sigma(\vec{w}^T\phi)
\end{equation}
which has a linear dependence on the number of parameters.Logistic regression model can also be written as
\begin{equation}
p(y|\vec{x},\vec{w})=\mathrm{Ber}(y|\mathrm{sigm}(\vec{w}^T\vec{\phi}))
\end{equation}
where $\vec{w}$ and $\vec{\phi}$ are extended vectors, i.e., $\vec{w}=(b, w_1, w_2,\cdots, w_D)$, $\vec{\phi}=(1, \phi_1, \phi_2,\cdots, \phi_D)$.

Note that the derivative of the logistic sigmoid function is
\begin{equation}
\dfrac{d\sigma}{da} = \sigma(1-\sigma)
\end{equation}

\subsubsection{evaluation}
For a data set $\{\phi_n,t_n \}$,where $t_n \in \{0,1\}$ and $\phi_n = \phi(\vec{x_n})$,with $n = 1,...,N$,the likelihood function can be written
\begin{equation}
p(\vec{t}|\vec{w}) = \prod_{n=1}^{N}y_n^{t_n}\{1-y_n \}^{1-t_n}
\end{equation}
where $\vec{t}=(t_1,...,t_N)^T$ and $y_n = p(\mathcal{C}_1|\phi_n)$.The negative logarithm of the likelihood,which gives the \textbf{cross-entropy} error function is in the form
\begin{equation}
E(\vec{w})=-\ln p(\vec{t}|\vec{w})=-\sum_{n=1}^{N}\{t_n\ln y_n +(1-t_n)\ln (1-y_n) \}
\end{equation}
Taking the gradient of the error function with respect to $\vec{w}$,we obtain
\begin{align}
\nabla E(\vec{w})
&=-\sum\limits_{n=1}^{N}\{t_n\dfrac{y_n(1-y_n)}{y_n}-
(1-t_n)\dfrac{y_n(q-y_n)}{1-y_n} \}\phi_n\\
&=-\sum_{n=1}^{N}\{t_n(1-y_n)-(1-t_n)y_n \}\phi_n\\
&=\sum_{n=1}^{N}(y_n-t_n)\phi_n
\end{align}
We see that the factor involving the derivative of the logistic sigmoid has cancelled.In particular,the contribution to the gradient from data point $n$ is given by the 'error' $y_n-t_n$ between the target value and the prediction of the model,times the basis function vector $\phi_n$,taking the same form as the gradient of the sum-of-squares error function for the linear regression model.

Note that maximum likelihood can exhibit severe \textbf{over-fitting} for data sets that are linearly separable.The singularity can be avoided by inclusion of a prior and finding a \textbf{MAP(maximum a posterior)} solution for $\vec{w}$,or equivalently by adding a \textbf{regularization} term to the error function.

\subsubsection{optimization}
See the following contents.
\subsection{Iterative reweighted least squares}
\subsubsection{Newton-Raphson update}
For logistic regression,there is no longer a closed-form solution,due to the nonlinearity of the logistic sigmoid function.The error function is concave,and can be minimized by an efficient iterative technique based on the \textbf{Newton-Raphson} iterative optimization scheme,which  uses a local quadratic approximation to the log likelihood function.
\begin{equation}
\vec{w}^{(new)} = \vec{w}^{(old)}-\vec{H}^{-1}\nabla E(\vec{w})
\end{equation}
where $\vec{H}$ is the Hessian matrix whose elements comprise the second derivatives of $E(\vec{w})$ with respect to the components of $\vec{w}$.

Now apply the Newton-Raphson method ot the linear regression model with sum-of-squares error function.The gradient and Hessian of this error function are given by
\begin{align}
\nabla E(\vec{w}) &= \sum_{n=1}^{N}(\vec{w}^T\phi_n-t_n)\phi_n = \vec{\Phi}^T\vec{\Phi}\vec{w} - \vec{\Phi}^T\vec{t} \\
\vec{H} = \nabla\nabla E(\vec{w}) &= \sum_{n=1}^{N}\phi_n\phi_n^T=\vec{\phi}^T\vec{\Phi}
\end{align}
where $\vec{\Phi}$ is the $N\times M$ design matrix,whose $n^{th}$ row is given by $\phi_n^T$.The Newton-Raphson update then takes the form
\begin{align}
\vec{w}^{(new)} &= \vec{w}^{(old)} -(\vec{\Phi}^T\vec{\Phi})^{-1}\{\vec{\Phi}^T\vec{\Phi}\vec{w}^{(old)} -\vec{\Phi}^T\vec{t} \} \\
&=(\vec{\Phi}^T\vec{\Phi})^{-1}\vec{\Phi}^T\vec{t}
\end{align}
which we recognize as the standard least-squares solution.

Now apply the Newton-Raphson update to the cross-entropy  error function for logistic regression model.Gradient and Hessian of this error function are given by
\begin{align}
\nabla E(\vec{w})&=\sum_{n=1}^{N}(y_n-t_n)\phi_n
=\vec{\Phi}^T(\vec{y}-\vec{t}) \\
\vec{H} &= \nabla\nabla E(\vec{w})= \sum_{n=1}^{N}y_n(1-y_n)\phi_n\phi_n^T=\vec{\Phi}^T\vec{R}\vec{\Phi}
\end{align}
where we have introduced the $N\times N$ diagonal matrix $\vec{R}$ with elements
\begin{equation}
R_{nn}=y_n(1-y_n)
\end{equation}
We see that the Hessian depends on $\vec{w}$ and is positive definite.

The Newton-Raphson update formula for the logistic regression model then becomes
\begin{align}
\vec{w}^{(new)} &= \vec{w}^{(old)} -(\vec{\Phi}^T\vec{R}\vec{\Phi})^{-1}\vec{\Phi}^T(\vec{y}-\vec{t}) \\
&=(\vec{\Phi}^T\vec{R}\vec{\Phi})^{-1}\{\vec{\Phi}^T\vec{R}\vec{\Phi}\vec{w}^{(old)} -\vec{\Phi}^T(\vec{y}-\vec{t}) \} \\
&=(\vec{\Phi}^T\vec{R}\vec{\Phi})^{-1}\vec{\Phi}^T\vec{R}\vec{z}
\end{align}
where $\vec{z}$ is an $N$-dimensional vector with elements
\begin{equation}
\vec{z} = \vec{\Phi}\vec{w}^{(old)}-\vec{R}^{-1}(\vec{y}-\vec{t})
\end{equation}
We see that the update formula takes the form of a set of normal equations for a weighted least-squares problem.Because the weighting matrix $\vec{R}$ depends on the parameter vector $\vec{w}$,we must apply the normal equations iteratively,each time using the new weight vector $\vec{w}$ to compute a revised weighing matrix $\vec{R}$.For this reason,the algorithm is known as \textbf{iterative reweighted least squares},or IRLS.The elements of the diagonal weighing matrix $\vec{R}$ can be interpreted as variances because the mean and variance of $t$ in the logistic regression model are given by
\begin{align}
\mathbb{E}[t] &=\sigma(\vec{x}) = y\\
var[t]&=\mathbb{E}[t^2] -\mathbb{E}[t]^2 = \sigma(\vec{x})-\sigma(\vec{x})^2 = y(1-y)
\end{align}
where we have used the property $t^2=t$ for $t\in \{0,1\}$.

\subsubsection{Linearization}
In fact,we can interpret IRLS as the solution to a linearized problem in the space of the variable $a=\vec{w}^T\phi$.Apply Taylor expansion and then we obtain
\begin{align}
a_n(\vec{w}) &\simeq a_n(\vec{w}^{(old)})+\dfrac{da_n}{dy_n}\bigg|_{w^{(old)}}(t_n-y_n) \\
&=\vec{\phi}_n^T\vec{w}^{(old)}-\dfrac{(y_n-t_n)}{y_n(1-y_n)} =z_n
\end{align}

\subsection{Multiclass logistic regression}
In our discussion of generative models for multiclass classification,the posterior probabilities are given by a \textbf{softmax} transformation of linear functions of the feature variables,so that
\begin{align}
p(\mathcal{C}_k|\phi) = y_k(\phi) &= \dfrac{\exp(a_k)}{\sum_j\exp(a_j)}\\
&=\dfrac{\exp(a_k)}{sum} \\
\end{align}
\textbf{Softmax function} is
\begin{align}
\delta(\vec{z}_j) &= \dfrac{\exp^{\vec{z}_j}}{\sum_{k=1}^{K}\exp^{\vec{z}_k}},\text{for } j=1,...,K.  \\
&=\dfrac{\exp^{\vec{z}_j-\max{\vec{z}}}}{\sum_{k=1}^{K}\exp^{\vec{z}_k - \max{\vec{z}}}} \\
\log \delta(\vec{z}_j) &= \vec{z}_j-\max{\vec{z}}-\log\sum_{k=1}^{K}\exp^{\vec{z}_k - \max{\vec{z}}}
\end{align}

where $sum$ denotes $\sum_j\exp(a_j)$ the 'activations' $a_k$ are given by
\begin{equation}
a_k = \vec{w}_k^T\phi
\end{equation}
$\vec{w}_k$ is a $M\times 1$ column vector parameter, and $\phi$ is the $M \times 1$ data column vector.

Here we consider the use of maximum likelihood to determine the parameters $\{\vec{w}_k\}$ of this model directly.The derivatives of $y_k$ with respect to all of the activations $a_j$ are given by
\begin{align}
\dfrac{\partial y_k}{\partial a_j} = \dfrac{\partial p_k}{\partial a_j}
&=\dfrac{e^{a_k}\vec{I}_{kj}sum-e^{a_k}e^{a_j}}{sum^2} \\
%&=\dfrac{\partial\dfrac{1}{1+\sum_{i\neq k}\exp(a_i-a_k)}}{\partial a_j}\\
%&=-\dfrac{e^{a_i-a_k}}{(1+\sum_{i\neq k}e^{a_i-a_k})^2}\\
&=y_k(I_{kj}-y_j)
\end{align}
where $\vec{I}_{kj}$ are the elements of the identity matrix.

Next we write down the likelihood function,using the $1-of-K$ coding scheme.
\begin{equation}
y = p(\vec{T}|\vec{w}_1,...,\vec{w}_K)=
\prod_{n=1}^{N}\prod_{k=1}^{K}p(\mathcal{C}_k|\phi_n)^{t_{nk}}
=\prod_{n=1}^{N}\prod_{k=1}^{K}y_{nk}^{t_{nk}}
\end{equation}
where $y_{nk}=y_k(\phi_n)$,and $\vec{T}$ is an $N\times N$ matrix of target variables with elements $t_{nk}$.Taking the negative logarithm then gives
\begin{equation}
E(\vec{w_1},...,\vec{w_K})=-\ln p(\vec{T}|\vec{w}_1,...,\vec{w}_K)
=-\sum_{n=1}^{N}\sum_{k=1}^{K}t_{nk}\ln y_{nk}
\end{equation}
which is known as the $cross-entropy$ error function for the multiclass classification problem.

Take the gradient of the error function,making use of the derivatives of the softmax function,we obtain
\begin{align}
\because
\frac{\partial E_n}{\partial a_j}
&=-\sum_kt_k\frac{\partial \log p_k}{\partial a_j}\\
&=-\sum_kt_k\frac{1}{p_k}\frac{\partial p_k}{\partial a_j}\\
&=-t_j(1-p_j)-\sum_{k\neq j}t_k\frac{1}{p_k}({\color{red}{-p_kp_j}})\\
&=-t_j(1-p_j)+\sum_{k\neq j}t_k({\color{red}{p_j}})\\
&=-t_j+\color{blue}{t_jp_j+\sum_{k\neq j}t_k({p_j})}\\
&=\color{blue}{p_j\left(\sum_kt_k\right)}-t_j \\
&=p_j-t_j \\
&=y_j-t_j \\
\\
\therefore
\nabla_{w_j}E(\vec{w_1},...,\vec{w_K}) &= \sum_{n=1}^{N}(y_{nj}-t_{nj})\phi_n \\
\sum_{n=1}^{N}\phi_n(y_{nj}-t_{nj})
\therefore
\nabla_{\vec{W}}E(\vec{w}_1,...\vec{w}_k) &= \vec{\Phi}^T(\vec{Y}-\vec{T})
\end{align}
where $\phi_n = \phi(\vec{x}_n)$ is basis function of $\vec{x}_n$ of dimension $M \times 1$, $\vec{Y}$ is the $N\times K$ prediction matrix, $\vec{T}$ is the $N\times K$ target label matrix, $\vec{\Phi}$ is the $N\times M$ design matrix, $\vec{W}$ is $M \times K$ weight matrix of parameters.
\begin{equation}
\end{equation}

To find a batch algorithm,we evaluate the Hessian matrix
\begin{equation}
\nabla_{w_k}\nabla_{w_j}E(\vec{w_1},...,\vec{w_K})
=-\sum_{n=1}^{N}y_{nk}(I_{kj}-y_{nj})\phi_n\phi_n^T
\end{equation}
The Hessian matrix is also positive definite and so the error function has a unique minimum.

\subsection{Probit regression}
A broad range of class-conditional distributions,described by the exponential family,the resulting posterior class probabilities are given by a logistic(or softmax) transformation acting on a linear function of the feature variables,but not all choices of class-conditional density give rise to such a simple form for the posterior probabilities.

The cumulative distribution function is given by
\begin{equation}
	\Phi(a) = \int_{-\infty}^{a}\mathcal{N}(\theta|0,1)d\theta
\end{equation}
which is known as the \textbf{probit function}.It has a sigmoidal shape.A evaluation of a closely related function is
\begin{equation}
erf(a) = \dfrac{2}{\sqrt{\pi}}\int_{0}^{a}\exp(-\theta^2/2)d\theta
\end{equation}
known as the \textbf{erf function} or \textbf{error function}.The generalized linear model based on a probit function is known as \textbf{probit regression}.

One issue that occur in practical applications is that of \textbf{outliers}.Because probit activation function they decay like $\exp(-x^2)$,so probit model can be significantly more sensitive to outliers.However,the effect of mislabelling is easily incoporated into a probabilistic model by introducing a probability $\epsilon$ that the target value $t$ has been flipped to the wrong value,
\begin{align}
p(t|\vec{x}) &= (1-\epsilon)\sigma(\vec{x})+\epsilon(1-\sigma(\vec{x}))\\
&=\epsilon + (1-2\epsilon)\sigma(\vec{x})
\end{align}
where $\sigma(\vec{x})$ is the activation function with input vector $\vec{x}$.Here $\epsilon$ may be set in advance,or it may be treated as a \textbf{hyperparameter} whose value is inferred from data.

\subsection{Canonical link functions}
We now show that there is a general result of assuming a conditional distribution for the target variable from the exponential family,along with a corresponding choice for the activation function known as the \textbf{canonical link function}.

Making use of the restricted form of exponential family distribution \ref{eqn:restricted exponential family} for target variable $t$
\begin{equation}
p(t|\eta,s)=\dfrac{1}{s}h(\dfrac{t}{s})g(\eta)\exp\{\dfrac{\eta t}{s}\}
\end{equation}
The conditional mean of $t$,which denoted by $y$,is given by
\begin{equation}
y \equiv \mathbb{E}[t|\eta]=-s\dfrac{d}{d\eta}\ln g(\eta)
\end{equation}
Thus $y$ and $\eta$ must be related,and we denote this relation through $\eta = \psi(y)$.

We define a \textbf{generalized linear model} to be one for which $y$ is a nonlinear function of a linear combination of the input(or feature) variables so that
\begin{equation}
y=f(\vec{w}^T\vec{\phi})
\end{equation}
where $f(\cdot)$ is known as the \textbf{activation function},and $f^{-1}(\cdot)$ is known as the \textbf{link function} in statistics.

The log likelihood function for this model,which,as a function of $\eta$ ,is given by
\begin{equation}
\ln p(\vec{t}|\eta,s) = \sum_{n=1}^{N}\ln p(t_n|\eta,s) = \sum_{n=1}^{N}\{\ln g(\eta_n)+\dfrac{\eta_n t_n}{s} \} + const
\end{equation}
where we assume that all observations share a common scale parameter(which corresponds to the noise variance for a Gaussian distribution for instance) and so $s$ is independent of $\eta$.The derivative with respect to the model parameter $\vec{w}$ is then given by
\begin{align}
\nabla_w\ln p(\vec{t}|\eta,s)
&=\sum_{n=1}^{N}\{\dfrac{d}{d\eta_n}\ln g(\eta_n)+\dfrac{t_n}{s} \} \dfrac{d\eta_n}{dy_n}\dfrac{dy_n}{da_n}\nabla a_n \\
&= \sum_{n=1}^{N}\dfrac{1}{s}\{t_n-y_n \}\psi'(y_n)f'(a_n)\phi_n
\end{align}
where $a_n=\vec{w}^T\vec{\phi_n}$,and we have used $y_n=f(a_n)$ together.There is a considerable simplification if we choose a particular form of the link function $f^{-1}(y)$ given by
\begin{equation}
f^{-1}(y) = \psi(y)
\end{equation}
In this case,the gradient of error function reduces to
\begin{equation}
\nabla\ln E(\vec{w}) = \dfrac{1}{s}\sum_{n=1}^{N}\{y_n-t_n\}\phi_n
\end{equation}
For the Gaussian $s=\beta^{-1}$,whereas for logistic model $s=1$.

\subsection{The Laplace Approximation}
\subsection{Model comparison and BIC}

\section{Bayesian Logistic Regression}
\subsection{Laplace approximation}
\subsection{Predictive distribution}

max_stars_repo_path: prml/Linear Models for Classification.tex

max_stars_repo_name: Alexoner/Statistical-formula

max_issues_repo_name: Alexoner/Statistical-formula

max_forks_repo_name: Alexoner/Statistical-formula


text:
\documentclass{article}

%% preamble
\usepackage{hyperref}
\usepackage{verbatim}
\usepackage{color}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{cancel}

\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 9.in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in
\newcommand{\myhrule}{ \begin{center}\rule{.9\linewidth}{.25mm}\end{center} }
\definecolor{darkgray}{rgb}{0.95,0.95,0.95}
\definecolor{heavygray}{rgb}{0.05,0.05,0.05}
\definecolor{foo}{rgb}{.8,0,.8}
\newcommand{\pad}{\vspace{8pt}\noindent}
\newcommand{\red}[1]{{\color{red}#1\color{black}}}
\newcommand{\myhref}[2]{\href{#1}{\color{foo}\underline{#2}\color{black}}}

\newcommand{\Mod}[1]{\ \mathrm{mod}\ #1}

\begin{document}

\title{CMDA 3634 Spring 2018 Homework 04}

%% change this to your name
\author{Pim Silpacharn}
\vspace{-64pt}\maketitle
\begin{center}\underline{You must complete the following task by 5pm on Tuesday 04/10/18.}\end{center}
Your write up for this homework should be presented in a {\LaTeX} formatted PDF document. You may copy the \LaTeX{} used to prepare this report as follows

\begin{enumerate}
\item Click on this  \myhref{https://www.sharelatex.com/read/twcrgdqfcvvr}{link}
\item Click on Menu/Copy Project.
\item Modify the HW04.tex document to respond to the following questions.
\item Remember: click the Recompile button to rebuild the document when you have made edits.
\item Remember: Change the author.

\end{enumerate}

\pad \emph{Each student} must individually upload the following files to the CMDA 3634 Canvas page at \myhref{https://canvas.vt.edu}{https://canvas.vt.edu}

\begin{enumerate}
\item \verb|firstnameLastnameHW04.tex| {\LaTeX} file.
\item Any figure files to be included by \verb|firstnameLastnameHW04.tex| file.
\item \verb|firstnameLastnameHW04.pdf| PDF file.
\end{enumerate}

\noindent In addition, all source code must be submitted to an online git repository as follows:
\begin{enumerate}
    \item While on the webpage for you git repository, go to Settings $\rightarrow$ Collaborators.
    \item Add nchalmer@vt.edu and zjiaqi@vt.edu as collaborators.
    \item Make a folder named HW04 in you repository and store all relevant source files to this assignment in this folder. Ensure your assignment files compile with {\verb make }.
\end{enumerate}

\pad You must complete this assignment on your own.

\vspace{16pt}
\begin{center}
\underline{\bf 130 points will be awarded for a successful completion.}
\vspace{8pt}\underline{\bf Extra credit will be awarded as appropriate.}
\end{center}

\newpage

\vspace*{1em}
\noindent{\bf Q0.1}(30 points) 30 points is awarded for attending class on April 3rd and submitting ``Adding CUDA to a Mandelbrot set generator'' code to GitHub. All code and any plotting results should be saved to your GitHub repository. Any additional comments can be included here in this homework.

\vspace*{1em}

\section*{ElGamal Public-key Cyptography}

Using the work done in the previous assignments, we now can perform the three main operations of an ElGamal cryptographic system; we can select a suitable prime and generator of $\mathbb{Z}_p$, encrypt elements of $\mathbb{Z}_p$, and also decrypt them. We then used MPI to crack an ElGamal encryption in parallel. In this assignment we will add some additional functionality to our ElGamal codes, and use parallelism given by OpenMP to improve performance.

\pad {\bf Encrypting Strings}:
In this assignment we'll deal with the question we have been avoiding since Homework 1: if we can only encrypt messages $m$ when $m\in \mathbb{Z}_p$, how do we encrypt `real' messages, i.e. strings?

In our programs we have been using the unsigned integer data type, {\verb unsigned } {\verb int }, which holds an integer number in memory as its 32-bit binary representation. The usual {\verb int } type is similar, but instead uses the first bit to hold the sign of the number (0 for positive, 1 for negative). Somewhat surprisingly, the 8-bit {\verb char } data type also has a `sign' bit. This means if you were to cast a single {\verb char } into an integer it could take values from $[-127,127]$. It will therefore be more convenient for us to use the {\verb unsigned } {\verb char } data type to hold our strings, as casting a single {\verb unsigned } {\verb char } to an {\verb unsigned } {\verb int } can take values from $[0,255]$.

The goal of our program will be to consider a string as an array of {\verb unsigned } {\verb char } variables and cast this string into an array of {\verb unsigned } {\verb int } variables so that each entry is an element of $\mathbb{Z}_p$. Once this is done, we should be able to perform our usual encryption and decryption functions and then cast back to an {\verb unsigned } {\verb char } array to recover our original message.

\vspace*{1em}
\noindent{\bf Q1.1}(5 points) To begin, we need to decide how many characters we can encode with a single element of $\mathbb{Z}_p$. Notice that, in order to encode a single character $\mathbb{Z}_p$ must contain at least $256$ elements, and hence $p$ must be at least 9 bits in size. At the beginning of the {\verb main } function, use the user's input bit length {\verb n } and your knowledge of the size of the prime $p$ to set the variable {\verb charsPerInt } to be the number of {\verb unsigned } {\verb char } variables we can represent with a single element of $\mathbb{Z}_p$.

\vspace*{1em}
\noindent{\bf Q1.2}(15 points) In the {\verb main } function we have set up an {\verb unsigned } {\verb char } string message. To cast this string into an array of elements of $\mathbb{Z}_p$ we need to divide it into pieces of {\verb charsPerInt } characters. However, if the length of the string isn't divisible by {\verb charsPerInt } we will run into issues. We therefore first pad the end of the string with {\verb ' } {\verb ' } characters so that its length is divisible by {\verb charsPerInt }. We have called this function {\verb padString }. Complete the implementation of this function in {\verb functions.c }. (Hint: Remember that a string's termination character is {\verb '\0' }. A useful function may be {\verb strlen } which will return the length of an input string.).

\vspace*{1em}

\vspace*{1em}
\noindent{\bf Q1.3}(15 points) Once the string is padded we must cast the string into an array of {\verb unsigned } {\verb int } values.  We have called this function {\verb convertCharToZ }. Complete the implementation of this function in {\verb functions.c }.

\vspace*{1em}
\noindent{\bf Q1.4}(15 points) The conversion of the integer array back to a string is done in the function {\verb convertZToChar }. Complete the implementation of this function in {\verb functions.c }.

\pad {\bf Acceleration with OpenMP}:
The strength of OpenMP is the ease at which we can parallelize sections of our program. Here we quickly identify and parallelize some embarrassingly parallel sections.

\vspace*{1em}
\noindent{\bf Q2.1}(5 points) The process of converting a string to an integer array and vise versa would appear to be embarrassingly parallel. Use OpenMP to parallelize these functions.

\vspace*{1em}
\noindent{\bf Q2.2}(5 points) The encryption and decryption of the array of integers in $\mathbb{Z}_p$ is also embarrassingly parallel. Use OpenMP to parallelize these functions.

\vspace*{1em}
\noindent{\bf Q2.3}(20 points) Searching for the secret key at the end of the {\verb main } function is certainly embarrassingly parallel. Use OpenMP to parallelize this loop. Afterwards, alter the loop so that if one of the OpenMP threads finds the secret key, all the threads will exit the loop. (Hint: With MPI, you did this using an Allreduce, but with OpenMP we don't need a reduction. We just need a way to share the information among the threads is a safe way).

\vspace*{1em}
\noindent{\bf Q2.4}(20 points) Using New River, run your program several times for $1, 2, 4, 8, 12, 16,$ and 20 OpenMP threads. Recreate the runtime and throughput plots from HW03 to investigate the performance of your program with OpenMP threads. Is OpenMP's parallel performance comparable to MPI? Why/why not?

\vspace*{1em}
\textbf{Answer:}
Each MPI process allocates the same amount of memory, whereas with OpenMP uses a shared memory model. I personally think that OpenMP is the easiest one to get started with. In the beginning you'll only need very few directives to parallelize some of your for-loops. It is hard to understand the full OpenMP standard and it is also hard to get the best performance out of this. In some cases, OpenMP is comparable to MPI but it really just depends on the scenario.

\vspace*{1em}
\noindent{\bf Bonus:}(20 points) It is sometimes more convenient to use only strings when doing cryptographic operations, only casting the characters to integer when we need to do operations in $\mathbb{Z}_p$. While you could rewrite our encryption and decryption functions to do something like this, keep them for now and instead write another function which can convert all of the cyphertext pairs $(\hat{m},a)$ into a single string, where the first several characters of the string correspond to the first cyphertext $(\hat{m},a)$ and so on. You will have to determine how many characters are required to encode a pair $(\hat{m},a)$, and be careful to remember that both $\hat{m}$ and $a$ can be any element of $\mathbb{Z}_p$ now. Print out the cyphertext as a string once it is converted.

%\myhrule

%\begin{thebibliography}{9}
%\bibitem{gol} Wikipedia -- description of the Game of Life
%\url{https://en.wikipedia.org/wiki/Conway's_Game_of_Life}

%\end{thebibliography}

\end{document}


text:
% Copyright 2008 V. Bos, T. van Deursen, and S. Mauw
% This file is part of the MSC Macro Package.
%
\documentclass[12pt,a4paper]{article}
\usepackage{a4wide}
\usepackage{url}
\usepackage{moreverb}
\usepackage{multicol}
\usepackage{msc}
\usepackage{lipsum}
\usepackage{paralist}
\usepackage{verbatim}
\usepackage{fancyvrb}
\usepackage{relsize}
% \usepackage{adjustbox,lipsum}
\usepackage{array}
\usepackage{listings}
\usepackage{minipage-marginpar}
\usepackage{graphicx}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{subfig}
\usepackage{array}
\usepackage{booktabs}

% we allow a ragged right
\setlength{\rightskip}{0pt plus 0.05\linewidth minus 0pt}
\newlength{\rpwidth}
\setlength{\rpwidth}{.5cm}
\newlength{\rpheight}
\setlength{\rpheight}{0.5\levelheight}
\newcommand{\rpN}{%
  \psframe(-0.5\rpwidth,-\rpheight)(0.5\rpwidth,0\rpheight)%
  \rput[B](0\rpwidth,-0.8\rpheight){\tiny \textsc{n}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpNE}{%
  \psframe(-\rpwidth,-\rpheight)(0\rpwidth,0\rpheight)%
  \rput[B](-.5\rpwidth,-0.8\rpheight){\tiny \textsc{ne}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpE}{%
  \psframe(-\rpwidth,-.5\rpheight)(0\rpwidth,.5\rpheight)%
  \rput[B](-.5\rpwidth,-0.3\rpheight){\tiny \textsc{e}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpSE}{%
  \psframe(-\rpwidth,0\rpheight)(0\rpwidth,\rpheight)%
  \rput[B](-.5\rpwidth,0.2\rpheight){\tiny \textsc{se}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpS}{%
  \psframe(-.5\rpwidth,\rpheight)(.5\rpwidth,0\rpheight)%
  \rput[t](0\rpwidth,0.8\rpheight){\tiny \textsc{s}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpSW}{%
  \psframe(0\rpwidth,0\rpheight)(\rpwidth,\rpheight)%
  \rput[B](.5\rpwidth,0.2\rpheight){\tiny \textsc{sw}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpW}{%
  \psframe(0\rpwidth,-.5\rpheight)(\rpwidth,.5\rpheight)%
  \rput[B](.5\rpwidth,-0.3\rpheight){\tiny \textsc{w}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}
\newcommand{\rpNW}{%
  \psframe(0\rpwidth,-\rpheight)(\rpwidth,0\rpheight)%
  \rput[B](.5\rpwidth,-0.8\rpheight){\tiny \textsc{nw}}%
  \pscircle[fillstyle=solid,fillcolor=black](0\rpwidth,0\rpheight){.5\labeldist}%
}

% The following code is taken from the doc package. It defines a global
% macro \bslash that produces a bslash (if present in the current font).
\makeatletter
{\catcode`\|=\z@ \catcode`\\=12 |gdef|bslash{\}}
\makeatother
\newcommand{\cmd}[1]{\texttt{\bslash #1}}

\newcommand{\acro}[1]{{#1}}

\newcommand{\MSC}{\acro{MSC}}
\newcommand{\HMSC}{\acro{HMSC}}
\newcommand{\MSCdoc}{\MSC{}doc}
\newcommand{\mscpack}{\MSC{} macro package}

\newcommand{\env}[1]{\texttt{#1}}
\newcommand{\opt}[1]{[#1]}
\newcommand{\cmdarg}[1]{\{\emph{#1}\}}
\newcommand{\coordarg}[1]{\emph{#1}}
\newcommand{\coordargs}[2]{(\coordarg{#1},\coordarg{#2})}
\newcommand{\lnsvalue}[3]{large/normal/small value #1/#2/#3}

\newenvironment{defs}{%
  \begin{list}{}%
              {\setlength{\labelwidth}{0pt}%
               \setlength{\labelsep}{1em}%
               \setlength{\leftmargin}{1em}%
               \setlength{\parsep}{1ex}%
               \setlength{\listparindent}{0pt}%
               \setlength{\rightmargin}{0pt}%
               \renewcommand{\makelabel}[1]{##1}%
               \raggedright%
              }%
  }{%
  \end{list}}

\title{
  VULFI: An Instruction-level Fault Injection Framework\\{\large User Manual}
}

\author{
 \begin{tabular}{c}
  \begin{tabular}{cc}
   Vishal Sharma &
   Ganesh Gopalakrishan  \\
   %
    \multicolumn{2}{c}{\scriptsize School of Computing, University of Utah} \\
   %
   \multicolumn{2}{c}{\scriptsize \texttt{\{vcsharma,ganesh\}@cs.utah.edu}} \\
  \end{tabular}\\
 \end{tabular}
}

\date{\small Version 1.0, last updated \today\\
%       Describing \mscpack{} version \mscversion
      }

\begin{document}

\maketitle

\newpage

\tableofcontents

\newpage

\section{Getting Started}
\label{started}
VULFI is an instruction-level fault injection framework which is developed using LLVM compiler infrastructure\cite{lattner2004,llvmwebref}.
%
VULFI targets LLVM's intermediate representation (IR) level instructions for fault injections.
%
VULFI is capable of targeting vector instructions (including architecture specific intrinsics) in addition to scalar instructions
for fault injections.
%
Visit VULFI's home page (\url{http://formalverification.cs.utah.edu/fmr/vulfi/}) and fill out the request form to obtain a copy
of VULFI.

\section{Requirements}
\label{prereq}

\subsection{Supported Languages}
\label{lang}
VULFI is capable of targeting any high-level program which can be compiled into LLVM IR. Currently, VULFI has been extensively tested
with C, C++, and ISPC programs.

\subsection{Software Requirements}
\label{softreq}
\begin{compactitem}
\item \textbf{OS} : Ubuntu 12.04 LTS.

 \item \textbf{LLVM} : VULFI has been developed using LLVM version 3.2. LLVM documentation:\\ \url{http://llvm.org/releases/3.2/docs/GettingStarted.html#getting-started}.

 \item \textbf{ISPC}: In case you are targeting ISPC programs for the fault injection, please use ISPC
 version 1.8.1 (download link: \url{http://sourceforge.net/projects/ispcmirror/files/})

 \item \textbf{Python}: Install Python version 2.7 or later.

 \item \textbf{OpenCV} (optional) : Install OpenCV version 2.3.1. Installation steps: \\
 \url{http://www.samontab.com/web/2012/06/installing-opencv-2-4-1-ubuntu-12-04-lts/}.
\end{compactitem}

%
% \section{Code Documentation}
% \label{docu}
% VULFI's code documentation is maintained at: \url{http://formalverification.cs.utah.edu/fmr/vulfi/doc/}.

\section{License \& Copyright Information}
\label{license}
Before using VULFI, please refer to the \texttt{LICENSE} file located in the topmost directory of the
distribution.
%
Please note that the VULFI distribution contains 3rd party programs which may have different license
and copyright requirements than VULFI.
%
You must always adhere to the individual license and copyright requirements of the 3rd party
programs in addition to the VULFI license and copyright requirements.
%
%
The license and copyright information for the 3rd party programs is provided in the \texttt{LICENSE} file.

% \subsection{VULFI License}
% \label{vulfi-license}
%
% \subsection{ISPC License}
% \label{ispc-license}
%
% \subsection{ISPC Benchmark License}
% \label{ispc-bench-license}
%
% \subsection{PARVEC Benchmark License}
% \label{parvec-license}

% \subsection{SCL Benchmark License}
% \label{scl-license}

\section{Installation}
\label{install}

\begin{enumerate}
 \item Change directory to the LLVM lib directory and copy VULFI source code to
 the current directory.
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd [LLVM_SRC_ROOT_DIR]/lib/
 \end{Verbatim}

 \item Change directory to the VULFI lib directory.
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd vulfi/lib
 \end{Verbatim}

 \item Run below commands to install a shared library (\texttt{relief.so}) containing VULFI implementation .
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
  $ make
  $ sudo make install
 \end{Verbatim}
\end{enumerate}

% \caption{An example C++ function \texttt{foo()}}
% \label{fig:foo}
% \vspace{-0.5cm}
% \end{figure}

% \subsection{OpenCV Installation}
% \label{opencv-install}
%
% \subsection{VULFI Installation}
% \label{vulfi-install}

\section{Usage}
\label{usage}
Performing fault injection using VULFI is a two-step process.
%
First, we instrument a target program using VULFI.
%
This step involves compiling the target program into an LLVM
bitcode file.
%
The target bitcode file is then linked with the VULFI's \texttt{Corrupt.bc}
bitcode file containing callable runtime fault injection APIs.
%
The resultant bitcode is then instrumented using the VULFI LLVM pass.
%
Second, we execute the resultant instrumented binary using python script
for performing fault injection.
%
The automation script provides a rich set of command line options
to effectively analyze the result of a fault injection run, and for generating a
detailed fault injection report.

\subsection{Step 1: Instrumenting a Target Program for Fault Injection}
Given that VULFI is implemented as an LLVM pass, a target program bitcode
is instrumented by running VULFI on it with \texttt{opt}
(\url{http://llvm.org/docs/WritingAnLLVMPass.html#running-a-pass-with-opt}).
Below are the command line options supported by VULFI:

\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm,commandchars=\\\{\}]
 -fn                                 Name of the functions (separated by a single
                                     white space) to be targeted for fault injections.
                                     Example: -fn "foo1 foo2 foo3"

 -fsa \{data,dint,dflo,ctrl,addr\}     Fault site selection algorithm; valid options:
                                     data: Target all pure-data fault sites for fault
                                     injection. Specifically, target all fault sites
                                     with its: (1) Lvalues of Integer or FloatingPoint
                                     types, and (2) forward program slice of the
                                     fault sites does contain an instruction
                                     with Lvalue of Pointer or Control type.

                                     dint: Same as "data" option except that target
                                     fault sites of only Integer types.

                                     dflo: Same as "data" option except that target
                                     fault sites of only FloatingPoint types.

                                     Note:- For Integer type, the supported
                                     bit-widths are 8, 16, 32, and 64. For
                                     FloatingPoint type, the supported
                                     formats are single (32-bit) and double
                                     precision (64-bit) floating-point values.

                                     ctrl: Target all fault sites of Integer or
                                     FloatingPoint type such that its forward program
                                     slice must contain atleast one contol instruction
                                     (e.g., "cmp" instruction).

                                     addr: Target all fault sites which either
                                     have its Lvalue of Pointer type or its
                                     forward program slice must contain atleast one
                                     instruction of Pointer type.
                                     (e.g., "getElementPtr" instruction).

 -rf \{0,1\}                           Range flag used to indicate if fault injection
                                     has to be done in a specific code section (using
                                     the source code line numbers). It is disabled
                                     (set to 0) by deafult.

 -flr                                Name of the csv file containing range information.
                                     This is an optional argument used only when -rf
                                     is enabled.

 -lang \{C,C++,ISPC\}                  The language option is used by VULFI to check for
                                     language specific information which may be useful
                                     for fault injection. For example, language
                                     specific metadata information.

 -arch \{x86,neon,nvvm,mips,spu\}      This architecture information is used when
                                     instrumenting architecture specific intrinsics
                                     for fault injection. This includes information
                                     such as whether an intrinsic uses mask. Currently,
                                     x86 is the only supported architecture; other
                                     architectures are just a placeholder and may
                                     be supported in future.

 -dbgf                               Name of the CSV file to which the list of static
                                     fault sites (instrumented for runtime fault
                                     injection) is written.
\end{Verbatim}

\subsection{An Example Instrumentation}
For a target C++ program file \texttt{foo.cpp}
having function definitions \texttt{foo()}
and \texttt{main()}, the first mandatory step is to
add following function prototypes at the beginning of the file:
%
%
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
extern int printFaultSitesData(void);
extern int printFaultInjectionData(void);
\end{Verbatim}
%
%
Next, add the below functions calls in the function \texttt{main()}
just before every return statements.
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
printFaultSitesData();
printFaultInjectionData();
\end{Verbatim}
%
%
%
The modified C++ program file \texttt{foo.cpp} is then
compiled into a bitcode \texttt{foo.bc} as shown below:
%
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ clang++ -emit-llvm foo.cpp -c -o foo.bc
 \end{Verbatim}
%
%
Similarly, the program file \texttt{Corrupt.C} located in the
directory \texttt{runtime/}, is compiled into a bitcode \texttt{Corrupt.bc}.
When using clang++ to compile the Corrupt.C file, add a macro definition
"\_\_CLANGPP" as described below. However, you must not use this macro when
using clang to compile the Corrupt.C file.

%
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
$ clang++ -emit-llvm -D__CLANGPP Corrupt.C -c -o Corrupt.bc
\end{Verbatim}
 %
 %
Next, both the bitcode files are linked using \texttt{llvm-link}:
%
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
$ llvm-link Corrupt.bc foo.bc -o foo-corrupt.bc
\end{Verbatim}
%
The resultant bitcode file is then instrumented by running VULFI pass
using \texttt{opt} as shown below:
%
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
$ opt -load /usr/local/lib/relief.so -vulfi  -fn "foo" -fsa "data" -lang "C++"
-dbgf "instruction_list_foo.csv" < foo-corrupt.bc > foo-inject.bc
\end{Verbatim}
%
%
Finally, the instrumented and the uninstrumented binaries (\texttt{foo-inject} and \texttt{foo}
respectively) are generated.
%
%
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
$ clang++  foo-inject.bc -o foo-inject

$ clang++  foo.bc -o foo
\end{Verbatim}
%
Both the binaries are executed under a given program input.
%
The execution output of \texttt{foo-inject} is compared with that of
\texttt{foo} to detect the occurrence of a silent data corruption in
the former execution output.
%

\subsection{Step 2: Automating Fault Injection using Python Script}
\label{cmdopt}

A binary instrumented using VULFI reads fault site selection and
fault injection related information from a configuration file.
%
The python script \texttt{driver.py} (located at \texttt{vulfi/scripts/})
automates the creation of the configuration file, running the instrumented
binary, and generating the fault injection report.
%
Therefore, we must always using the python script \texttt{driver.py} for
fault injection runs.
%
The command line options supported by the python script are listed
below:

 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm,commandchars=\\\{\}]
Usage: python driver.py [-h] [-e] [-s] [-c] [--ov] [--pv] [--iter EXEC_COUNT]
                 [--fcp {bex,img,num}] [--th TH] [--exec1 EXEC1] [--cmd1 CMD1]
                 [--out1 OUT1] [--cho] [--cfs] [--exec2 EXEC2] [--cmd2 CMD2]
                 [--out2 OUT2] [--rslt RSLT] [--ficsv FICSV] [--fia {cbr,abr}]
                 [--fib {1,2,3,4,5,6,7,8}] [--fid {msb,lsb}] [--fbu FBU]
                 [--fbl FBL] [--np NUM_PROB] [--dp DEN_PROB]
                 [--ff {eql,max,min,nlm}] [--fc FICOUNT]

Descriptions of the arguments:
-h, --help                Shows this help message and exits.

-e                        Execution mode.

-s                        Standalone mode, always used in combination with
                          -e option.

-c                        Comparison mode, always used in combination with
                          -e option.

--exec1 EXEC1             Location of the 1st executable, EXEC1, to be
                          executed with -s or -c mode.

--cmd1 CMD1               Command line args (if any) for the exe EXEC1.

--out1 OUT1               Name (including complete path) of the output file
                          generated by EXEC1.

--exec2 EXEC2             Location of the instrumented executable, EXEC2,
                          to be executed only with -c.

--cmd2 CMD2               Command line args (if any) for the exe EXEC2.

--out2 OUT2               Name (including complete path) of the output file
                          generated by EXEC2.

--ov                      Calculate overhead.

--pv                      Override user provided fault injection probability
                          and perform fault injection using a fault injection
                          probability of 1/N where N is the total number of
                          dynamic instructions.

--iter EXEC_COUNT         Number of iterations.

--fcp \{bex,img,num\}       Method used for comparing outputs of a fault-free
                          execution and a fault execution. Valid options for
                          comparing the outputs:
                          1. bex = Compare and match the outputs bit by bit.
                          2. img = If the outputs are image files, perform
                          PSNR based comparison.
                          3. num = If the outputs are of numeric type then
                          express the difference between the outputs
                          as L2Norm.

--th TH                   Threshold value; only valid with --fcp=img or num.
                          Report SDC if the difference in the outputs calculated
                          using "img" or "num" based method exceeds the threshold TH.

--cho                     This performance flag indicates that use cached output
                          for exec1 if available. This is disabled by default.
                          WARNING: Enable this option ONLY if you are sure that the
                          execution output of exec1 is expected to be the same
                          during every execution iteration.

--cfs                     This performance flag indicates that use cached fault
                          site count. Use it only when --cho flag is enabled.

--rslt RSLT               Name of the output file where result will be written.

--ficsv FICSV             Name of the input csv file where fault injection report
                           information will be written.

--fia \{cbr,abr\}            Fault injection algorithm, valid options:
                           1. cbr - cumulative byte ordering; target one or more
                           contiguous bytes for fault injection.
                           2. abr - absolute bit range;target a bit-range for fault
                           injection.

--fib \{1,2,3,4,5,6,7,8\}    Number of contiguous bytes to be considered for fault
                           injection.
                           Note: this option must be provided when using -fia=cbr

--fid \{msb,lsb\}            Direction from which no. of contiguous bytes to be
                           considered for fault injection.
                           Note:this option must be provided when using -fia=cbr.

--fbu FBU                  Upper bound for the fault injection bit-range.
                           Note:this option must be provided when using -fia=abr.

--fbl FBL                  Lower bound for the fault injection bit-range.
                           Note:this option must be provided when using -fia=abr.

--np NUM_PROB              Numerator of fault injection probability expressed as a
                           fraction.

--dp DEN_PROB              Denominator of fault injection probability expressed as
                           a fraction.

--ff\{eql,max,min,nlm\}      Puts restriction on the fault injection count;
                           valid options:
                           1. eql - continue injecting faults until fault injection
                           count becomes equal to the value passed using --fc.
                           2. max - continue injecting faults until fault injection
                           count is less than or equal to the value passed using
                           --fc.
                           3. nlm - inject faults in all eligible fault sites.

--fc FICOUNT               Fault injection countm this value is used in combination
                           with --ff.
 \end{Verbatim}

%  \subsection{An Example of Fault Injection Run using Python Script}

%   \subsection{Fault Injection Algorithm}
%  \label{fia}
%  \subsubsection{Absolute Bit Range (ABR)}
%  \subsubsection{Cumulative Byte Range (CBR)}

\section{Examples}
\label{examples}
\subsection{Example 1 : Scalar Benchmark - Quicksort }
\label{ex1-qsort}

\subsubsection{Step 1: Compile}
\begin{enumerate}
%
%
\item Change the directory to the source folder:
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd vulfi/benchmarks/scalar_benchmarks/src
 \end{Verbatim}
 %
 %
 \item Make sure that the paths provided in \texttt{common.mk} file
 are correct.
 %
 %
 \item Change directory to qsort folder and compile the program:
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd qsort/
 $ cd make
 \end{Verbatim}
 %
 %
 \end{enumerate}

\subsubsection{Step 2: Run}
\begin{enumerate}
 \item Change the directory to the run folder:
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd vulfi/benchmarks/scalar_benchmarks/run
 \end{Verbatim}
 %
 %
 \item Make sure that the paths provided in \texttt{common.mk} file
 are correct.
 %
 %
 \item Change the directory to the qsort folder. The \texttt{common\_qsort.mk} file
 contains the python command line used for performing the fault injection.
 The script file \texttt{run\_qsort.sh} launches the execution accepting
 a single argument to denote the number of fault injection runs to
 be carried out.

 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd qsort/
 $ ./run_qsort.sh 100
 \end{Verbatim}
 \end{enumerate}

 \subsubsection{Step 3: Result Analysis}

 \begin{enumerate}
  \item \textbf{Fault Sites:} Refer to the files \texttt{addr\_dbgData\_qsort\_addr.csv} and \\
  \texttt{ctrl\_dbgData\_qsort\_ctrl.csv} located at \texttt{\textbackslash{vulfi}\textbackslash{benchmarks}\textbackslash{scalar\_benchmarks}
  \textbackslash{src}\textbackslash{qsort}}. These files list the fault sites chosen for runtime fault injection
  for \texttt{address} and \texttt{control} fault site categories respectively. You may notice that there are no \texttt{data} fault
  sites found for the qsort program.

  \item \textbf{Fault Injection Result:} The fault injection result for each fault site category is saved in the directory
  \texttt{\textbackslash{vulfi}\textbackslash{benchmarks}\textbackslash{scalar\_benchmarks}\textbackslash{run}\textbackslash{qsort}\textbackslash{serial}}.
  For example, each row in the file \texttt{result\_qsort\_ctrl.csv} correspond to a unique fault injection run during which a fault was injected.
  The result file provides useful information such as instruction and bit position chosen at runtime for the fault injection. It also categorizes
  the result of the fault injection into one of the categories - \texttt{sdc, benign, crash,}  and \texttt{error}.

 \end{enumerate}

\subsection{Example 2 : ISPC Benchmark - Raytracing }
\label{ex2-rt}

The steps to run \texttt{Raytracing} benchmark is similar to the previous example. However, this benchmark requires
the ISPC compiler\cite{ispc12,ispcwebref}
to be present on the machine.
%
%
A precompiled version of the ISPC compiler which is known to work with the current version of VULFI could be downloaded from:
\url{https://www.dropbox.com/s/5qjilxsqho14ell/ispc.tar.gz?dl=0}.
%
%
Please make sure that you follow the license requirement of the ISPC compiler (\url{https://ispc.github.io/})
before downloading the precompiled binaries.
%
%

\subsubsection{Step 1: Compile}
\begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd vulfi/benchmarks/ispc_benchmarks/src/rt/
 $ make
 \end{Verbatim}

\subsubsection{Step 2: Run}
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
 $ cd vulfi/benchmarks/ispc_benchmarks/run/rt/
 $ ./run_rt.sh 5
 \end{Verbatim}

\subsubsection{Step 3: Result Analysis}
Figures~\ref{fig:original} and~\ref{fig:corrupted} are the outputs
generated by the \texttt{Raytracing} benchmark for a fault-free
and a faulty runs respectively.

\begin{enumerate}
 \item \textbf{Files listing fault sites: }
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
addr_dbgData_rt_avx2.csv
addr_dbgData_rt_sse4.csv
addr_dbgData_rt_ser.csv
ctrl_dbgData_rt_avx2.csv
ctrl_dbgData_rt_sse4.csv
ctrl_dbgData_rt_ser.csv
data_dbgData_rt_avx2.csv
data_dbgData_rt_sse4.csv
data_dbgData_rt_ser.csv
 \end{Verbatim}
 \item \textbf{Files with the fault injection results: }
 \begin{Verbatim}[fontsize=\relsize{-1},frame=single,framerule=0.1mm]
result_rt_avx2_addr.csv
 result_rt_sse4_addr.csv
 result_rt_ser_addr.csv
 result_rt_avx2_ctrl.csv
 result_rt_sse4_ctrl.csv
 result_rt_ser_ctrl.csv
 result_rt_avx2_data.csv
 result_rt_sse4_data.csv
 result_rt_ser_data.csv
 \end{Verbatim}

\end{enumerate}

\begin{figure}
  \begin{minipage}{0.45\textwidth}
  \centerline{\includegraphics[width=7.2cm,height=6.2cm]{image/teapot.png}}
  \caption{Golden Image}
  \label{fig:original}
  \end{minipage}
  \begin{minipage}{0.45\textwidth}
  \centerline{\includegraphics[width=7.2cm,height=6.2cm]{image/teapot-corrupt_fip_1_1000.png}}
  \caption{Corrupted Image}
  \label{fig:corrupted}
  \end{minipage}
\end{figure}

% \section{Extending VULFI}
% \label{extend}
%
% \subsection{The \texttt{SiteSelector} Class}
% \label{siteselector}
%
% \subsection{The \texttt{FaultInjector} Class}
% \label{faultinjector}
%
% \subsection{Runtime Library file \texttt{Corrupt.C}}
% \label{corruptlib}
%
% \subsection{Updating Intrinsic definitions}
% \label{intdef}

\bibliographystyle{plain}
\bibliography{biblio}

\end{document}


text:
\chapter*{Epilogue}
\addcontentsline{toc}{chapter}{Epilogue}
\thispagestyle{plain}

%Here is the place where you may want to tell a little story or a fairy tale
%which has some relevance for your thesis, such as ``Once upon a time, \dots''.
%The Epilogue is optional.

max_stars_repo_name: MatCast/master-thesis

max_issues_repo_name: MatCast/master-thesis

max_forks_repo_name: MatCast/master-thesis