stat-learning.tex

\documentclass[english,course,oneside]{Notesext}
%\documentclass[english,course,kindle]{Notesext}

% now for lecture 2022
% supervised by Carsten, Philipp, Tingying (thx!)

% lecture 2021 -> Carsten

%now for lecture 20202
%inverted classroom, so show one (sub)-chapter each or so per week


%Proposed new lecture series (2019)
% see https://docs.google.com/document/d/1kM1MlXbLuql9sTR7G7oHG7s1dtHrDfmlUwU5N8wUkyg/edit
%Lecture 01 (24.04): Introduction to Statistical Learning:
%Lecture 02 (08.05->30.04): (chapter 2) Basic Concepts of statistical learning, including linear regression, nearest neighbour, parametric vs non-parametric, bayesian classifier, curse of dimensionality, model accuracy, bias-variance trade-off
%Lecture 03 (08.05 -> Fabian remotely): (chapter 3) High-dimensional regression and shrinkage methods, including Ridge regression and PCA, Bayesian linear regression, Lasso regression, Bayesian interpretation of Lasso, Elastic Net, regularization parameter selection criteria.
%Lecture 04 (15.05  ->Benjamin): chap 3 continued
%Lecture 05 (22.05) Linear classifier: linear regression for classification, linear discriminant and quadratic discriminant analysis
%Lecture 06 (29.05): logistic regression + regularizations, nonlinear optimization and multi-class ext
%Lecture 07 (05.06 -> Tingying): Nonlinear classifier: Decision tree/Random Forest, Boosting?
%Lecture 08 (12.06 -> Tingying): Introduction to neural networks, separating hyperplane and deep learning, Tutorial on Keras.
%Lecture 09 (19.06): Convolutional neural network (Carsten’s lecture slides for TiCB, which includes explanations of neural network architecture, till mid May), still need to add loss function, optimization, tricks needed to get things to run, a word on frameworks (Keras or so), debugging: bias vs variance. Exercise: we can use Felix Nature methods paper, exercise prepared
%Lecture 10 (26.06): Unsupervised learning: Gaussian mixture model & k-means
%Lecture 11 (03.07): Dimension reduction, PCA, nonlinear dimred
%Lecture 12 (10.07): auto-encoder, (app denoising auto-encoder, Lukas slides, mid June), variational inference (Murphy), variational autoencoder, exercise from Lukas
%Lecture 13 (17.07): outlook “newer topics” without much theory: Generative Adversarial Networks: GAN, image-to-image translation, cycle GAN. A good online exercise; Recurrent Neural Networks (RNN/LSTM). Graph Convolutional Networks (GCN).


% see http://poisson.phc.unipi.it/~maggiolo/index.php/2008/12/latex-class-for-lecture-notes/
%course, seminar or talk: a course is a medium-length document structured in sections, subsections and paragraphs; a seminar is a short document without structure or with subsections; a talk is the sheets you take with you when you give a seminar.
%
%The second option defines the style of the document; if you want to tweak some parameters, there are also these fine-tuning options:
%headertitle, headersection, headersubsection or headerno: what to display on the top of the pages;
%theoremnosection, theoremsection or theoremsubsection: numbering of theorems, definitions, etc.;
%cleardoublepage or nocleardoublepage: whether we want empty double pages after a section ending;
%oneside or twoside: margins and headers optimized for one-sided or two-sided printing;
%onecolumn or twocolumn: how many columns per page.
%
%The class defines also some commands:
%for course, the command \lecture{duration}{dd}{mm}{yyyy} which writes on the margin information about the lesson?s number, duration and date;
%for talk, the commands \separator, which simply draws a horizontal line, and \tosay{message} which prints the message inside a box (useful to remember things you don?t have to write on the blackboard);
%\mymarginpar{message}, which behaves like a regular \marginpar but with a custom style

\usepackage[nomytheorems,pagesaving]{mymacros}
%\newcommand{\Rp}{\R_{\geq0}} % nervt
\newcommand{\Rp}{\R_+} % nervt
\newtheorem{question}[theorem]{Question}

% Graphical models chapter (macros not in mymacros.sty)
\newcommand{\bE}{\mathbb{E}}
\newcommand{\cG}{\mathcal{G}}
\newcommand{\cE}{\mathcal{E}}
\newcommand{\cM}{\mathcal{M}}
\newcommand{\cS}{\mathcal{S}}
\newcommand{\cN}{\mathcal{N}}
\newcommand{\indep}{\perp\!\!\!\perp}
\newcommand{\dep}{\not\indep}
\newcommand{\inlineQ}[1]{\par\noindent\textit{#1}\par}
\newcommand{\pa}{\mathrm{pa}}
\newcommand{\paG}[1]{\mathrm{pa}_{\cG}(#1)}
\newcommand{\chG}[1]{\mathrm{ch}_{\cG}(#1)}
\newcommand{\deG}[1]{\mathrm{de}_{\cG}(#1)}
\newcommand{\ndG}[1]{\mathrm{nd}_{\cG}(#1)}

\usepackage{verbatim}

% coole arrows, adapted from http://tex.stackexchange.com/questions/38670/broken-arrow-symbol
\newcommand*{\DashedArrow}[1][]{\mathbin{\tikz [baseline=-0.15ex,-latex, densely dashed,#1] \draw [#1] (0pt,0.5ex) -- (1.3em,0.5ex);}}%
\newcommand*{\InhibArrow}{\DashedArrow[-|]}
\newcommand{\Data}{\mathcal{D}}
\newcommand{\dims}[1]{{n_{#1}}}
%\newtheorem{algorithm}[theorem]{Algorithm}
%\usepackage[figure]{algorithm2e}
\usepackage[linesnumbered,ruled]{algorithm2e}



% equation environment
\newcommand{\eq}[1]{\begin{align} #1 \end{align}}

\usepackage{cite}
\usepackage[T1]{fontenc}

% Skip figures missing from the public source export instead of failing the build.
\let\origincludegraphics\includegraphics
\newcommand{\maybeincludegraphics}[2][]{%
  \IfFileExists{#2.pdf}{\origincludegraphics[#1]{#2}}{%
  \IfFileExists{#2.png}{\origincludegraphics[#1]{#2}}{%
  \IfFileExists{#2.jpg}{\origincludegraphics[#1]{#2}}{%
  \IfFileExists{#2.eps}{\origincludegraphics[#1]{#2}}{%
  \IfFileExists{#2}{\origincludegraphics[#1]{#2}}{%
  \fbox{\parbox[c][2cm][c]{0.9\linewidth}{\centering\scriptsize Figure omitted from source release.}}%
  }}}}}%
}
\let\includegraphics\maybeincludegraphics



\title{MA4802: Statistical learning}
\subject{MA4802: Statistical learning}
%Statistical Learning - A 2+2SWS course for the Data Science Master program
% Fabian Theis
%Target audience: master students from math or computer science
%~14 slots
% https://www-m12.ma.tum.de/teaching/summer-semester-2017/statistical-learning-ma4802/
\author{Fabian J. Theis}
\email{fabian.theis@helmholtz-muenchen.de, {www.helmholtz-munich.de/en/icb/research-groups/theis-lab}}
%\speaker{SPEAKER}
\date{17}{4}{2024}
\dateend{17}{7}{2024}
\place{Department of Mathematics, TU Munich}

\begin{document}

\section{Introduction}
%\section*{Abstract}
Statistical learning aims to find a relationship between a set of variables denoted as features, realized by a series of samples. In the supervised learning case these samples are labeled, and the goal is to predict the label an unlabeled feature sample correctly. Statistical learning extends classical concepts in statistics such as regression by introducing extensions for a series of models ranging from local approximations over sparse coding to graphical modeling, density estimation and clustering methods. In this course we want to review a series of basic models from a probabilistic perspective. The lecture draws heavily on two textbooks by Hastie \& Tibshirani \& Friedman \cite{Hastie:2009wp} and Murphy \cite{Murphy:2012uq}, both great reads.

Please note, these lecture notes are for participants of this course only, since some material may be be protected by copy-right and may not distributed without asking permission.

For an example-based introduction on how statistical learning is generally important, in particular in the growing field of data sciences, see presentations slides online (\url{http://moodle.tum.de}) and the 2015 Science special section on machine learning/AI \cite{Stajic:2015jo}.

%\section{A short primer in molecular cell biology}
%Planned to do this, but maybe too detailed. Please read up on it online (\url{http://moodle.tum.de}) if you are interested.

%\newpage

\section{Basic concepts in statistical learning} % 1-2 slots

\subsection{Types of statistical/machine learning} % supervised/unsupervised, some examples

Statistical learning aims to learn a relationship between a multivariate set of random variables $X=(X_1,\ldots, X_p)$, typically from $\R$.

In \defn{unsupervised learning}, given a set of samples $D=\{x_i\in\R^p|i=1,\ldots, N\}$ of $X$, one seeks to estimate `interesting structure´ in $X$ from $D$ such as clusters or geometrical relationships.

In the more well-defined \defn{supervised learning},
we additionally observe a \defn{response} or (classically) \defn{dependent variable} $Y$.
The goal is to explain it by $X$, then also known as \defn{predictors}, \defn{features} or (classically) \defn{independent variables}.

$Y$ may be `quantitative´, i.e. taking real values, or `qualitative´ i.e. ordinal or categorical (Q: examples); most important example of the latter is the binary case $\{0,1\}$.
Depending on the output type, supervised learning leads to \defn{regression} (quantitative $Y$) and \defn{classification} or \defn{pattern recognition} (qualitative $Y$)\footnote{of course there is logistic regression and more generally ordinal regression and other models not precisely adhering to this simplified notion}.

Probabilistically, the input-output relation is measured by $P(X,Y)$.
The algorithms are called `learning´ because we want to estimate the relation of $Y$ to $X$ from so called \defn{training data} i.e. iid\footnote{in the simplest case; there exists extensions to training data with temporal or spatial structure} samples $D=\{(x_i,y_i)\in\R^p\times\R|i=1,\ldots, N\}$ of $(X,Y)$.

There are other types that mix the distinction of training on labels versus blind pattern mining such as reinforcement learning or semi-supervised learning; these however are beyond the scope of the lecture.



\subsection{Examples of supervised learning} % supervised/

Assume an unknown functional relationship
\begin{equation}\label{eqn:YfX}
Y=f(X)+\epsilon
\end{equation}
between input and output, with the random variable (RV) $\epsilon$ being independent of $X$ and having zero-mean. The goal of learning is to estimate $f$ from $D$ (\defn{function approximation}) and then make predictions $\hat y=\hat f(x)$ for some new input $x$.

\subsubsection{Need for probabilistic predictions}
To quantify uncertainty in both the model and the estimation, we determine the probability distribution over the output
$P(Y=y|X=x,D)$ for a given input $X=x\in\R^p$, where we implicitly condition on the used model (constraints on $f$).

To get to a deterministic output, the mode
\begin{equation}\label{eqn:MAP}
\hat y=\hat f(x)=\argmax_c P(Y=c|X=x,D)
\end{equation}
also known as the \defn{MAP (maximum a posteriori)} estimate is often used. Alternatively the \defn{posterior mean} $E_{Y|X}[P(Y|X=x,D)]$
may be determined as a summary statistics.
In the following, we will usually abbreviate $P(Y=y|X=x,D)\equiv P(y|x,D)$ and drop dependency on the data $D$.

%Given an estimate $\hat f$ together with a prediction $\hat y=\hat f(x)$, we can ask about its accuracy. Assuming fixed $\hat f$ and $x$, we can decompose the mean square error as follows:
%$$
%MSE(\hat y)=
%$$


\subsubsection{Linear regression} %or move up
%More explicitly we can write down the model as
%$$p(Y|X,\theta)=N(y|\mu(X), \sigma^2(X))$$
Consider the setting of quantitative supervised learning, where lables $y$ have actual values.
Assume that $\epsilon$ in the above model \eqnref{eqn:YfX} is a
Gaussian with zero-mean and fixed variance $\sigma^2$, then
\eq{\label{eqn:gaussassump}
P(y|x,\beta)=\NNN(y|f(x), \sigma^2),
}
where $\NNN(y|f(x), \sigma^2)$ denotes the normal distribution with mean $f(x)$
and variance $\sigma^2$. Hence the MAP estimate becomes $\hat y = f(x)$, which in this case also coincides with the posterior mean.

\begin{exercise}Why?\end{exercise}

In the simplest case, we assume $f$ is affine linear, so
\eq{
Y=f(X)+\epsilon= X^\top \beta+\epsilon , \quad \beta \in \mathbb{R}^{p+1},
}
for a vector of weights $\beta$ and $x_0:=1$.
Or without vector notation
\eq{
Y=\beta_0 + \sum_{i=1}^p X_i\beta_i +\epsilon .
}
%So
%\eq{
%P(y|x,\beta)=\NNN(y|x^\top \beta,\sigma^2).
%}

The parameters $\beta$ are estimated from $D$ for example by least squares: minimize the residual sum of squares
\eq{
\RSS(\beta)=\sum_{i=1}^N (y_i-x_i^\top \beta)^2.
}
Note that this is independent of the magnitude of the fixed variance $\sigma^2$. A probabilistic motivation of RSS is given in section \ref{sec:statdec}.
With a $(N\times p)$-dimensional sample matrix $\X$ (i.e., observations $x$
of $X$ are stored as rows) and a $N$-dimensional sample vector $\y$, rewrite this as
\eq{
\RSS(\beta) = ||\y - \X \beta||^2 = (\y-\X \beta)^\top(\y-\X \beta).
}
Differentiation gives the \defn{normal equations}
$$\X^\top(\y-\X \beta)=0$$
and for nonsingular $\X^\top\X$ the unique solution
$$\beta=(\X^\top\X)^{-1}\X^\top\y.
$$
See example Figure \ref{fig:linclass}, more details in the exercises.

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/intro-linclass}
\caption{Linear regression used for classification of samples with Blue=0, Orange=1, separated at $x^\top\beta=0.5$.
From \cite{Hastie:2009wp}, Fig. 2.1}
\label{fig:linclass}
\end{figure}

%TODO: or slides? also in terms of classification, with $P(y|x,\theta)>0.5$ for class selection?
%how to show in blackboard lecture?

\subsubsection{Logistic regression} %or move up
For qualitative learning, in particular binary classification, we replace the Gaussian distribution for $Y$ in equation \eqnref{eqn:gaussassump} by a discrete distribution, the Bernoulli distribution:
\eq{
P(y|x)=\Ber(y|\mu(x))
}
with the mean being a function $\mu(X=x):R^p\rightarrow [0,1]$.
For the mean of a Bernoulli distribution holds $\mu(x)=E(Y|X=x)=P(Y=1|X=x)$.
To ensure that $\mu(x)\in[0,1]$, but still use a linear model, we feed the linear model into the sigmoidal function
$$\sigm(x):=(1+\exp(-x))^{-1}.$$
The model then reads
$$P(y|x,\beta)=\Ber(y|\sigm(x^\top \beta)).$$
with the $p$ dimensional parameter vector $\beta$.
The MAP estimator then simply corresponds to
$$\hat y(x)=1 \qquad \Leftrightarrow \qquad P(y=1|x,\beta)>1/2$$
We will discuss estimation and regularization in logistic regression later.

%\subsubsection{$k$-nearest neighbor regression/classification} %or move up
\subsubsection{Nearest neighbor methods} %or move up
Instead of the above global, linear model, consider the following local  method: Determine $P(Y=c|x,D)$ from the average of close-by samples in your data $D$. The \defn{$k$-nearest neighbor} (kNN) fit is then defined as
$$\hat y(x) = \frac{1}{k} \sum_{x_i \in N_k(x)} y_i$$
with $N_k(x)$ being the $k$-closest samples of $D$ to $x$ in some metric,
or as probabilistic model
$$P(Y=c|x,D,k)= \frac{1}{k} \sum_{x_i \in N_k(x)} \unitm_{y_i}(c),$$
with the indicator function $\unitm_{y_i}(c)=1$ if ${y_i}=c$ and $0$ otherwise.
Again take decision for classification at $P(y|x,D,k)>0.5$, which is equivalent to majority voting in the class.

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/intro-15nearestneigbor}
\caption{$k$-nearest neighbor classification with $k=15$.
From \cite{Hastie:2009wp}, Fig 2.2.}
\label{fig:nnclass}
\end{figure}

For an example of this method, see Figure \ref{fig:nnclass}. In the exercises, we will see the effect of varying $k$'s, as well as a comparison to linear classification.

Many models in statistical learning are a mixture of the above two models, e.g. kernel methods replacing the $0/1$-weights of $k$-nearest neighbor, local regression, linear models with basis expansion for higher flexibility or neural networks as iterations of sigmoidal models.

\subsection{Statistical decision theory}\label{sec:statdec}
%TODO: a bit double to prob formulation in 2.2, remove redundancy

We want to predict $Y$ from $X$, i.e. seek a function $f(X)$ as predictor of $Y$. To quantify prediction preformance between predicted labels $y$ and true labels $y'$, we use a loss function $L(y,y')$, which is positive definite ($L(y,y')\geq 0$), e.g.~the square loss $L(y,y'):=(y-y')^2$.
We can then choose $f$ by minimizing the (squared) \defn{expected prediction error}
$$EPE(f):= E( L (Y,f(X))) = E(Y-f(X))^2 = \int (y-f(x))^2 P(dx,dy).
$$
By conditioning we get
$$EPE(f) = E_X E_{Y|X} (L(Y,f(X))|X),$$
and % because of definiteness
it suffices to minimize $EPE(f)$ pointwise (i.e.~for a single point $X=x$):
\eq{
f(x)=\argmin_c E_{Y|X} (L(Y,c)|X=x).
}

For the square loss function, differentiation yields the conditional expectation
$$f(x)=E_{Y|X} (Y|X=x)=E(Y|X=x),$$
also known as regression function.

In words: the best prediction of $Y$ at $X=x$ is the conditional mean if we measure error by $L_2$.

Its empirical estimate is realized by nearest-neighbor methods: instead of averaging over all points at $x$ with samples $x_i=x$, we approximate this in a neighborhood and use averaging instead of expectation:
$$\hat f(x)=\frac{1}{k} \sum \{y_i|x_i \in N_k(x)\}$$

Note: Under mild regularity conditions on $P(X,Y)$, $\hat f(x)\to E(Y|X=x)$ for $N,k\to\infty$ and $k/N\to0$.

However, this is not a general approximator, because often $N$ is too small, in particular in high dimensions.

Let's instead take additional assumptions (regularization) e.g. by a model based approach in a linear approximation:
$$f(x)= x^\top \beta$$
Putting this into the $EPE$ definition we get
$$\beta=(E(XX^\top))^{-1}E(XY)$$

Note: The least-squares estimator introduced above replaces these expectations by averages.

So least-squares $f(x)$ is approximated by globally linear function, whereas in nearest-neighbor by a locally constant function. Extensions are e.g.
$$f(X)=\sum f_j(X_j)$$ with univariate $k$-nearest neighbor.


Other loss functions:
Replacing $L_2$ loss by $L_1$ loss $L(y,y')=|y-y'|$, leading to
$$\hat f(x)=\median (Y|X=x)$$
For discrete $Y$, one may choose the $0-1$ loss function
$$L_{01}(y,y'):=\left\{\begin{array}{cc}0 & y=y'\\1&y\neq y'\end{array}\right.$$
Then minimal $EPE$ leads to
$$f(x)=\argmin_c E_{Y|X} (L(Y,c)|X=x)=
\argmin_c \sum_{y\in \Dom(Y)} L(y,c) P(Y=y|X=x)
$$
which equals
\eq{\label{eqn:BayesClassifier}
f(x)=
\argmin_c (1- P(Y=c|X=x))=\argmax_c P(Y=c|X=x)
}
for the $0-1$ loss function, which is known as the \defn{Bayes classifier}: choose the most probable class using the conditional discrete $P(Y|X)$.

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/intro-optimalbayes}
\caption{Optimal Bayes classifier, using the known generating densities (mixture of Gaussians).
From \cite{Hastie:2009wp}, Fig 2.2.}
\label{fig:optimbayes}
\end{figure}

Again this is approximated by kNN; for binary classification this implies simply thresholding at $0.5$, see figure \ref{fig:optimbayes}.


\subsection{Parametric vs non-parametric models}

A \defn{statistical model} $H$ is a set of probability distributions on a common set, e.g. $\R^p$. A \defn{parametric model} is parametrized by a finite number of parameters.

\begin{example}
Consider the parametric model $H=\{\NNN(\mu, 1) | \mu\in\R\}$. Given iid samples $x_1,\ldots,x_N$, its parameter $\mu$ is estimated by $\hat \mu=\frac{1}{N}\sum x_i$.
\end{example}

\begin{example}
Consider the nonparametric model $H=\{P | \var P<\infty\}$. Given iid samples $x_1,\ldots,x_N$, the mean again is estimated by $\hat \mu=\frac{1}{N}\sum x_i$.
\end{example}

Linear and logistic regression are parametric learning models, whereas $k$-nearest neighbor is non-parametric.


\subsection{The curse of dimensionality}

Although kNN classifiers may work well for many samples, they have serious issues in high dimensions, due to a fact known as the \defn{curse of dimensionality}:
In high dimensions, samples are really sparse (i.e. far from each other).

To be more precise:
Consider input samples $x_i$ uniformly distributed in the unit hypercube $[0,1]^p$. If we want to observe a fraction $r$ of all samples within a neighborhood, it has to be  (in average) of size $e_p(r):=r^{1/p}$.
While in one dimension $e_1(r)=r$, for 10 dimensions we have to cover 80\% of the total volume to cover 10\% of all samples, $e_{10}(0.1)=0.80$. Even if we only want to capture 1\%, the size is still $e_{10}(0.01)=0.63$. Hence the `neighborhood´ is not really local anymore (see Figure \ref{fig:cursedim}).

\begin{figure}\centering
\includegraphics[width=\columnwidth]{images/intro-cursedim}
\caption{The curse of dimensionality. In $p=10$ dimensions we need to cover 80\% of the range of each coordinate to capture 10\% of the data.
From \cite{Hastie:2009wp}, Fig 2.6.}
\label{fig:cursedim}
\end{figure}

Similarly, one can show that the all sample points are close to the edges of the sampling space in high dimensions:
Assume that $x_i$ are uniformly samples from the unit ball $\{x| \|x\|<1\}$ centered at $0$. The median distance from $0$ to the closest sample is
\eq{
d(p,N)=\left(1-\frac{1}{2}^{1/N}\right)^{1/p}
}
E.g. $d(10,500)\approx 0.52$, so most data points are closer to the boundary than to $0$. This implies that we mostly face the situation of estimation at the edges of the sampling space, where extrapolation is more difficult.

\subsection{Assessing model accuracy and the bias-variance trade-off}
%model selection and bias variance, see \cite{Hastie:2009wp} 2.9

All (most) models contain a complexity parameter, such as the number of basis functions used in a regression or the number $k$ of nearest neighbors in kNN. What is the optimal model complexity for the given data?

%Consider the $k$-nearest neighbor regression $\hat f_k(x)$ as function of $k$.
Assume that the data arise from model $Y=f(X)+\epsilon$ with independent noise $\epsilon$ with $E(\epsilon)=0$ and $\var(\epsilon)=\sigma^2$, and assume a priori fixed training samples $x_i$. The EPE or test error can then be decomposed into

\begin{eqnarray*}
\EPE(\hat f)(x)&=&E_{Y|X} ((Y-\hat f(x))^2|X=x) \\
&=& \sigma^2 + \underbrace{E_D\left(\hat f(x)-f(x)\right)^2}_ {\MSE(\hat f(x))}\\
&=&\sigma^2 + E_D(\hat f(x)^2)
-\left(E_D(\hat f(x))\right)^2
+\left(E_D(\hat f(x))\right)^2\\
&&\qquad-2E_D(\hat f(x))f(x)+f(x)^2\\
&=& \sigma^2 + \underbrace{\left( E_D(\hat f(x))-f(x)\right)^2}_{\bias(\hat f(x))^2} + \underbrace{E_D \left( \hat f(x) - E_D(\hat f(x))\right)^2}_{\var(\hat f(x))},
\end{eqnarray*}
where we used that the noise is independent and centered at 0 in the second line, and took expectations $E_D$ over the training set $D$.
The first term is called the \defn{irreducible error}, the variance of the target, which we cannot control. In statistical learning, we aim at minimizing the two right terms, the mean square error of the estimator $\hat f(x)$ of $f(x)$.

If we for simplicity assume that the values of $x_i$ in the samples have been fixed in advance i.e. are nonrandom, then for the $k$-nearest neighbor regression $\hat f_k(x)$ as function of $k$ we get
$$\EPE(\hat f)(x) =
\sigma^2 + \left(f(x) - \frac{1}{k} \sum_{x_i \in N_k(x)} f(x_i) \right)^2 + \frac{\sigma^2}{k}\\
$$
For increasing $k$ the bias term will likely increase due to samples being further off in larger neighborhoods, thereby the average also being more off. The variance will decrease with $k$. So as $k$ varies there is a trade off between bias and variance, see figure \ref{fig:biasvariance}.

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/intro-biasvariance}
\caption{Bias variance trade-off.
From \cite{Hastie:2009wp}, Fig 2.11.}
\label{fig:biasvariance}
\end{figure}

We will learn techniques, e.g. cross validation, to control for overfitting and achieve generalization.

\subsection{No free lunch}
G Box: All models are wrong, but some are useful.

Which model to choose? Cross validation may help empirically, but there is no universal best model, as stated in the so called \defn{no free lunch theorem}. This is based on the fact that while some assumptions work for one domain, they may work poorly in another, which can be exploited to construct a set of samples that a predictor minimizing $\EPE$ for another set of samples misses. For more details see \cite{ShalevShwartz:2014wp}

%Wolpert 1996, see also math basics of ML

%c.	some basics from prob theory/stat/optimization or so needed later?


\newpage
\part{Supervised learning}

\section{High-dimensional regression}

\subsection{Review of linear regression, estimation, testing and covariate selection} %1 slot

As seen before, linear regression is modeled by
\eq{\label{eqn:linreg}
P(y|x,\theta)=\NNN(y|x^\top\beta, \sigma^2),
}
Here as in the following we summarize all model parameters in the parameter vector $\theta$.

The model can describe nonlinear relationships by including nonlinear functions of the input random vector, e.g. basis expansions $X, X^2, X^3$ of a single variable $X$, interaction variables $X_1X_2$ or factor coding of a categorical variable. Multidimensional output $Y$ can similarly be handled but is often split into two models.

\subsubsection{Estimation}
The model parameters $\theta$ are most commonly estimated by \defn{maximum likelihood}, i.e. by
\eq{
\hat\theta := \argmax_\theta P(D|\theta)
}
Equivalently we may minimize the negative log likelihood
\eq{
l(\theta) := -\log P(D|\theta) = -\sum_{i=1}^N \log p(y_i|x_i,\theta)
}
for iid samples $D$.

MLE for the linear regression model \eqnref{eqn:linreg} yields
\eq{
l(\theta) &= -\sum_{i=1}^N \log\left(
(2\pi\sigma^2)^{-1/2} \exp\left(-(2\sigma)^{-1}(y_i-x_i^\top\beta)^2\right)
\right)\\
&=(2\sigma)^{-1} RSS(\beta)- (N/2) \log(2\pi\sigma^2)
}
Hence for minimal $l(\theta)$, $RSS(\beta)$ needs to be minimized, which implies
$$\hat\beta=(\X^\top\X)^{-1}\X^\top\y$$
for nonsingular $\X^\top\X$ as seen above.
The variance of the estimator $\hat\beta$ follows from the likelihood as
$$\var(\hat\beta)=(\X^\top\X)^{-1}\sigma^2.$$
The MLE estimate for the variance $\sigma^2$ results in
$$\hat\sigma^2=\frac{1}{N}\sum_{i=1}^N (y_i-\hat y_i)^2$$
with $\hat y_i:=x_i^\top\hat\beta$, which may be replaced by
$$\hat\sigma^2=\frac{1}{N-p-1}\sum_{i=1}^N (y_i-\hat y_i)^2$$
to guarantee an unbiased estimator.


\subsubsection{Testing}

To draw inference about parameters and model, we assume that model \eqnref{eqn:linreg} holds with true parameters $(\beta, \sigma)$.
So
$$\hat\beta \sim \NNN(\beta, (\X^\top\X)^{-1}\sigma^2)$$
and
$$(N-p-1)\hat\sigma^2 \sim \sigma^2 \chi^2_{N-p-1},$$
which are both independent. Hence we can test if a particular coefficient $\beta_j \neq 0$ by considering the Z-score
$$
z_j=\frac{\hat\beta_j}{\hat\sigma\sqrt{((\X^\top\X)^{-1})_{jj}}}
$$
It is distributed as $t_{N-p-1}$.

\subsubsection{Optimality of the estimator}
Theorem (Gauss-Markov):
The least-squares estimator has the smallest variance among all linear unbiased estimators.

\begin{exercise}
Proof? (idea:
$\tilde\beta=((\X^\top\X)^{-1}\X^\top+\A)\y$. No bias, hence $\A\X=0$, so $\var\tilde\beta=\var\hat\beta+\sigma^2\A\A^\top$, i.e. $\A=0$ if variance is to be minimal)
\end{exercise}

Consider the mean-square error of an estimator $\tilde\beta$ of $\beta$. If it is unbiased, its MSE equals the variance
$$\MSE(\tilde\beta)= E(\tilde\beta-\beta)^2=\var(\tilde\beta)$$
According to GM, $\hat\beta=(\X^\top\X)^{-1}\X^\top\y$ has the lowest variance and hence the lowest MSE.

However, there may be biased estimators with lower MSE, as we will see in the following.

%TODO: subset selection?? -> exercises


\subsection{Shrinkage methods and ridge regression}% ridge regression, Lasso, elastic net)} %1-2 slots

In this section, we describe `robust regression' methods, which trade-off bias for lower variance, which, in particular, in high dimensions or noisy data, may be preferable. These estimators may be derived by modifying the noise model (i.e. the likelihood) or putting prior constraints $P(\theta)$ on the parameters and doing MAP estimation
$$\hat\theta := \argmax_\theta P(D|\theta)P(\theta)
$$

Some common model assumptions are summarized in this table (see \cite{Murphy:2012uq}, p. 224):
\begin{center}
\begin{tabular}{ l | l | l }
Likelihood & Prior & Name\\
\hline
Gaussian & Flat & Least squares\\
Gaussian & Gaussian & Ridge regression\\
Gaussian & Laplace & Lasso regression\\
Laplace or Student & Flat & Robust regression
\end{tabular}
\end{center}

\subsubsection{Ridge regression}
In situations in which we face an ill conditioned regression problem, the variance of the ordinary least squares (OLS) estimates $\hat{\beta}$ might become very large - especially when we face {\it multicoliniarity} in the covariats, i.e., where one or more covariats are highly linearly related, or the regression is overparameterized. This can be easily seen by realizing that $\X^T\X$ is a matrix with $1$ on the diagonal and the pairwise Pearson correlation coefficient at the off diagonal:
\begin{equation*}
(X^TX) = \begin{pmatrix}
1 & \rho_{1,2} & \cdots &\rho_{1,m}\\
\rho_{2,1} &1 & \cdots & \rho_{2,m} \\
\vdots  & \vdots &\ddots & \vdots\\
\rho_{m,1} & \rho_{m,2} & \cdots & 1
\end{pmatrix}.
\end{equation*}
By increasing the off-diagonals (through highly correlated covariates) the inverse will become large, which has an immediate effect on the variance of the estimates since:
\begin{equation*}
\text{var}(\beta) = (\X^T\X)^{-^1}\sigma^2.
\end{equation*}

In nearly perfect multicolinearity, $\X^T\X$ might become singular, leading to a non-invertible matrix. In such a case the OLS estimate does not even exist!

Ridge regression was originally motivated by finding a {\it biased} estimator of $\beta$ for singular covariate matrices. It was derived by realizing that a nonsingular matrix can always be recovered by simply adding a sufficiently small positive value $\lambda > 0$ to the diagonal of the singular $\X^T\X$. This results in the following augmented closed form solution of $\hat{\beta}^{\mathtt{ridge}}$:
\begin{equation}
\hat{\beta}^{\mathtt{ridge}} = \left( \X^T\X + \lambda\mathbf{I}\right)^{-1}\X^T\y,
\end{equation}

which is a solution to the following unconstrained
\begin{equation}
\label{FRR}
	\hat{\beta}^{\mathtt{ridge}} =  \argmin_\beta  \left[  \sum_{i=1}^N (y_i - x_i^\top\beta )^2 + \lambda \underbrace{\sum_{j=1}^p \beta_j^2}_{||\beta||_2}   \right]
\end{equation}

The unconstrained optimization problem can be equivalently expressed as a constrained optimization problem, which will give us a better intuition of the effect the added term to the RSS has on $\beta$. The constrained optimization problem is stated as follows:
\begin{align}
\label{CO_RR}
\begin{split}
\hat{\beta}^{\mathtt{ridge}} =&  \argmin_\beta  \sum_{i=1}^N (y_i - x_i^\top\beta )^2\\
&\mathtt{s.t.} ~~||\beta||_2 \leq t
\end{split}
\end{align}
The two optimization problems are connected.
This nicely illustrates that ridge regression constraints all coefficients (except the intercept $\beta_0$) to an upper limit $t$. The  $\ell_2$ norm constraint defines a ball who's radius is constrained by $t$. The solution to the optimization problem will lie on the intersection of the constrained ball and the RSS estimate (Figure \ref{fig:est2b}).

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/est2}
\caption{Estimation picture for the Lasso (left) and ridge regression
(right). Shown are contours of the error and constraint functions. The solid blue
areas are the constraint regions $|\beta_1| + |\beta_2| \leq t$ and $\beta_1^2 + \beta_2^2 \leq t^2$ respectively, while the red ellipses are the contours of the least squares error function.
From \cite{Hastie:2009wp}, Fig 3.11.}
\label{fig:est2b}
\end{figure}

There is a one-to-one correspondence between the parameters $\lambda$ in \eqref{FRR} and $t$ in \eqref{CO_RR}. The complexity parameter $\lambda \geq 0$ controls the amount of shrinkage: the larger the value $\lambda$ (the smaller the value $t$) the greater the amount of shrinkage and hence bias-variance trade-off. The problem \eqref{FRR} is sometimes referred to as {\it Tikhonov regularization} or {\it penalized least square.}

Note that the ridge regression estimator is related to the classical OLS estimator, $\hat \beta^{\mathtt{ls}},$ in the following manner
\[
  \hat \beta^{\mathtt{ridge}}  = (\unitm + \lambda (\X^T \X)^{-1})^{-1} \hat \beta^{\mathtt{ls}},
\]
assuming that $\X^T \X$ is non-singular. This relationship can be verified by simply applying the definition of $\hat \beta^{\mathtt{ls}}$ and using the fact $\B^{-1} \A^{-1} = (\A\B)^{-1}$

In the case of orthonormal inputs, the ridge estimates are just a scaled version of the least squares estimates, i.e., $\hat \beta^{\mathtt{ridge}}  = \frac{1}{1+\lambda} \hat \beta^{\mathtt{ls}}.$

This is instructive, as it shows that, in this simple case, the ridge estimator is simply a downweighted
version of the least square estimator.


\subsubsection{Ridge Regression and PCA*}

The {\it singular value decomposition} (SVD) gives us some additional insight into the nature of ridge regression. Let $$\X = \U \D \V^T$$ be the SVD of $N \times p$ matrix $\X,$ where $\U$ and $\V$ are orthogonal matrices of order $N \times p$ and $p \times p$ respectively; and $\D$ is a diagonal matrix of order $p \times p$ containing the {\it singular values} of $\X$. The columns of $\U$ are called the {\it left singular} vectors, span the column space of $\X$; whereas the columns f $V$, which are called the {\it right singular} vectors, space the row space of $X$. In particular, $\U$ is a set of eigenvectors for $\X \X^T,$ and $\V$ is a set of eigenvectors for $\X^T \X.$ The non-zero singular values of $\X$ are the square roots  of the eigenvalues of both $\X \X^T$ and $\X^T \X.$

Using SVD we can write the least squares fitted vector as
\[
 \X \hat \beta^{\mathtt{ls}} = \X (\X^T \X)^{-1} \X^T \y = \U \U^T \y.
\]
Note that $\U^T \y$ are the coordinates of $\y$ with respect to the orthonormal basis $\U.$ The ridge solution is given as
\[
 \X \hat \beta^{\mathtt{ridge}}  =  \X ( (\X^T \X + \lambda \unitm)^{-1} \X^T \y),
\]
 we will define $\HH_{\lambda} =  \X (\X^T \X + \lambda \unitm)^{-1} \X^T.$  Using the SVD of $\X,$ observe that the so-called Gram matrix $\X^T \X$ can be decomposed as follows
 \[
 	\X^T \X = \V \D \U^T \U \D \V^T = \V \D^2 \V^T,
 \]
 which is the eigendecomposition of $\X^T \X.$ Next, we can apply this to the hat matrix,

 \begin{eqnarray}
 \label{ridge_SVD1}
 	  \HH_\lambda & = & \U \D \V^T (\V \D^2 \V^T +\lambda \unitm)^{-1} \V \D \U^T \nonumber \\
					     & = & \U \D (\D^2 + \lambda \unitm)^{-1} \D \U^T,
\end{eqnarray}
since $\V \D^2 \V$ and $\lambda \unitm$ commute and are hence simultaneously diagonalizable. Therefore, it also follows
that $\HH_\lambda$ is diagonalizable, with respect to $\U$ and with eigenvalues given by $ \D (\D^2 + \lambda \unitm)^{-1} \D,$ which is a diagonal matrix of order $p \times p.$
Thus, since the trace of a matrix is equal to the sum of its eigenvalues, it readily
follows that
\[
 \tr (\HH_\lambda) = \sum_{j=1}^p \frac{d_j^2}{d_j^2 +\lambda} ,
\]
 where $d_j$ is the $j'$th diagonal entry of $\D$ and since $\lambda \geq 0,$ we have  $\frac{d_j^2}{d_j^2 +\lambda}  \leq 1.$ Like linear regression, ridge regression computes the coordinates of $\y$ with respect to the orthonormal basis $\U$ and then it shrinks these coordinates by the factor $\frac{d_j^2}{d_j^2 +\lambda}.$ This means that a greater amount of shrinkage is applied to the coordinates of basis vectors with smaller $d_j^2.$


\begin{exercise}
 What does a small value of $d_j^2$ mean? Discuss in terms of  Bayesian or frequentist approaches.
 What is the behavior of $ \tr (\HH_\lambda)$ if $\lambda \rightarrow 0$ and $\lambda \rightarrow \infty?$
 %The following figures could be used for motivation and as hints.
\end{exercise}




\subsubsection{Bayesian linear regression}
Although ridge regression is a useful way to compute a point estimate, sometimes we want to compute the full posterior over $\beta$ for a standard linear regression problem: $y_i = x_i^T \beta + \varepsilon_i$ and $\varepsilon_i \sim \mathcal N(0, \sigma^2),$ and $i=1,\ldots, N.$


This corresponds to the following likelihood function:
\[
	P(\y | \X, \beta, \sigma^2) \propto \exp (\frac{1}{2\sigma^2} (\y - \X \beta)^T (\y - \X \beta) ).
\]
For simplicity, we will assume that a prior distribution for $\beta$ is known, i.e., $\beta \sim \mathcal N (0, \frac{\sigma^2}{\lambda}) \propto \exp(-\frac{\lambda}{2 \sigma^2} \| \beta \|^2).$ Thus,  we focus on computing the posterior distribution $P(\beta | \X, \y, \sigma^2).$

Using Bayes' rule for Gaussians
\begin{eqnarray*}
	P(\beta | \X, \y, \sigma^2) &=& \frac{P(\y| \beta, \X, \sigma^2) P(\beta | \X,
	\sigma^2)}{P(\y | \X, \sigma^2)} \\
	      &=& \frac{\prod_{i=1}^N  P(y_i | \beta_i, x_i, \sigma^2) P(\beta | \X, \sigma^2)}{P(\y | \X, \sigma^2)},
\end{eqnarray*}
and substituting ${P(\y | \X, \sigma^2)}$ with some rescaling factor $Z'$ the posterior is given by
\begin{eqnarray*}
P(\beta | \X, \y, \sigma^2)  & = & \frac{1}{Z'} \prod_{i=1}^N      \exp (-\frac{1}{2\sigma^2} (y_i - x_i^T \beta)^2  \exp(-\frac{\lambda}{2 \sigma^2} \| \beta \|^2) \\
					& = &   \frac{1}{Z'}  \exp (-\frac{1}{2\sigma^2} \sum_{i=1}^{N} (y_i - x_i^T \beta)^2)  \exp(-\frac{\lambda}{2 \sigma^2} \| \beta \|^2) \\
					& = & \frac{1}{Z'}  \exp (-\frac{1}{2\sigma^2} [(\y - \X \beta)^T (\y - \X \beta) + \lambda \beta^T \beta]) \\
					& = &   \frac{1}{Z'}  \exp( -\frac12 [ \frac{1}{\sigma^2} \y^T \y + \frac{1}{\sigma^2} \beta^T (\X^T \X + \lambda \unitm)\beta - \frac{2}{\sigma^2} \beta^T \X^T \y]) \\
					& = & \mathcal N (\beta | \mu, \Sigma),
\end{eqnarray*}
where $\Sigma = \sigma^2 (X^T X +\lambda \unitm)^{-1}$ and $\mu = \frac{1}{\sigma^2} \Sigma X^T \y = (X^T X + \lambda I)^{-1} X^T \y.$

%In this case, the posterior mean $\mu$ is exactly the classical ridge estimate,
%More generally, the most likely parameter $\argmax_\beta P(\beta | X, \y, \sigma^2)$ is also the least-cost parameter $\argmin_\beta J(\beta).$ In the Gaussian case, mean and most-likely coincide.

The Bayesian inference approach not only gives a mean / optimal $\mu,$
 but also a variance $\Sigma$ of that estimate. This is a core benefit of the Bayesian view since it naturally provides a probability distribution over predictions, not only a single prediction.


%% LASSO REGRESSION
\subsection{The Lasso}

We have seen that ridge regression essentially re-scales the least squares estimates. The Lasso (Tibshirani 1996), which is also the shrinkage method, by contrast, tries to produce a sparse solution, in the sense that several of the slope parameters will be set to zero.

As for ridge regression, the Lasso can be expressed as a constrained minimization problem,
\begin{eqnarray}
	\label{Lasso}
	\hat \beta^{\mathtt{lasso}} &=& \argmin_\beta \sum_{i=1}^N \left(y_i - \beta_0 - \sum_{j=1}^p x_{ij} \beta_j \right)^2 \\
	&& \mbox{subject to}  \sum_{i=1}^p |\beta_j | \leq t.
\end{eqnarray}
In the signal processing literature, the lasso is also known as {\it basis pursuit}. We can also write the lasso problem in the equivalent {\it Lagrangian form}
\begin{equation}
	\label{LassoLagrange}
	\hat \beta^{\mathtt{lasso}} = \argmin_\beta \underbrace{\sum_{i=1}^N \left(y_i - \beta_0 - \sum_{j=1}^p x_{ij} \beta_j \right)^2 }_{J(\beta)}+\lambda \sum_{i=1}^p |\beta_j |.
\end{equation}

Notice the similarity to the ridge regression problem \eqref{FRR} and \eqref{CO_RR}: the $\ell_2$ ridge penalty $\sum_{i=1}^p \beta_j^2$ is replaced by the $\ell_1$ Lasso penalty $\sum_{i=1}^p |\beta_j|.$ This
latter constraint makes the solutions nonlinear in the $y_i.$ Thus, contrary to ridge regression, the Lasso does not admit a closed-form solution. The functional \eqref{LassoLagrange} is convex but not {\it strictly} convex, so that the solution is not unique.

The above constrained minimization is a quadratic programming problem, whose solution can be efficiently approximated (e.g., coordinate descent, interior point methods, $\ldots$).  The historical efficiency of algorithms to fit lasso models can
be summarized as follows:

\vspace{0.4cm}

\begin{center}
\begin{table}[h]
\begin{tabular}{ c | c | c | c }
  \hline
   \mbox{Year} & \mbox{Algortithm} & \mbox{Operations} & \mbox{Practical Limit}  \\
     \hline
  	1996 & \mbox{Quadratic programming} & $O(N2^p)$ & $\sim 100$  \\
	2003 & \mbox{LARS} & $O(Np^2)$ & $\sim 10,000$ \\
	2008 & \mbox{Coordinate descent} & $O(Np)$ & $\sim 1, 000,000$ \\
  \hline
\end{tabular}
  \caption{The historical efficiency of algorithms to fit lasso models.}
  \end{table}
\end{center}

%However,  efficient algorithms exist for computing the entire path of solutions as $\lambda$ varied, with the same computational cost as for ridge regression.
%Because of the nature of the constraint, making $t$ sufficiently small will
%cause some of the coefficients to be exactly zero.
\subsubsection{Lasso parameter estimation by coordinate descent}

Coordinate descent belongs the the first-order optimization methods. It optimizes each coordinate $\beta_j$ at a time iteratively using subgradient information $\frac{\partial f}{\partial \beta_j}(\beta )$(Algorithm \ref{algo:coordinat_descent}).
\begin{algorithm}[h!]
\label{algo:coordinat_descent}
\SetAlgoLined
  \textbf{Initialize} $\beta^{(0)}$\;
  \Repeat{\it converged}{
     \For{$j=1,...,P$}{
		$\beta^{(i+1)}_j = \beta^{(i)}_j - \alpha \frac{\partial f}{\partial \beta_j}(\beta )$\;
     }
    }
    \caption{Generalized coordinat descent algorithm in pseudo code with step size $\alpha$ and subgradient $\frac{\partial f}{\partial \beta_j}(\beta )$}
\end{algorithm}

If a closed form solution to the single variable optimization problem exists, then we can update the parameter directly at each step.
It turns out, Lasso possesses a closed form solution for the single variable optimization problem!

To derive at the closed form solution for an individual coordinate of the Lasso, we first will re-write the gradient of $J(\beta)$ with respect to $\beta_j$:
\begin{align}
\begin{split}
\frac{\partial}{\partial \beta_j} J(\beta ) &= -2\sum_{i=1}^N X_{ij} \left[y_i - \beta_0 - \sum_{j=1}^P\beta_j X_{ij}\right]\\
&= -2\sum_{i=1}^NX_{ij}\left[y_i -\beta_0 - \sum_{k \neq j}^P\beta_k X_{ik} - \beta_j X_{ij}\right]\\
&= -2\underbrace{\sum_{i=1}^NX_{ij}\left[y_i - \beta_0 - \sum_{k \neq j}^P\beta_k X_{ik} \right]}_{\hat{=}\rho} + 2\beta_j\underbrace{\sum_{i=1}^NX_{ij}^2}_{\hat{=}z_j}\\
&\hat{=} -2\rho +2\beta_j z_j
\end{split}
\end{align}

The problem arises when taking the derivative of the $\ell_1$ norm $\lambda\frac{\partial}{\partial \beta_j}|\beta_j|$ since the derivate is only defined for $\beta_j \neq 0$ as :
\begin{equation}
\frac{\partial}{\partial \beta_j} |\beta_j| = \begin{cases}
-1 \quad &\beta_j <  0\\
1 \quad &\beta_j > 0
\end{cases}
\end{equation}

To address the problem point $|\beta_j| = 0$ we use the notion of a {\it subgradient}, which is an extension of differential/gradients to non-differentiable points.
Remember that for convex and differentiable $f$:
$$ f(y) \geq f(x) + \partial f(x) (y-x) \qquad \mbox{for all } x,y,$$
i.e. the linear approximation always underestimates $f$.
More generally, the subgradient of a convex $f$ at $x$ is defined as the family of linear forms $V\in \partial f(x)$ that lower bound $f(x)$:
\begin{equation}
f(y) \geq f(x) + V(y-x) \qquad \mbox{for all } y
\end{equation}
Note that this always exists and equals the differential $\partial f(x)$ if $f$ is differentiable at $x$ (and this definition also works for non-convex $f$ but does not always exist then).
Example: The subgradient with respect to $|\beta_j|$ are all linear forms that have a slope within the interval $[-1, 1]$.

Putting everything together we arrive at the subgradient of $|\beta_j|$ as:
\begin{equation}
\partial_{\beta_j}|\beta_j| = \begin{cases}
-1 \quad& \beta_j < 0\\
[-1,1] \quad& \beta_j = 0\\
1 \quad& \beta_j > 0\\
\end{cases}
\end{equation}

and the subgradient of the Lasso as:
\begin{align}
\begin{split}
\partial_{\beta_j}J(\beta)^{\mathtt{lasso}} &= 2z_j\beta_j -2\rho +\begin{cases}
-\lambda \quad& \beta_j < 0\\
[-\lambda,\lambda] \quad& \beta_j = 0\\
\lambda \quad& \beta_j > 0\\
\end{cases}\\
&=\begin{cases}
2z_j\beta_j -2\rho-\lambda \quad& \beta_j < 0\\
[-2\rho-\lambda,-2\rho+\lambda] \quad& \beta_j = 0\\
2z_j\beta_j -2\rho+\lambda \quad& \beta_j > 0\\
\end{cases}\\
\end{split}
\end{align}

To derive the closed form solution for $\beta_j$ we set the subgradient to $0$ for the three individual cases independently.
\\
\\
\textbf{Case 1} ($\beta_j < 0$):
\begin{align}
\begin{split}
&2z_j\beta_j -2\rho-\lambda = 0\\
&\beta_j = \frac{\rho_j+\frac{\lambda}{2}}{z_j}.
\end{split}
\end{align}
For $\beta_j$ satisfying the inequality the following constraint $\rho_j< \frac{-\lambda}{2}$ must hold.
\\
\\
\textbf{Case 2} ($\beta_j = 0$):\\
For  $\beta_j = 0$ being an optimum the interval $[-2\rho-\lambda, -2\rho_j+\lambda]$ has to contain contain $0$. This only holds true if $-\frac{\lambda}{2}\leq\rho_j\leq\frac{\lambda}{2}$.
\\
\\
\textbf{Case 3} ($\beta_j > 0$):\\
\begin{align}
\begin{split}
&2z_j\beta_j -2\rho+\lambda = 0\\
&\beta_j = \frac{\rho_j-\frac{\lambda}{2}}{z_j}.
\end{split}
\end{align}
For $\beta_j$ satisfying the inequality the  following constraint $\rho_j> \frac{\lambda}{2}$ must hold.

Altogether this gives us the closed form solution of the Lasso with respect to $\beta_j$ as:
\begin{equation}
\hat{\beta}_j^{\mathtt{lasso}} = \begin{cases}
\frac{\rho_j+\frac{\lambda}{2}}{z_j} \quad& \rho_j < -\frac{\lambda}{2}\\
0 \quad& \rho \in  \left[-\frac{\lambda}{2},\frac{\lambda}{2}\right]\\
\frac{\rho_j-\frac{\lambda}{2}}{z_j} \quad& \rho_j  > \frac{\lambda}{2}\\
\end{cases}\\
\end{equation}

which is known as {\it soft-thresholding} update rule. The coordinate descent algorithm for the Lasso is then defined as:

\begin{algorithm}[h!]

\SetAlgoLined
  \KwIn{Regularization parameter $\lambda$}
  \KwResult{Optimal coeffecients $\hat{\beta }^{\text{LS}}$}
 \textbf{Initialize} $\beta = (X^TX+\lambda I)^{-1}X^Ty$\;
 \textbf{Precompute}  $z_j =\sum_{i=1}^N X_{ij}^2$\;
  \Repeat{\it converged}{
     \For{$j=1,...,P$}{
     \begin{align*}
        \rho_j &=\sum_{i=1}^N X_{ij}\left [ y_i -\sum_{k\neq j}^P\beta_k X_{ik}  \right ]\\
\hat{\beta}_j^{\mathtt{lasso}} &= \begin{cases}
\frac{\rho_j+\frac{\lambda}{2}}{z_j} \quad& \rho_j < -\frac{\lambda}{2}\\
0 \quad& \rho \in  \left[-\frac{\lambda}{2},\frac{\lambda}{2}\right]\\
\frac{\rho_j-\frac{\lambda}{2}}{z_j} \quad& \rho_j > \frac{\lambda}{2}\\
\end{cases}\\
     \end{align*}
     }
    }
\end{algorithm}

\subsubsection{Bayesian interpretation of Lasso}

The solution of \eqref{LassoLagrange} can be interpreted as posterior mode estimates when
the regression parameters $\beta_j'$s have independent and identical Laplace priors. More precisely, we consider a prior of the form
\[
	P(\beta | \tau) = \frac{1}{2 \tau}\prod_{j=1}^p \exp(- |\beta_j| / \tau),
\]
with $\tau = \frac{1}{\lambda}.$
We will use a uniform prior on the offset term, $P(\beta_0) \propto 1.$

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/Laplacian}
\caption{The Laplacian prior assigns more weight to regions near
zero than the normal prior.}
\label{fig:estlaplcian}
\end{figure}


Then we perform MAP estimation with this prior. The penalized negative log likelihood has the form
\[
	f(\beta) = -\log P(\y | \beta) - \log P(\beta | \tau) = -\log P(\y | \beta)  + \tau \|\beta\|_1,
\]
where $\| \beta \|_{1} = \sum_{j=1}^p |\beta_j |$ is the $\ell_1$ norm of $\beta.$

%{\bf NOTE:} Some info on the minimization algorithms!!!!

%% Generalisation REGRESSION
\subsubsection{More on Lasso}

Since 2005 with the introduction of the LARS algorithm, there have been a lot of activity and fruitful research in developing algorithms for fitting regularization paths for a variety of different problems, incl. high-dimensional problems.  In addition, $\ell_1$ regularization has taken on a life of its own, leading to the development of the field {\it compressed sensing} in the signal processing literature \cite{Donoho2006a}. We  touch briefly upon some of the Lasso modifications.

Moreover, we recall that the solution of $\ell_1$ regularization is not unique, i.e., the various solution have the same
prediction properties but different selection properties.

\subsubsection{The Dantzig Selector*}
In many real-life problems the number of variables / parameters $p$ is much larger than the number of observations $N.$ Then when fitting the model $\y = \X^T \beta$ we end up with infinitely many solutions, since our system is underdetermined.

Candes and Tao (2007) proposed a new estimator, the so-called {\it Dantzig selector}, to estimate $\beta$ as the solutions to the $\ell_1$ regularization problem.
\begin{eqnarray*}
	&& \min_{\beta} \|\beta\|_{1} \\
	 \mbox{ subject to } && \| \X^T (\y - \X \beta) \|_{\infty} \leq \sqrt{2 \log p} \sigma,
\end{eqnarray*}
where $\sigma$ is a noise level, and $\| \cdot \|_{\infty}$ denotes the $\ell_\infty$ norm, the maximum absolute value of the components of the vector. In this form it resembles the lasso, replacing
squared error loss by the maximum absolute value of its gradient.

It is known from the work on the Dantzig selector that one can reconstruct $\hat \beta^{\texttt{DS}}$ with the following accuracy
\begin{equation}
	\label{recon_rate}
	\| \hat \beta^{\texttt{DS}} - \beta \|_2^2 \leq C  s \sigma^2 \log p,
\end{equation}
where $s=\# \mathtt{supp}(\beta)$ and $s \leq N/2.$ The constant $C$ depends on the matrix $\X$ and is different from infinity, if the relevant variables are not correlated. However, error rate \eqref{recon_rate} is achieved under the assumption that someone has provided us with the exact location of the non-zero components of $\beta,$ i.e., the support of $\beta$ is known.

\subsubsection{The Grouped Lasso*}
In some problems, the predictors belong to pre-defined groups; for example
genes that belong to the same biological pathway. In
this situation it may be desirable to shrink and select the members of a
group together.
In this case, the classical lasso is going to select just one of them (this can be bad for interpretability but maybe good for compression), whereas the {\it grouped lasso} provide such an option. Suppose that the $p$ predictors are divided into $L$ groups, with $p_l$ the number in group $l$. We use the notation for a matrix $\X_l$ to represent the predictors corresponding to the $l$th group, with corresponding coefficient vector $\beta_l$. The grouped lasso minimizes the convex criterion
\begin{equation}
\label{GroupedLasso}
	\min_{\beta \in \R^p} \left(  \| \y - \beta_0 \mathbf 1 - \sum_{l=1}^L \X_l \beta_l \|_2^2 + \lambda  \sum_{l=1}^L \sqrt{p_l} \| \beta_l \|_2 \right),
\end{equation}
where the $\sqrt{p_l}$ terms accounts for the varying group sizes, and $\| \cdot \|_2$ is the Euclidean norm (not squared!).  Since the Euclidean norm of a vector $\beta_l$ is zero only if all of its components are zero, this procedure encourages sparsity at both the group and individual levels.   That is, for some values of $\lambda,$ an entire  group of predictors may drop out of the model.


%% Generalisation REGRESSION
\subsubsection{Discussion and Generalization}

We discuss and compare two approaches for restricting the linear regression model: the ridge regression and the Lasso. It is worthwhile mentioning that the described techniques could be applied to estimating any regression function $f(\X),$ not only of linear type.

In the case of an orthonormal input matrix $\X$ the two procedures have explicit solution. Each method applies a simple transformation to the least squares estimate $\hat \beta_j$ as shown in the table below.
\vspace{0.4cm}

\begin{center}
\begin{table}[h]
\begin{tabular}{ c | c  }
  \hline
   \mbox{Estimator} & \mbox{Formula} \\
     \hline
  	\mbox{Ridge regression} & $\hat \beta_j / (1+\lambda)$ \\
	\mbox{Lasso} & $\sgn (\hat \beta_j) (|\hat \beta_j| - \lambda)_+$ \\
  \hline
\end{tabular}
  \caption{Estimators of $\beta_j$ in the case of orthonormal columns of $\X,$ $\lambda$ is the regularization parameter chosen by the corresponding technique, $\sgn$ denotes the sign of its argument $(\pm 1),$ and $x_+$ denotes 'positive part' of $x$.}
  \end{table}
\end{center}

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/est1}
\caption{Ridge regression and lasso estimators (broken red line) for the case of orthonormal columns of $X.$ The $45^\circ$ line in gray shows the unrestricted estimate
for reference. From \cite{Hastie:2009wp}, Fig 3.4.}
\label{fig:estridge}
\end{figure}

Ridge regression does a proportional shrinkage. Lasso translates each
coefficient by a constant factor $\lambda,$ truncating at zero. This is called {\it soft-thresholding}.  For the nonorthogonal case, let us consider an illustrative example in 2D, i.e., $p=2$ and compare Lasso / ridge regression (see Figure \ref{fig:est2}). The residual sum of squares has elliptical
contours, centered at the full least squares estimate. The constraint region for ridge regression is the disk $\beta_1^2 + \beta_2^2 \leq t^2,$ while that for Lasso is the diamond $|\beta_1| + |\beta_2| \leq t.$ Both methods find the first point where the elliptical contours hit the constraint region. Unlike the disk, the diamond
has corners; if the solution occurs at a corner, then it has one parameter $\beta_j$ equal to zero. When $p>2,$ the diamond becomes a rhomboid, and has
many corners, flat edges and faces; there are many more opportunities for
the estimated parameters to be zero.


\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/est2}
\caption{Estimation picture for the Lasso (left) and ridge regression
(right). Shown are contours of the error and constraint functions. The solid blue
areas are the constraint regions $|\beta_1| + |\beta_2| \leq t$ and $\beta_1^2 + \beta_2^2 \leq t^2$ respectively, while the red ellipses are the contours of the least squares error function.
From \cite{Hastie:2009wp}, Fig 3.11.}
\label{fig:est2}
\end{figure}


We will now further generalize ridge regression and the lasso, and view them both in frequentist and bayesian formulations. Consider the criterion
\begin{equation}
\label{RegressionGeneral}
\bar \beta = \argmin_\beta \left[ \sum_{i=1}^N (y_i - \beta_0 - \sum_{j=1}^p x_{ij} \beta_j)^2 + \lambda \sum_{j=1}^p |\beta_j|^q    \right],
\end{equation}
for $q \geq 0.$
The contours of constant value of $\sum_j |\beta_j|^q$ are shown in Figure \ref{fig:contours}, for the case of two inputs.

\begin{figure}\centering
\includegraphics[width=.9\columnwidth]{images/RR_Shapes}
\caption{Contours of constant value of $\sum_j |\beta_j|^q$ for given values of $q$
From \cite{Hastie:2009wp}, Fig 3.12.}
\label{fig:contours}
\end{figure}

\subsubsection{Bayesian formulation}
Thinking of $|\beta_j|^q$ as the log-prior density for $\beta_j,$ these are also the equicontours
of the prior distribution of the parameters. The value $q=0$ corresponds
to variable subset selection, as the penalty simply counts the number
of nonzero parameters; $q=1$ corresponds to the Lasso; whereas $q =2$ to ridge regression. For $q \geq 1$ the prior is not uniform in direction, but concentrates more mass in the coordinate directions. The prior corresponding
to the $q=1$ case is an independent double exponential (or Laplace) distribution for each input.
For $q <1$ the constraint region is non-convex, which leads to non-convex optimization problem.

In this view, the Lasso and ridge regression are Bayes estimates with different priors. Note, however, that they are derived
as posterior modes, that is, maximizers of the posterior. It is more common
to use the mean of the posterior as the Bayes estimate. Ridge regression is
also the posterior mean, but Lasso  is not.
%% Generalisation REGRESSION
\subsection{Elastic Net}

Could  other values of $q$ rather than $q=1$ and $q =2$ bring better performance? Should one estimate $q$ from the data? The experience show that it is not worth the effort and actually  the values of $q \in (1,2)$ suggest the compromise between the Lasso and ridge regression (see Figure \ref{fig:contours}). Although this is the case, with $q>1,$ the penalty term $|\beta_j|^q$ is differentiable at $0$, and so does not share the ability of Lasso for setting coefficients exactly to $0.$ Partially for this reason as well as for computational tractability, \cite{ZouHastie05} introduced the {\it elastic-net} penalty
\begin{equation}
\label{ELNET}
	\lambda \sum_{j=1}^p (\alpha |\beta_j|^2 + (1-\alpha) |\beta_j|),
\end{equation}
which represents another compromise between ridge regression and the Lasso. In particular, the elastic-net selects variables like Lasso, and shrinks together the coefficients of correlated predictors like ridge function. It also has considerable computational advantages over the $\ell_q$ penalties.


\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/EN}
\caption{Contours of constant value of $\sum_j |\beta_j|^q$ for $q=1.2$ (left plot), and the elastic-net penalty \eqref{ELNET} with $\alpha =0.2$ (right plot). Although visually very similar, the elastic-net has sharp (non-differentiable) corners, while the $q=1.2$ penalty does not.
From \cite{Hastie:2009wp}, Fig 3.13.}
\label{fig:ELN}
\end{figure}

\subsection{Parameter and model selection}

The choice of the regularization parameter $\lambda$ in both ridge regression and in the lasso (as well as in the elastic net) is more of an art than a science. This parameter can be constructed as a {\it complexity parameter,} since as $\lambda$ increases, less and less effective parameters are likely to be included in both ridge and lasso regression. Therefore, one usually adopts a model selection perspective and compare different choices of $\lambda$ using, for instance, {\it cross-validation} or an {\it information criterion}. In the following we review briefly some of the existing approaches for the parameter selection problem as well as provide recipes for a few (but most well-known) of them.

 To motivate these approaches, we assume that we are only provided with a finite set of training examples, usually smaller than we wanted, and we need to adjust regularization parameter(s), so that our predictor has an ability to generalize, i.e., perform well on a new previously unseen data. A general idea is to split the given training data into disjoint subsets: {\it training set} and {\it testing set.}

In  the {\bf hold-out method} the given data set is split into two groups: {\it training set} (used to train the classifier) and {\it testing set} (used to estimate the error rate of the trained classifier).  It is usually applied to determine a stopping point for iterative methods. However, this method has two main drawbacks: (i) In problems where we have a sparse dataset we may not be able to afford the  'luxury' of setting aside a portion of the dataset for testing; (ii) Since it is a single train-and-test experiment, the hold-out estimate of error rate  will be misleading if we happen to get an 'unfortunate' split.

\begin{figure}\centering
\includegraphics[width=.6\columnwidth]{images/HOM}
\caption{Hold-out method}
\label{fig:HOM}
\end{figure}


The limitations of the hold-out method can be overcome with a family of resampling methods (such as cross validation, bootstrap) at the expense of more computations.

{\bf Random Sampling} performs $K$ data splits on the dataset.  Each split randomly selects a (fixed) number examples without replacement. Then for each data split we retrain the classifier from scratch with the training  examples and compute the error $RS_i$ with the test examples. The true error estimate is obtained as the average of the separate estimates $RS_i = \frac 1K \sum_{i=1}^K RS_i,$ which is significantly better than the hold-out estimate.

{\bf $K-$Fold cross validation} is a common technique for estimating the performance of a predictor. A general idea is  for each of $K$ experiments to use $K-1$ folds for training and the remaining one for testing. It is similar to Random Sampling. However, the advantage of $K-$fold cross validation is that all the examples in the dataset are eventually used for both training and testing.

\begin{figure}\centering
\includegraphics[width=.45\columnwidth]{images/RS}
\includegraphics[width=.45\columnwidth]{images/CV}
\caption{Random sampling (left panel) and 4-fold cross validation (right panel)}
\label{fig:RS}
\end{figure}


\begin{algorithm}
\SetAlgoLined
  \KwData{training data $D$, fold size $K$ (often 5 or 10)}
  \KwResult{regularization parameter $\lambda$}

Partition the training data $D$ into $K$ separate sets of equal size $D = (D_1, D_2, \ldots, D_K).$

\For{multiple $\lambda$}{
\For{$k=1,\ldots,K$}{
Fit the model with parameter $\lambda$ to the other $K-1$ parts (excluding the $k$th fold $D_k$). Denote the resulting estimates by $\hat \beta^\lambda_{-k}.$

Compute the prediction error in predicting $T_k$ part:
\[
	RSS_k (\lambda) = \sum_{(x_i, y_i) \in D_k} (y_i - x_i^T \hat  \beta^\lambda_{-k})^2.
\]
}

Compute the overall cross-validation error
\[
	CV(\lambda) = \frac1K \sum_{k=1}^K RSS_k (\lambda).
\]
}
Choose $\lambda^*$ that minimizes $CV(\lambda).$

\emph{Test Phase:}\\
Compute the estimates $\hat \beta^{\lambda^*}$ on the entire training set $D = (D_1, \ldots, D_K).$\\
Apply $\hat \beta^{\lambda^*}$ to the test set to assess the test error.
\caption{$K$-fold cross validation for choosing the regularization parameter $\lambda$}
\end{algorithm}


\begin{exercise}
	How many separate sets (folds) are needed? Is it better to have a small number or a large number of folds?
	What is leave-one-out cross validation?
\end{exercise}
%% REGULArization Parameter choice
%\subsection{Choice of regularisation parameter}

In the above approaches, although one model does not directly 'learn' from testing examples, it has 'seen' the testing examples and use their information to adjust its hyperparameters. Therefore, it is recommended to hold out a 'clean' testing set which the model does not 'touch' at all and make the final evaluation of model performance on the 'clean' testing set. So the entire dataset is split into three disjoint sessions: a training set used to train the model, a validation set used to provide an unbiased estimate of the model fit on the training set while tuning model hyperparameters, and a testing set used to evaluate the final model performance. This became a common practice in the neural network training.

\begin{figure}\centering
\includegraphics[width=.6\columnwidth]{images/TrainValTestSplit}
\caption{Data Split of training, validation and testing sets.}
\label{fig:TrainValTestSplit}
\end{figure}


\subsubsection{Information criteria based parameter choice}
Alternatively, one might use the Akaike information criterion (AIC) and the Bayes Information criterion (BIC) to select the regularization parameter(s). It is a computationally cheaper alternative to find the optimal value of $\lambda$ as the regularization path is computed only once instead of $K+1$ times when using $K$-fold cross validation. However, such criteria needs a proper estimation of the degrees of freedom of the solution, are derived for large samples (asymptotic results) and assume the model is correct, i.e. that the data are actually generated by this model. They also tend to break when the problem is badly conditioned (more features than samples).

{\bf Akaike information criterion (AIC)}  for choosing $\lambda$ is
\[
	AIC(\lambda) = \log \left[ \frac 1N \| \y - \X \hat \beta(\lambda) \|^2 \right] + 2 \frac{df (\lambda)}{N},
\]
where for the Lasso estimator, $df(\lambda) = \#$ of nonzero coefficients in the model fitted with $\lambda.$ In general case, $df(\lambda)$ is the {\it effective degrees of freedom} or {\it effective number of parameters.}


{\bf Bayesian information criterion} is more aggressive than AIC in seeking a sparse model. In particular,
\[
	BIC(\lambda) = \log \left[ \frac 1N \| \y - \X \hat \beta(\lambda) \|^2 \right] + 2 \log(N) \frac{df (\lambda)}{N}.
\]



%ARD? maybe later, next version

%todo: a bit more disc on graph lasso etc?
\newpage
\section{Classification}


A \defn{probabilistic classifier} predicts discrete class labels $y$ from observations $x$, i.e. is described by $P(Y=y|X=x)$ as discussed in the learning theory section.
To build a \defn{probabilistic classifier} $P(Y|X)$, one may
\begin{itemize}
\item create joint model $P(Y,X)$ and then condition on $X$
(\defn{generative approach})
\item directly estimate $P(Y|X)$ (\defn{discriminative approach})
\end{itemize}

\begin{exercise}
What are advantages and disadvantages of a generative model?
\end{exercise}
% just summarize some of discriminative vs egnerative models, nice disc in Murphy, chap 8.6 or so


\subsection{Linear regression for classification}\label{sec:linregclass}
\begin{exercise}
How can we use regression for classification? For $2$ classes, for $C>2$ classes?
\end{exercise}

Given $C$ classes of $Y$, replace the univariate $Y$ by $C$ indicators
$$Y_c := \left\{\begin{array}{ll}1 & Y=c\\0 & \mbox{otherwise}\end{array}
\right.$$
We take a discriminative approach by approximating the discrete
$P(Y_c=y|x,\theta)$ by the continuous model $\NNN(y|x^\top\beta_c, \sigma^2)$. The latter implies solving $C$ univariate regression problems.

In practice we can do this by considering the $N\times C$ indicator response matrix $\Y$ corresponding to concatenated samples of the $Y_c$'s, which yields the $(p+1)\times C$ coefficient matrix
$$\hat\beta:=(\hat\beta_c)_{c=1,\ldots,C}=(\X^\top\X)^{-1}\X^\top\Y.$$

The minimal expected prediction error (EPE) estimate then classifies new input $x$ (implicitly augmented by 1 to include intercept) as
\eq{\label{eqn:EPElinregclass}
\hat y(x) = \argmax_c E(Y_c|X=x,\hat\beta)= \argmax_c
x^\top\hat\beta_c
}

\begin{figure}
\includegraphics[width=\columnwidth]{images/class-REGvsLDA}\\
\includegraphics[width=.48\columnwidth]{images/class-REGvsLDA1d}
\caption{Linear regression used for classification of samples in a $C=3$ class problem in $\R^2$. Black lines indiciate class boundaries found by equation \eqnref{eqn:EPElinregclass}.
Linear regression (left) cannot separate the blue class, which is masked by the others, whereas LDA (right) can. Reason is that the regression for the blue class indicator variable cannot put the boundary to the left or right, as illustrated when projecting along the line joining the centers (bottom).
From \cite{Hastie:2009wp}, Fig 4.3 and 4.4}
\label{fig:class-REGvsLDA}
\end{figure}

This approach works well for $C=2$ classes but easily fails for $C>2$ due to the rigid nature of the linear model (see Figure \ref{fig:class-REGvsLDA}). This effect is reduced in higher dimensions, but we can take it into account with modified linear models as seen in the following.


%\subsection{Linear classification}
\subsection{Discriminant analysis}
Let us now take a generative approach. Given prior class probabilities $\pi_c:=P(Y=c)$ and using Bayes theorem, the probabilistic classifier $P(Y|X)$ can be constructed as
\eq{\label{eqn:Bayesclass}
P(Y=c|X=x,\theta)=\frac{P(X=x|Y=c,\theta)\pi_c}{\sum_{d=1}^C P(X=x|Y=d,\theta)\pi_d}.
}
Hence for classification it suffices to model the class conditional densities $P(X|Y)$.

In the funnily enough called \defn{(Gaussian) discriminant analysis}, we consider the generative Gaussian model
\eq{\label{eqn:GDA}
P(X=x|Y=c,\theta)=\NNN(x|\mu_c,\Sigma_c).
}

%Then with uniform prior on $Y$
So
$$
P(Y=c|x,\theta)\propto\NNN(x|\mu_c,\Sigma_c)\pi_c.
$$
%Assuming $\Sigma_c$ is diagonal, maximizing this yields the naive Bayes classifier \eqnref{eqn:BayesClassifier} from the introduction i.e.
%$$
%\hat y(x)= \argmax_c \left(\log P(x|\theta_c) + \log\pi_c\right).
%$$
For example, choosing a uniform prior gives
$$
\hat y(x)= \argmin_c (x-\mu_c)^\top\Sigma_c^{-1}(x-\mu_c).
$$

\subsubsection{Quadratic discriminant analysis}
With the Gaussian model of the class densities from above, we can determine decision boundaries between classes $c$ and $c'$ as
$$\{ x\in\R^p | P(Y=c|x,\theta)=P(Y=c'|x,\theta)\}$$
which is equivalent to requiring the joint likelihoods $P(x,Y=c|\theta)$ to be the same, or
$$P(x|Y=c,\theta)\pi_c=P(x|Y=c',\theta)\pi_{c'}.$$
Plugging in the model yields decision boundaries of the form
$\delta_c(x)=\delta_{c'}(x)$ with
\eq{\label{eq:qda}
\delta_c(x):=%&\log P(x,Y=c|\theta)\\
%\nonumber =&
-\frac{1}{2} \log \det \Sigma_c - \frac{1}{2} (x-\mu_c)^\top\Sigma_c^{-1}(x-\mu_c) + \log \pi_c,
}
i.e. quadratic functions. Hence this model is also known as \defn{quadratic discriminant analysis (QDA)}. See figure \ref{fig:class-QDA} for two examples.

\begin{figure}\centering
\includegraphics[width=\columnwidth]{images/class-QDA}\\
\caption{Quadratic discrimant analysis for $C=2$ (left) and $C=3$ (right).
From \cite{Murphy:2012uq}, Fig 4.3.}
\label{fig:class-QDA}
\end{figure}

\begin{exercise}
How to estimate the $\mu_c$ and $\Sigma_c$?
\end{exercise}

\subsubsection{Parameter estimation for discrimant analysis}
To determine the parameters, the simplest approach is to do MLE for the training data i.e. to maximize
$$P(D|\theta)=\prod_{i=1}^N P(x_i,y_i|\theta)=
\prod_{i=1}^N P(x_i|y_i,\theta)P(y_i|\theta)$$
with corresponding negative log likelihood
\begin{eqnarray*}
l(\theta)&=&
\sum_{i=1}^N\sum_{c=1}^C 1_{y_i=c} (\log \NNN(x_i|\mu_c,\Sigma_c) + \log\pi_c)\\
&=&\sum_{c=1}^C \left(\sum_{i | y_i=c} \log \NNN(x_i|\mu_c,\Sigma_c)\right) + \sum_{c=1}^C N_c\log\pi_{c}
%+ \sum_{i=1}^N\log\pi_{y_i}
\end{eqnarray*}
The minimization separates into optimization of the $\pi_c$ (with constraint $\sum_c \pi_c=1$ as in naive Bayes) as well as $C$ separate Gaussian MLE problems, which with $N_c:=|\{i|y_i=c\}|$ altogether yield
\begin{eqnarray*}
\hat\pi_c&=&\frac{N_c}{N}\\
\hat\mu_c&=&\frac{1}{N_c} \sum_{i | y_i=c}x_i\\
\hat\Sigma_c&=&\frac{1}{N_c} \sum_{i | y_i=c}(x_i-\hat\mu_c)(x_i-\hat\mu_c)^\top\\
\end{eqnarray*}
or reduced denominator for unbiased estimates of covariance.

MLE is prone to overfitting in particular in large dimensions, in particular estimates of $\hat\Sigma_c$ are singular for $N_c<p$ or at least ill-conditioned. As in regression, there are various ways of regularization, as discussed in the following.


\subsubsection{Linear discriminant analysis}
In \defn{linear discriminant analysis (LDA)} we assume that each class shares the same covariance structure $\Sigma=\Sigma_c$, thereby strongly reducing dimensionality of the parameter space.

From QDA, we know that the class boundaries are determined by $\delta_c(x)=\delta_{c'}(x)$. The $\delta_c(x)$ from equation  \ref{eq:qda} can now be simplified as follow
\eq{
\delta_c(x):=-\frac{1}{2} \log \det \Sigma - \frac{1}{2} x^\top\Sigma^{-1}x + x^\top\Sigma^{-1}\mu_c
- \frac{1}{2} \mu_c^\top\Sigma^{-1}\mu_c
+ \log \pi_c
}
The first two terms are independent of $c$, hence the class boundaries are linear and determined by $\tilde\delta_c(x)=\tilde\delta_{c'}(x)$ with \defn{linear discriminant functions}
$$\tilde\delta_c(x)= x^\top\Sigma^{-1}\mu_c - \frac{1}{2} \mu_c^\top\Sigma^{-1}\mu_c + \log \pi_c.
$$

Estimation is again be via MLE, which yields
$$\hat\Sigma=\frac{1}{N} \sum_{c=1}^C \sum_{i | y_i=c}(x_i-\hat\mu_c)(x_i-\hat\mu_c)^\top$$
or normalized by $N-K$ for an unbiased covariance estimate.

Alternatively, estimation can also be done by data sphering (or whitening, i.e. transforming $X$ to have $\Sigma=\unitm$) followed by classification to closest class mean in the transformed space (modulo effect of class priors).

See figure \ref{fig:class-REGvsLDA} for how LDA solves the 3 class problem, where linear regression struggled with masking. However for two classes, the LDA decision boundary is very similar to classification by linear regression (section \ref{sec:linregclass}), equal if the classes are of same size.

Note that none of the above methods can solve linearly not separable problems such as data coming from a fuzzy XOR function.

%TODO: Alex, exercise on this? exercise 4.2 from Hastie, and integration with log reg, see Murphy 4.2.3

\begin{exercise}
How to approximate QDA by LDA?
\end{exercise} % by
\vspace{-1em}\noindent$\to$ A common way for such approximation is to extended the covariate space with $X_iX_j$ as new covariates ($i\leq j$).


\subsubsection{Regularizations}

Selecting equal covariances in QDA is a strong model simplification to reduce overfitting. Other ideas have been proposed as well:
\begin{itemize}
\item diagonal covariances, which leads to the naive Bayes classifier
\item diagonal covariances in LDA, which leads to \defn{diagonal LDA}:
then according to \eqnref{eq:qda},
$$\bar\delta_c(x)=
-\sum_{i=1}^p\frac{(x_i-\mu_{ci})^2}{2\sigma_i^2} + \log \pi_c,
$$
where the determinant was left out since it is independent of $c$, and the $\sigma_i$'s are the diagonal entries of $\Sigma$. In high dimensions this simple model may outperform QDA and LDA considerably.
\item priors on the QDA model, followed by MAP estimation, which is known as \defn{regularized LDA}: as trade-off between QDA and LDA, consider
$$\hat\Sigma_c(\lambda):=\lambda \hat\Sigma_c + (1-\lambda)\hat\Sigma,\qquad \lambda\in[0,1],$$
with class-specific and pooled covariance estimates $\hat\Sigma_c$ and $\hat\Sigma$ from above. In Bayesian terms this can be formulated as an inverse Wishart prior. As in Lasso, $\lambda$ can be chosen by design or by CV.

Similarly starting from LDA, one may shrink its covariance to diagonal LDA by
$$\hat\Sigma(\lambda):=\lambda \hat\Sigma + (1-\lambda)\sigma^2\unitm,\qquad \lambda\in[0,1].$$
After sphering, these shrinkage approaches can be implemented efficiently, similar to ridge regression.

\item low-dimensional projection followed by QDA and LDA:
QDA and LDA implicitly reduce dimension, since distances only in the $(C-1)$-dimensional subspace spanned by the $\mu_c$'s matter. So one may simply project onto that, interesting for visualization for $K=3$ for example. For larger $K$ dimension reduction of the centroid space by e.g. PCA is common and leads to reduced-rank LDA/QDA.
\end{itemize}

\begin{exercise}
Which model is the best?
\end{exercise}
%again mention bias variance trade off -> no free lunch = no unique answer, choose by CV?

\subsubsection{Nonlinear extensions}
As seen with LDA one cannot learn nonlinearly separable problem such as an XOR function. Instead more complex model assumptions are necessary:

\begin{itemize}
\item One popular extension of QDA is to replace the Gaussian distribution for $P(x|Y=c,\theta)$ by a \defn{Gaussian mixture model}:
$$P(x|Y=c,\theta)=\sum_{k=1}^{K_c} w_{c,k} \NNN(x|\mu_{c,k}, \Sigma_{c,k})$$
with $w_{c,k}>0, \sum_{k=1}^{K_c} w_{c,k}=1$. The added flexibility of course results in higher model complexity, and the inference of the class density parameters cannot be done in closed form anymore, indeed is not convex. We will learn more about mixture models later, and for more details now see \cite{Hastie:2009wp}, chapter 6.8.
\item Instead of a non-normal class density model, we may also choose a fully non parametric density estimate in each class, e.g. by kernel density estimation, which implies
$$P(x|Y=c,\theta)=\frac{1}{N_c}\sum_{i=1}^{N} 1_{y_i=c}\NNN(x|x_i,, \sigma^2 \unitm)$$
or more generally
$$P(x|Y=c,\theta)=\frac{1}{N_c}\sum_{i=1}^{N} 1_{y_i=c}\kappa_h(|x-x_i|)$$
with some kernel $\kappa_h$ parametrized by kernel width $h$. Much can be said about approximation by kernels, but this is out of scope of this lecture.
\end{itemize}
% \begin{comment}

\subsection{Logistic regression}\label{sec:logreg}

Now take a discriminative approach, i.e.~we model the posterior $P(Y|X)$. Let us consider the two-class case $c\in\{0,1\}$ first.

In logistic regression, we use a linear model fed into a sigmoidal to ensure mean in $[0,1]$ of a Bernoulli distribution: The logistic model reads
\eq{\label{eqn:logregfull}
P(Y=c|x,\theta)=\Ber(c|\sigm(x^\top \beta))
}
with sigmoidal or logistic function $\sigm(x)=(1+\exp(-x))^{-1}$.
Thresholding the probabilities at $0.5$ gives a linear prediction boundary, with a normal vector specified by $\beta$.

The inverse of the logistic function
$$L(p):=\log\left(\frac{p}{1-p}\right)$$ % because 1/(1+e^-x)(1+e^-x-1)/(1+e^-x)=e^x
is called the \defn{logit function}.
%The linear prediction boundary is then given by $1/2=\sigm(x^\top \beta)(1-\sigm(x^\top\beta)$.
%s(x)(1-s(x))=s(x)(1+e^-x-1)/(1+e^-x)=...
With mean $\mu=P(Y=1|x,\beta)$ of the Bernoulli distribution, the model can be written as
$$E(Y|X,\beta) = \sigm(x^\top \beta)$$
or equivalently
$$L(E(Y|X,\beta))=L(\mu)=x^\top \beta.$$
The term
$$L(\mu)=\log\left(\frac{P(Y=1|x,\beta)}{P(Y=0|x,\beta)}\right)$$
is also called the \defn{log-odds}, which in logistic regression, we model as linear function of the input.

\subsubsection{Likelihood and estimation}
With $\mu_i=\sigm(x_i^\top \beta)$,
the negative log likelihood of \eqnref{eqn:logregfull} is given by
\begin{eqnarray*}
l(\theta) &=& - \sum_{i=1}^N \log\left(\mu_i^{1_{y_i=1}}\left(1-\mu_i\right)^{1_{y_i=0}}\right)\\
&=& - \sum_{i=1}^N \left(y_i\log\mu_i+(1-y_i)\log(1-\mu_i)\right).
\end{eqnarray*}
The latter is also known as cross-entropy.

\begin{exercise}
How can we estimate parameters? MLE.\\
How can we minimize $l(\theta)$?
\end{exercise}

We infer the parameters as before by maximum likelihood. But in contrast to all previous models, we cannot solve this in closed form but instead have to numerically minimize it.
As we will see in the exercises, the derivatives of $l(\theta)$ are given by
\begin{eqnarray}\label{eqn:logreggrad}
\nabla_\theta l(\theta) &=& \left(d/d\theta l(\theta)\right)^\top = \sum_{i=1}^N x_i(\mu_i-y_i) = \X^\top(\mu-\y)\\
\nabla^2_\theta l(\theta) &=& H(\theta) = \sum_{i=1}^N \mu_i(1-\mu_i)x_ix_i^\top = \X^\top\S\X
\end{eqnarray}
with gradient $\nabla_\theta l(\theta)$ and Hessian $H(\theta) = \nabla^2_\theta l(\theta)$ of $l$ and diagonal matrix $\S:=\diag(\mu_i(1-\mu_i))$. Since $\mu_i\in(0,1)$, the latter is positive definite (see exercises) and then so is $H$. Hence the problem is strictly convex and has a unique solution.

Stay tuned on how to find it numerically.


\subsubsection{Unconstrained local optimisation in a nutshell}
In unconstrained optimization, we want to solve
\begin{equation*}
\min_\theta  l(\theta)
\end{equation*}
The basic idea to solve this is to construct a sequence or optimizer path $\theta^{k}, k=0,1,2,\ldots$
%$\theta^{0}, \theta^{1}, \ldots, \theta^{k}, \theta^{k+1}, \ldots$
along which $l(\theta)$ decreases.
%$l(\theta^{0}) > l(\theta^{1}) > \ldots > l(\theta^{k}) > l(\theta^{k+1}) > \ldots, $
%using local properties of the objective function.
%
%\begin{exercise}
%Why do we want to use only local properties?
%\end{exercise} % not true for line search
%
The optimizer path is commonly constructed as
\begin{equation*}
\theta^{k+1} = \theta^{k} + t^k \Delta_\theta^k
\end{equation*}
in which $\Delta_\theta^k \in \mathbb{R}^{\dims{t}}$ is search direction and $t^k \in \mathbb{R}_+$ the step length.

\begin{exercise}
How to determine $\Delta_\theta^k$?
\end{exercise}

We select the search direction as steepest descent direction for some vector norm $||.||$:
$$\Delta_\theta^k = \arg \min_{v} \frac{\nabla_\theta l(\theta^k)^\top v}{||v||^2}$$
%$\Delta_\theta^k$ in this norm with strongest decrease.
For example:
\begin{itemize}
\item Euclidean norm $||v|| := v^T v$: $\Delta_\theta^k = - \nabla_\theta l(\theta^k)$ $\rightarrow$ gradient descent
\item General quadratic norm $||v|| := v^T P v$ with positive definite $P$:  $\Delta_\theta^k = - P^{-1} \nabla_\theta l(\theta^k)$
\end{itemize}
Note that local norms are possible, $P = P(\theta^{k})$.

\begin{exercise}
Draw example: quadratic objective + direction using euclidian norm and best quadratic norm
\end{exercise}

\begin{exercise}
How to determine $t^k$?
\end{exercise}

We may select step length either by
exact line search $$t^k:=\argmin_{t>0} l(\theta^k+t \Delta_\theta^k)$$ or by approximations such as heuristic choice of learning rate $t^k:=\eta(k)$ e.g. $\eta(k):=\eta/k$, or
backtracking line search, i.e. successive geometric reduction of $t^k$:

\begin{algorithm}
\KwIn{initial guess for $t^k$ (e.g., $t = 1$), parameters $\alpha\in(0,1/2),\beta\in(0,1)$}
\Repeat{$l(\theta^k+t^k \Delta_\theta^k)<l(\theta^k)+\alpha t^k \nabla l(\theta^k)^T \Delta_\theta^k$}{
$t^k\leftarrow \beta t^k$
}
\end{algorithm}

\includegraphics[width=.8\columnwidth]{images/backtrackingLineSearch}

A local optimization algorithm can therefore be summarized as shown above.

\begin{algorithm}
%\SetLine
\KwIn{starting point $\theta^0$, $k = 0$}
\Repeat{stopping criterion holds}{
determine steepest descent direction $ \Delta_\theta^k$\\
\emph{line search}: choose step size $t^k$\\
\emph{update}: $\theta^{k+1} \leftarrow \theta^{k} + t^k \Delta_\theta^k$, $k \leftarrow k + 1$
}
\end{algorithm}

Common stopping criteria are
\begin{itemize}
\item small change of objective function: $|l(\theta^{k+1}) - l(\theta^{k})|<\epsilon_l$
\item small change of parameters: $||\theta^{k} - \theta^{k-1}||_2<\epsilon_\theta$
\item number of function evaluations: $k > N_k$
\end{itemize}


\subsubsection{Newton's method}
An alternative to gradients and line search is to minimize the second-order Taylor approximation
$$l(\theta^k + v) \approx l(\theta^k) + \nabla_\theta l(\theta^k)^\top v + \frac{1}{2} v^\top H(\theta^k) v.$$
Minimizing the right-hand side i.e. the quadratic form approximation yields
$$0=\nabla_\theta l(\theta^k) +H(\theta^k)v$$
and therefore the combined step-size and update
$$t \Delta_\theta^k = - H^{-1}(\theta^k) \nabla_\theta l(\theta^k),$$
%$$t \Delta_\theta^k = \arg \min_{v}l(\theta^k) + \nabla_\theta l(\theta^k)^\top v + \frac{1}{2} v^\top H(\theta^k) v$$

Newton's method has local quadratic convergence (instead of linear of steepest descent) but fails if $H(\theta^k)$ is not positive definite. Robust extensions are trust region methods, where $||v||$ is constrained or regularizations of $H(\theta^k)$.

\subsubsection{Application to logistic regression}
Applying Newton's method to logistic regressions with gradient and Hessian from \eqnref{eqn:logreggrad}, we get
\begin{eqnarray*}
\theta^{k+1}&=&\theta^k - H^{-1}(\theta^k) \nabla_\theta l(\theta^k)\\
&=& \theta^k +  (\X^\top\S_k\X)^{-1} \X^\top(\y-\mu)\\
&=& (\X^\top\S_k\X)^{-1} \X^\top\S_k (\X\theta^k +\S_k^{-1}(\y-\mu))\\
&=& (\X^\top\S_k\X)^{-1} \X^\top\S_k \z_k
\end{eqnarray*}
with working response $\z_k:=(\X\theta^k +\S_k^{-1}(\y-\mu))$.
This rewritten Newton method thus solves the weighted least-squares problem
$$\argmin_\beta (\z_k-\X\beta)^\top\S_k(\z_k-\X\beta).$$
Since $\S_k=\S_k(\beta)$ changes in each update, the algorithm is also called \defn{iteratively reweighted least-squares}.

\subsubsection{Regularizations}
Regularizations (or Bayesian priors) in form of L1-norm and L2-norm are straight forward and implemented as in the linear regression case. They are common and very useful for large dimensions or correlated covariates.

But also for the large sample setting regularization is important ---  in a linearly separable setting, MLE is obtained if $\|\beta\|\to\infty$ (why?). This corresponds to an infinitely steep sigmoidal at the sample boundary
$$\{x| \sigm \beta^\top x=1/2\}= \{x| \beta^\top x=0\}.$$
The scale of $\beta$ does not change the decision so in the optimization it is tuned to assign maximal probability mass to the training data thus not taking into account distance to this boundary. This may not generalize well. Here L2-regularization (ridge) can help a lot. For implementation the logistic regression negative log likelihood $l(\theta)$ is extended by
$$l_\lambda(\theta) = l(\theta) + \lambda \theta^\top\theta,$$
or more generally
$$l_\lambda(\theta) = l(\theta) + \lambda \|\theta\|_p^p,$$
which is easily propagated into gradient, Hessian and the IRLS learning rule. We proceed similarly for L1-, elastic net or other regularizations, if we e.g. want to do feature selection.

\subsubsection{More than two classes}
\begin{exercise}
How to setup a model for $C>2$ classes?
\end{exercise}

The two-class logistic regression model is defined as
$$P(Y=1|x,\beta)= \frac{\exp(\beta^\top\x)}{1+\exp(\beta^\top\x)}.$$
With $\beta_0:=0$ and $\beta_1:=\beta$, we can equivalently write
$$P(Y=c|x,\beta)= \frac{\exp(\beta_c^\top\x)}{\exp(\beta_0^\top\x)+\exp(\beta_1^\top\x)}.$$

The \defn{multinomial logistic regression} model generalizes this to $C>2$ by defining
$$P(Y=c|x,\beta)= \frac{\exp(\beta_c^\top\x)}{\sum_{d=1}^C\exp(\beta_d^\top\x)}$$
with $\beta_C:=0$ to ensure identifiability. This is equivalent to modeling the log-odds versus a reference category, i.e. the $C-1$ linear models
$$\log\left(\frac{P(Y=c|x,\beta)}{P(Y=C|x,\beta)}\right)=\beta_c^\top\x
$$
for $c=1,\ldots,C-1$ --- the equation for $c=C$ actually implies $\beta_C=0$ if it is to hold for all $\x$.

Likelihood and inference are done similar as to the 2-class case, although due to higher dimensionality quasi-Newton optimization may need to be used.

Note that in contrast to linear regression for classification, masking does not occur due to the added flexibility of the log odds comparisons, also see exercises.





%\subsection{Extensions of logistic regression}
%maybe nonparam log regr. (or quantile regression etc? hm, no)
%bayesian log reg?
\subsection{Classification and regression trees (CART)}
%murphy 16.2, essentially merges hastie anyway
%also nice: James chap 8

We will introduce tree-based methods for regression and classification. The are based on the idea of segmenting the predictor space into a number of simple regions, in which a constant function is assumed. The segmentation can be summarized by a tree, hence the name. See figure \ref{fig:CARTexample} for an example.

\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/CART-example}
\caption{Example of input space partitions (top, right one from recursive splitting) and corresponding tree and function defined on that. From \cite{Hastie:2009wp}, Fig 9.2.}
\label{fig:CARTexample}
\end{figure}

\subsubsection{Regression trees}
Suppose first that a partition of the $p$-dimensional input space into $M$ distinct regions $R_1,\ldots, R_M\subset\R^p$ is given. The \defn{regression tree} predictor is given by
$$\hat f(x)=\sum_{m=1}^Mc_m 1_{R_m}(x)$$
with indicator function $1$ as used in nearest neighbor or prototype methods. Clearly this is also an adaptive basis-function model (ABM).

\begin{exercise}
How to estimate $c_m$?
\end{exercise}

If as before we estimate $f$ by minimizing the (squared) expected prediction error under least-squares criterion (corresponding to MLE for a normal error model), it is easy to see that for fixed $\{R_m\}$
$$\hat c_m=\ave \{y_i | x_i\in R_m\}$$
with average
$$\ave S = \frac{1}{|S|}\sum_{s\in S} s$$
for some set $S\subset\R^p$.

\subsubsection{Classification trees}

\begin{exercise}
How to model a classification tree?
\end{exercise}

In the classification setting with classes $c=1,\ldots, C$, let
$$\hat p_{mc}:=\ave \{1_{y_i=c} | x_i\in R_m\}
= \frac{1}{|\{x | x\in R_m\}|} \sum_{x_i\in R_m}1_{y_i=c}
$$
denote the percentage of samples of class $c$ in region $R_m$.
The \defn{classification tree} or \defn{decision tree} predictor is given by
$$\hat f(x)=\sum_{m=1}^Mc_m 1_{R_m}(x)$$
with
$$c_m:=\argmax_c p_{mc},$$
which means simply selecting the most probable class in the region. In practice this means storing the class frequencies per region (in each tree leaf). See example in figure \ref{fig:CARTexampleClass}.


\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/class-picExampleMurphy}\\
\vspace{1em}
\includegraphics[width=.35\columnwidth]{images/class-picExampleMurphy-CART}
\caption{Classification trees. Top: Simple data example from \cite{Murphy:2012uq}, Fig 1.1. Bottom: Classification/decision tree for the data, with leafs indicating class membership (class yes, class no). The tree contains mostly pure class leafs.}
\label{fig:CARTexampleClass}
\end{figure}


\subsubsection{Growing the tree}
\begin{exercise}
How to choose best partition?
\end{exercise}

Finding the best partition in terms of least squares is in general NP-hard, hence often greedy algorithms are chosen (in methods then called CART, C4.5 or ID3).
% \cite{Hyafil:1976cx} was cited here, that does not exist. Maybe
% @article{hyafil1976bounds,
% title={Bounds for selection},
%  author={Hyafil, Laurent},
%  journal={SIAM Journal on Computing},
%  volume={5},
%  number={1},
%  pages={109--114},
%  year={1976},
%  publisher={SIAM}
% }

Our goal is to iteratively split the input data set along a splitting variable $j$ and a splitting value $s$. Define the half planes
$$R_-(j,s):=\{x|x_j\leq s\} \qquad R_+(j,s):=\{x|x_j> s\}.$$
Least-squares implies that the optimal (greedy) split is given by
$$(j^\ast,s^\ast)=\argmin_{j,s} \min_{c_-} \sum_{x_i\in R_-(j,s)} (y_i-c_-)^2 + \min_{c_+} \sum_{x_i\in R_+(j,s)} (y_i-c_+)^2.$$
Clearly, the inner minimization is solved by
$$\hat c_-=\ave \{y_i | x_i\in R_-(j,s)\}
\qquad
\hat c_+=\ave \{y_i | x_i\in R_+(j,s)\}$$
This procedure is repeated recursively in each region (i.e. for samples only from that region) to build the tree, until some depth has been reached.

If we want to fit a classification/decision tree instead, we need to replace the least squares cost by:
%As in logistic regression, we first fit a multinoulli model to the data in each leaf for the test $X_j<s$ by estimating class-conditional probabilities
\begin{itemize}
\item misclassification rate (or 0-1-loss):
$$
\frac{1}{|\{x_i\in R_m\}|}\sum_{x_i\in R_m}1_{y_i\neq c_m}
$$
with $c_m$ the majority class as defined above. This equals $1-p_{mc_m}$.
\item entropy:
$$
-\sum_{c=1}^C p_{mc}\log p_{mc}
$$
which avoids homogeneous sample distributions and is minimal if one $p_{mc}=1$. %This is equivalent to maxinimizing the information gain
\item Gini index:
$$\sum_{c=1}^C p_{mc}(1-p_{mc})=1-\sum_{c=1}^C p_{mc}^2$$
which is the expected error rate.
\end{itemize}

\begin{exercise}
How large should a tree be? When to stop growing?
\end{exercise}

If the tree is too small, we will not get good classification rates, if it is too large, we get overfitting. There are various ways to approach this bias variance tradeoff for trees. A common one is to grow a full tree and then prune it back.

\begin{exercise}
Why not do a simple stopping criterion in terms of cost function improve? (hint: XOR example)
\end{exercise}

Pruning the tree is commonly done by removing branches that give the minimal increase in classification error.  The optimal pruning is then selected by cross validation, see \cite{Hastie:2009wp}, section 9.2.2 and references.


\subsection{Bagging}
Tree-based methods are simple and useful for interpretation, e.g. they follow decision rules common for disease classification by physicians.
However, they simply do not classify as accurately as some of the other methods --- in parts because of the greedy estimation algorithm, and in part because of the rigid model. Also they have stability issues with small data changes possibly resulting in quite different trees.

To reduce the variance of the estimator, one can average over many estimates. In \defn{bagging} (= bootstrap aggregation), $M$ so called weak learners $f_m$ trained on boostrap samples (i.e. randomly chosen samples with replacement) are aggregated as
$$
f(x)=\ave \{f_m(x)\}.
$$
or by majority voting in the case of classification.

Combining a large number of trees can often result in strong improvements in prediction accuracy --- at the expense of loss in interpretability.

\subsection{Random forests}
However rerunning the same method on different data subsets results in highly correlated predictors, which limits the variance reduction

To improve this, random forests are ensemble of decision/regression trees learned from randomly chosen subsets of input variables as well as on randomly chosen data subsets.

To be more precise: As in bagging, a random forest is an ensemble of decision trees learned from bootstrapped training samples. But every time a split is considered, a random subsample of $m$ (often $m\approx\sqrt{p}$) predictors is chosen as split candidates from the full set of $p$ predictors. A new sample of $m$ predictors is taken at each split.

\begin{exercise}
Why does this improve the classification and not make it worse?
\end{exercise}

This works well, because if for example a strong predictor is present, it is chosen on the top level for all trees in bagging, which results in highly correlated predictions in bagging. In RFs, other predictors get a chance too, which can be seen as `decorrelating´ the trees.

Another reason for robustness of RFs is the use of \defn{ouf of the bag (OOB) samples}:
For each observation $(x_i,y_i)$, its random forest  predictor is constructed by averaging only those trees corresponding to bootstrap samples, in which $(x_i,y_i)$ did not appear. Hence the OOB error estimate is very similar to cross validation error estimates.

Random forests are easy to implement and inherently robust, hence wildly popular.

\subsection{Boosting}
\subsubsection{Ensemble learning}
Ensemble learning refers to learning a weighted combination of base models of the form
$$f(y|x,\pi)=\sum_{m\in\mathcal{M}}\omega_mf_m(y|x)$$
where the $f_m$ is a base model and ensemble learning can be interpreted as a \textbf{committee decision} where each base model $f_m$ gets a weighted "vote" (weighted by $\omega_m$). As shown in the schematic plot of Fig. \ref{fig:ensemblelearning}, \textbf{bagging} and \textbf{boosting} are two most popular ensemble methods.
\begin{itemize}
\item \textbf{Bagging}: Each base model $f_m$ is trained by a random subset of the data in a parallel fashion. One representative machine learning method that uses bagging is random forest, in which individual base models are deep decision trees which usually have low bias and high variance. Yet the committee, i.e. random forest, reduces model variance with majority voting of predictions made by individual decision trees.
\item \textbf{Boosting}: Base models are trained in a sequential fashion. Each base model learns from mistakes made by the previous base model. An representative machine learning method is Adaboost, which consists of shallow trees (also known as weak learners) with high bias and low variance.  Adaboost reduces bias by combining a sequence of weak learners into a strong classifier.
\end{itemize}

\begin{figure}\centering
\includegraphics[width=.99\columnwidth]{images/ensemblemethods.png}
\caption{Two popular ensemble methods: bagging vs boosting, taken from \cite{towardsdatascience}}
\label{fig:ensemblelearning}
\end{figure}

\begin{exercise}
Random forest and Adaboost both belong to the ensemble family. Please list three key differences between these two methods.
\begin{itemize}
\item bagging (trained in a parallel fashion) vs boosting (trained in a sequential fashion);
\item base model: deep decision trees (low bias \& high variance) vs shallow trees (decision stumps, high bias \& low variance);
\item stochastic vs deterministic.
\end{itemize}
\end{exercise}


\subsubsection{AdaBoost}
AdaBoost, short for Adaptive Boosting, is proposed by Yoav Freund and Robert Schapire \cite{Freund1997}, who won the 2003 Goedel Prize for their work. AdaBoost is ensemble algorithm which makes predictions based on a weighted sum of outputs of a sequence of base models in the form:
$$F_T(x)=\sum_{t=1}^T{\alpha_th_t(x)}$$
where the base model $h_t$ is also called weak learner, as it usually performs poorly but slightly better than random guessing. One typical weak learner used in AdaBoost is a decision stump which is a one-level decision tree on just one single input feature. AdaBoost constructs the strong classifier $F_t(x)$ in a iterative fashion and at each iteration $t$, a new weak learner $h_t$ and its associated coefficient $\alpha_t$ are chosen such as the empirical loss is minimized.
\begin{eqnarray*}
\mathcal{L}(y,F_t(x))&=&\sum_{i=1}^N\mathcal{L}\left(y_i,F_{t-1}(x_i)+\alpha_th_t(x_i)\right)\\
&=&\sum_{i=1}^N\text{exp}\left[-y_i(F_{t-1}(x_i)+\alpha_th_t(x_i))\right]
\end{eqnarray*}
where $F_{t-1}(x)$ is the boosted classifier that has been built up to the previous stage of training. Instead to directly minimise empirical expected prediction error $\sum_{i=1}^N\unitm(y_i\neq{F_t(x_i)})$, AdaBoost chooses to minimize its convex upbound, an exponential loss function  $\sum_{i=1}^Ne^{-y_iF_t(x_i)}$ (shown in Fig. \ref{fig:adaboostloss}).

\begin{figure}\centering
\includegraphics[width=.49\columnwidth]{images/adaboostloss.png}
\caption{Illustration of various loss functions for binary classification. AdaBoost uses the exponential loss function, which is a convex upbound of the empirical 0-1 loss.
From \cite{Murphy:2012uq}, Fig 16.9}
\label{fig:adaboostloss}
\end{figure}

\subsubsection{Training AdaBoost}
With $N$ training data points $(x_1,y_1),...,(x_N,y_N),y_i\in\{-1,1\}$, we can train AdaBoost with the following iterative updating algorithm \cite{Freund1997}:

\begin{algorithm}
    \SetKwInOut{Input}{Input}
    \SetKwInOut{Output}{Output}
    \textbf{Initialization}\;
    Weight $\omega_1(i)=1/n$ for $i=1,...,N$ (equal weights for each sample)\;
    Initiate the strong classifier $F_0(x)=0$ (we do not yet know anything)\;
    \For{$t=1:T$} {
    Fit a weak learner $h_t(x)$ to minimize weighted error $\epsilon_t=\frac{\sum_{i=1}^N\omega_t(i)\unitm(h_t(x_i)\neq{y_i})}{\sum_{i=1}^N\omega_t(i)}$\;
    Update the strong classifier $F_t=F_{t-1}+\alpha_th_t$, where $\alpha_t=\frac{1}{2}\log{\frac{1-\epsilon_t}{\epsilon_t}}$\;
    Update weight $\omega_{t+1}(i)=\frac{\omega_{t}(i)}{Z_t}\exp\left(-\alpha_ty_ih_t(x_i)\right)$, where $Z_t=2\sqrt{\epsilon_t(1-\epsilon_t)}$\;
    }
    \textbf{Return} $F_T(x)=\sgn\left[\sum_{t=1}^T \alpha_t{h_t(x)}\right]$
    \caption{AdaBoost for binary classification with exponential loss}
\end{algorithm}

An example of this algorithm in action, using decision stumps as the weak learner, is given in Figure \ref{fig:adaboosttraining}. We see that after many iterations, we can carve out a complex decision boundary.

\begin{figure}\centering
\includegraphics[width=.99\columnwidth]{images/AdaBoosttraining.png}
\caption{Example of adaboost using a decision stump as a weak learner. The degree of blackness represents the confidence in the red class. The degree of whiteness represents the confidence in the blue class. The size of the datapoints represents their weight. Decision boundary is in yellow. (a) After 1 round. (b) After 3 rounds. (c) After 120 rounds. Figure generated by boostingDemo, written by Richard Stapenhurst. From \cite{Murphy:2012uq}, Fig 16.10}
\label{fig:adaboosttraining}
\end{figure}

\begin{exercise}
Show the iterative update of AdaBoost is to minimize the exponential error $\sum_{i=1}^Ne^{-y_iF_t(x_i)}$ by fitting an additive model $F_t=F_{t-1}+\alpha_th_t$
\begin{itemize}
\item we first prove the claim $\omega_t(i){\propto}e^{-y_iF_{t-1}(x_i)}$:\\
\begin{eqnarray*}
\omega_t(i)&=&\frac{\omega_{t-1}(i)}{Z_{t-1}}\exp\left(-\alpha_{t-1}y_ih_{t-1}(x_i)\right) \\
&=&\frac{\omega_{t-2}(i)}{Z_{t-1}Z_{t-2}}\exp\left(-(\alpha_{t-1}y_ih_{t-1}(x_i)+\alpha_{t-2}y_ih_{t-2}(x_i))\right) \\
&=&\frac{1}{N}\frac{\exp\left(-y_i\sum_{j=1}^{t-1}\alpha_jh_j(x_i)\right)}{\prod_{j=1}^{t-1}Z_j}{\propto}e^{-y_iF_{t-1}(x_i)}
\end{eqnarray*}
\item given $\alpha_t$, we need to choose a weak learner $h_t$ to minimise the exponential loss, i.e.
\begin{eqnarray*}
\mathcal{L}(y,F_t)&=&\sum_{i=1}^N\exp\left[-y_i(F_{t-1}(x_i)+\alpha_th_t(x_i))\right] \\
&=&\sum_{i:h_t(x_i)=y_i}e^{-y_iF_{t-1}(x_i)}e^{-\alpha_t}+\sum_{i:h_t(x_i){\neq}y_i}e^{-y_iF_{t-1}(x_i)}e^{\alpha_t}\\
&=&\sum_ie^{-\alpha_t}e^{-y_iF_{t-1}(x_i)}+(e^{\alpha_t}-e^{-\alpha_t})\sum_{i:h_t(x_i){\neq}y_i}e^{-y_iF_{t-1}(x_i)}
\end{eqnarray*}
We note that the first term in fact does not depend on $h_t$. Hence, it suffices to minimize the second term, which is proportional to $\sum_{i=1}^N\omega_t(i)\unitm(h_t(x_i)\neq{y_i})$.
\item then we fix $h_t$, we optimize $\alpha_t$:
 \begin{eqnarray*}
\frac{\partial\mathcal{L}(y,F_t)}{\partial\alpha_t}&=&\frac{\partial}{\partial\alpha_t}\sum_{i=1}^N\exp\left[-y_i(F_{t-1}(x_i)+\alpha_th_t(x_i))\right] \\
&=&\frac{\partial}{\partial\alpha_t}\sum_{i:h_t(x_i)=y_i}e^{-y_iF_{t-1}(x_i)}e^{-\alpha_t}+\frac{\partial}{\partial\alpha_t}\sum_{i:h_t(x_i){\neq}y_i}e^{-y_iF_{t-1}(x_i)}e^{\alpha_t}\\
&=&-\sum_{i:h_t(x_i)=y_i}e^{-y_iF_{t-1}(x_i)}e^{-\alpha_t}+\sum_{i:h_t(x_i){\neq}y_i}e^{-y_iF_{t-1}(x_i)}e^{\alpha_t}=0
\end{eqnarray*}
so $e^{\alpha_t}\epsilon_t=e^{-\alpha_t}(1-\epsilon_t)\Rightarrow\alpha_t=\frac{1}{2}\log{\frac{1-\epsilon_t}{\epsilon_t}}$
\end{itemize}
\end{exercise}

\subsubsection{*A statistical view of AdaBoost}
In 2000, Friedman et al. \cite{Friedman2000} developed a statistical view of the AdaBoost algorithm. They interpreted AdaBoost as stagewise estimation procedures for fitting an additive logistic regression model, as the exponential loss function $\mathcal{L}(y,F(x))=\exp(-yF(x))$ is minimized at:
\begin{eqnarray*}
\frac{\partial}{\partial{F(x)}}\mathbb{E}\left[e^{-yF(x)}|x)\right]
&=&\frac{\partial}{\partial{F(x)}}\left[p(y=1|x)e^{-F(x)}+ p(y=-1|x)e^{F(x)}\right]\\
&=&-p(y=1|x)e^{-F(x)}+p(y=-1|x)e^{F(x)}\\
&=&0{\Rightarrow}\frac{p(y=1|x)}{p(y=-1|x)}=e^{2F(x)}
\end{eqnarray*}
So we can see that AdaBoost try to approximate (half) the log-odds ratio. AdaBoost trackles the optimal $\hat{F(x)}$ in a sequential fashion with each iteration $t$ to minimize
$$(\alpha_t,\gamma_t)=\argmin_{\alpha,\gamma}\sum_{i=1}^N\mathcal{L}(y_i,F_{t-1}(x_i)+\alpha{h(x_i;\gamma})) $$
where $\gamma$ specifies the parameters of the weaklearner $h$. When we use decision stumps as the $h$ structures, $\gamma$ will be which feature to choose and the threshold to split classes. We continue this for a fixed number of iterations T. In fact T is the main tuning parameter of AdaBoost. Often we pick it by monitoring the performance on a separate validation set, and then stopping once performance starts to decrease; this is called early stopping. Alternatively, we can use model selection criteria such as AIC or BIC (see previous section of Lasso regression for details).

\subsection{Optimal Separating Hyperplanes and Support Vector Machine}

The optimal separating hyperplane as introduced by Vapnik separates the two classes and maximises
the distance to the closest point from either class. The  optimal separating hyperplane provides a unique solution to the separating hyperplane problem. In addition,  by maximising the margin between the two classes on the training data, this leads to better classification performance on test data. We define the Margin as $M$ and obtain the following optimisation problem:
\begin{eqnarray*}
	&&\max_{\beta, \beta_0} M \\
	 \text{subject to } && \beta^T x_i + \beta _0\geq1\quad\text{if}{\quad}y_i=+1\\
          &&\beta^T x_i + \beta _0\leq-1\quad\text{if}{\quad}y_i=-1
\end{eqnarray*}

The set of conditions ensure that all the points are at least a signed distance $M$ from the decision boundary defined by $\beta$ and $\beta_0$, and we seek the largest such $M$ and associated parameters.
As the margin $M=\frac{1}{\|\beta\|}$, the above optimisation problem becomes equivalent to
\begin{eqnarray*}
	&&\min_{\beta, \beta_0} \frac12 \| \beta \|^2 \\
	 \text{subject to } && y_i (\beta^T x_i + \beta _0)>= 1 \quad i =1,\ldots, N.
\end{eqnarray*}
Here the constraints define an empty slab or margin around the linear decision boundary of thickness $1/\| \beta \|.$ Hence we choose $\beta$ and $\beta_0$ to maximise its thickness (see Figure \ref{fig:OSH}). This is a convex optimisation problem and we can used Lagrange function to solve it. For further details and the derivation we refer to \cite{Hastie:2009wp}. Samples at the margin, i.e. $\beta^T x_i + \beta _0=\pm1$ are called the \defn{support vectors}. So finding the optimal separating hyperplanes that maximize the margin is also called the \defn{support vector machine (SVM \cite{Cortes1995}) }.

To extend SVM to cases in which the data are not linearly separable, we introduce the \defn{hinge loss} function:
$$\max(0,1-y_i(\beta^Tx_i+\beta_0))$$
This function is \defn{zero if a sample satisfies $y_i(\beta^Tx_i+\beta_0)\geq1$}, i.e. at the right side of the margin.
For samples on the wrong side of the margin, the function's value is \defn{proportional to the distance from the margin}. Hence the new optimization goal is:
$$\min_{\beta,\beta_0}\left[\frac{1}{n}\sum_{i=1}^N\max(0,1-y_i(\beta^Tx_i+\beta_0))\right]+\lambda\| \beta \|^2$$
where $\lambda$ is a hyper-parameter that controls the trade-off between increasing the margin size and ensuring that the $x_i$ lies on the correct side of the margin.

So far, SVM is only applied to linear problems, but it can also be extended for nonlinear classification by using the 'kernel trick'. For further details about the 'kernel trick' we refer to Chapter 14 in \cite{Murphy:2012uq}.


\begin{figure}\centering
\includegraphics[width=.6\columnwidth]{images/OSH}
\caption{The shaded region delineates
the maximum margin separating the two classes. There are three support points
indicated, which lie on the boundary of the margin, and the optimal separating
hyperplane (blue line) bisects the slab. Included in the figure is the boundary found
using logistic regression (red line), which is very close to the optimal separating
hyperplane. From \cite{Hastie:2009wp}, Fig 4.16}
\label{fig:OSH}
\end{figure}



\subsection{Separating hyperplanes and neural networks}
%{\bf Our goal for today is to understand the basic concepts of one of the oldest machine learning methods, the perceptron learning algorithm, for the problem of separating hyperplanes.}

%We have seen that LDA and logistic regression both estimate linear decision boundaries in similar but slightly different ways. For the rest of this section we will discuss separating hyperplane classifiers.
So the separating hyperplane algorithm tries to separate the data using linear decision boundaries into different classes as well as possible. When the classes overlap, it can be generalized to support vector machine, which constructs nonlinear boundaries by constructing a linear boundary in an enlarged and transformed feature space.

%To classify points with least squares we take the sign of a linear combination of data points and assign a label equivalent to $+1$ or $-1.$

In order to understand the difference between LDA and separating hyperplane, let us consider a toy example (see Figure \ref{fig:toy_ex}) with 20 data points belonging to two classes in $\R^2.$ These data can be separated by a linear boundary. Two blue lines on the figure depicts two of infinitely many possible {\it separating hyperplanes,} we are going to consider.

\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/toy_ex_SH}
\caption{ A toy example with two classes separable by a hyperplane. The orange line is the least squares solution, which misclassifies one of the training points; whereas blue separating hyperplanes found by the perceptron learning algorithm with different random starts.
From \cite{Hastie:2009wp}, Fig 4.14.}
\label{fig:toy_ex}
\end{figure}

The orange line is the least squares solution to the problem, obtained by regressing the $-1/1$ response $Y$ on $X$ (with intercept); the line is given by
\begin{equation}
	\label{LDA_LS}
	\{ x: \hat \beta_0 + \hat \beta_1 x_1 + \hat \beta_2 x_2 = 0 \}.
\end{equation}
This least squares solution is not able to separate the points and makes one error. Recall that the use of least squares regression as a classifier, shown to be identical to LDA.
As we can see from this toy example, even when the training data can be perfectly separated by hyperplanes, LDA or other linear methods developed under a statistical framework may not achieve perfect separation. Instead of using a probabilistic argument to estimate parameters, here we consider a geometric argument.

Classifiers of the type \eqref{LDA_LS}, that compute a linear combination of the input features and return the sign, were called {\it perceptrons} in the engineering literature in the 1950s (Rosenblatt).
Linear methods all determine slightly different decision boundaries, Rosenblatt's algorithm seeks to minimize the distance between the decision boundary and the misclassified points. The concept of a perceptron was fundamental to the later development of the neural network models. The perceptron is a simple type of neural network which models the electrical signals of biological neurons. In fact, it was the first neural network to be algorithmically described.


\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/Simpleperceptron}
\caption{Diagram of a linear perceptron.
From Wikipedia.}
\label{fig:perceptron}
\end{figure}

{\bf Vector Algebra Interlude:}  A {\it hyperplane} or {\it affine set} $L$ is defined by the equation $f(x) = \beta_0 + \beta^T x = 0;$ in $\R^2$ this is a line (see Figure \ref{fig:hyperplane}).
Below we list some properties:

\begin{enumerate}
\item For any two points $x_1$ and $x_2$ lying in $L, \beta^T (x_1 - x_2) = 0,$ and hence $\beta^* = \beta / \| \beta \|$ is the vector normal to the surface of $L.$
\item For any point $x_0$ in $L$, $\beta^T x_0 = -\beta_0.$
\item The signed distance of any point $x$ to $L$ is given by
\[
	\beta^{*^T} (x-x_0) = \frac{1}{\| \beta \|} (\beta^T x + \beta_0) = \frac{1}{ \| f'(x) \| } f(x).
\]
\end{enumerate}


\begin{figure}\centering
\includegraphics[width=.7\columnwidth]{images/hyperplane}
\caption{The linear algebra of a hyperplane (affine set).
From \cite{Hastie:2009wp}, Fig 4.15.}
\label{fig:hyperplane}
\end{figure}

Hence $f(x)$ is proportional to the signed distance from $x$ to the hyperplane defined by $f(x) = 0.$

\subsubsection{Rosenblatt's Perceptron Learning Algorithm}
The {\it perceptron learning algorithm} seeks a linear boundary between two classes. A linear decision boundary can be represented by ${\beta}^T {x}+\beta_{0}$.  The algorithm begins with an arbitrary hyperplane $\underline {\beta}^T x+\beta_{0}$  (initial guess). Its goal is to minimize the distance between the decision boundary and the misclassified data points. This is illustrated in Figure \ref{fig:percept}. It attempts to find the optimal  $\beta$ by iteratively adjusting the decision boundary until all points are on the correct side of the boundary. It terminates when there are no misclassified points.

\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/Misclass}
\caption{This figure shows a misclassified point and the movement of the decision boundary.
From Wikipedia.}
\label{fig:percept}
\end{figure}

\begin{exercise}
Discuss derivation of the algorithm using the properties from vector algebra.
\end{exercise}

Consider as before $Y = \{ -1, +1 \}.$ If a response $y_i =1$ is misclassified, then $\beta^T x_i + \beta_0 <0,$ and the opposite for a misclassified response with $y_i = -1.$ The goal is to minimise
\begin{equation}
\label{perceptron}
	D(\beta, \beta_0)=-\sum_{i \in \mathcal M} y_i (\beta^T x_i + \beta_0),
\end{equation}
where
$$\mathcal M:=\left\{ i \left| \begin{array}{ll}\beta^T x_i + \beta_0 <0 & y_i=1\\ \beta^T x_i + \beta_0 >0 & y_i=-1\end{array}\right. \right\}$$
indexes the set of misclassified points. The quantity is non-negative and proportional to the distance of the misclassified points to the decision boundary defined by $\beta^T x + \beta_0 = 0.$
Note that often we simply include the intercept $\beta_0$ in $\beta$ by defining $x_0:=1$.

In order to minimize $D(\beta, \beta_0),$ we will use a {\it gradient descent}, which however my get trapped in a local minimum. To continue, the following derivatives are needed (assuming $\mathcal M$ is fixed):
\begin{eqnarray*}
	\partial \frac{D(\beta,\beta_0)}{\partial \beta} &=& - \sum_{i \in \mathcal M} y_i x_i, \\
	\partial \frac{D(\beta, \beta_0)}{\partial \beta_0} &=& - \sum_{i \in \mathcal M} y_i.
\end{eqnarray*}

The algorithm in fact uses {\it stochastic gradient descent} to minimise this piecewise linear criterion. This means that rather than computing the sum of the gradient contributions of each observation followed by a step in the negative gradient direction, a step is taken after each observation is visited. Hence, the misclassified observations are visited in some sequence, and the parameters $\beta$ are updated via:
\[
\begin{pmatrix}  \beta \\ \beta_0 \end{pmatrix} \leftarrow \begin{pmatrix}  \beta \\ \beta_0 \end{pmatrix} + \rho \begin{pmatrix}  y_i x_i \\ y_i \end{pmatrix}.
\]
Here $\rho$ is  the magnitude of each step called the 'learning rate' or the 'convergence rate'. It can be taken to be 1 without loss of generality. The algorithm continues until
$\begin{pmatrix}
  {\beta}^{\mathrm{new}}\\
  {\beta_0}^{\mathrm{new}}
\end{pmatrix}
=
\begin{pmatrix}
  {\beta}^{\mathrm{old}}\\
  {\beta_0}^{\mathrm{old}}
\end{pmatrix}$  or until it has iterated a specified number of times. If the classes are linearly separable, it can be shown that the algorithm converges to a separating hyperplane in a finite number of steps, i.e., there are no misclassified points. Figure \ref{fig:toy_ex} shows two solutions a toy problem, each started at a different random guess.

\begin{algorithm}
%\SetLine
\caption{Stochastic Gradient Descent}
{\bf Initialise} $\beta$ and $\rho;$\\
{\bf Define} a cost function $D(\beta);$\\
\Repeat{convergence criterion is fulfilled}{
Randomly permute data; \\
\For{$i=1\ldots N$}{
$\beta \leftarrow \beta - \rho \grad_\beta D(\beta);$\\
Update $\rho$
}}
\end{algorithm}

%One of the restrictions on using gradient descent is that the cost function $D(\beta)$ must be differentiable (and hence continuous). This means mis-classification cannot be used as the cost function for gradient descent schemes. Gradient descent is not usually guaranteed to decrease the cost function

{\bf Comment on gradient descent algorithm.} Consider yourself on the peak and you want to get to the land as fast as possible. So which direction should you step? Intuitively it should be the direction in which the height decreases fastest, which is given by the gradient. However, if the mountain has a saddle shape and you initially stand in the middle, then you will finally arrive at the saddle point (local minimum) and get stuck there.
In addition, note that in the final form of our gradient descent algorithm, we get rid of the summation over $i$ (all data points). Actually, this is an alternative of the original gradient descent algorithm (sometimes called batch gradient descent) known as Stochastic gradient descent, where we approximate the true gradient by only evaluating on a single training example. This means that $\beta$ gets improved by computation of only one sample. When there is a large data set, say, population database, it's very time-consuming to do summation over millions of samples. By Stochastic gradient descent, we can treat the problem sample by sample and still get decent result in practice.\\
\begin{exercise}
What are the drawbacks of the stochastic gradient descent? (e.g. same step size for all parameters, adjustment of $\rho.$) And how to choose $\rho?$
\end{exercise}

\begin{algorithm}
%\SetLine
\caption{Perceptron Algorithm}
\KwIn{linearly separable data set $x_i \in \R^d, y_i \in \{-1,+1\}$ for $i=1,\ldots,N;$}
{\bf Initialise} the vector of parameters $\beta^0 = (\beta^0, \beta_0^0);$\\
$k\leftarrow 0;$\\
\Repeat{converged}{
$k\leftarrow k+1;$ \\
$i \leftarrow k \mod N;$ \\
\eIf{$\hat{y}_i \neq y_i$}{$\beta^{k+1} \leftarrow \beta^k + y_i x_i$}{no update}
}
\end{algorithm}

{\bf There are a number of problems with the algorithm:}
\begin{enumerate}
\item  The output values of a perceptron can take on only one of two values $\{+1, -1\}$; that is, it only can be used for two-class classification.
\item When the data is separable, there are infinitely many solutions, and which one is found depends on the starting values.
\item Convergence rates depend on the size of the gap between classes. If the gap is large, then the algorithm converges quickly. However, if the gap is small, the algorithm converges slowly. This problem can be eliminated by using basis expansions technique. To be specific, we try to find a hyperplane not in the original space, but in the enlarged space obtained by using some basis functions.
\item If the data is not separable, then the perceptron algorithm will not converge since it cannot find a linear classifier that classifies all of the points correctly.
\item The speed of convergence of the algorithm is also dependent on the value of $\rho,$ the learning rate. A larger value of $\rho$ could yield quicker convergence, but if this value is too large, it may also result in 'skipping over' the minimum that the algorithm is trying to find and possibly oscillating forever between the last two points, before and after the min.
\item A perfect separation is not always available even desirable. If observations comes from different classes sharing the same input, the classification model seems to be overfitting and will generally have poor predictive performance.
\item The {\it perceptron convergence theorem} states that if there exists an exact solution (in other words, if the training data set is linearly separable), then the perceptron learning algorithm is guaranteed to find an exact solution in a finite number of steps. Note, however, that the number of steps required to achieve convergence could still be substantial, and in practice, until convergence is achieved we will not be able to distinguish between a nonseparable problem and one that is simply slow to converge
\end{enumerate}

\subsection{Multi-layer Perceptron (MLPs)}
\subsubsection{Multi-layer Perceptron}

The perceptron model presented above is very limited: it is theoretically applicable only to linearly separable
data.  In order to handle non-linearly separable data, perceptron is extended to a more complex structure,
namely multi-layer perceptron (MLP).  A MLP is a neural network in which neuron layers are stacked such
that  the  output  of  a  neuron  in  a  layer  is  only  allowed  to  be  an  input  to  neurons  in  the  upper  layer.

A neural network is a two-stage regression or classification model, typically
represented by a network diagram as in Figure \ref{fig:NN}. This network
applies both to regression or classification.

 A neural network resembles the brain in two respects:
\begin{enumerate}
\item Knowledge is acquired by the network from its environment through a learning process.
\item Interneuron connection strengths, known as synaptic weights, are used to store the acquired knowledge.
\end{enumerate}

It is a multistage regression or classification model represented by a network. It network applies both to regression or classification.
\begin{enumerate}
\item A regression problem typically has only one unit $y_1$  in the output layer, but these networks can handle multiple quantitative responses in a seamless fashion.
\item In a $k-$class classification problem, there are usually $k$ target measurements units $y_1,...,y_k$  in the output layer that each represent the probability of class $k$ and each $y_k$ is coded $(0,1)$.
\end{enumerate}

\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/NN}
\caption{Schematic of a single hidden layer, feed-forward neural network.
From \cite{Hastie:2009wp}, Fig 11.2}
\label{fig:NN}
\end{figure}

\begin{figure}\centering
\includegraphics[width=\columnwidth]{images/ANNhistory}\caption{A brief history of artificial neural networks, source: Di et al. (2018)}
\label{fig:ANNhistory}
\end{figure}



{\bf Activation Function} is a term that is frequently used in classification by NN.
In perceptron, we have a $\sign$ function that takes the sign of a weighted sum of input features. Since the sign function is not continuous at 0 and we cannot take derivative of it, we replace it by a smooth function $\sigma$ and call it the activation function.  The choice of this function $\sigma$  is determined by the properties of the data and the assumed distribution of target variables, but for multiple binary classification problems the logistic function, also known as inverse-logit (sigmoid function), is often used: $\sigm(z) = \frac{1}{1+e^{-z}}.$

An MLP is capable to capture high non-linearity of data:  one can approximate any continuous function at an arbitrary small error by using complex-enough MLPs. Clearly MLPs are highly non-identifiable.

\begin{theorem}
A two-layer MLP is a 'universal approximator', i.e. it can approximate any $\C^\infty$-function $f:\R\to\R$ of compact support to given precision for sufficiently high number of hidden neurons.
\end{theorem}
Note: easy to generalize to $f:\R^n\to\R^m$

Proof-idea: construct $f$ by approximations in local 'patches‘, see figure \ref{fig:MLPuniversal}.
WLOG, assume $f$ has support in $[0,1]$.
Choose linear output, so model reads
$y = w^\top \sigma(v x + v_0)$
with $v,v_0,w \in \R^{2n}$
Set $v:=(\lambda,\ldots,\lambda)$ sufficiently high, $v_{0,2i}:=v_{0,2i-1}:=-\lambda i/n$, and $w_{2i}:= f(i/n), w_{2i+1}:= -w_{2i}$.
Then the two hidden neurons $(2i,2i+1)$ are (approximately, $\lambda>>1$) constant equal to $f(x)$ in $[i/n, (i+1)/n]$ and $0$ elsewhere.
So $y(x)$ being the sum of all (weighted) hidden neurons approximates $f(x)$ as step function.

\begin{figure}\centering
\includegraphics[width=.55\columnwidth]{images/ANNuniversal}
\caption{Illustration of how to approximate arbitrary $\C^\infty$-functions by sufficiently complex MLPs.}
\label{fig:MLPuniversal}
\end{figure}



\subsubsection{Fitting MLPs}

Learning parameters of an MLP can again be done by MLE.
For illustration we derive the algorithm for a two-layer regression neural network. Extensions to more general cases are straight-forward, see exercises.

An MLP is simply a series of stacked logistic regression models, with the final layer being another logistic regression or a linear regression model depending on the output.
The model can be written as
\begin{eqnarray*}
P(Y=y|x,\theta) &=& \NNN(y | w^\top z(x), \sigma^2)\\
z(x) &=& g(Vx) = (g(v_1^\top x), \ldots, g(v_H^\top x))
\end{eqnarray*}
with $H$ neurons in the hidden layer, and one-dimensional activation function $g$, often taken as logistic or soft-max function. If $g$ is linear, the model simply collapses to a (possibly over-parametrized) linear/logistic regression model.
Class output as in logistic regression is modeled by
$$P(Y=c|x,\theta) = \Ber(c | \sigm(w^\top z(x))).$$

In summary, we can write the model as
$$
x\stackrel{V}{\to}a\stackrel{g}{\to}z\stackrel{W}{\to}b\stackrel{h}{\to}\hat y
$$
with intermediate (‚pre-synaptic‘) activities $a$ in the hidden and $b$ in the output layer, with $g$ and $h$ being componentwise nonlinearities i.e. $\partial g_i/\partial x_j=0$ if $i\neq j$, same for $h$. Often $h$ is chosen as identity for regression or as concatenation of logistic functions for (classification).
Altogether we can write
$$
\hat y = h(Wg(Vx)).
$$


Learning parameters in an MLP is not as easy as in (logistic) regression, since the cost function/likelihood is not convex anymore.

Given $k$-dimensional outputs, up to scaling the negative log-likelihood is then given by
$$l(\theta)= \sum_{i=1}^N \sum_{k=1}^K (\hat y_{ik}(\theta)-y_{ik})^2 =: \sum_{i=1}^N l_i(\theta) $$
or
$$l(\theta)= -\sum_{i=1}^N \sum_{k=1}^K y_{ik}\log\hat y_{ik}(\theta)=: \sum_{i=1}^N l_i(\theta) $$
%Denote by index $i$ the $i$-th sample.

Maximum likelihood is carried out numerically by gradient descent on $l$ i.e. by
$$\theta^{(k+1)}:=\theta^{(k)}-\eta^{(k)}\nabla l(\theta^{(k)})$$
Often a stochastic gradient descent approach is used i.e. gradient for a single sample $y_i$ is calculated and the parameter is updated directly. One may also simply add up the gradients for batch update.

Derivation of the single-sample gradients for the update is straight-forward but yields some interesting insights, in the following for the squared-loss:
\begin{eqnarray*}
\frac{\partial l_i}{\partial w_{km}}&=&
\frac{\partial l_i}{\partial \hat y_{ik}}
\frac{\partial \hat y_{ik}}{\partial b_{ik}}
\frac{\partial b_{ik}}{\partial w_{km}}=
-2(y_{ik}-\hat y_{ik}) h'_k(w_k^\top z_i)z_{im}\\
\frac{\partial l_i}{\partial v_{mj}}&=&
\frac{\partial l_i}{\partial \hat y_{i}}
\frac{\partial \hat y_i}{\partial b_{i}}
\frac{\partial b_i}{\partial z_{im}}
\frac{\partial z_{im}}{\partial a_{im}}
\frac{\partial a_{im}}{\partial v_{mj}}=
-\sum_{k=1}^K 2(y_{ik}-\hat y_{ik}) h'_k(w_k^\top z_i)w_{km}g'(v_m^\top x_i)x_{ij}
\end{eqnarray*}
We update the weights of the network by (possibly stochastic) gradient descent, with learning rate as discussed above.

The gradient can be rewritten as
\begin{eqnarray*}
\frac{\partial l_i}{\partial w_{km}}&=&\delta_{ki}z_{im}\\
\frac{\partial l_i}{\partial v_{mj}}&=&\gamma_{mi}x_{ij}
\end{eqnarray*}
with `errors´ $\delta_{ki}$ and $\gamma_{mi}$, the latter determined by
\begin{equation}\label{eqn:backprop}
\gamma_{mi}=g'(v_m^\top x_i)\sum_{k=1}^K \delta_{ki}w_{km}.
\end{equation}
This \defn{backpropagation equation} shows that the gradient in an MLP can be implemented by a two-pass algorithm of local (node-wise) operations, namely a forward pass with fixed weights to compute the network output $\hat y_{ik}$, and a backward pass, where first the errors $\delta_{ki}$ at the last layer are computed and then back-propagated via \eqnref{eqn:backprop}. This is also known as the delta rule and generalizes the Rosenblatt algorithm for perceptrons.

To prevent overfitting, regularization in MLPs is crucial. Commonly, a standard $\NNN(0,\alpha\unitm)$ prior i.e. $l_2$ regularization is used as in ridge regression, which is known as \defn{weight decay} in the field. The gradient can easily be adapted to this. See figure \ref{fig:MLPexample} for an example.

\begin{figure}\centering
\includegraphics[width=.49\columnwidth]{images/NN-example-nodecay}
\includegraphics[width=.49\columnwidth]{images/NN-example}
\caption{Using a two-layer MLP with $H=10$ without and with weight decay on the example from the introduction.
From \cite{Hastie:2009wp}, Fig 11.4.}
\label{fig:MLPexample}
\end{figure}

Another popular regularization method in MLPs is \defn{dropout}, which consists simply of zeroing out some fraction (typically 50\%) of the nodes in each layer before calculating the subsequent layer, as shown in figure \ref{fig:dropout}. Intuitively, \defn{dropout} forces a network\'s output not to depend too precariously on the exact activation pathway through the network, hence preventing the co-adaptation of feature detectors. Classic \defn{dropout} is only used in the model training phase but not in the testing phase. Yet recent research use dropout at test time as a heuristic approach for estimating the confidence of neural network predictions: if the predictions agree across many different dropout masks, then we might say that the network is more confident. Please refer to \cite{Gal2015Dropout} for more details.

\begin{figure}\centering
\includegraphics[width=.9\columnwidth]{images/dropout}
\caption{Dropout MLPs. Left: A standard neural net with 2 hidden layers. Right: An example of a thinned net produced by applying dropout to the network on the left. Crossed units have been dropped. From \cite{JMLR:v15:srivastava14a}, Fig 1.}
\label{fig:dropout}
\end{figure}


\subsubsection{Deep learning}
In practice neural networks have excellent classification performance if many samples are given and a strongly nonlinear structure is needed, or if special data structure properties can be used --- for example convolutional neural networks are applied to automatically learn features and then classify imaging and also sequencing data.

The key issue in learning deep MLPs was that back-propagation did not work well due to (i) the vanishing gradient problem (in early layers gradient becomes too small) and (ii) plateaus in the learning surface. This has recently been overcome by unsupervised parameter initialization known as generative pre-training as well as more recently adapted architectures (e.g. ResNet, ReLu activation function), and has resulted in an explosion of the deep learning field, see \cite{LeCun:2015dt} for a review.

We also refer to Fig \ref{fig:ANNmonitoring} for illustration on batch updates and monitoring during training of deep networks, and to Fig \ref{fig:ANNcode} for demonstration of how easy it is to setup deep NNs with current code frameworks.


\begin{figure}\centering
\includegraphics[width=\columnwidth]{images/ANNmonitoring}
\caption{Training neural networks for supervised learning, and monitoring of learning process. From \cite{Eraslan:2019bt}.}
\label{fig:ANNmonitoring}
\end{figure}

\begin{figure}\centering
\includegraphics[width=\columnwidth]{images/ANNcode}
\caption{Example code for training neural networks using Keras. From \cite{Eraslan:2019bt}.}
\label{fig:ANNcode}
\end{figure}



\subsection{Convolutional neural network (CNN)*}
% a lot of the course notes are based on Stanford CS231 course, Andrew Ng's deeplearning.ai course and d2l.ai course

In deep learning, many different network architectures, and in particular layer schemes have been proposed, see Fig {fig:ANNarchitectures} for an overview.

\begin{figure}\centering
\includegraphics[width=\columnwidth]{images/ANNarchitectures}
\caption{Neural network architectures, layers and their parameter-sharing schemes. From \cite{Eraslan:2019bt}.}
\label{fig:ANNarchitectures}
\end{figure}

Here, we focus on convolutional neural networks (CNNs), a powerful family of neural networks that were designed to process sequence data, often images. As an image consists of pixels arranged in a two-dimensional (2D) grid, applying a MLP directly on a image need to flatten the image into a 1D vector, resulting in information loss of the rich image structures. Moreover, we will get qualitatively identical results out of a MLP whether we preserve the original order of image pixels or if we permute the columns of our image pixels before learning the parameters (assume we flatten the image along columns), which certainly conflicts with our human perception of an image. Last but not least, once we start dealing with large-scale image data, the parameter size of these structure-less MLP networks can grow unwieldy. For example, using a simple three-layer MLP with 100 neurons in the hidden layer on an RGB image of size $224\times224\times3$, would lead to neurons that have $224\times224\times3\times100 = 15,052,800$ weights, which is even more than the total weights of the modern Resnet18 network with 68 layers ($11,689,512$ weights).

CNN-based network architectures now dominate the field of computer vision to such an extent that hardly anyone these days would develop a commercial application or enter a competition related to image recognition, object detection, or semantic segmentation, without basing their approach on CNNs. As shown in Fig. \ref{fig:cnn}, a typical CNN network for image classification is made up of several cascading layers, namely, convolutional layer, activation layer (typically ReLu), pooling layer (typically Max pooling) and fully-connected layer. Normalization and dropout (introduced in MLPs) are optional but also used frequently in CNN.

\begin{figure}\centering
\includegraphics[width=.99\columnwidth]{images/CNN/cnnondigital}
\caption{A typical CNN to classify handwritten digits, which includes convolutional layer, max-pooling layer, ReLU activation layer and fully-connected layer. Images are taken from \cite{towardsdatascience_Saha}.}
\label{fig:cnn}
\end{figure}


\subsubsection{From dense layers to convolutions}
Let us consider what a MLP would look like with  $I_h{\times}I_w$ images as inputs, represented as as $I_h{\times}I_w$ matrices. Let us assume the hidden representation is similarly organized as $I_h{\times}I_w$ matrices and use $x[i,j]$ and $h[i,j]$ to denote pixel location $(i,j)$ in an image and hidden representation, as:
$$h[i,j]=\sum_{k,l}W[i,j,k,l]{\cdot}x[k,l]=\sum_{a,b}V[i,j,a,b]{\cdot}x[i+a,j+b]$$
where we re-index the subscripts $(k,l)$ such that $k=i+a$ and $l=j+b$ (i.e. set $V[i,j,a,b]=W[i,j,i+a,j+b]$ ). The indices $a,b$ run over both positive and negative offsets that covers the entire image. For any given location $(i,j)$ in the hidden layer $h[i,j]$, we compute its value by summing over pixels in the input layer, weighted by $V[i,j,a,b]$.

Since we are processing image data, there are a few key principles for building neural networks, namely:
\begin{itemize}
\item \textbf{Translation Invariance:} our network should, to a certain extent, respond similarly to the same object regardless of where it appears in the image.
\item \textbf{Locality:} our network should, to a certain extent, focus on local regions, without regard for what else is happening in the image at greater distances.
\end{itemize}

The first principle \textbf{translation invariance} implies the weight $V[i,j,a,b]$ should not depend on the position $(i,j)$. The second principle \textbf{locality} suggests that we should limit our receptive field in a local region, namely, $|a|,|b|\leq\Delta$, hence we have:
$$h[i,j]=\sum_{a,b}V[a,b]{\cdot}x[i+a,j+b]=\sum_{a=-\Delta}^{\Delta}\sum_{b=-\Delta}^{\Delta}V[a,b]{\cdot}x[i+a,j+b],$$
which is a \defn{cross-correlation}. If we flip the kernel $V'[-a,-b]=V[a,b]$, we get:
$$h[i,j]=\sum_{a=-\Delta}^{\Delta}\sum_{b=-\Delta}^{\Delta}V'[a,b]{\cdot}x[i-a,j-b],$$
which is a \defn{convolution}. For a CNN, correlation is equivalent to convolution, but simply with a 'flipped' kernel, named as 'weights' in the framework of CNN. Fig. \ref{fig:crosscorrelation} illustrates a 2d cross-correlation.

\begin{figure}\centering
\includegraphics[width=.6\columnwidth]{images/crosscorrelation}
\caption{Two-dimensional cross-correlation operation. The shaded portions are the first output element and the input and kernel array elements used in its computation:  $0\times0+1\times1+3\times2+4\times3=19$. Images are taken from \cite{d2l}}
\label{fig:crosscorrelation}
\end{figure}

\begin{exercise}
Suppose we have an input of size $I_h\times{I_w}$ and kernel size of $k_h{\times}k_w$, what is the computational complexity of their convolution? \\
Answer: $(I_h-k_h+1)\times(I_w-k_w+1)\times{k_h}\times{k_w}$ multiplications, can be approximated into $I_h\times{I_w}\times{k_h}\times{k_w}$ if the kernel size is small.
\end{exercise}

\subsubsection{Multi-channel convolutions}
An RGB image consists of three channels and hence is a three dimensional tensor, represented as $x[i,j,k]$. Hence the convolutional mask has to adapt accordingly by adding a third dimension, as $V[a,b,k]$.
Moreover, we would like to use multiple convolution operations in our hidden representation, e.g. some filters are responsible to filter out noise, e.g. via Gaussian filtering, whilst others to recognize edges, hence adding a further dimension to $V$, yielding:
$$h[i,j,c]= \sum_{a=-\Delta}^{\Delta}\sum_{b=-\Delta}^{\Delta}\sum_kV[a,b,k,c]{\cdot}x[i-a,j-b,k]$$
where $k$ and $c$ correspond to the number of input and output channels, respectively. Fig. \ref{fig:depthwiseconvolution} (left) shows a diagram of a multi-channel convolution that leads to a single feature map in the hidden representation.

\begin{exercise}
Suppose we have an input 3d tensor of size $I_h\times{I_w}\times{c_\text{in}}$, we use a multi-channel convolution with $c_\text{out}$ filters, each of size of $k_h{\times}k_w$
\begin{itemize}
\item What is the number of weights in this convolution to learn? \\
Answer: $k_h{\times}k_w{\times}c_\text{in}{\times}c_\text{out}$
\item What is the computational complexity?
Answer: $I_h\times{I_w}\times{c_\text{in}}{\times}k_h{\times}k_w{\times}c_\text{out}$
\item How to reduce the number of weights for convolution? \\
Answer: use a small filter size, e.g. $3{\times}3$ instead of $11{\times}11$, separable convolution
\end{itemize}
\end{exercise}

\subsubsection{Depthwise separable convolutions}
Depthwise separable convolution was proposed in Xception and MobileNet. As shown in Fig. \ref{fig:depthwiseconvolution} (right), a depth-wise convolution consists of two essential steps:
\begin{itemize}
\item \textbf{Depthwise convolution} a 2d spatial convolution that performs independently over each channel of an input
\item \textbf{Pointwise convolution} a $1\times1$ convolution that mixes the outputs of individual channels.
\end{itemize}

\begin{exercise}
Again, suppose we have an input 3d tensor of size $I_h\times{I_w}\times{c_\text{in}}$, we use a depthwise separable convolution with $c_\text{out}$ filters:
\item What is the number of weights in this convolution to learn? \\
Answer: $k_h{\times}k_w{\times}c_\text{in}+c_\text{in}{\times}c_\text{out}$
\item What is the computational complexity?\\
Answer: $I_h\times{I_w}{\times}k_h{\times}k_w\times{c_\text{in}}+I_h{\times}I_w{\times}c_\text{in}{\times}c_\text{out}$
\item Suppose our input is RGB image of a $224{\times}224{\times}3$, we use $16$ filters of size $3\times3$ for convolution. What is the weights for regular convolution and depthwise separable convolution? How do you compare their computational cost?\\
Answer: Parameters: $3\times3\times3\times16=432$ vs. $3\times3\times3+3\times16=75$;\\
Computation cost: $224\times224\times3\times3\times3\times16\approx21.7e7$ vs. $224\times224\times3\times3\times3+224\times224\times3\times16\approx3.8e7$
\end{exercise}

\begin{figure}\centering
\includegraphics[width=.99\columnwidth]{images/depthwiseconvolution}
\caption{Standard multi-channel convolution vs depthwise separable convolution.}
\label{fig:depthwiseconvolution}
\end{figure}

\subsubsection{Padding, stride and dilation}
As shown in Fig. \ref{fig:crosscorrelation}, when we apply convolutional layers, we will lose pixels on the border of an image, and this accumulates as we apply many successive convolutional layers. One solution to this problem is to add extra pixels to fill around the boundary of our input image, known as padding, thus increasing the effective size of the image. We typically set the padded values to be $0$. If we add a total of  $p_h$  rows of padding (roughly half on top and half on bottom) and a total of  $p_w$  columns of padding (roughly half on the left and half on the right), the output shape will be $(I_h+p_h-k_h+1)\times(I_w+p_w-k_w+1)$. When we use a filter size of $3\times3$, we typically choose $p_h=p_w=2$, i.e. padding one pixel at each border location, hence maintaining the output size to be the same as the input.

In previous example of Fig. \ref{fig:crosscorrelation}, we compute the convolution with default stride of one pixel at a time. However, sometimes, either for computational efficiency or because we wish to downsample, we move our window more than one pixel at a time, skipping the intermediate locations. Fig. \ref{fig:stride} plots an example of cross-correlation with strides of 3 and 2 for height and width respectively. In general, when the stride for the height is  $s_h$ and the stride for the width is  $s_w$, the output shape is $(I_h+p_h-k_h+1)/s_h{\times}(I_w+p_w-k_w+1)/s_w.$

\begin{figure}\centering
\includegraphics[width=.69\columnwidth]{images/stride}
\caption{Cross-correlation with strides of 3 and 2 for height and width respectively. The shaded portions are the output element and the input and core array elements used in its computation:  $0\times0+0\times1+1\times2+2\times3=8$,  $0\times0+6\times1+0\times2+0\time3=6$. Images are taken from \cite{d2l}}
\label{fig:stride}
\end{figure}

Fisher et al. \cite{YuKoltun2016} introduces one more hyperparameter to the convolution layer called the \defn{dilation}. Dilation is to insert zeros within each filter. As shown in Fig. \ref{fig:dilation}, a normal convolution is equivalent to 1-dilated convolution that the filter is compact, whilst $n$-dilated convolution is to insert $n$ zeros between two non-zero filter element. The advantage of using dilation convolution is that the receptive field grows without increasing the size of the filter weights, so that large context information can be effectively captured.

\begin{figure}\centering
\includegraphics[width=.99\columnwidth]{images/dilation}
\caption{Systematic dilation supports exponential expansion of the receptive field without increasing parameter size of the filter. (1) normal convolution with a $3\times3$ filter, equivalent to 1-dilated convolution,each element in the output has a receptive field of  $3\times3$. (2) 2-dilated convolution, each output element has a receptive field of $7\times7$. (3) 3-dilated convolution, each output element has a receptive field of $15\times15$. Images are cited from \cite{YuKoltun2016}}
\label{fig:dilation}
\end{figure}

\subsubsection{Activation layer}
Just as the the \defn{Activation Function} used in MLPs, CNNs also need an activation function that transforms the output of convolution layer. A popular activation function used in CNNs is the \defn{Recitfied Linear Unit (ReLu)}, which is defined as:
$$\text{ReLu}(z)=\max(0,z).$$
Just like Sigmoid function, \defn{ReLu} is also nonlinear, as shown in Fig. \ref{fig:activation}. Additionally, \defn{ReLu} has a few advantages: i) it is faster to compute than Sigmoid, ii) it reduces problem of 'vanishing gradient' and iii) it leads to sparse neuron activations. Yet \defn{ReLu} is non-differentiable and zero gradient for the negative $z$ could causing the 'dying' \defn{ReLu} problem that it no longer respond to the variations in error/input. A solution is that we replace the horizontal line for the negative $z$ with a slightly inclined line, e.g. $f(z)=0.01z$ for $z<0$ which is called \defn{Leaky ReLu}.

\begin{figure}\centering
\includegraphics[width=.69\columnwidth]{images/activation}
\caption{Logistic Sigmoid vs ReLu activation function.}
\label{fig:activation}
\end{figure}


\subsubsection{Pooling layer}
It is common to periodically insert a pooling layer in-between successive convolution layers in a CNN architecture. The purpose of pooling is to progressively reduce the spatial size of the representation to reduce the amount of parameters and computation in the network, and hence to also control overfitting. As shown in Fig. \ref{fig:maxpooling}, one common pooling layer is called \defn{MaxPooling}, in which a filter of size $2\times2$ applied with a stride of $2$ downsamples every depth slice in the input by $2$ along both height and width. The depth dimension remains unchanged. It should also be noted that pooling layer does not contain any learnable parameters.

\begin{figure}\centering
\includegraphics[width=.69\columnwidth]{images/maxpooling}
\caption{Maximum pooling with a pooling window shape of  $2\times2$ . The shaded portions represent the first output element and the input element used for its computation:  $\max(0,1,4,5)=5$. Images are taken from \cite{d2l}}
\label{fig:maxpooling}
\end{figure}

\subsubsection{Fully connected layer}
Fully connected layer like MLPs is also used in CNN structure such as AlexNet and VGG. But more recent CNN structures, such as Resnet, do not contain fully connected layer any more (mostly because of the large size of the parameters involved in the fully connected layer).

\subsubsection{Normalization layer}
Many machine learning methods require some preprocessing to the input data so that it resembles a normal distribution, i.e. zero mean and unit standard derivation. We usually achieve it by subtracting each data point with the mean of the training data and dividing it by the standard derivation. For CNNs, the same problem appears in an intermediate layer because the input of that layer, i.e. the activation from pervious layers, is constantly changing during training. This change in the input distribution to internal layers of a CNN is known as \defn{internal covariate shift}, which could slow down the network training. The possible solution is to use a batch normalization layer in the training, which does the following:
\begin{itemize}
\item Calculate the mean and variance of the layers input on each mini batch:
\begin{eqnarray*}
\mu_B&=&\frac{1}{m}\sum_{i=1}^mx_i \quad\quad  \text{Batch mean} \\
\sigma_B^2 &=& \frac{1}{m}\sum_{i=1}^m(x_i-\mu_B)^2 \quad\quad \text{Batch variance} \\
\end{eqnarray*}
\item Normalize the layer inputs using the previously calculated batch statistics:\\
$$
\bar{x_i}=\frac{x_i-\mu_B}{\sqrt{\sigma_B^2+\epsilon}}
$$
where a small $\epsilon$ is introduced to avoid zero division.\\
\item Scale and shift to obtain the layer outputs: \\
$$y_i=\gamma\bar{x_i}+\beta.$$
Note that $\gamma$ and $\beta$ are \textbf{learned during training} along as other parameters of the network. In other words, batch normalization does the denormalization by changing only these two weights for each activation, instead of losing the stability of the network by changing all the weights.
\end{itemize}

As shown in Fig. \ref{fig:batchnormalization}, using batch normalization normally allows faster training and improves training accuracy. Additionally, as the normalization is calculated per mini-batch, this adds a bit of noisy to the network and hence batch-normalization also provides some extra regularization, which could lead to better generalization.

\begin{figure}\centering
\includegraphics[width=.49\columnwidth]{images/training_nobn}
\includegraphics[width=.49\columnwidth]{images/training_bn}
\caption{Training a CNN with ReLU activation without and with using batch normalization.
From \cite{towardsdatascience_Jansma}.}
\label{fig:batchnormalization}
\end{figure}

\subsubsection{Multitask models, multimodal models and transfer learning}
So far, we have talked only about single-task model, often simply two-class learning or univariate output $y$.  Being highly focused on our single task, we could ignore information that comes from the training other tasks. Specifically, this information comes from the training signals of related tasks might help us do even better on the original task. By sharing representations between related tasks, we can enable our model to generalize better. This approach is called \defn{Multitask Learning (MTL)}. As shown in figure \ref{fig:extension} b, multitask model has a general structure with several shared hidden layers between all tasks, while keeping output layers to be task-specific.

\begin{figure}\centering
\includegraphics[width=.9\columnwidth]{images/ANNextensions1}
\includegraphics[width=.9\columnwidth]{images/ANNextensions2}
\caption{ From \cite{Eraslan2019review}.}
\label{fig:extension}
\end{figure}

In contrast to multitask models (which have multiple outputs),\defn{multimodal models} aim to deal with different types of inputs, often different classes of multivariate inputs that should be processed together (similar to group lasso). To address this, neural networks are designed to have input and intermediate layers that are modality-specific and then only share output layers (figure \ref{fig:extension} c).

In many biological and medical applications, we have a limited amount of training data. So training a CNN from scratch may not be ideal. A popular solution is to use \defn{transfer learning}, that is, we first train a base network on a base dataset and task where the data is abundant, and then we repurpose the learned features, or transfer them, to a second target network to be trained on a target dataset and task where the data is usually sparse. This process tends to work better than training from scratch on the target dataset and task. Since the low-level features are usually more general than the high-level features, meaning suitable to multiple tasks, instead of specific to a single task. It is a common practice to freeze shallow layers and only retrain/fine tune deep layers in transfer learning. %A trained ImageNet



\subsubsection{Visualizing deep neural network}
When training a CNN to identify objects in images, we want to be able to interpret how the model make its predictions. For example, we want to explain why the network classifies a particular image as a cat. Here visualizations can uncover what a network is learning.The prevalent visualization methods can be broadly classified into two categories:
\begin{itemize}
\item Perturbation based visualizations
\item Backpropagation based visualizations
\end{itemize}
Fig. \ref{fig:featureimportance} shows a schematic of feature importance for genomics.
\begin{figure}\centering
\includegraphics[width=.9\columnwidth]{images/ANNinterpretation}
\caption{Feature importance scores highlight parts of the input most predictive for the output. From \cite{Eraslan2019review}.}
\label{fig:featureimportance}
\end{figure}
Using images, we briefly introduce several prevalent visualization methods below:
\begin{itemize}
\item \textbf{Deconvolution network:} A deconvolution network, proposed in \cite{zeiler2011}, can be thought of as a CNN model that uses the same components (filtering, pooling) but in reverse, so instead of mapping pixels to features does the opposite, i.e. map features back to pixels. Using deconvolution network, it demonstrates that a CNN learns a hierarchical feature representation of an image (figure \ref{fig:hierarchicalfeatures}): the filters in the lower layers often look like Gabor-filters which focus on local features (partially due to a small receptive field) whilst the filters in the higher layers turn to learn a more abstract representation of the entire image (thanks to an enlarged receptive field).

\begin{figure}\centering
\includegraphics[width=.9\columnwidth]{images/CNN/lowlevelfeature}
\includegraphics[width=.9\columnwidth]{images/CNN/highlevelfeature}
\caption{Visualizations of the filters learned in each layer of a CNN model with zoom-ins showing the variability of select features. From \cite{zeiler2011}.}
\label{fig:hierarchicalfeatures}
\end{figure}

\item \textbf{Occlusion} \cite{zeiler2011} also explored the technique of occluding patches of the network and monitoring the prediction and activations of feature map in the last layer. Through occlusion of different parts of the image (figure \ref{fig:occlusion}, it helps to understand importance of different regions of that image in the final prediction.

\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/CNN/occlusion}
\caption{Occlusion of different parts of images for model prediction. From \cite{zeiler2013stochastic}.}
\label{fig:occlusion}
\end{figure}

\item \textbf{Class model visualization: } Class model visualization was proposed in \cite{Simonyan2013}.The task is to generate an image which best represent a category. Mathematically, it aims to find an image $I$ which maximizes the score $S_c(I)$ that the network assigns to a category $c$. Class model visualization method starts off with some image $I$, compute the gradient with respect to $I$ using back-propagation, and then use the gradient ascent to find a better $I$, repeating the process until it finds a locally optimal $I$.This is very similar to the optimization process for training a network, but instead of optimizing the weights, it optimizes the image, keeping the weights fixed.

\begin{figure}\centering
\includegraphics[width=.9\columnwidth]{images/CNN/classmodelvisualization}
\caption{Class model visualizations for 3 different classes. From \cite{Simonyan2013}.}
\label{fig:classmodelvisualization}
\end{figure}

\item \textbf{Class saliency map:} Saliency map was also proposed in \cite{Simonyan2013}. The goal here is to find the pixels of an image which contribute most towards a particular classification. Mathematically, we take the derivative of the class score $S_c(I)$ with respect to the input image space $I$, evaluated at our input image $I_0$, i.e. $\frac{\partial{S_c}}{\partial{I}}|_{I_0}$. By contrast to the iterative updating of $I$ in the previous method, we only need to back-propagate the gradient once for each image and hence the computation time is short. As shown in Fig. \ref{fig:saliencymap}, the image specific saliency map can be used for localizing an object of interest (i.e. the dog here) without the training of any object locations. It is only given (image, category) class labels, but learns to localize: this is called weakly-supervised object localization.

\begin{figure}\centering
\includegraphics[width=0.5\columnwidth]{images/CNN/saliencymap}
\caption{Image-specific Saliency Map shows pixels that are most important for the image being classified as a dog. From \cite{Simonyan2013}.}
\label{fig:saliencymap}
\end{figure}


\item \textbf{Class activation map:} Class activation map was proposed in \cite{Zhou2016}. The idea is very similar to saliency map to locate the object of interest leading to a specific classification. While saliency map looks at the output class score with respect to the input using back-propagation, class activation map modifies the network architecture so that the forward propagation can perform both the classification and localization, as shown in Figure \ref{fig:activationmap}.

\begin{figure}\centering
\includegraphics[width=0.99\columnwidth]{images/CNN/activationmap}
\caption{The predicted class score is mapped back to the previous convolutional layer to generate the class activation maps, which highlights the class-specific discriminative region. From \cite{Zhou2016}.}
\label{fig:activationmap}
\end{figure}

\item \textbf{Visualizing the embedding of the learnt features:} We can also visualize the embedding of the learnt features using dimension reduction methods, e.g. PCA (principle component analysis) and tSNE (t Distributed Stochastic Neighbor Embedding). For example, figure \ref{fig:tsne} shows a tSNE visualization of the last-layer features of a CNN trained with image flow cytometry, which well reconstructs the cell-cycle.

\begin{figure}\centering
\includegraphics[width=0.8\columnwidth]{images/CNN/Fig3_HTML}
\caption{Cell-cyle reconstruction by tSNE visualization of the last-layer features. From \cite{Eulenberg2017}.}
\label{fig:tsne}
\end{figure}
\end{itemize}





% \begin{comment} % no GLMs here for now
\begin{comment}
\subsection{Generalized linear models (GLMs)}
Let us generalize the distribution models used so far in normal and logistic regression.

\subsubsection{Exponential family}
A pdf
$$p(x|\theta)=\frac{1}{Z(\theta)}h(x)\exp\left(\theta^\top\phi(x)\right)$$
for $x\in \mathcal{X}^m$ for some measure space $\mathcal{X}$ and $\theta\in\R^d$ with partition function
$$Z(\theta):=\int_{\mathcal{X}^m} h(x)\exp\left(\theta^\top\phi(x)\right)dx$$
is said to be of the \defn{exponential family}. $\phi(x)\in\R^d$ is called sufficient statistics, and $h(x)$ a scaling constant, often 1.
If $\phi(x)=x$, then the above defines a natural exponential family.

We can rewrite the definition as
$$P(x|\theta)=h(x)\exp\left(\theta^\top\phi(x) - A(\theta)\right)$$
with $A(\theta)=\log Z(\theta)$.

\begin{exercise}
Give some examples!
\end{exercise}

Examples are Bernoulli:
$$\Ber(x|\mu)=\mu^x(1-\mu)^{1-x}=
\exp\left(x\log\mu+(1-x)\log(1-\mu)\right)=
\exp\left(\theta^\top\phi(x)\right)
$$
with a linear however overcomplete (two-dimensional) parametrization; the minimal one is
$$\Ber(x|\mu)=\mu^x\mu^{1-x}=(1-\mu)\exp\left(x\log\frac{\mu}{1-\mu}\right)
$$
but also multinoulli distribution or uni/multivariate Gaussian (with quadratic $\phi$), however not the uniform distribution (finite, parametric support).
%$$\Cat(x|\mu)=\prod_{k=1}^K\mu_k^{x_k}=\exp\left(\sum_{k=1}^K}
%$$

Note: derivatives of the log partition function generate cumulants of the sufficient statistics. E.g. for univariate $\theta$, get
\begin{eqnarray*}
\frac{dA}{d\theta}&=&\frac{d}{d\theta} \log Z(\theta)\\
&=&\frac{\frac{d}{d\theta} \int \exp\left(\theta\phi(x)\right)h(x)dx}{Z(\theta)}\\
&=&\frac{\int \phi(x)\exp\left(\theta\phi(x)\right)h(x)dx}{\exp A(\theta)}\\
&=&\int \phi(x)\exp\left(\theta\phi(x)-A(\theta)\right)h(x)dx\\
&=&\int \phi(x)p(x|\theta)dx = E(\phi(x)),
\end{eqnarray*}
and
\begin{eqnarray*}
\frac{d^2A}{d\theta^2}
&=&\int \phi(x)\exp\left(\theta\phi(x)-A(\theta)\right)h(x)
(\phi(x)-A'(\theta))dx\\
&=&\int \phi(x)p(x|\theta)(\phi(x)-A'(\theta))dx\\
&=&E(\phi(x)^2)-E(\phi(x))^2 = \var(\phi(x)),
\end{eqnarray*}
using $A'=E(\phi(x))$;
similarly for higher order cumulants such as variance and for the multivariate setting.

Inference is done via MLE, which is unique and results in moment matching.


\subsubsection{GLMs}
Our goal is to generalize the distribution/noise models used before to the exponential family. Let us first consider the univariate case:

$$p(y_i|\theta,\sigma^2)=\exp\left(\frac{y_i\theta - A(\theta)}{\sigma^2} + c(y_i,\sigma^2)\right)$$
with dispersion parameter $\sigma$, often set to 1; here we have defined the normalizing constant by $h(y_i)=:\exp c(y_i,\sigma^2)$.

E.g. for logistic regression $\theta=\log\frac{\mu}{1-\mu}$ with $\mu:=E(y)=p(y=1)$.

Then define GLM by $\theta=\beta^\top x_i$ i.e.
$$\mu_i=g^{-1}(\beta^\top x_i)$$
with invertible link function $g$. Examples are normal distribution (link $g=\id$), binomial distribution (link equals the log odds $g=\log(\mu/(1-\mu))$) and poission distribution ($g=\log$).

So
$$\log p(y_i|x_i,\beta,\sigma^2)=\exp\left(\frac{y_i\beta^\top x_i - A(\beta^\top x_i)}{\sigma^2} + c(y_i,\sigma^2)\right)$$

Linear and binomial, in particular logistic regression are examples of GLMs, but so are also Poisson regression (with $\log p=y_i\log\mu_i - \mu_i - \log(y_i!)$).

\noindent Examples:\\
\noindent 1. linear regression can be written as
$$p(y_i|x_i,\beta,\sigma^2)=\frac{y_i\mu_i-\mu_i^2/2}{\sigma^2}- \frac{1}{2}\left(\frac{y_i^2}{\sigma^2} + \log(2\pi\sigma^2)\right)$$
with $\mu_i=\beta^\top x_i$. So $A(\theta)=\theta^2/2$, and $E(y_i)=\mu_i$. % and $\var(y_i)=\sigma^2$.

\noindent 2. poisson regression: for $y\in\N_0$, the poisson distribution is defined by
$$\log p = y\log \mu  - \mu - \log y!$$
which is also from the exponential family. Poisson regression can therefore be interpreted as a GLM with
$$\mu=\exp\beta^\top x$$
or $\theta=\log\mu= \beta^\top x$, so $A=\exp\theta$. Hence its derivates are the same and $E(y)=\var y= \mu$.

Determining the gradients is straight forward. Inference can again be done by MLE, for example in an IRLS algorithm.






%outlook towards poisson regression (one exercise to that please?)

%Multi-task learning (murphy p296)


%\subsection{Resampling: cross validation, bootstrap}
%\section{Nonlinear classification}
%\subsection{K-nearest neighbor classification}
%intro {K-nearest neighbor classification}  - briefly has been done in intro already, no just refer





\subsection{Review of nonlinear models}

\begin{itemize}
\item covered in lecture
\begin{itemize}
\item polynomial regression - just extend covariates by products
\item quadratic discriminant analysis
\item k-nearest neighbor regression/classification (prototype methods)
\item multi-layer perceptrons/neural networks
\item decision trees and random forests
\end{itemize}
\item others, not discussed:
\begin{itemize}
\item regression and smoothing splines
\item kernel methods, support vector regression/classification
\item generalized additive models (GAMs): include in additive (generalized) linear regression model smooth nonlinear functions $f_j$ of the covariates $X_j$, which are to be estimated as well. E.g. for regression
$$P(y|x,\theta)=\NNN(y|\beta_0+\sum_{j=1}^pf_j(x_j),\sigma^2),
$$
or
$$E(y|x)=\beta_0+\sum_{j=1}^pf_j(x_j),$$
similar for classification. The $f_j$'s may be modeled by any combination of the above nonlinear models (often splines are used), see \cite{Hastie:2009wp}, section 9.1. Murphy calls similar models
$$\hat f(x)=\beta_0+\sum_{j=1}^p\beta_jf_j(x_j)$$
also adaptive basis-function models (ABMs), see \cite{Murphy:2012uq}, section 16.1. Estimation is done by the so-called backfitting algorithm.
\end{itemize}
\end{itemize}

\newpage
\end{comment} 

\part{Unsupervised Learning}

\section{Mixture models}

In the previous lectures we were concerned with predicting the values of one or more outputs or response variables $Y = (Y_1, \ldots, Y_m)$ for a given set of input or predictor variables $X^T = (X_1, \ldots, X_p).$  The predictors, we have considered, are based on the training data $D = \{ (x_i, y_i) \in \R^p \times \R | i =1,\ldots,N\}$ of previously solved cases, where the joint values of all of the variables are known. We called this situation as {\it supervised learning} or {\it learning with a teacher.}

In this lecture we start with {\it unsupervised learning} or {\it learning without a teacher.} In this case, one has a set of $N$ observations ${x_1,\ldots,x_N}$ of a random $p-$vector $X$ having joint density $P(X).$ The goal is to directly infer the properties of this probability density without the help of a supervisor or teacher providing correct answers or degree-of-error for each observation. The dimension of $X$ is sometimes much higher than in  supervised learning and the properties of interest are often more complicated.

In low-dimensional case $(p \leq 3)$ there are a variety of effective nonparametric methods for directly estimating the density $P(X)$ itself at all $X-$values, and representing it graphically. However, these methods fail badly in high dimensions due to the curse of dimensionality. One must settle for estimating rather crude global methods, such as Gaussian mixtures or various simple descriptive statistics that characterise $P(X)$. The latter tries to characterise $X-$values or collections of such values, where $P(X)$ is relatively large.

Other methods, such as principal components, self-organising maps, for instance, attempt to identify low-dimensional manifolds within
the $X$-space that represent high data density. This provides information
about the associations among the variables and whether or not they can be
considered as functions of a smaller set of 'latent' variables. Cluster analysis
attempts to find multiple convex regions of the $X-$space that contain
modes of $P(X)$. This can tell whether or not $P(X)$ can be represented by
a mixture of simpler densities representing distinct types or classes of observations.
Mixture modeling has a similar goal.

Contrary to supervised learning, where a clear measure of success is provided to be used to judge adequacy in particular situation and compare the effectiveness of different methods, there is no such measure in unsupervised learning. It is difficult to ascertain the validity of inferences drawn from the output of most unsupervised learning algorithms.

\subsection{Mixture models and the EM algorithm}

{\it Latent variable models.} Unsupervised learning framework can be recasted in the context of {\it latent variable models (LVM),} i.e.,  models with hidden variables. In general, there are $L$ latent variables $z_{1},\ldots, z_{L},$ and $p$ visible variables $x_{1},\ldots,x_{p},$ where usually $p>>L.$ Figure \ref{fig:LVM} illustrates some generic LVM structures. In addition, depending on the form of likelihood $P(x|z)$ and the prior $P(z),$ one can generate a variety of different models. In the following we discuss in details the simplest form of LVM, the so-called {\it mixture models}, when $z_i \in \{1,\ldots, K\},$ representing a discrete latent state.


\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/LVM}\\
\caption{A latent variable model represented as directed graph model \cite{Murphy:2012uq}, Fig 11.2. (a) Many-to-many (b) One-to-many (c) Many-to-one (d) One-to-one}
\label{fig:LVM}
\end{figure}

 Let
\[
P(z) = \pi, P(x|z=k) = p_k(x),
\]
where $p_k$ is the $k'$th base distribution for the observations, which can be of any type. The {\it mixture model} is given as a convex combination of $p_k'$s
\begin{equation}
	\label{mix_model}
	P(x|\theta) =\sum_{k=1}^K \omega_k p_k(x | \theta_k),
\end{equation}
where $0\leq \omega_k \leq 1$ and $\sum_{k=1}^K \omega_k =1.$

Examples are:
\begin{enumerate}
\item {\it Histogram} with piecewise constant
$$p_k(x| \theta_k) = c_k 1_{R_k} $$
a histogram is simply a sum of unit functions localized in their discrete bins $R_k$.

\item {\it Mixtures of Gaussians or Gaussian mixture models (GMM)} is the most widely used mixture model. In this model, each base distribution in the mixture is a multivariate
Gaussian with mean $\mu_k$ and covariance matrix $\Sigma_k,$ i.e.
$$p_k(x| \theta_k) = \NNN(x | \mu_k, \Sigma_k)$$
Then
\begin{equation}
	\label{GMM}
	P(x|\theta) =\sum_{k=1}^K \omega_k \NNN(x | \mu_k, \Sigma_k),
\end{equation}

Figure \ref{fig:GMM} shows a mixture of 3 Gaussians in 2d/ Each mixture component is represented by a different set of elliptical contours. Given a sufficiently large number of mixture components. A GMM can be used to approximate any density defined on $\R^N.$

\begin{figure}\centering
\includegraphics[width=.8\columnwidth]{images/GMM}\\
\caption{A mixture of 3 Gaussians in 2d from \cite{Murphy:2012uq}, Fig 11.3. (a) Contours of constant probability for each component in the mixture (b) A surface plot of the overall density.}
\label{fig:GMM}
\end{figure}

\item {\it Mixtures of Multinoullis} used to define density models on many kinds of data. In this case, the model is given as a product of Bernoulli densities:
\begin{equation}
	\label{MM}
	P(x|z, \theta) =\prod_{j=1}^p \Ber(x_{j}| \mu_{jk}) = \prod_{j=1}^p \mu_{jk}^{x_{j}} (1-\mu_{jk})^{1-x_{j}},
\end{equation}
where $x_i \in [0,1]^p$ is a bit vector, $\mu_{jk}$ is the probability that bit $j$ turns on in cluster $k.$
\end{enumerate}

\subsubsection{Application of mixture models}

There are two main application of mixture models. The first is to use them as flexible {\it black-box} density model, $P(x).$ This can be useful for variety of tasks such as data compression, outlier detection, nonlinear LDA/QDA etc. The second, more common, application is to use them for clustering.

Shortly, we first fit the mixture model and then compute the posterior probability of the point to belong to cluster $k,$ i.e., $P(z = k|x,\theta).$ This is also known as {\it responsibility} of cluster $k$ for the point and is computed using Bayesian rule:
\begin{equation}
	\label{SC}
r_{k} = P(z=k | x, \theta) = \frac{P(z = k | \theta) P(x | z = k, \theta)}{\sum_{k'=1}^K P(z=k'| \theta) P(x | z=k', \theta)}.
\end{equation}
In the above mixture model, $P(z = k | \theta)=\omega_k$.
The  procedure is also called {\it soft clustering} and is identical to the computations performed when using a generative classifier.

We can represent the amount of uncertainty in the cluster assignment by using $1- \max_k r_{k}.$ Assuming this is small, it may be reasonable to compute a {\it hard clustering} using the MAP estimate, given by
\begin{equation}
	\label{HC}
z^* = \argmax_{k} r_k = \argmax_k \log P(x | z = k, \theta) + \log P (z = k | \theta).
\end{equation}

However, it is known that mixture models can produce quite a lot of errors in the clustering task. Therefore, one should be quite careful in interpreting the results and, if possible, add little bit of supervision, or using informative priors.

\subsubsection{The EM Algorithm}

As we have seen for many models computing the ML or MAP parameter estimate is easy provided all the values of all the relevant random variables are observed, i.e., if we have complete data. However, if some data is missing and/or contains latent variables, then computing the ML/MAP estimate becomes hard.

One approach is to use a generic gradient-based optimizer to find a local minimum of the
negative log likelihood or NLL, given by
\begin{equation}
	\label{NLL}
NLL (\theta) = - \frac{1}{N} \log P(D | \theta).
\end{equation}
However, we often have to enforce constraints, such as the fact that covariance matrices must be positive definite, mixing weights must sum to one, etc., which can be tricky. In
such cases, it is often much simpler (but not always faster) to use an algorithm called {\it expectation maximization (EM)}. This is a simple iterative algorithm, often with closed-form updates at each step. Furthermore, the algorithm automatically enforces the required constraints.

EM exploits the fact that if the data were fully observed, then the ML/ MAP estimate would be easy to compute. In particular, EM is an iterative algorithm which alternates between inferring the missing values given the parameters (E step), and then optimising the parameters given the 'filled in' data (M step).

{\bf EM Basic idea.} We are interested in maximising the log likelihood of the observed data
\[
	l (\theta) =  \sum_{i=1}^N \log P(x_i | \theta) = \sum_{i=1}^N \log \left[  \sum_{z_i} P(x_i, z_i | \theta)\right].
\]
Due to the optimisation difficulties, i.e., $\log$ cannot be pushed inside the sum, we define the complete data log likelihood
\[
	l_c (\theta) = \sum_{i=1}^N \log P(x_i, z_i | \theta).
\]
Since $l_c$ cannot be computed ($z_i$ is unknown), we define the expected complete data log likelihood as follows:
\[
	Q(\theta, \theta^{t-1}) = E [l_c(\theta) | D, \theta^{t-1}],
\]
where $t$ is the current iteration number, $Q$ is called the auxiliary function. The goal of the E step is to compute $Q(\theta, \theta^{t-1})$ with $z_i$ estimated in previous Iteration $t-1$, or rather, the terms inside of which the MLE depends on it; these are known as the expected sufficient statistics (ESS). In the M step, the $Q$ function is optimised wrt $\theta$:
\[
	\theta^t = \argmax_\theta Q(\theta, \theta^{t-1})
\]
or
\[
	\theta^t = \argmax_\theta Q(\theta, \theta^{t-1}) + \log P(\theta).
\]
The EM algorithm monotonically increases the log likelihood of the observed data (plus the log prior, if doing MAP estimation), or it stays the same. So if the objective ever goes down, there must be a bug in a code. (This is a surprisingly useful debugging tool!).

Below we explain how to fit a mixture of Gaussians using EM.

{\bf EM for GMMs.} We assume the number of mixture components, $K,$ is known. The expected complete data log likelihood is given by
\begin{eqnarray*}
Q(\theta, \theta^{t-1}) & = & E \left[   \sum_i \log P(x_i, z_i | \theta)\right] \\
& = & \sum_i E \left[  \log \left(\prod_{k=1}^K (\omega_k P(x_i | \theta_k))^{\mathbb I(z_i = k)}   \right) \right] \\
& = & \sum_i \sum_k E [\mathbb I (z_i = k)] \log(\omega_k P(x_i | \theta_k)) \\
& = & \sum_i \sum_k P(z_i = k | x_i, \theta^{t-1})  \log(\omega_k P(x_i | \theta_k)) \\
& = & \sum_i \sum_k r_{ik} \log \omega_k + \sum_i \sum_k r_{ik} \log P(x_i | \theta_k),
\end{eqnarray*}
where $r_{ik} = P(z_i = k | x_i, \theta^{t-1})$ is the responsibility that cluster $k$ takes for data point $i$. This is computed in the E step, described below.

{\it The E Step} has the following simple form, which is the same for any mixture model:
\[
 	r_{ik} = \frac{\omega_k^{t-1} P(x_i | \theta_k^{t-1})}{\sum_{k'} \omega_{k'}^{t-1} P(x_i | \theta_{k'}^{t-1})}.
\]
In the {\it M Step}, we optimise $Q$ wrt $\omega$ and the $\theta_k.$ For $\omega$ we have:
\[
	\omega_k = \frac1N \sum_i r_{ik} = \frac{r_k}{N},
\]
where $r_k = \sum_i r_{ik}$ is the weighted number of points assigned to cluster $k.$ Now we derive the M step for the $\mu_k$ and $\Sigma_k:$
\begin{eqnarray*}
	l(\mu_k, \Sigma_k) & = & \sum_k \sum_i r_{ik} \log P(x_i | \theta_k) \\
		& = & -\frac12 \sum_{i} r_{ik} [ \log |\Sigma_k| + (x_i -\mu_k)^T \Sigma_k^{-1} (x_i - \mu_k)].
\end{eqnarray*}

One can show that
\begin{eqnarray*}
\mu_k & = & \frac{\sum_i r_{ik} x_i}{r_k} \\
\Sigma_k & = & \frac{\sum_i r_{ik} (x_i - \mu_k) (x_i -\mu_k)^T}{r_k} = \frac{\sum_i r_{ik} x_i x_i^T}{r_k} - \mu_k \mu_k^T.
\end{eqnarray*}

The intuition behind these equations is the following: the mean of cluster $k$ is just the weighted average of all points assigned to cluster $k,$ and the covariance is proportional to the weighted empirical scatter matrix.

After computing the new estimates, we set $\theta^t = (\omega_k, \mu_k, \Sigma_k)$ for $k=1,\ldots,K$ and return to the step E.


{\bf $K-$means algorithm} is a popular modification of the EM algorithm for GMMs. Let $\Sigma_k = \sigma^2 \mathbf I_D$ and $\omega_k = \frac1K$ is fixed, so only the cluster centres, $\mu_k \in \R^p$ have to be estimated. Now let's search for 'hard clusters' i.e. consider the following delta-function approximation to the posterior computed during the E step:
\[
P(z_i = k | x_i, \theta) \approx \mathbb I(k = z_i^*),
\]
where $z_i^* = \argmax_k P(z_i = k | x_i, \theta).$ This is sometimes also called hard EM,  since we are making a hard assignment of points to clusters. Since we assumed an equal  covariance matrix for each cluster, the most probable cluster for $x_i$ can be computed as:
\[
z_i^* = \argmin_k \| x_i - \mu_k \|^2_2.
\]

Hence in each E step, we must find the Euclidean distance between $N$ data points and $K$ cluster centres, leading to $O(NKD)$ time. This can be speed up using various techniques (see Elkan 2003 for more details). Given the hard cluster assignments, the M step updates each cluster center by computing the mean of all points assigned to it:
\[
	\mu_k = \frac{1}{N_k} \sum_{i: z_i = k} x_i.
\]


\begin{algorithm}
\text{Initialize} $\mu_k;$ \\
\Repeat{converged}{
$z_i = \argmin_k \| x_i -\mu_k \|^2_2; $\\
$\mu_k = \frac{1}{N_k} \sum_{i: z_i = k} x_i;$
}
\caption{$K-$means algorithm}
\end{algorithm}


% possible extensions:
% k-medoids, k-median, projective k-means etc
% vector quantization (murphy 11.4.2.6)
% initialization
% nice: proof of conv/theoretical basis for EM 11.4.7

\subsection{Convergence of EM}
One can show that EM monotonically increases the observed data likelihood $P(\Data|\theta)$ until a local maximum (or saddle point, which are commonly unstable) is reached; similarly for $k$-means. For a brief proof we refer to \cite{Murphy:2012uq}, section 11.4.7.

\subsection{General comments on clustering}
$k$-means is a method for partitional clustering. A very popular and simpler method for clustering is hierarchical clustering, which aggregated samples by similarity, building a similarity tree, which may be cut off a certain depth to generate clusters.


\section{Dimension reduction}


Goal: represent $X$ by lower dimensional `representative´ random vector. Since this representation is a prior unknown, this is also called a \defn{latent variable model}. Here we will initially discuss the linear case only.

\subsection{Factor analysis}
Mixture models can be interpreted as models with a single discrete (or multivariate binary) latent variable to generate observations, selected from $k$ prototypes.
Instead we now assume $k$-dimensional continuous prototypes  encoded by the latent random vector $Z\in\R^n$ with normal prior
$$p(z)=\NNN(z|\mu_0,\Sigma_0),$$
and usually $n<p$.

For continuous observations, as in linear regression (and mixture modeling), we represent the data $X$ by Gaussian, where the mean is modeled as linear superposition of the latent inputs
$$p(x|z,\theta)=\NNN(x|\W z+\mu,\Sigma).$$
This model is called \defn{factor analysis (FA) model}, the $p\times n$-matrix $\W$ is called \defn{factor loading matrix}. The noise correlation $\Sigma$ is assumed to be diagonal to make the model explain the observed correlation  by $Z$. If $\Sigma=\sigma^2\unitm$, the model is called \defn{probabilistic principal component analysis (PPCA)}, see later and illustration in figure\ref{fig:PPCA_illustration}. In other words
$$X=\W Z + \epsilon$$
with $\epsilon\sim\NNN(\mu,\sigma^2\unitm)$.

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/PPCA_illustration}
\caption{Illustration of PPCA data generation.
From \cite{Murphy:2012uq}, Fig 12.1.}
\label{fig:PPCA_illustration}
\end{figure}

\begin{exercise}
How does the full observation density $p(x|\theta)$ look like?
\end{exercise}

FA models the density by a lower number of parameters, in particular
\begin{eqnarray*}
p(x|\theta) &=& \int \NNN(x|\W z+\mu,\Sigma)\NNN(z|\mu_0,\Sigma_0) dz\\
&=& \NNN(x|\W\mu_0+\mu, \Sigma + \W\Sigma_0\W^\top)
\end{eqnarray*}
i.e. the marginal is again Gaussian.

\begin{exercise}
What does this imply for identifiability?
\end{exercise}

Hence, WLOG assume $\mu_0=0$ or $\mu=0$ and $\Sigma_0=\unitm$ --- otherwise rescale $\tilde\mu=\mu-\W\mu_0$, and similarly
$\tilde\W=\W\Sigma_0^{-1/2}$ since then
$$\Cov(\tilde x|\theta)=\tilde\W\Cov(z|\theta)\tilde\W^\top + \Cov(\epsilon|\theta)=
\W\Sigma_0^{-1/2} \Sigma_0 \Sigma_0^{-1/2}\W^\top+\Sigma =
\W\W^\top+\Sigma $$

\begin{exercise}
That's it?
\end{exercise}

There is another unidentifiability: for any $\M\in O(n)$, i.e. $\M\M^\top=\unitm$, and $\tilde\W=\W\M$
$$\Cov(\tilde x)=\tilde\W\tilde\W^\top + \Sigma
=\W\M\M^\top\W^\top + \Sigma = \W\W^\top + \Sigma = \Cov(x)
$$

Finally also the diagonal part of $\Sigma$ may be absorbed in $\W$.

Identifiability may be enforced by
\begin{itemize}
\item (i) Assuming $\W\in O(p)$ or $\W^\top\W=\unitm$ (i.e. from Stiefel manifold) and $\Sigma=\unitm$, with order by decreasing variance of latent components (so $\Sigma_0$ is only assumed to be diagonal not one), see PCA later.
\item Regularization: sparse priors on $\W$ or other structural assumptions.
\item non-Gaussian prior for latent factors: see independent component analysis later.
\end{itemize}

\begin{exercise}
How to do inference?
\end{exercise}

Depending on the regularization/additional assumption, (frequentist) inference differs. The MLE estimate for the normal distribution boils down to estimating $\W$ and $\Sigma$ from the covariance estimate
$$\hat\C=\W\Sigma_0\W^\top+\Sigma,$$
which for assumption (i) from above can be done by eigenvalue decomposition, see also PCA.

Also EM algorithm can be derived for the estimation, which amongst others allows efficient coping with missing values.

Posterior inference i.e. $p(s|x)$ can be done by Bayes rules, which for Gaussians is nice since it stays in the class.


\subsection{Principal component analysis}
FA with $\Sigma=\sigma^2\unitm$ leads to the probabilistic PCA model for $\sigma^2>0$. For vanishing noise, it is simply known as PCA.

PCA, also known as Karhunen-Loève transformation is a common multivariate data analysis tools; it aims to transform data to a feature space, where few `main features' (principal components) explain most of the data in terms of covariance.

WLOG assume centered data i.e. vanishing mean of $X$ and hence of $Z$.
Classically (`synthesis view'), PCA is introduced as a dimension reduction technique minimizing the expected reconstruction error
$$J(\W,Z)=E(\|X-\hat X\|^2),$$
with $\hat X=\W Z$ subject to (left)-orthonormal $\W\in\R^{p\times n}$ i.e. $\W^\top\W=\unitm$.

In figure \ref{fig:PCA_digits}, the reconstruction process is visualized for three examples.

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/PCA/digits_pca}
\caption{Reconstruction $\hat X$ of 3 images from NIST image data base using PCs trained from 1000 samples each  after centering, with varying number $m$ (should be $n$ for this text, sorry) of principal components.
% see pca_digits.m
}
\label{fig:PCA_digits}
\end{figure}

This idea of compression may be used for filter bank learning in natural images and image compression, see figure \ref{fig:PCA_jakob}.

\begin{figure}\centering
\includegraphics[width=.5\textwidth]{images/PCA/filters_from_jakob}\\
\includegraphics[width=.49\textwidth]{images/PCA/jakobrec_1pcs}
\includegraphics[width=.49\textwidth]{images/PCA/jakobrec_2pcs}\\
\includegraphics[width=.49\textwidth]{images/PCA/jakobrec_4pcs}
\includegraphics[width=.49\textwidth]{images/PCA/jakobrec_32pcs}
\caption{PCA filters estimated from 10.000 random 32x32  image patches from the original image (centered).
Reconstructions $\hat X$ shown for $n=1,2,4,32$  principal components.
% see jakob.m
}
\label{fig:PCA_jakob}
\end{figure}

In terms of samples this is equivalent to
$$J(\W,\Z)=\frac{1}{N}\|\X-\Z\W^\top\|_F^2$$
with low-dimensional encoding or principal components $\Z\in\R^{N\times n}$, and squared Frobenius norm (i.e. sum of squared matrix entries).

\begin{exercise}
How to minimize this?
\end{exercise}

\begin{theorem}
$J(\W,\Z)$ is minimal at $(\hat\W, \X\hat\W)$ with $\hat\W$ containing the first $n$ eigenvectors of the empirical covariance $\hat\Sigma$ of $X$ corresponding to the $n$ largest eigenvalues (with multiplicities).
\end{theorem}

\begin{proof}
Either from linear algebra, or explicitly:

Let us start with $n=1$ or the first column $w_1$ of $\W$ (analysis view of PCA).
In figure \ref{fig:firstpca_visualize}, we visualize reconstruction of the first PC in two dimensions.

The reconstruction error is given by (with $x_i$ $i$-th row of $X$)
$$J(w_1,z_1)=\frac{1}{N}\sum_{i=1}^N \|x_i-z_{i1}w_1^\top\|^2
= \frac{1}{N}\sum_{i=1}^N (x_ix_i^\top-2z_{i1}w_1^\top x_i^\top + z_{i1}^2 )
$$
since $\|w_1\|=1$. This is optimal at the orthogonal projection $z_{i1}(w_1)=w_1^\top x_i^\top$. Hence
$$J(w_1)=\frac{1}{N}\sum_{i=1}^N (x_i x_i^\top-z_{i1}^2 )=const-\var(z_1)
$$
In particular minimizing the reconstruction error is equivalent to maximizing the variance of the projection i.e. the measure for meaningful dimension reduction in PCA is maximization of variation.

The quadratic form
$$\var(z_1)=w_1^\top\hat\Sigma w_1$$
with symmetric $\hat\Sigma$ is maximized, subject to unit norm $w_1$, at eigenvector corresponding to largest eigenvalue.
(why? see directly e.g. by Lagrange method)

For $n=2$ (or any $n$ by induction), we pick $w_1$ as above, and now search for $w_2$ such that $w_j^\top w_k=\delta_{jk}$. The error
$$J(w_1,z_1,w_2,z_2)=\frac{1}{N}\sum_{i=1}^N \|x_i-z_{i1}w_1^\top-z_{i2}w_2^\top\|^2$$
gives same solution for $w_1$ as before, and for $w_2$ as before the projection $z_{i2}=w_2^\top x_i^\top$.
Hence
$$J(w_2)=const - \var(z_2)= const - w_2^\top\hat\Sigma w_2$$
The latter quadratic form subject to orthogonality on the first eigenvector is maximized by the second one.
\end{proof}

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/PCA/firstpca_visualize}
\caption{Visualization of the reconstruction of the first PC for a two-dimensional correlated Gaussian (left), with cost function $J$ (right).}
\label{fig:firstpca_visualize}
\end{figure}


\subsubsection{Link to singular value decomposition}
Instead of eigenvalue decomposition of the covariance matrix, we can instead perform an $n$-dimensional singular value decomposition (SVD) of $\X$:
$$\X=\U\S\V^\top$$
with diagonal $\S\in\R^{n\times n}$ and (left)-orthonormal $\U$ and $\V$. In particular this is an `economy-sized SVD´ i.e. zero-diagonal elements in a full SVD would have been pruned away. This makes this decomposition quite efficient to calculate, $O(Nn \min(n,N))$. Also no potentially high-dimensional square matrices need to be stored (in contrast to covariance calculation). The covariance decomposition can then be easily constructed by
$$N\hat\C= \X^\top\X=\V\S^\top\U^\top\U\S\V^\top=\V\D\V^\top$$
with $\D=\S\S$, so $\hat\W=\V$.

In particular
$$\hat\Z=\X\hat\W=\U\S\V^\top\V=\U\S$$
and
$$\hat\X=\Z\W^\top=\U\S\V^\top$$
so the SVD is the LS-optimal approximation of $\X$.

\subsubsection{Probabilistic PCA}
Probabilistic PCA assumes the simplified FA model with $\Sigma=\sigma^2\unitm$ and orthogonal $\W$ i.e.
$$p(x|z,\theta)=\NNN(x|\W z,\sigma^2\unitm)$$
and
$$p(z)=\NNN(z|0,\unitm).$$
The log-likelihood is given by
$$
\log p(\Data|\theta)=
-\frac{N}{2}\log\det\C-\frac{1}{2}\sum_{i=1}^N x_i^\top\C x_i
$$
with $\C=\W\W^\top+\sigma^2\unitm$. Similar to PCA one can show that its maxima are given by
$$\hat\W=\V(\D-\sigma^2\unitm)^{1/2}\M$$
with $\V$ containing the first $n$ eigenvectors of the empirical covariance matrix of $X$ corresponding to the largest $n$ eigenvalues stored in $D$, and $\M\in O(n)$.
The MLE of $\sigma$ is given by
$$\hat\sigma^2=\frac{1}{p-n}\sum_{i=n+1}^p d_i,$$
i.e. the average variance associated with the discarded dimensions.

Hence this coincides with PCA for $\sigma\to0$.

Note that extension of this model to categorical data $x$ is done as for logistic regression by using the softmax function and modeling
$$p(x|z,\theta)=\prod_{c=1}^C \Ber\left(x_c\left|
\frac{\exp(\beta_c^\top\z)}{\sum_{d=1}^C\exp(\beta_d^\top\z)}
\right.\right)$$
with normal prior on $z$ and potential intercept $(\beta_c)_0$.
Inference is usually done by EM.

\subsection{Choosing the number of dimensions}
Just as in mixture modeling (in particular clustering) or regression, one may ask how to determine the complexity parameter number of components $n$.

As before the answer should be to deal with the variance-bias-tradeoff by evaluation on a test set (or cross validation etc).

However PCA directly does not lead to the expected increase in accuracy on a test set, see figure \ref{fig:PCAoverfit}, top row. Here the root mean square error $\sqrt{J}$ is plotted --- similarly often the sum of the last $d-n+1$ eigenvalues is shown. The RMSE does not increase on the test data set since the PCA cost function has not been derived as a probabilistic problem, it just compresses the data, and this always gets better also for training data not perfectly in line with the original data.

If instead PPCA with MLE is analyzed, we get the expected U-shaped curves on test data, set \ref{fig:PCAoverfit}, bottom row.

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/PCA/PCAoverfit}\\
\includegraphics[width=\textwidth]{images/PCA/PCAoverfitMLE}
\caption{Selecting the number of PCs $n$ for the MNIST data set. Top row: PCA root mean-square error $\sqrt{J}$ for training and test set.
Bottom row: negative log likelihood of PPCA on training and test data set.}
\label{fig:PCAoverfit}
\end{figure}


There are other extensions of PCA such as partial least squares or canonical correlation analysis, please see \cite{Murphy:2012uq}, sec. 12.5.1ff.
%\subsection{PLS}
%and latent factor regression?
%CCA and even
%murphy p406ff

%\bibliographystyle{plain}\bibliography{paperslibrary}
%\end{document}




\subsection{Autoencoders}
An \defn{autoencoder} is a kind of unsupervised neural network that is used for dimension reduction and feature representation. As in synthesis view of PCA, it is a processing unit trained to match its output to the given input, however now with the added flexibility of a nonlinear function. As shown in Figure \ref{fig:autoencoder}, an autoencoder consists of an encoder function $h=f(x)$ maps the inputs into hidden representation whilst an decoder $y=g(h)$ tries to reconstruct the input so that $g(f(x)){\approx}x$. We use $\approx$ rather than $=$ here as autoencoders are designed to be unable to learn to copy perfectly, i.e. prevent autoencoders from learning the trivial identity function. Usually we constrain the hidden layer $h$ in the middle to be a narrow so-called \defn{bottleneck layer} to allow autoencoder to copy only approximately, and the reconstruction error is quantified by some loss $L(x,g(f(x)))$.

\subsubsection{Training autoencoders}
Training is done as usual in neural networks by (mini-batch) gradient descent on the loss. When using deep multi-layer perceptrons as encoder and decoder, one could potentially face the \textbf{"vanishing gradient"} problem. As explained before, we can circumvent this by using pretraining (as explained in Figure \ref{fig:autoencoder}), ReLu activation, batch normalization, or replacing dense connection with convolutional layers (e.g. for images) to build deep an autoencoder.

\begin{exercise}
Prove one hidden layer autoencoder with linear activation is equivalent to PCA, more precisely, the weights of a linear autoencoder span the same subspace as the principal components found by PCA?
\end{exercise}
Let x be the input, we have our one-layer hidden representation of linear autoencoder as:
\begin{eqnarray}
h&=&W_1x+b_1\\
y&=&W_2h+b_2
\end{eqnarray}
To train the autoencoder, we minimize our reconstruction error:
$$L(x,y)=\min_{W_1,b_1,W_2,b_2}\|x-(W_2(W_1x+b_1)+b_2)\|^2$$
Setting partial derivatives with respect to $b_i$, $i=1,2$ respectively, to zero and inserting the solution into the above equation, simplifies the optimization as follows:
$$\min_{W_1,W_2}\|x-W_2W_1x\|^2$$
This includes a classical linear inverse problem, and can be solved as follows: Setting again the partial derivative with respect to $W_1$ to zero, $W_1$ is the left Moore-Penrose pseudoinverse of $W_2$ (assuming WLOG full rank of $W_2$), yielding:
$$W_1=W_2^{\dagger}=(W_2^TW_2)^{-1}W_2^T$$
and the optimization:
$$\min_{W_2}\|x-W_2W_2^{\dagger}x\|^2$$
We can see the second term $W_2W_2^{\dagger}x$ is the orthogonal projection of $x$ onto the column space of $W_2$, where its columns are not necessarily orthonormal. This problem is very similar to PCA and its synthesis view, but without the orthonormality constraint. In particular, this implies that although the weights of a linear autoencoder span the same subspace as the principal components found by PCA, the vectors at the hidden layers are not necessarily the same as the PCA eigenvectors (the later ones are orthonormal by definition).

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/Deep_autoencoder}\\
\caption{Training a deep autoencoder. (a) First we greedily train some random Boltzmann machines (RBMs). (b) Then we construct the auto-encoder by replicating the weights. (c) Finally we fine-tune the weights using back-propagation. From Figure 28.3 in \cite{Murphy:2012uq}.}
\label{fig:autoencoder}
\end{figure}

\subsubsection{Sparse autoencoder*}
A \defn{sparse autoencoder} is an autoencoder whose training criterion involves a sparsity penalty $\Omega(h)$ on the hidden layer $h$ in addition to the normal reconstruction error.
$$L(x,g(f(x)))+\Omega(h),$$
The sparse regularization term forces only a small number of the hidden units to be active at once (shown in Figure \ref{fig:sparseautoencoder}).

The sparse autoencoder framework can be considered as approximating maximum likelihood training of a generative model with visible variables $x$ and latent variables $h$, as:
$${\log}p_\text{model}(x)=\log\sum_hp_\text{model}(h,x)$$
We can think of the autoencoder as approximating this sum with a point estimate for just one highly likely value for $h$. With this chosen $h$, we are maximizing
$${\log}p_\text{model}(h,x)={\log}p_\text{model}(h)+{\log}p_\text{model}(x|h)$$
we further make a spare regularization on $h$, e.g. using the Laplace prior,
$p_\text{model}(h_i)=\frac{\lambda}{2}e^{-\lambda\|h_i\|}$ corresponds to the lasso sparsity penalty, as
\begin{eqnarray}
\Omega(h)&=&\lambda\sum_i\|h_i\|, \\
-{\log}p_\text{model}(h)&=&\sum_i\left(\lambda|h_i|-\log\frac{\lambda}{2}\right)=\Omega(h)+\text{const}
\end{eqnarray}

Due to the penalty function, the dimension of the bottleneck layer can also be higher than the input. A very successful method for one-hidden-layer sparse coding for image denoising is \texttt{dictionary learning}, with an example given in Figure \ref{fig:dictionarylearning}.


\begin{figure}\centering
\includegraphics[width=0.5\textwidth]{images/autoencoder/sparseautoencoder}\\
\caption{Only a few hidden units are activated in sparse autoencoder.}
\label{fig:sparseautoencoder}
\end{figure}

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/dictionarylearning}\\
\caption{Image denoising with dictionary learning}
\label{fig:dictionarylearning}
\end{figure}

\subsubsection{Denoising autoencoder*}
The objective of denoising autoencoders is to clean the corrupted input, i.e., denoising. To achieve this aim, \defn{denoising autoencoders} (DAE) intentionally add noise to the raw input before providing it to the network. The training process of a denoising autoencoder works as follows:
\begin{itemize}
\item The initial input $x$ is corrupted into $\tilde {x}$ through a corruption process $C(\tilde{x}|x)$
\item The corrupted input $\tilde {x}$ is then mapped to a hidden representation with the same process of the standard autoencoder, $h=f(\tilde {x})$
\item From the hidden representation the model reconstructs the uncorrupted input $x=g(h)$.
\end{itemize}
Using the learnt hidden representation as a manifold, we can visualize the process of corruption and denoising (figure \ref{fig:denoisingmanifold} top).
\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/denoising_manifold}\\
\includegraphics[width=\textwidth]{images/autoencoder/denoising_network2}\\
\caption{RNA sequence denoising, from \cite{Eraslan2019scRNA}}
\label{fig:denoisingmanifold}
\end{figure}

Additionally, we can also explicitly change the reconstruction criterion in denoising autoencoders to achieve a good representation. For example, for RNA-sequence denoising, we can explicitly model the influence of dropout in the reconstruction process, as shown in figure \ref{fig:denoisingmanifold} bottom.

\subsubsection{Variational autoencoder*}
%TODO: this needs more careful introduction; also full prob formulation as before; ref to generative model and discuss this. talk about inference, ELBO in more detail
The basic idea behind a \defn{variational autoencoder (VAE)} is that instead of mapping an input to fixed vector, input is mapped to a distribution, see figure \ref{fig:vae}. The only difference to a normal autoencoder is that bottleneck vectors are now learning distributions, e.g. in the common normal case they are now learning two different vectors: one representing the mean of the distribution and the other representing the standard deviation of the distribution.

The loss function for a variational autoencoder is defined as:
\begin{equation}
\mathcal{L}(\phi,\theta)=-\mathbb{E}_{h{\sim}q_{\phi}(h|x)}[{\log}p_\theta(x|h)]+\mathbb{KL}(q_{\phi}(h|x)\|p_\theta(h))
\end{equation}
The first term represents the reconstruction error and the second term is a regularizer where KL stands for Kullback-Leibler divergence between the encoder's distribution $q_{\phi}(h|x)$ and true distribution of latent variables $p_\theta(h)$. This divergence measures how much information is lost when using $q$ to represent $p$.

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/normal_autoencoder}\\
\includegraphics[width=\textwidth]{images/autoencoder/variational_autoencoder}\\
\caption{Normal autoencoder vs. variational autoencoder. From \url{http://kvfrans.com/variational-autoencoders-explained/}.}
\label{fig:vae}
\end{figure}

\subsubsection{The probability perspective of a variational autoencoder*}
In the probability model framework, a variational autoencoder contains a specific probability model of data $x$ and latent variables $h$.  The goal of an variational autoencoder is to maximize the log likelihood of the data in our model:
$$\max_\theta\log{p_\theta(x)}$$
where $\theta$ is the model parameter, which is the network weights in a variational autoencoder. We can compute the above likelihood by marginalizing out the latent variables:
$$p_\theta(x)={\int}p_\theta(x,h)dh={\int}p_\theta(x|h)p(h)dh.$$
Unfortunately, this integral requires exponential time to compute as it needs to be evaluated over all configurations of latent variables. We therefore choose to approximate it lower bound, constructed as:
\begin{eqnarray*}
{\log}p_\theta(x)&\leq&{\log}p_\theta(x)-\mathbb{KL}(q_{\phi}(h|x)\|p_\theta(h|x)) \\
&=&\mathbb{E}_{q_{\phi}(h|x)}\left[{\log}p_\theta(x,h)-{\log}q_\phi(h|x)\right] \\
&=&\text{ELBO}(\phi,\theta).
\end{eqnarray*}
The above $\leq$ holds as any $\mathbb{KL}(\cdot\|\cdot)\geq0$. This lower bound is usually called \defn{Evidence Lower Bound (ELBO)}.

\begin{exercise}
Prove $\text{ELBO}(\phi,\theta)=-\mathcal{L}(\phi,\theta)$
\end{exercise}

\begin{eqnarray*}
\text{ELBO}(\phi,\theta)&=&\mathbb{E}_{q_{\phi}(h|x)}\left[{\log}p_\theta(x,h)-{\log}q_\phi(h|x)\right] \\
&=&\mathbb{E}_{q_{\phi}(h|x)}\left[{\log}p_\theta(x|h)+{\log}p_\theta(h)-{\log}q_\phi(h|x)\right] \\
&=&\mathbb{E}_{q_{\phi}(h|x)}{\log}p_\theta(x|h)-\mathbb{KL}(q_{\phi}(h|x)\|p_\theta(h)) \\
&=&-\mathcal{L}(\phi,\theta)
\end{eqnarray*}
The approach taken by variational autoencoder is to introduce an inference model $q_\phi(h|x)$ that learns to approximate the intractable true distribution of hidden parameters $p_\theta(h)$. The variational posterior $q_\phi(h|x)$ is parameterized by a encoder neural network, $f(x)$, that accepts input $x$, and outputs the mean and variance of a Normal distribution, shown in figure \ref{fig:vae}:
$$\mu_h(x,\phi),\log \sigma_h(x,\phi)=f(x)$$
$$q_\phi(h|x)=\mathcal{N}(h|\mu_h(x,\phi),\sigma_h^2(x,\phi))$$
Furthermore, the authors of \cite{Diederik2014vae} proposed a \defn{reparameterization trick}, i.e., they reparameterize the
hidden parameter $h{\sim}q_\phi(h|x)=\mathcal{N}(h|\mu_h(x,\phi),\sigma_h^2(x,\phi))$, as:
$$h=h(\varepsilon, x, \phi)=\varepsilon\sigma_h(x,\phi)+\mu_h(x,\phi),\quad\varepsilon\in\mathcal{N}(0,I)$$
As shown in Fig. \ref{fig:reparameterization}, \defn{reparameterization trick} allows gradient to flow through the sampling node.

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/vae_reparameterization}
\caption{A training-time variational autoencoder implemented as a feed-forward neural network, where $p(x|z)$ is Gaussian. The network on the left is without the \defn{reparameterization trick}, and network on the right is with it. Image from \url{https://arxiv.org/abs/1606.05908}.}
\label{fig:reparameterization}
\end{figure}

\textbf{Generation of images with VAE}
Variational autoencoder can be used by simulating dataset by passing samples on the latent representation through the trained decoder, as shown in figure \ref{fig:vae_imgen}.
\begin{figure}\centering
\includegraphics[width=0.6\textwidth]{images/autoencoder/vae_imgen}
\caption{Synthetic digits generated by an variational autoencoder trained with MNIST dataset.}
\label{fig:vae_imgen}
\end{figure}


%In this way, the objective function can be rewritten with respect to $\varepsilon$:
%\begin{eqnarray*}
%\mathcal{L}(\phi,\theta)&=&\mathbb{E}_{h\sim{q_{\phi}(h|x)}}\left[{\log}p_\theta(x,h)-{\log}q_\phi(h|x)\right]\\
%&=&\mathbb{E}_{\varepsilon\sim{\mathcal{N}(0,I)}}\left[{\log}p_\theta(x,h(\varepsilon,x,\phi))-{\log}q_\phi(h(\varepsilon,x,\phi)|x)\right]
%\end{eqnarray*}
%Thus the gradients with variational parameters $\phi$ can be directly moved into the expectation, enabling an unbiased low-variance Monte Carlo estimator:
%$$\nabla_\phi\mathcal{L}(\phi,\theta)\approx\frac{1}{k}\sum_{i=1}^{k}\nabla_\phi\left[{\log}p_\theta(x,h(\varepsilon_i,x,\phi))-{\log}q_\phi(h(\varepsilon_i,x,\phi)|x)\right] $$

\subsection{Generative adversarial networks*}
Although VAE can synthesize images, the synthetic images using VAE are usually blurred and not sufficiently realistic, as shown in Figure \ref{fig:vae_recon}. In 2014, a breakthrough paper introduced Generative adversarial networks (GANs, \cite{Goodfellow2014gan}) to generate fake photorealistic images, based on an ingenious way to leverage the power of discriminative models to get good generative models. The key idea behind GANs is that a data generator is good if we cannot tell fake data apart from real data with a discriminator, as illustrated in figure \ref{fig:gan}.

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/vae_recon}
\caption{Random samples from training set compared to their VAE reconstruction. From \url{https://towardsdatascience.com/gans-vs-autoencoders-comparison-of-deep-generative-models-985cf15936ea}}
\label{fig:vae_recon}
\end{figure}

\begin{figure}\centering
\includegraphics[width=0.5\textwidth]{images/autoencoder/gan}
\caption{Architecture of a Generative Adversarial Network.}
\label{fig:gan}
\end{figure}

More specifically, a GAN consists of two networks:
\begin{itemize}
\item a generator network $G(.)$ that takes a random input $z$ and generates an "fake" output $x_g=G(z)$ that should follow the "true" probability distribution (whether it matches the "true" distribution or not is judged by a discriminative network)
\item a discriminative network $D(.)$ that takes two inputs: one input from the "true" data $x$ and the other input from the "fake" one from the generator $G(z)$ and returns the probability of both inputs to be a "true" data.
\end{itemize}
Both networks are in competition with each other in such a way that the generator network attempts to fool the discriminator network and in the mean time, the discriminator network adapts to the new fake data, written as pesudo-code:

\begin{algorithm}
\For{\text{number of training iterations}}{
Sample minibatch of $m$ noise samples $\{z^{(1)},...,z^{(m)}\}$ from noise prior $p_g(z)$; \\
Sample minibatch of $m$ samples $\{x^{(1)},...,x^{(m)}\}$ from data generating distribution $p_\text{data}(x)$;\\
Update the discriminator by ascending its stochastic gradient:
$$\nabla_{\theta_d}\frac{1}{m}\sum_{i=1}^m\left[\log D(x^{(i)})+\log (1-D(G(z^{(i)}))    \right]$$ \\
Sample minibatch of $m$ noise samples $\{z^{(1)},...,z^{(m)}\}$ from noise prior $p_g(z)$; \\
Update the generator by descending its stochastic gradient:
$$\nabla_{\theta_g}\frac{1}{m}\sum_{i=1}^m\log (1-D(G(z^{(i)}))) $$ \\
}
\caption{Minibatch stochastic gradient descent training of generative adversarial network. From \cite{Goodfellow2014gan}.}
\end{algorithm}

\textbf{Exemplary GAN applications}
GAN is a very popular topic in deep learning and here we list some common applications:
\begin{itemize}
\item \textbf{Data Augmentation} - Aiming to increase training data size by synthesis.
\item \textbf{Style Transfer and Manipulation} - The input to the generator is no-longer a noise sample $z$, but an image, and the target is the same image of the transferred style in another domain. This can be trained either using the paired training data or un-paired ones (cycle-GAN). Some examples are shown in Figure \ref{fig:gan_styletransfer}.
\item \textbf{Signal Super Resolution} - Artificially increasing the resolution of images. Methodology-wise similar to Style Transfer.
\end{itemize}

\begin{figure}\centering
\includegraphics[width=\textwidth]{images/autoencoder/gan_styletransfer}
\caption{Style transfer using Cycle-GAN, from \cite{Zhu2017cyclegan}.}
\label{fig:gan_styletransfer}
\end{figure}




%\bibliographystyle{plain}\bibliography{paperslibrary}
%\end{document}

\subsection{Independent component analysis}

For a given random vector independent component analysis (ICA) tries to find its statistically independent components, usually after an unknown linear transformation. This idea can be used also to solve the blind source separation (BSS) problem which is, given only the mixtures of some underlying independent source signals, to separate the mixed signals thus recovering the original sources. Figure \ref{im:icaimages} shows how to apply ICA to separate two mixed natural images.
%Here neither the sources nor the mixing process is known, hence the term \emph{blind} source separation.
In contrast to correlation-based transformations such as PCA, ICA renders the output signals as statistically independent as possible by evaluating higher-order statistics, see figure \ref{fig:PCA_vs_ICA}.

%The idea of ICA was first expressed by Jutten and Herault \cite{HeJu86} \cite{JuHeCoSo91} while the term 'ICA' was later coined by Comon in \cite{comon}. However the field became popular only with the seminal paper by Bell and Sejnowski \cite{bell1995} who elaborated upon the Infomax-principle, which was first advocated by Linsker \cite{linsker1989} \cite{linsker1992}.  Cardoso and Laheld \cite{CaLa96} as well as  Amari \cite{AmariNatGrad98} later simplified the Infomax learning rule introducing the concept of a natural gradient which accounts for the non-Euclidean Riemannian structure of the space of weight matrices.

Often ICA is performed after pre-processing with PCA, which implies that instead of searching for any linear transformation only orthogonal transformations have to be identified. Many ICA algorithms have been proposed with the FastICA-algorithm \cite{HyvaFastICA99} being the one of the most efficient and popular ones.

\begin{figure}
\centering
\ifthenelse{\isundefined{\colorimages}}{
\includegraphics[width=0.5\columnwidth]{images/icabasics/michafabisource}
\includegraphics[width=0.5\columnwidth]{images/icabasics/michafabimixed}
\includegraphics[width=\columnwidth]{images/icabasics/michafabirecovered}}{
\includegraphics[width=0.5\columnwidth]{images/icabasics/michafabisource_color}
\includegraphics[width=0.5\columnwidth]{images/icabasics/michafabimixed_color}
\includegraphics[width=\columnwidth]{images/icabasics/michafabirecovered_color}}
\caption[Use of ICA for performing BSS]{
Use of ICA for performing BSS. The top row shows the source images. These were linearly mixed by $\matA=\left(\begin{array}{cc} 0.95 & 0.61\\ 0.23 & 0.49\end{array}\right)$. The second line shows those two mixture images. A self-adaptive version of the natural gradient ICA algorithm, section \ref{sec:MLalgos}, is used to separate these two images. When comparing the recovery matrix with the mixture matrix, we get a low crosstalking error of $0.053$. The resulting images are plotted in the third line; here each recovered source was plotted together with its multiplication by $-1$. This is done for visual purposes only; in general the sign cannot be recovered, it is part of the ICA indeterminacies. %, see section \ref{sec:ICAprops}.
\label{im:icaimages}
}
\end{figure}


\begin{figure}\centering
\includegraphics[width=.7\textwidth]{images/icabasics/PCA_vs_ICA}
\caption{PCA versus ICA: Both applied to linearly transformed 2d unit cube. PCA (red) identifies directions of strongest/lowest variance (orthogonal), whereas ICA (green) identifies true mixing directions (parallel to cube edges).}
\label{fig:PCA_vs_ICA}
\end{figure}

\subsubsection{Linear BSS and ICA}
Consider the noiseless linear instantaneous blind source separation (BSS) model with as many sources as observation ($p=n$, otherwise dimension reduction by PCA is common):
\begin{equation}
\label{eq:BSSmodel}
\vecX=\matA\vecS
\end{equation}
with an independent $(p=n)$-dimensional random vector $\vecS$ and invertible matrix $\matA$. (Sorry for the notation clash, this was $Z$ and $\W$ before).
The task of linear BSS is to find $\matA$ and $\vecS$ given only $\vecX$.


\begin{exercise}
When does a solution exist? And it is unique?
\end{exercise}

\subsubsection{Linear separability}
An obvious indeterminacy of this problem is that $\matA$ can be found only up to scaling and permutation because for scaling $\matL$ and permutation matrix $\matP$
$$\vecX=\matA\matL\matP \matP^{-1}\matL^{-1}\vecS$$
and $\matP^{-1}\matL^{-1}\vecS$ is also independent.
Here, an invertible matrix $\matL$ is said to be a \defn{scaling matrix}, if it is diagonal.
We will show that under mild assumptions to $\vecS$ there are no further indeterminacies of linear BSS:

\begin{theorem}[Separability of linear BSS]
\label{thm:uniqlin}
Let $\matA$ be an invertible matrix and $\vecS$ be an independent random vector.
Assume that $\vecS$ has at most one Gaussian or deterministic component and that the covariance of $\vecS$ exists.
Then if $\matA\vecS$ is again independent, $\matA$ is the product of a scaling and a permutation matrix.
\end{theorem}

This theorem has been used to show separability by Comon\cite{comon}
%and extended by Erikkson and Koivunen \cite{ErKo03}
based on the Darmois-Skitovitch theorem\cite{darmois,skitovitch53}.
We will give a much easier proof without having to use the Darmois-Skitovitch theorem.

The theorem indeed shows separability, because if $\vecX=\matA\vecS=\matA'\vecS'$ are two BSS solutions of $\vecX$, then $\vecS$ and $\vecS'=\matA'^{-1}\matA\vecS$ are independent. Hence according to theorem \ref{thm:uniqlin} $\matA'^{-1}\matA$ is the identity except for permutation and scaling indeterminacy, which implies that $\matA'$ and $\matA$ differ only by right multiplication of a scaling and a permutation matrix.

\subsubsection{2d positive $C^2$-density proof}

For illustrative purpose we will prove separability for a two-dimensional random vector $\vecS$ with positive twice continuously differentiable density $p_\vecS$. For the full proof, we refer to \cite{Theis:2004jea}.
Let $\matA$ be an invertible $2\times2$ matrix. It is enough to show that if $\vecS$ and $\matA \vecS$ are independent then either $\matA$ is the product of a scaling and a permutation matrix or $\vecS$ is Gaussian.

$\vecS$ is assumed to be independent, so its density factorizes:
$p_{\vecS}(\vecs) = g_1(s_1) g_2(s_2),$
for $\vecs\in\R^2$.
First note that the density of $\matA\vecS$ is given by
\begin{align*}
p_{\matA\vecS}(\vecx)&=|\det \matA|^{-1} p_{\vecS}(\matA^{-1}\vecx)\\
&=c g_1(b_{11}x_1+b_{12}x_2)g_2(b_{21}x_1+b_{22}x_2)
\end{align*}
for $\vecx\in\R^2$, $c\neq 0$ fixed. Here $\matB = (b_{ij}) = \matA^{-1}$.
$p_\vecS$ was assumed to be positive, then so is $p_{\matA\vecS}$. Using independence of $\matA\vecS$, $\ln p_{\matA\vecS}(\vecx)$ is a sum of functions depending only on either $x_1$ or $x_2$, so
\begin{align*}
&\frac{\partial^2}{\partial x_1 \partial x_2} \ln p_{\matA\vecS}(\vecx)\\
&=b_{11}b_{12}h_1''(b_{11}x_1+b_{12}x_2) + b_{21}b_{22} h_2''(b_{21}x_1+b_{22}x_2)
=0
\end{align*}
for all $\vecx\in\R^2$, where $h_i:=\ln g_i$. By setting $\vecy := \matB\vecx$, we therefore have
\begin{equation}
\label{eq:twodseq}
b_{11}b_{12}h_1''(y_1) + b_{21}b_{22} h_2''(y_2)=0
\end{equation}
for all $\vecy\in\R^2$, because $\matB$ is invertible.

Now, if $\matA$ (and therefore also $\matB$) is the product of a scaling and a permutation matrix, then equation \ref{eq:twodseq} holds. If not, then $\matA$ and hence also $\matB$ have at least three nonzero entries. By equation \ref{eq:twodseq} also the fourth entry has to be nonzero, because the $h_i''$ are not zero (otherwise $g_i(y_i)=\exp(ay_i+b)$, which is not integrable). Furthermore,
$$b_{11}b_{12}h_1''(y_1) = -b_{21}b_{22} h_2''(y_2)$$
for all $\vecy\in\R^2$, so the $h_i''$ are constant, say $h_i''\equiv c_i$, and $c_i\neq 0$ as noted above.
Therefore the $h_i$ are polynomials of degree $2$, and the $g_i = \exp h_i$ are Gaussians ($c_i<0$ because of the integrability of the $g_i$) as was to be shown.



\subsubsection{ICA using maximum-likelihood estimation}
\label{sec:MLica}
There exist many different algorithms for solving linear ICA. Here we present an early one based on MLE.

Consider the noiseless square ($n=p$) ICA model
$$\vecX=\matA\vecS.$$
With $\matB=\matA^{-1}$ as before and the transformational properties of densities, we can write
$$p_\vecX(\matA\vecs) = |\det \matB|p_\vecS(\vecs)= |\det \matB|\prod_{i=1}^n p_i(\vecs)$$
with $p_i:=p_{S_i}$ the source component densities.

%Setting $\vecx:=\matA\vecs$ yields $\vecs=\matB\vecx$. If we denote with $\vecb_i^\top$ the rows of $\matB$ that is
%$$\matB=(\vecb_1 | \ldots | \vecb_n)^\top$$
%then $s_i=\vecb_i^\top\vecx$ and therefore
Let $\vecb_i$ denote the rows of $\B$. Then $$p_\vecX(\vecx)=|\det\matB|\prod_{i=1}^n p_i(\vecb_i^\top\vecx).$$
and the likelihood reads
$$p(\mathcal{D}|\theta)= \prod_{j=1}^N p_{\vecX|\matB}(\vecx_j|\matB)
= \prod_{j=1}^N |\det\matB| \prod_{i=1}^n p_i(\vecb_i^\top\vecx_j)
$$
And log-likelihood, with sample mean $E$,
$$\frac{1}{N}\ln p(\mathcal{D}|\theta) = E\left(\sum_{i=1}^n \ln p_i(\vecb_i^\top\vecx(t))\right) + \ln|\det\matB|.$$

The main problem we are facing now is that in addition to the parametric model --- the estimation of $\matB$ --- also the unknown source densities have to be estimated, which cannot be directly described by a finite set of parameters. So we are dealing with so-called \defn{semiparametric estimation}.

If we still want to use maximum likelihood estimation in order to find $\matB$, two different solutions can be taken depending on prior information:
\begin{itemize}
\item
Either due to prior information, the source densities $p_i$ are known. Then the likelihood of the whole model is indeed described only by $p(\mathcal{D}|\matB)$ because $\matB$ is the only unknown parameter.
\item
If no additional information is given, the source densities $p_i$ will have to be approximated using some sort of parameterized density families.
\end{itemize}

Indeed, the second route can be taken without too much hassle as is shown by the following theorem. It claims that for ICA estimation it is enough to locally describe each $p_i$ by a simple binary density family meaning a family with only two elements --- this is quite astonishing as the space of density families is obviously very large.

\begin{theorem}
\label{thm:loccons}
Let $\tilde p_i$ be the estimated IC densities, assume $\tilde p_i>0$. Let
$$g_i(s):= \frac{d}{ds} \ln \tilde p_i(s)=\frac{\tilde p_i'}{\tilde p_i} (s)$$
be the (negative) \defn{score functions} and let the $y_i:=\vecb_i^\top\vecx$ be whitened. Then the maximum likelihood estimator is locally consistent if
\begin{equation}
\label{eqn:locconscondition}
E(S_ig_i(S_i)-g_i'(S_i))>0
\end{equation}
for $i=1,\ldots,n$.
\end{theorem}

Here \defn{locally consistent} means that locally the estimated matrix $\tilde\matB$ converges to $\matB$ in probability for $T\rightarrow\infty$.

For a proof of this theorem, see for example \cite{HKO2001}, theorem 9.1.

Note that equation \ref{eqn:locconscondition} is invariant under small perturbations to the estimated densities $\tilde p_i$, because this equation only depends on the sign of $sg_i(s)-g_i'(s)$, so the local consistency of the maximum likelihood estimator is stable under small perturbations.

After choice of a binary density family, the likelihood can be minimized by standard gradient descent (or more elaborate methods) with adaptive choice of super or sub-gaussian approximate source. 
% For more details we refer to \cite{HyvarinenKarhunenOja2001} and \cite{Theis_Thesis02,Theis_Thesis03}.  % not known




%This idea enables us to use a simple binary density family. Define densities
%\begin{eqnarray}
%\label{eqn:pplus}
%\tilde p^+(s)&:=&\frac{c_+}{\cosh^2(s)} \\
%\label{eqn:pminus}
%\tilde p^-(s)&:=&\frac{c_-\cosh^2(s)}{\exp\left(s^2/2\right)}
%\end{eqnarray}
%with constant $c_\pm$ such that $\int \tilde p^\pm =1$. Calculation shows that $c_+ = 0.5$ and $c_-\approx 0.0951$.
%
%Taking logarithms we note that
%\begin{eqnarray*}
%\ln \tilde p^+(s)&=&\ln c_+ - 2\ln\cosh(s) \\
%\ln \tilde p^-(s)&=&\ln c_- - \left(\frac{s^2}{2} - \ln \cosh^2(s)\right)
%\end{eqnarray*}
%so $\tilde p^+$ is supergaussian and $\tilde p^-$ is subgaussian. This can also be seen in figure \ref{im:ML_tildepm}.
%
%\begin{figure}
%\centering
%\includegraphics[width=0.49\columnwidth]{images/icabasics/ML_tildepplus}
%\includegraphics[width=0.49\columnwidth]{images/icabasics/ML_tildepminus}
%\caption[A binary density family]{
%A binary density family.
%The left density is given by $\tilde p^+(s):=0.5\cosh^{-2}(s)$ (equation \ref{eqn:pplus}) and the right one by $\tilde p^-(s):=0.0951\cosh^2(s)\exp\left(-s^2/2\right)$ (equation \ref{eqn:pplus}).
%\label{im:ML_tildepm}
%}
%\end{figure}
%
%The score functions $g^\pm$ of these two densities are easily calculated as
%$$g^+(s)=(-2\ln\cosh s)'=-2\tanh s$$
%for $\tilde p^+$ and
%$$g^-(s)=(-\frac{s^2}{2}+\ln\cosh s)'=-s+\tanh s$$
%for $\tilde p^-$.
%Putting the score functions into equation \ref{eqn:locconscondition} then yields condition
%$$E(-S_i\tanh S_i + (1-\tanh S_i)^2) >0$$
%respectively
%$$E(S_i\tanh S_i - (1-\tanh S_i)^2) >0$$
%(because $E(S_i^2)=1$)
%for local consistency of the maximum likelihood estimator.
%
%If we assume that the source components fulfill $E(S_i\tanh S_i - (1-\tanh S_i)^2)\neq$ (similar to the assumption $\kurt(S_i)\neq 0$ in the kurtosis maximization algorithms), we have shown that either $\tilde p^+_i$ or $\tilde p^-_i$ fulfills equation \ref{eqn:locconscondition}. So, in order to guarantee local consistency of the estimator, for choosing the density of each source component, we simply have to choose the correct $\tilde p^+_i$. Then theorem \ref{thm:loccons} guarantees that the maximum likelihood estimator with this approximated source density still gives the correct unmixing matrix $\matB$ (as long as the mixtures have been whitened).
%
%Note that if we put $g(s)=-s^3$ into equation \ref{eqn:locconscondition}, we get the condition $\kurt(S_i)>0$ for local consistency. So in some sense, the choice of $\tilde p^\pm_i$ corresponds to whether we minimize or maximize kurtosis as we did in section \ref{sec:fastica}.
%
%\subsection{Algorithms}
%\label{sec:MLalgos}
%
%\subsubsection{Euclidean and natural gradient}
%
%In the next section, we want to maximize the likelihood from above using gradient ascent.
%For this we have to calculate the gradient of a function defined on a manifold of matrices. The gradient of a function is defined as the dual of the differential of the function with respect to the scalar product. As the standard scalar product on $\R^n$ is $\vecx^\top\vecy$, the ordinary gradient is simply the transpose of the derivative of the function. Here, we are interested in the gradient of a function defined on the open submanifold $\Gl(n)$ of $\R^{n^2}$. On $\Gl(n)$ we can either use the standard (Euclidean) scalar product (standard Riemannian metric) to get the \defn{Euclidean gradient}
%$$\grad^{\eukl} f(\matW) := \grad f(\matW) := (\matD f(\matW))^\top$$
%or we take a metric that is invariant under the group structure (multiplication) of $\Gl(n)$ to get the \defn{natural gradient}
%$$\grad^{\nat} f(\matW) := (\grad^{\eukl} f(\matW)) \matW^\top \matW.$$
%More details are given for example in \cite{TheisDiplom_00}, chapter 2.
%
%We also write for the Euclidean gradient
%$$\frac{\partial}{\partial \matW} f(\matW):= \grad^{\eukl} f(\matW).$$
%
%\begin{lemma}
%\label{lem:LnDet}
%$$\frac{\partial}{\partial \matW} \ln \det \matW = \matW^{-\top}$$
%for $\matW \in \Gl(n)$.
%\end{lemma}
%
%\begin{proof}
%We have to show that
%$$\frac{\partial}{\partial w_{ij}} \ln \det \matW = (\matW^{-1})_{ji}$$
%holds for $i,j=1,\ldots,n$. Using the chain rule we get
%$$\frac{\partial}{\partial w_{ij}} \ln \det \matW = \frac{1}{\det \matW} \frac{\partial}{\partial w_{ij}} \det \matW.$$
%According to the Cramer rule for the inverse, we have
%$$(\matW^{-1})_{ji} = (-1)^{i+j} \frac{1}{\det \matW} \det \matW^{(ij)},$$
%where $\matW^{(ij)} \in \Mat((n-1) \times (n-1); \R)$ denotes the matrix which comes from $\matW$ by leaving out the $i$-th row and the $j$-th column. The proof is finished if we show
%$$\frac{\partial}{\partial w_{ij}} \det \matW  = (-1)^{i+j} \det \matW^{(ij)},$$
%For this develop $\det \matW$ by the $i$-th row to get
%$$\det \matW = \sum_{k=1}^n (-1)^{i+k} w_{ik} \det \matW^{(ik)}.$$
%Then, taking derivative by $w_{ij}$ shows the claim.
%\end{proof}
%
%\begin{lemma}
%\label{lem:GPhi}
%For $\matW \in \Mat(n \times n;\R)$ and $p_i \in \C^{\infty}(\R, \R)$, $i=1,\ldots,k$
%$$\frac{\partial}{\partial \matW} \sum_{i=1}^n \ln p'_i( \matW\vecx)_i = \vecg(\matW\vecx) \vecx^\top,$$
%for $\vecx\in\R^n$, where for  $\vecy \in \R^n$
%$$\vecg(y) := \left( \frac{p''_i(y_i)}{p'_i(y_i)} \right)_{i=1}^n \in \R^n.$$
%\end{lemma}
%
%\begin{proof}
%We have to show that
%$$\frac{\partial}{\partial w_{ij}} \sum_{k=1}^n \ln p'_k( \matW\vecx)_k = \frac{p''_i(y_i)}{p'_i(y_i)} x_j$$
%But this follows directly from the chain rule.
%\end{proof}
%
%
%
%\subsubsection{Bell-Sejnowski algorithm}
%
%With the following algorithm, Bell and Sejnowski gave one of the first easily applicable ICA algorithms \cite{bell1995}. It maximizes the likelihood from above by using gradient ascent.
%
%The goal is to maximize the likelihood (or equivalently the log likelihood) of the parametric ICA model. If we assume that the source densities are differentiable, we can locally do this using gradient ascent. The Euclidean gradient of the log likelihood can be calculated using lemmata \ref{lem:LnDet} and \ref{lem:GPhi} to be
%$$\frac{1}{T} \frac{\partial \ln L(\matB)}{\partial \matB} = \matB^{-\top}+\vecE(\vecg(\matB\vecX)\vecX^\top)$$
%with the $n$-dimensional score function $\vecg=g_1\times\ldots\times g_n$. So the local update algorithm goes as follows:
%
%\emph{Algorithm:}(\defn{gradient ascent maximum likelihood}) Choose $\eta>0$ and $\matB(0)\in \Gl(n)$. Then iterate for whitened mixtures $\vecX$
%\begin{eqnarray*}
%\Delta\matB(t) &:=& \matB(t)^{-\top}+\vecE(\vecg(\matB(t)\vecX)\vecX^\top)\\
%\matB(t+1) &:=& \matB(t) + \eta\Delta\matB(t)
%\end{eqnarray*}
%
%Instead of using this batch update, we can use a stochastical version by substituting expectation by samples to get
%$$\Delta\matB(t) := \matB(t)^{-\top}+\vecg(\matB(t)\vecx(t))\vecx(t)^\top)$$
%for a sample $\vecx(t)\in\R^n$.
%
%This algorithm was quite a revolution in its early days, however it faces problems such as convergence speed and the numerically problematic matrix inversion in each update step.
%
%\subsubsection{Natural gradient algorithm}
%
%These problems were mostly fixed by Amari \cite{AmariNatGrad98}, who used the natural instead of the Euclidean gradient:
%$$\frac{1}{T} \grad^{\nat} L(\matB) = \frac{1}{T} \left(\grad^{\eukl} L(\matB)\right) \matB^\top\matB = (\unitm +\vecE(\vecg(\vecY)\vecY^\top))\matB$$
%with $\vecY:=\matB\vecX$. Using
%$$\Delta\matB(t):=(\unitm +\vecE(\vecg(\vecY)\vecY^\top))\matB$$
%gives both better convergence and numerical stability, as simulations confirm.
%
%
%\subsubsection{Score functions}
%
%Still, it is not clear which score functions are to be used.
%As we saw before, the score functions of the binary density family $\tilde p^\pm$ are
%\begin{eqnarray*}
%g^+(s)&=& -2\tanh s\\
%g^-(s)&=& \tanh s - s.
%\end{eqnarray*}
%For the above two algorithms, the componentwise nonlinearities $g_i$ are then chosen online according to equation \ref{eqn:locconscondition}: If
%$$E(-S_i \tanh S_i + (1-\tanh^2 S_i))>0$$
%then we use $g^+$ for the $i$-th component, if not $g^-$. As said before, this is done online after prewhitening.
%

\section{Graphical Models \& Causality}
\subsection{Motivation}

So far, our starting point has always been a \emph{joint distribution} $P({\X})$ over a finite set of variables $\X = \{X_1, \ldots, X_p\}$, $p \in \N$.

\paragraph{Examples}
\begin{itemize}
  \item Supervised: $X_1, \ldots, X_{p-1}$ are real-valued features ($p$ in the tens or hundreds for tabular data, $p$ in the thousands or millions for images, etc.) and $Y:= X_p$ is a categorical (classification) or real-valued (regression) target variable and we are interested in the conditional $P(Y \, | \, X_1, \ldots, X_{p-1})$.
  Most of the time, we have deterministic predictions, i.e., we build models for $\bE[Y \, |\, X_1, \ldots, X_{p-1}]$.
  \item Unsupervised: We want to identify modes, or disconnected high-density regions of $P({\X})$, or write it as a mixture of simple densities; we want to write the joint distribution as a factorized one.
\end{itemize}

We could do all this only assuming that we have i.i.d.\ samples from a joint $P(\X)$, which amounts to purely \emph{observational} data.
In particular, we didn't mind if we were using \emph{how} variables relate to each other.
However, in most settings data is generated in some structured fashion.
There are natural laws, \emph{mechanisms or processes} governing how systems evolve and thus what we get to observe.
The resulting joint distribution ``hides'' these processes and thereby removes some information that may be relevant for our inferences.
For example, consider measurements of smoking $X$ and lung cancer $Y$.
As always, we assume that these two variables follow some joint distribution $X, Y \sim P$.
In real data, we find a strong correlation between the two, for example by fitting a model to $P(Y\, | \, X)$ from a finite sample.
However, the joint distribution $P$ cannot provide any insights into \emph{why} these two variables are correlated.
Is it because ``smoking causes cancer'' $X \to Y$?
It appears implausible that a cancer directly causes someone to smoke $Y \to X$ (years before even getting the cancer).
But perhaps there are other factors $Z$ (socio-economic status, or genetic mutations), that make people more susceptible to nicotine/smoking and simultaneously increase the risk of getting lung cancer $X \gets Z \to Y$?
All three scenarios could give rise to the same observed joint $P(X, Y)$, thus no method depending only on i.i.d.\ samples from the joint could distinguish between them.
How can we ever come to know that ``smoking causes lung cancer'' without measuring \emph{everything}?
Assuming we do know that smoking increases the risk of lung cancer (and not vice versa), can we quantify how much of the observed association is due to direct causation?
Mostly, we do not know the processes and mechanisms involved (otherwise we may not even need statistical learning), but perhaps we can improve our inference by making \emph{some assumptions about them}.

\paragraph{Takeaway}
Sometimes it appears useful to \emph{look beyond the joint distribution} and its associations by assuming/imposing/inferring more structure about it.
One useful way to talk about and represent such structure, is via \textbf{graphical models} and \textbf{structural equation models}.
Causality (causal inference) is a large field with many subfields. A broad classification by task is
\begin{itemize}
    \item \textbf{causal discovery:} given only $P$ infer how variables causally relate to each other. Who listens to whom?
    \item \textbf{cause-effect estimation:} given $P$ and a causal structure, quantify the strength of causal influences.
\end{itemize}

\paragraph{Further reading}
These two lectures are only going to be a small teaser of what constitutes full lecture series. For further reading, please see:
\begin{itemize}
    \item \href{https://djsaunde.github.io/read/books/pdfs/probabilistic%20graphical%20models.pdf}{Probabilistic Graphical Models} by Daphne Koller and Nir Friedmann
    \item \href{http://bayes.cs.ucla.edu/BOOK-2K/}{Causality} by Judea Pearl
    \item \href{https://mitpress.mit.edu/books/elements-causal-inference}{Elements of Causal Inference} by Jonas Peters, Dominik Janzing, and Bernhard Sch\"olkopf (click on ``Open Access'' to download the book for free)
    \item \href{https://www.cs.cmu.edu/afs/cs.cmu.edu/project/learn-43/lib/photoz/.g/web/.g/scottd/fullbook.pdf}{Causation, Prediction and Search} by Peter Spirtes, Clark Glymour, and Richard Scheines
    \item \href{http://www.statslab.cam.ac.uk/~qz280/teaching/causal-2019/reading/Lauritzen_1996_Graphical_Models.pdf}{Graphical models} by Steffen Lauritzen
\end{itemize}

We took a lot of inspiration (with many parts copied essentially verbatim) from Jonas Peters' \href{http://web.math.ku.dk/~peters/jonas_files/scriptChapter1-4.pdf}{causality script}.


\subsection{Graphs}

Let $\V = \{1, \ldots, p\}$ and $\X = \{X_1, \ldots, X_p\}$ have joint distribution $P({\X})$ and density $p(\x)$.

\begin{definition}
  A \textbf{graph} $\cG = (\V, \cE)$ consists of \textbf{nodes} or \textbf{vertices} $\V$ and \textbf{edges} $\cE \subset \V \times \V$ with $(v, v) \notin \cE$ for any $v \in \V$.
\end{definition}


\begin{definition}
We introduce some standard graph terminology.
\begin{enumerate}
    \item A node $i$ is called a \textbf{parent} of $j$ if $(i, j) \in \cE$ and $(j,i) \notin \cE$ and a \textbf{child} if $(j, i) \in \cE$ and $(i,j) \notin \cE$.
    The set of parents of $j$ is $\paG{j}$ and the set of children is $\chG{j}$.
    \item Nodes $i, j$ are \textbf{adjacent} if $(i,j) \in \cE$ or $(j,i) \in \cE$.
    \item $\cG$ is \emph{fully connected} if all pairs of nodes are adjacent.
    \item There is an \textbf{undirected edge} between $i,j$ if $(i,j), (j,i) \in \cE$.
    An edge $(i,j) \in \cE$ is \textbf{directed} if it is not undirected.
    We then write $i \to j$ for $(i, j) \in \cE$.
    \item Three nodes are called \textbf{v-structure} if one is a child of the other two, which are not adjacent ($1 \to 3 \gets 2$).
    \item  The \textbf{skeleton} of $\cG$ is the graph $\cG$ with undirected edges between all adjacent nodes.
    \item A \textbf{path} in $\cG$ is a sequence of distinct vertices $i_1, \ldots, i_n$ ($n \ge 2$) such that $i_k, i_{k+1}$ are adjacent for $k \in \{1, \ldots, n-1\}$.
    A path is \textbf{directed} if $i_k \to i_{k+1}$ for all $k$ and $i_n$ is a \textbf{descendant} of $i_1$.
    A node $i$ is neither a descendant nor non-descendant of $i$.
    The set of all descendants of $j$ is $\deG{j}$ and non-descendants $\ndG{j}$.
    \item If $i_{k-1} \to i_k \gets i_{k+1}$ in a path, $i_k$ is a \textbf{collider relative to this path} ($1 \to 3 \gets 2$).
    \item $\cG$ is a \textbf{directed acyclic graph (DAG)} if all edges are directed and there is no directed cycle, i.e., no pair $i,j$ with directed paths from $i$ to $j$ and from $j$ to $i$.
    \item In a DAG, a path between $i_1$ and $i_n$ is \textbf{blocked by a set $\S$} with $i_1, i_n \notin \S$ if there is a node $i_k$ such that either
    \begin{enumerate}
        \item $i_k \in \S$ and $i_{k-1} \to i_k \to i_{k+1}$ or $i_{k-1} \gets i_k \gets i_{k+1}$ or $i_{k-1} \gets i_k \to i_{k+1}$
        \item $i_{k-1} \to i_k \gets i_{k+1}$ and neither $i_k$ nor any of its descendants is in $\S$.
    \end{enumerate}
    A path that is not blocked by a set $\S$ is \textbf{active}.
    \item Two subsets $\A, \B \subset \V$ are \textbf{d-separated} by $\S \subset \X$, if all paths between $i\in \A$ and $j\in \B$ are blocked by $\S$. We write $\A \perp \B | \S$.
\end{enumerate}
\end{definition}

D-separation encodes information about which variables can ``influence'' each other (not just association).

\inlineQ{Example!}

\subsection{Probabilistic Graphical Models}

\fcolorbox{black}{white}{How do graphs and distributions relate to each other: $P({\X}) \leftrightarrow \cG$?}


\paragraph{Factorizations.}
A joint distribution $P(\X)$ can be \textbf{factored}, for example:
\begin{equation}
  p(\X) = p(X_1) \, p(X_2 \, | \, X_1)\, p(X_3 \, | \, X_1, X_2) \, \ldots \, p(X_p\, |\, X_1, \ldots, X_{p-1}) \:.
\end{equation}
Here $p$ is the density of $P$ and we always assume it exists.
We can build a DAG $\cG$ by taking all factors $p(X_j \, | \, \X')$ and add edges $X_k \to X_j$ for all $X_k \in \X' \subset \X$.
\inlineQ{Why is this a DAG?}
We say that $p(\X)$ \textbf{factors according to the DAG} $\cG$.
There are $p!$ factorizations of $P(\X)$ according to a fully connected DAG, one for each permutation $\sigma \in S_n$ and DAG.
\inlineQ{Why?}

What about the remaining DAGs? \inlineQ{How does the number of DAGs grow as a function of the number of vertices?}
Let $\cG$  be a fully connected DAG and $p(X_j \, | \, \X')$ the factor for $X_j$ in the corresponding factorization.
A set of variables $\pa_j := \pa(X_j) \subset \X$ is the \textbf{Markovian parents of $X_j$} if it is a minimal subset of $\X'$ such that $p(X_j \, | \, \pa_j) = p(X_j \, |\, \X')$.
We then have the \textbf{Markov factorization property}
\begin{equation}
  p(\X) = \prod_{j=1}^p p(X_j \, | \, \pa_j) \:.
\end{equation}

\begin{example}
  Observing samples from a joint $P$, we may have reason to believe that there is a factorization according to some DAG with ``small factors'', i.e., $\pa_j$ contains only few variables for each $j$.
  If we then want to estimate the joint $P$, we only need to estimate the simpler conditionals with few inputs.
  In case of a chain $X_1 \to X_2 \to \ldots \to X_p$, we only need to estimate $p(X_{k+1} \, |\, X_{k})$ for $k \in \{1, \ldots, p-1\}$ and $p(X_1)$.
  For a general joint distribution over $p$ binary variables we need $2^p-1$ parameters, for a chain only $2p-1$. \inlineQ{Why?}
\end{example}
There is one potential factorization of $P(\X)$ for each DAG and a distribution can factorize according to multiple DAGs.

\begin{definition}
  A \textbf{Bayesian network (or DAG model)} $(\cG, P)$ consists of a DAG $\cG$ and a joint distribution $P$ that factorizes according to $\cG$.
\end{definition}


\fcolorbox{black}{white}{Bayesian networks give a crude connection $P(\X) \leftrightarrow \cG$. Can we do more?}

\begin{example}[First order Markov models]
If variables influence each other in a chain $X_1 \to X_2 \to \ldots \to X_p$, the ``future is independent of the past given the present'' $X_{t+1} \indep \{X_1, \ldots, X_t-1\} \, | \, X_t$, known as the Markov property.
\end{example}

We have similar properties for Bayesian networks.

\begin{definition}
Given $\cG$ and a joint distribution $P$, $P$ satisfies
\begin{enumerate}
    \item the \textbf{local Markov property} w.r.t.\ $\cG$ if
    \begin{equation*}
        X_j \indep \ndG{X_j} \,|\, \paG{X_j} \text{ for all } j \in \V.
    \end{equation*}
    \item the \textbf{global Markov property} w.r.t.\ $\cG$ if
    \begin{equation*}
        \A \perp \B \, |\, \C \implies \A \indep \B \, |\, \C.
    \end{equation*}
    for all disjoint sets $\A, \B, \C \subset \X$.
    \item the \textbf{Markov factorization property} w.r.t.\ $\cG$ if
    \begin{equation*}
      p(\X) = \prod_{j=1}^p p(X_j \, | \, \paG{X_j}) \:,
    \end{equation*}
\end{enumerate}
\end{definition}

\begin{theorem}
  All three Markov properties are equivalent.
\end{theorem}
For a proof, see \cite{lauritzen1996graphical}.

When we have a graph $\cG$ of ``how the data came about'', it tells us something about conditional independencies in the joint distribution $P$.

\begin{definition}[Markov equivalence]
We define
\begin{equation*}
    \cM(\cG) := \{P \, |\, P \text{ satisfies the global Markov property w.r.t.\ } \cG\}.
\end{equation*}
Two DAGs $\cG_1, \cG_2$ are \textbf{Markov equivalent} if $\cM(\cG_1) = \cM(\cG_2)$.
\end{definition}

\begin{theorem}
  Two DAGs are Markov equivalent if and only if they have the same skeleton and the same v-structures.
\end{theorem}

Before being able to apply this to make statements about to what extent causal discovery is possible, we need one more assumption.

\begin{definition}
  The joint $P$ is said to be \textbf{faithful to the DAG $\cG$} if
  \begin{equation*}
      \A \perp \B \, |\, \C \impliedby \A \indep \B \, |\, \C.
  \end{equation*}
  for all disjoint sets $\A, \B, \C \subset \X$. \inlineQ{Compare to the global Markov property.}
\end{definition}
Faithfulness excludes the possibility of ``path cancellations''.


\fcolorbox{black}{white}{\parbox{\textwidth}{If $P$ satisfies the global Markov property \emph{and} is faithful w.r.t.\ $\cG$, d-separations in the graph map one-to-one to conditional independencies.}}

This also allows us to make statements about whether a causal graph can be inferred from the joint distribution.
Even if we could test for conditional independencies perfectly in observed data and there are no unmeasured confounders, we can generally only identify a causal graph up to its equivalence class (without further assumptions).

\paragraph{Summary.}
We don't have time to go into details, but probabilistic models can be useful for
\begin{itemize}
  \item Estimating joint densities from easier conditionals.
  \item Reading off conditional independencies from graphs.
  \item Expert systems and probabilistic modeling (belief propagation, message passing).
  \item Next up: \textbf{Causal inference}
\end{itemize}


\section{Causal inference}

Causality critically relies on structured, modular mechanisms underlying the joint distribution $P$.
So far, the conditional factors did not necessarily encode directionality.
We would now like to think of a factor $p(X_i \,|\, X_j)$ as a natural mechanism or process that tells us how we can get $X_i$ from $X_j$, or how $X_j$ causes $X_i$.
Ideally, such a mechanism is ``invariant'', i.e., it doesn't depend on the marginal ``input'' distribution $p(X_j)$, e.g.,
\begin{equation}
    p(X, Y) = \underset{\text{invariant mechanism}}{\underbrace{p(Y \, | \, X)}} p(X).
\end{equation}
If this is the causal factorization, When we change $p(X)$, we still get the right joint.
We will now introduce a mathematical model to express such mechanisms and how they give rise to the joint $P$.

\subsection{Structural equation models}

\fcolorbox{black}{white}{\parbox{\textwidth}{But \textbf{what does $X_i \to X_j$ mean formally} in terms of $P$?}}

\begin{definition}
  A \textbf{structural equation model} (SEM) is a tuple $\cS = (\S, P_{\N})$, where $\S$ is a collection of $p$ equations
  \begin{equation}
      S_j \,:\qquad X_j = f_j(\paG{X_j}, N_j), \text{ for} j \in \{1, \ldots, p\},
  \end{equation}
  where $\paG{X_j}$ are called \textbf{parents} of $X_j$ (assuming $\cG$ is a DAG) and $P_{\N}$ is the joint distribution of the \textbf{noise variables} $N_j$, which are jointly independent.
\end{definition}

\begin{proposition}
  An SEM defines a unique joint distribution $P(\X)$ such that $X_j \overset{d}{=} f_j(\paG{X_j}, N_j)$ for all $j \in \{1, \ldots, p\}$
\end{proposition}
\begin{proof}
  This follows from the fact that $\cG$ is acyclic and the existence of a topological ordering of $\cG$.
\end{proof}

\begin{definition}
  A permutation $\pi \in S_p$ is called a topological ordering of a the variables in a DAG $\cG$, if
  \begin{equation}
    j \in \deG{i} \implies \pi(i) < \pi(j).
  \end{equation}
\end{definition}

\begin{proposition}
  For each DAG there exists a topological ordering.
\end{proposition}
\begin{proof}
  Exercise! (Look it up in Jonas Peters script.)
\end{proof}

\paragraph{Program viewpoint of SEMs}

Assume from now on that the nodes of $\cG$ are numbered according to a topological ordering.

We can view an SEM as a program for generating samples from the implied joint distribution $P(\X)$:

\begin{align*}
 n_1 &\sim P_{\N}(N_1) \\
 &\ldots \\
 n_p &\sim P_{\N}(N_p) \\
 x_1 &\gets f_1(n_1) \\
 x_2 &\gets f_2(\paG{x_2}, n_2) \\
 &\ldots \\
 x_p &\gets f_p(\paG{x_p}, n_p) \\
 \textbf{return} &\quad (x_1, \ldots, x_p)
\end{align*}

This viewpoint makes a few interesting aspects of causal modeling explicit:
\begin{enumerate}
    \item We have the choose which varaibles to include. All others are ``abstracted away'' in the noise distribution $P_{\N}$.
    \item Causal relations are given by fixed (invariant), \emph{deterministic} mathematical functions $f_1, \ldots, f_p$. All randomness (variation) comes from $P_{\N}$.
    \item The distribution (value) of each variable is fully determined by the distribution (values) of its parents and noise.
    \item We could change (\emph{intervene upon}) every ``line in the program'' without changing the causal relationships.
\end{enumerate}

Such \emph{interventions} will provide a useful way to think about causal effects.

\subsection{Interventions}

\begin{definition}
  Given an SEM $\cS = (\S, P_{\N})$, we call the replacement of one (or more) structural equations (without generating cycles) an \textbf{intervention} resulting in a new SEM $\tilde{\cS}$ with a new joint distribution, denoted by
  \begin{equation}
      \tilde{P}(\X) = P\biggl(\X \,|\, do(X_j = \tilde{f}(\tilde{\pa}(X_j), \tilde{N}_j))\biggr)
  \end{equation}
\end{definition}

Mostly, we will work with \textbf{perfect} interventions, where $\tilde{f}(\tilde{\pa}(X_j), \tilde{N}_j)$ puts a point mass on some $a \in \R$ resulting in $P(\X \, | \, do(X_j = a))$.
(Also called: ideal, structural, surgical, independent, deterministic)
In the \emph{program viewpoint} we get

\begin{align*}
 n_1 &\sim P_{\N}(N_1) \\
 &\ldots \\
 n_p &\sim P_{\N}(N_p) \\
 x_1 &\gets f_1(n_1) \\
 x_2 &\gets f_2(\paG{x_2}, n_2) \\
 &\ldots \\
 {\color{blue}x_j} &{\color{blue}\gets a} \\
 &\ldots \\
 x_p &\gets f_p(\paG{x_p}, n_p) \\
 \textbf{return} &\quad (x_1, \ldots, x_p)
\end{align*}

\begin{example}
  Consider the SEM
  \begin{align*}
      X &= N_X \\
      Y &= 4X + N_Y
  \end{align*}
  with $N_X, N_Y \sim \cN(0,1)$.
  Let's compare the marginal distribution of $Y$ in the original and an intervened model $do(X=a)$:
  \begin{align*}
      \text{original} \qquad P(Y) &= \cN(0, 17) \\
      \text{intervened} \qquad P(Y\,|\, do(X=a)) &= \cN(4a, 1) = P(Y\,|\, X=a)
  \end{align*}
  Now, let's compare the marginal distribution of $X$ in the original and an intervened model $do(Y=a)$:
  \begin{align*}
      \text{original} \qquad P(X) &= \cN(0, 1) \\
      \text{intervened} \qquad P(X\,|\, do(Y=a)) &= \cN(0, 1) \ne P(X\,|\, Y=a)
  \end{align*}
  While $X$ and $Y$ are clearly correlated, intervening on $X$ changes the distribution of $Y$ but not vice-versa.
\end{example}

The ``symmetry'' in the joint $P$ is broken when looking at interventions in SEMs.
We relied on additional information via a graph to \emph{move beyond the joint}.

\fcolorbox{black}{white}{\parbox{\textwidth}{SEMs provide a \emph{directional} mathematical model for $X_i \to X_j$.}}

\begin{definition}
  Given an SEM $\cS$, there is a \textbf{(total) causal effect} from $X_i$ to $X_j$ if and only if
  \begin{equation*}
      X_i \dep X_j \text{ in } P(\X \, |\, do(X_i=\tilde{N}))
  \end{equation*}
  for some variable $\tilde{N}$.
\end{definition}

\begin{proposition}
  Given an SEM $\cS$, the following are equivalent
  \begin{enumerate}
      \item There is a causal effect from $X$ to $Y$.
      \item There are $a, b$ such that $P(Y \,|\, do(X=a)) \ne P(Y \,|\, do(X=b))$.
      \item There is an $a$, such that $P(Y \,|\, do(X=a)) \ne P(Y)$.
      \item $X \dep Y$ in $P(X,Y \, |\,do(X=\tilde{N}))$ for any $\tilde{N}$ whose distribution has full support.
  \end{enumerate}
\end{proposition}
\begin{proof}
  $(1 \Rightarrow 2)$: If there is a causal effect, then for some $\tilde{N}$ we have $X \dep Y$ in $P(\X \,|\, do(X_i=\tilde{N}))$. Hence there exist values $a,b$ with $P(Y \,|\, do(X=a)) \ne P(Y \,|\, do(X=b))$.

  $(2 \Rightarrow 3)$: Take $a$ from the pair with different interventional distributions; then $P(Y \,|\, do(X=a)) \ne P(Y)$ unless $P(Y \,|\, do(X=a)) = P(Y)$ for all $a$, which contradicts $(2)$.

  $(3 \Rightarrow 4)$: If $P(Y \,|\, do(X=a)) \ne P(Y)$ for some $a$, then $X$ and $Y$ are dependent in the interventional joint $P(X,Y \,|\, do(X=\tilde{N}))$ for any $\tilde{N}$ with full support (dependence is preserved when the intervened variable retains positive mass on the relevant values).

  $(4 \Rightarrow 1)$: Dependence in $P(X,Y \,|\, do(X=\tilde{N}))$ implies $X \dep Y$ in $P(\X \,|\, do(X_i=\tilde{N}))$, i.e.\ a causal effect from $X$ to $Y$.
\end{proof}

\begin{proposition}
  If there is no directed path from $X$ to $Y$, then thre is no causal effect. There can be no causal effect despite a directed path. (Faithfulness violation)
\end{proposition}
\begin{proof}
  Follows from the Markov property of the interventional SEM. In a perfect intervention, all incoming arrows are being removed, causing $X$ and $Y$ to be d-separated if there is no direct path from $X$ to $Y$.
  For the second part, consider $X= N_X, Z = 2X + N_Z, Y = 4X-2Z+N_Y$.
\end{proof}

\subsection{Counterfactuals}

\begin{definition}
  Consider an SEM $\cS=(\S, P_{\N})$. Given some observations $\x$, the \textbf{counterfactual SEM} is given by replacing the distribution of noise variables
  \begin{equation*}
      \cS_{\X=\x} := (\S, P(\N \, |\, \X=\x))
  \end{equation*}
  The new set of noise variables need not be mutually independent.
  A \textbf{counterfactual statement} is a do-statement (intervention) in $\cS_{\X=\x}$.
  (This also works when we only observe a subset of the variables.
\end{definition}

 There are many confusing notations for counterfactuals.
 One possibility:
 \begin{equation*}
  P_{\underset{\text{intervention}}{\underbrace{X_j \gets a}}}(\underset{\text{var of interest}}{\underbrace{Y}} \,|\, \underset{\text{observations}}{\underbrace{\X = \x}})
 \end{equation*}

 This is also often thought of in a three step procedure:
 \begin{enumerate}
     \item \textbf{Abduction}: Update the noise distribution given the observations.
     \item \textbf{Action:} Intervene on some variable(s).
     \item \textbf{Prediction:} In the modified model, predict variable(s) of interest.
 \end{enumerate}


\begin{example}
  You decide to skip topic A when studying for your StatsLearning exam and go drinking instead ($D=1$) and then fail ($F=1$).
  Assume the correct SEM is
  \begin{align*}
      D &= N_D \\
      F &= D \cdot N_F + (1 - D) \cdot (1 - N_F)
  \end{align*}
  with $N_F \sim \mathrm{Ber}(0.05)$ (5\% chance of topic A being in the exam). ``What would have happened had you not gone drinking instead?''
  The observation is $D=1, F=1$. This fully determines the noise variables in the abduction step: $N_D = 1, N_F = 1$. This means that topic A was in the exam.
  Therefore
  \begin{equation*}
      P_{D \gets 0}(F = 0 \, |\, D=1, F=1) = 1,
  \end{equation*}
 you would not have failed had you not gone out to drink.
\end{example}


There is a ``hierarchy of understanding'' (often called the ladder of causal inference) in the following sense:
Many SEMs can entail the same observed joint distribution, but different interventional distributions.
Many SEMs can entail the same observed joint and set of interventional distributions, but different counterfactual distributions.

\begin{example}
  Let $N_1, N_2 \sim \mathrm{Ber}(0.5)$ and $N_3 \sim U(\{1, 2, 3\})$ be jointly independent.
  Compute the joint as well as all interventional distributions for the following two SEMs:
  \begin{align*}
     X_1 &= N_1 \\
     X_2 &= N_2 \\
     X_3 &= (\mybold{1}[N_3 > 0] \cdot X_1 + \mybold{1}[N_3 = 0]\cdot X_2) \cdot \mybold{1}[X_1 \ne X_2] + N_3 \cdot \mybold{1}[X_1 = X_2]
  \end{align*}
  and
  \begin{align*}
     X_1 &= N_1 \\
     X_2 &= N_2 \\
     X_3 &= (\mybold{1}[N_3 > 0] \cdot X_1 + \mybold{1}[N_3 = 0]\cdot X_2) \cdot \mybold{1}[X_1 \ne X_2] + (2 - N_3) \cdot \mybold{1}[X_1 = X_2]
  \end{align*}
  Finally, compute the counterfactual distribution
  \begin{equation*}
      P_{X_1 \gets 0}(X_3 \, |\, X_1 = 1, X_2 = 0, X_3 = 0)
  \end{equation*}
  for both SEMs.
\end{example}

For many more intriguing properties and explanations of interventions and counterfactuals keep on reading, it's super fascinating.
For an easy entry, I suggest the \emph{Book of why} by Judea Pearl.
For a concise, but similarly technical intro, look at Jonas Peters script and then move on to some of the books mentioned in the intro.

% CONTINUE HERE!


\section{Summary}
Topics covered: We tried under the framework of statistical inference cover some supervised and unsupervised learning methods, starting from linear and logistic regression to discriminant analysis, decision trees and random forests, mixture models and clustering and finally dimension reduction and latent linear models.

`Text-book´/common topics left-out:
\begin{itemize}
\item kernels, reproducing kernel Hilbert spaces
\item support vector machines/regression
\item Graphical models, DAGs (Bayes networks) and undirected ones (in particular Markov random fields/Gaussian Graphical models)
\item more on learning theory, sample boundaries of PAC learning
\item ...
\end{itemize}

Topics for a possible outlook/second-term course:
\begin{itemize}
\item deep learning (very popular)
\item data science and big data machine learning
\item applied topic: machine learning in computational biology
\item ...
\end{itemize}

% 

\bibliographystyle{plain}\bibliography{paperslibrary}
\end{document}


%\bibliographystyle{plain}\bibliography{paperslibrary,ICA,pubtheis}
%\end{document}



\subsection{Outlook to nonnegative matrix factorization \& independent component analysis} 1-2 slots, optional


\subsection{Boosting}
%maybe leave out, but it's nice, maybe shortly
Valeriya on 20th? or start directly with dim red/PCA, latent models a la murphy?


\subsection{interpretion of black box models}
murphy 16.8

\section{Graphical models (undirected, Gaussian) optional} % 2-3 slots
use Murphy chapter 10

\subsection{Modeling conditional dependencies}
\subsection{Gaussian and discrete BBN's}
\subsection{Structure learning algorithms (constraint /score based)}


\bibliographystyle{plain}\bibliography{paperslibrary}
\end{document}